Reliable uncertainty estimates in deep learning with efficient Metropolis-Hastings algorithms

Domain: biology
Published: 2026-03-17
Editorial outcome: major_revision

Abstract
Approaching problems with data-driven models often requires reliable uncertainty estimates. Bayesian neural networks can offer these for deep learning models. Without the knowledge to set informative prior distributions, sampling methods such as Hamiltonian Monte Carlo are a robust choice. However, these come with prohibitive computational costs. We study two ways to incorporate computationally light-weight Metropolis-Hastings acceptance steps into deep neural networks and stochastic gradient Hamiltonian Monte Carlo. The first method proposes noisy acceptance steps computed on batched training samples rather than the entire set during the simulation of the stochastic dynamics accepting a small minima preserving bias. The second method sacrifices bias-free sampling of Hamiltonian Monte Carlo in favor of stochastic gradient driven trajectories. While the first is analytically plausible, the second is inspired by the Hamiltonian ensemble concept. Prediction accuracy is improved by up to 5.8% over deterministic and by up to 4.3% over Bayesian approaches while still guaranteeing calibration of the predictions. We observe that sampling methods facilitate model predictions with merely a third of the ensemble while maintaining prediction accuracy. In conclusion, the methods combine efficiency and regularization of stochastic gradients, showing strong performance despite the sampling bias.

Introduction
Deep Learning has become a working horse of modern machine learning. In the presence of large datasets, state-of-the-art results are achieved in image classification, segmentation, reinforcement learning, language processing and many more1. Throughout these applications the models match or outperform humans on the same tasks. However, deep neural networks lack uncertainty awareness2,3. For applications in medicine4, autonomous driving5 and other dependable applications6 this is a serious drawback. A solution has the potential to unfold the field of high performing deep learning for tasks, that require information on reliability. In contrast to deep learning1, Bayesian learning is well known for robust modeling in presence of noise and for uncertainty estimation based on parameter distributions7,8. It is used in control or traditional statistical modeling and builds upon a well grounded theory8. The concept of Bayesian neural networks combines both methods to extract features with deep architectures and still acquire uncertainty estimations on the predictions. Even with simple uninformed prior distributions this enables complicated parameter posteriors9 in contrast to simple post-processing3 or Bayesian classifiers on top of a deep neural network6,10–12. The key to Bayesian neural networks remains the inference. Since the posterior distribution is analytically intractable, an approximation is required to estimate the posterior probability density and increased attention is required for big data applications due to practical computational limits2,6,13. The approximation is typically realized by either variational inference or Markov Chain Monte Carlo sampling. Variational inference uses a fixed parametric posterior that is optimized to minimize the Kullback-Leibler divergence to the true parameter distribution2,14,15. This approach was first published as a Bayesian neural network and uses variational formulations of backpropagation. The main limitation remains the fixed parameter posterior, which limits the diversity in the Bayesian evaluation process. Markov Chain Monte Carlo methods16 sample from the parameter posterior, such that the collected samples allow to solve the predictive distribution with Monte Carlo integration using gradient sampling,

p(y*|x*;S) = \int p(y*|x*,\theta)p(\theta|S)d\theta: (1)

Here, the predictive distribution on a new data point (y*, x*) is computed by marginalizing over the model parameters \theta that are inferred based on the dataset S = (y_j, x_j)_J. While deep ensembles17 can be considered a primitive collection of optima approximation to p(\theta|S), Bayesian sampling methods converge towards the true posterior. While early Markov Chain Monte Carlo methods are very inefficient on high-dimensional problems, Hamiltonian Monte Carlo overcame these problems by utilization of a momentum variable. Its Stan package implementation is the current state-of-the-art18–20. The problem for deep neural networks is the required gradient of the potential energy term

U(\theta) = - \sum_J log p(y_j|x_j, \theta) - log p(\theta): (2)

Every step requires gradients on the whole dataset S21. The natural solution is to turn to stochastic gradient driven dynamical systems that allow the inclusion of step-wise additional noise terms. These simulate first or second order Langevin-dynamics, which have the posterior as the stationary distribution. Typically, the dynamics

d\theta = M^{-1}m dt

dm = -\gamma m dt - \nabla_\theta U~ (\theta) dt + M^{1/2} \sqrt{2 \gamma T} dW (3)

with M \in R^{d \times d} invertible and W a standard Wiener process with variance I are computed with greed sub-sampling

\nabla_\theta U~ (\theta) = - \frac{|J|}{|B|} \sum_{j \in B} \nabla_\theta log p(y_j|x_j, \theta) - \nabla_\theta log p(\theta): (4)

This state of the dynamical system (3) consists of the model parameters and a momentum state m. We orientate our notation on22–24 that covers typical stochastic gradient Langevin dynamics25 and stochastic gradient Hamiltonian Monte Carlo26. With T=0, equation (3) becomes the time continuous version of stochastic gradient descent with momentum23. This similarity allows to reuse prior knowledge on damping \gamma or the discrete step-size \epsilon from classical deep learning trainings. The advantages of Hamiltonian Monte Carlo mainly depend on its Metropolis-Hastings acceptance step27. This eliminates a bias caused by using the proposal points computed with the dynamical system as is. To use the advantages of acceptance steps in stochastic gradient driven sampling, we consider different implementations of the discrete dynamics of (3). Therefore, three methods are chosen and Table 1 gives an overview. For other discretizations we refer to Leimkuhler et al.28 for a full overview. The algorithms differ in symmetry24,29,30 as it is required for acceptance steps. The only method to apply acceptance steps in the stochastic gradient setting was AMAGOLD proposed in ref. 29. This paper approaches to overcome the convergence problems on deep learning tasks with a new Metropolis-Hastings acceptance proposal. These acceptance steps stabilize stochastic dynamics and as they are essential for the success of Hamiltonian Monte Carlo we evaluate their success in Bayesian neural network inference. Acceptance steps are used in two methods that allow for efficient Metropolis-Hastings acceptance steps in stochastic gradient driven sampling of neural networks to make use of the robust and knowledge independent sampling qualities known from Hamiltonian Monte Carlo. Our key findings:

• We propose generalized SGHMC, a new scalable Metropolis-Hastings corrected sampling algorithm that accepts simulation steps based on batch-evaluations21,31. This method allows to run efficient step-wise acceptance while introducing an extrema preserving bias.
• We propose Hamiltonian Trajectory Ensemble (HTE), an ensemble method based on stochastic gradient trajectories that does not converge to the true posterior, but uses acceptance steps in an ensemble setting.
• We grid search hyperparameters for our algorithms as well as all competitor SGHMC methods demonstrating the critical effect of step size \epsilon in the simulation.
• We show that the Bayesian ensemble is not decisive for its accuracy, but its calibration. Therefore, we run a greedy sub-sampling experiment tried to eliminate samples from the Bayesian ensembles.

In essence, our methodology provides the full pipeline from fine-tuning the sampler to improving efficiency of test inference, demonstrating the potential of Bayesian ensembles in deep learning. Prediction accuracy and calibration scores are pushed beyond the limits of deterministic and variational Bayesian methods without the costs of classical Hamiltonian Monte Carlo.

Results
The results section is put ahead of the details on methods4. While our evaluation demonstrate the behavior of stochastic gradient sampling algorithms, we sample multivariate Gaussian distribution with diagonal covariance matrix. The algorithms run on the exact gradients of the negative log-likelihood disturbed by standard normally distributed noise n ~ N(0,1). For the set of Gaussian distributions we sample the means from an exponential distribution Exp(0.1) and the variances from a uniform distribution Unif[0.1, 1]. This is close to the expected mean-field model in many variational approaches on Bayesian neural networks6. Similar to the following experiments on image datasets, the number of samples was limited to 1000 and the tuning was done by structural variations of the step size. The results are shown in Fig. 1. We observe increasing Kullback-Leibler divergence between the true and estimated distribution, with increasing number of dimensions. While AMAGOLD outperforms all competitors at dimension 10, the other acceptance methods start to outperform AMAGOLD and SGHMC from dimension 100 on. From dimension 1000 the Kullback-Leibler divergence stays constant for the proposed algorithms.

Exploration of trajectory parameters

The experiments on image datasets focus on two different neural network architectures, which cover aspects of modern deep learning architectures, but remain reasonable small to allow ablation study experiments. A total number of 50,000 epochs per model for the first experiment only made this necessary. We divide the ablation experiments in phases showing the considered parameter dependencies. An overview is given in Table 2. The first model is a LeNet32 architecture. It represents purely feed-forward architectures based on a feature extractor part with convolutional layers and a following dense classifier. Although larger architectures from this group of models can outperform the small models, its size allows us to compute the large-scale experiments in a reasonable amount of time. Since the effect of hyperparameters in deep learning tasks can be crutial for performance, we present this study for stochastic gradient Hamiltonian Monte Carlo. The costs made us choose the second model too, which represents models with residual skip connections. We use a Wide-ResNet-10-2, which we simply refer to as ResNet throughout this section,33. While the LeNet model has about 270,000 parameters, the ResNet has about 300,000. We train the LeNet model on the balanced EMNIST dataset34 including 47 classes of MNIST-like images of different characters. The ResNet is trained on the very popular CIFAR-10 dataset35. As in Izmailov et al.36, the batch normalization layers are replaced with feature response normalization to keep the samples independent of each other37. In both cases, the default train-test split is 95%/5% training/validation data leaving 107,160 train images for EMNIST and 47,500 for CIFAR-10.

In absence of the true predictive distributions, uncertainty quantification is usually compared on performance and calibration metrics, which allow the comparison of the predicted mean and therefore with deterministic neural networks as well. We choose accuracy metric for the performance and expected calibration error (ECE) for calibration assessment as the most common ones for balanced classification tasks. The metrics are all computed based on test sets using the datasets default train/test split (CIFAR-10: 10,000; EMNIST: 18,800). Section “Bayesian neural network comparison” contains further comparison with other Bayesian inference algorithms applied to these neural networks.

With the similarity of the Euler-Maruyama discretization (cp. Eq. (3)) to stochastic gradient descent (SGD), we apply prior knowledge on momentum decay \alpha, which is usually chosen from [0.9, 1.0). We restrict the choice of \epsilon and \gamma to satisfy the constraint

(1 - \epsilon \gamma) = \alpha: (5)

For the equations in Table 1 only one free parameter remains: The step size \epsilon. For very small step sizes we propose a cosine schedule as it was referenced in Table 2. Following the constraint, the step size schedule keeps the momentum variable constant and adapts both step size \epsilon and damping parameter \gamma accordingly to increase \epsilon. This method speeds up the convergence of burn-in, before the sampling phase starts.

While we concentrate on the algorithmic parameters from phases (2) on, we start with the weight decay parameter and get an overview on loss training curves, when we apply the Langevin dynamics (cp. Eq. (3)). The weight decay regularization appears similar to the regularization through the prior in the Gaussian case differing only in a scaling factor of the training set size. In the context of deep learning we will use the term weight decay from here on. In Fig. 2 we show an subset of the generated trajectories. For different step sizes, we observe the following influence of varying weight decay: Up to weight decays 2.0 and 0.02 for the LeNet and ResNet respectively, the trajectories are very alike. Above this threshold the loss remains considerably larger and we observe occasional instabilities as in the upper left corner (cp. Fig. 2). We choose the weight decay to be at least above the stability threshold to be able to converge to useful minima for the given classification problem. Since weight decay gives considerable performance gains on deterministic methods, we argue to not choose it too small either. In the neural network architectures analyzed in this paper, the instabilities observed with weight decays did not occur in the deterministic case, what sheds light on the importance choosing this parameter for the sampling process. During our experiments we shared the results of ref. 36 that cold posteriors23 are not necessary to get competitive results. Considering T as a parameter, we fix it to T=1 throughout the following experiments (cp. Eq. (3)).

For the following simulations we run 5 parallel sampling chains with a burn-in of 1000 epochs and collect 1000 models during sampling phase after 10 epochs each, starting at a random network initialization.

Figure 3 shows the results of phases (2) and (3). The figure shows the performance measures over the step size, while the momentum variable is fixed to 0.9. The red and blue lines show the results for baseline SGHMC and SGHMC including a curvature estimate similar to RMSprop38. The two models are separated by the different line-types. We set the range of step sizes from small step sizes that show progress during burn-in until large step sizes, which caused divergence of the trajectory. Along the baseline methods, we compare both methods with our constant momentum learning rate schedule, observing top calibration scores in these simulations. For the accuracy metric, the baseline methods scores slightly better on the LeNet problem (0.1%) and a lot better on the ResNet experiment (2.8%). Most importantly, the plot shows the behavior of accuracy scores with increasing step size. In multiple cases, the top performance is achieved with the largest step size close to the stability threshold. In all cases, the step above the last shown point to increasing step size caused the dynamics to diverge, what were presented with nan. In most cases, the top performance in accuracy was also the best calibrated Bayesian ensemble.

The results of Phases (4/5) have been moved to Appendix A. We ran the top experiment from the step size variations (cp. Fig. 3) checking whether it remains the top simulation in both curvature estimation scenarios. First, we changed the momentum value to 0.99 testing its influence on the optimized sampling. Second, in the last column the batch size is changed from 80 to 400. We observe small improvements towards the larger batch size, but the changes are negligible in comparison to the effect of changing the step size of the simulation. Still, for the following simulations, we keep the larger batch size.

The experiments so far demonstrate parameter dependencies in stochastic Langevin-dynamics applied to typical deep learning tasks. The following experiments include different integrators, comparing all main methods applicable on deep learning tasks. This includes the comparison of acceptance steps methods for more accurate Markov Chain Monte Carlo trajectories.

Ablation study

To understand the effects of the different components of the algorithms used throughout this paper, this section summarizes results on SGHMC chains with different integrators as well as different acceptance step algorithms. We observe multiple interesting results shown in Table 3. First, we observe the AMAGOLD algorithm to perform poorly regarding accuracy. The need to get an acceptance rate high enough to finish the simulation occasionally limit the possible step size range. While we also tried to speed up the burn-in using our initial learning rate schedule, we report on the best results we achieved. Second, our proposed method for OBABO acceptance performs similar to its non-acceptance counterpart. While it considerably improves over AMAGOLD acceptance simulations, the question arises how much influence the unbiased trajectory has on the ensemble performance. The following section evaluates the difference by repeating simulations with different random seeds, but keeping the simulation parameters fixed. The effects of the stochastic proposal algorithm using the Hamiltonian acceptance step, results in performance boosts on EMNIST and CIFAR-10 classification. For the LeNet on EMNIST it reaches a plus of 0.7%, for the ResNet it reaches 5.2% over the Euler ensemble. At the same time, the calibration measures decreases in comparison to the best competitors, but remains in the range of the Euler simulations.

As we already discussed, the acceptance rate was one crutial factor to influence the performance metrics reported here. While there is often an optimal acceptance rate to balance exploration and exploitation, we faced some issues deciding on hyperparameters with the acceptance rate. Mostly, because the changes in step-size, that were required to change the acceptance rate where much larger than the changes we observed to change the accuracy by multiple percentage points. We decided to leave the acceptance rate aside for the tuning phase, although it was still important to finish simulations. The typical dynamics of the acceptance rate was starting at almost 100% at the beginning of burn-in and settling at a stable value for the rest of the simulation. We ended up with acceptance rates of about 50% for both our algorithms with the most promising hyperparameter configuration regarding the evaluation metrics. While this important quality of the sampling algorithms is mainly left aside for the proposed evaluation, we extend the discussion in Appendix C.

Bayesian neural network comparison

This section compares the results on sampling algorithms in this paper to different Bayesian neural network approaches. In particular, we compare against regular deep ensembles and optimization based variational inference approaches as well as deterministic trainings. We include two different methods for variational inference on Bayesian neural networks: The Bayes by backpropagation approach2,39 extended by the flipout method40 and the variational online Gauss Newton (VOGN) method proposed in ref. 41 and proven to be competitive in ref. 42. We compare based on accuracy and expected calibration error again, but also report maximum calibration error and the number of epochs to account for the computational costs independent of the hardware setup. For an improved overview on the field MC-dropout, stochastic weight averaged Gaussian (SWAG) and Cyclical learning rate schedulers are included10,12,43. In this section, we repeat the simulations 5 times with different random seeds to report on statistical variation, which was not possible for all ablation experiments.

In addition to the hyperparameter tuning study, we add a regression experiment of time series forcasting for this final comparison. We experiment with the chaotic Lorenz’96 differential equation system, which is a common benchmark for time series forcasting problems. The system is simulated with 5 states, of which 50 steps are used for intializing the model and 50 are predicted for computing the loss. The evaluation predicts 150 time steps ahead, on which we compare maximum absolute error (MAE), mean squared error (MSE) and the R2 score. The model is a one layer LSTM with ReLU activation with a state-space of 100. Although this problem was not included in the hyperparameter study, we observe similar results, what demonstrates the generality of the results. Comparing the results in Table 4, we see the top results are achieved with the Bayesian sampling algorithms. Still, we observe differences between the problems: While the ensemble created from different model initializations does consistently outperform the deterministic baseline, it sets the strongest baseline in the time series forcasting case. The HTE algorithm does only match this baseline, while the other Bayesian inference algorithms could not reach this performance. Also, the two classification problems appear to be quite different. While the EMNIST classification benefits from all sampling algorithms equally, the CIFAR-10 classification mainly benefits from the HTE mechanism. We approach an explanation in the discussion. The results confirm a performance increase of Markov Chain Monte Carlo algorithms over variational inference based methods on the considered problems including the Bayesian neural networks. The costs for this performance remains the increased computational burden in comparison to variational inference methods. Although the proposed sampling methods considerably reduce the costs in comparison to Hamiltonian Monte Carlo experiments as in ref. 36, the number of epochs against the deterministic training is increased considerably, as we show in the column reporting the number of epochs these algorithms ran on average. These costs need to be balanced against the performance and calibration increase from Table 4.

We also report on out of distribution detection on two additional datasets. For the EMNIST classification we compute the predictive entropy for the train dataset of Fashion MNIST44. For the CIFAR-10 classification we use the test split of the Street View House Numbers45 dataset. The results are visualized in Fig. 4 by showing the sorted entropy scores over the out of distribution dataset. The shaded area show the standard deviation over the five random seed repetitions from Table 4. The optimal curve would be returning the maximum entropy for all data points, so the higher the curve, the better the out-of-distribution detection. We observe HTE to outperform its competitors, although not being comparably calibrated. The GSGHMC curve is at least partly above the SGHMC baseline, but does not meet the HTE line.

Ensemble size reduction

The ensembles of 1000 members make the evaluation and prediction process considerably more costly than running the inference of a single model. That raises the question how many models need to be kept to maintain the performance achieved with the large ensemble, but reducing the cost of evaluation. Therefore, we run a simple greedy selection on the ensembles, keeping only samples that appear to be essential for a given metric not to be impaired. We visualize the results in Fig. 5. The number of remaining ensemble members is displayed on the y-axis after the reduction of ensembles from the prior section. The variance is computed both over the different random seeds and 5 greedy repetitions with different ordering each. We can see that the Accuracy and Calibration have quite different requirements with regard to the ensemble size. While preserving the accuracy requires less than 50% of the ensemble, the calibration requires almost all of it. While there are methods to compress Bayesian ensembles in a principled manner46, we wanted to demonstrate, that the Bayesian sampling is not necessarily the last step before having a model ready for production. The advances in accuracy are less costly to preserve than the diversity needed for calibration. That applying a reduction technique can still be favorable is further discussed in Appendix F.

Methods
To introduce the concepts of Metropolis-Hastings acceptance we review the key concepts in the following and introduce the background to the two sampling algorithms, which have been evaluated in this paper. Afterwards, we introduce the generalized stochastic gradient Hamiltonian Monte Carlo algorithm and the Hamiltonian trajectory ensemble.

Errors in Markov Chain Monte Carlo estimation

There are multiple relevant errors that occur during the numeric approximation of the posterior probability distribution with Markov Chain Monte Carlo28. However, a moment matching at time step t does not exactly gather the error the user of a sampling method is interested in. For Markov Chains the long term error is most relevant, to compute correct Monte Carlo expectations. Therefore, the representation of p(\theta|S) through the ensemble is what matters. It consists of an error caused by discretization (cp. Eq. (3)), which is proportional to the step size. Also, there is an error due to the limited number of samples collected during sampling phase. These errors are often referred to as sampling bias and sampling error, respectively. The error or bias caused by the discretization can not be reduced by running the simulation for a longer time. On the other hand, a reduction of the step size reduces the mixing of stored Markov Chain Monte Carlo samples. Similar to a bias-variance trade-off, these two errors balance each other. Simulations with a constant number of steps do either precisely simulate the dynamics or gather more samples from a larger parameter space to represent the posterior accurately.

Previous acceptance step approaches

The sampling bias of Markov Chain Monte Carlo simulation can be reduced by considering the endpoints of trajectories proposal points z_p, which are accepted or rejected in a Metropolis-Hastings manner27. We introduce z as an auxiliary state, which summarizes the concatenation of \theta and m, in case the sampler utilizes a momentum state. The transition probability from initial value z_0 = (\theta_0, m_0) to z_p = (\theta_p, m_p) is required to fulfill t(z_0 \to z_p) = t(z_p \to z_0) for an unique stationary distribution, which is called detailed balance condition. The Metropolis-Hastings acceptance probability takes the form

r(z_0, z_p) = min(1, \frac{\pi(z_p)t(z_p \to z_0)}{\pi(z_0)t(z_0 \to z_p)}): (6)

We observe that non reversible Markov Chain Monte Carlo chains directly cause acceptance rates of zero. This applies in particular to Euler-Maruyama29. If the backward transition probability is not zero, the acceptance step guarantees reversibility by design, if the momentum is completely resampled. An advantage of acceptance steps is the stabilizing effect on the dynamics. Since stochastic gradients can prevent convergence for integrators of any order63, the application of acceptance steps allow to run more robust sampling chains. In the following, we present the direct competitor approach of stochastic gradient acceptance steps for sampling and its limitations in the context of large models. Zhang et al.29 propose AMAGOLD, an accumulating strategy for r(z_0, z_p) using the observation that we can accept proposal points based on the product of probability densities from all state transitions between z_0 and z_p. With the log-probability this becomes a sum, which we can compute step-by-step in the minibatch-setting with little extra costs. They use the skew-symmetric integrator from Table 1 and build a skew-reversible chain, which has non-zero acceptance rates (cp. Alg. 2). The algorithm allows to explore the posterior based on stochastic gradients, but requires U(\theta) for accepting the proposal state. We run the algorithm on our test cases in all configurations tested for any other method as in Table 3 with the results in Fig. 3. In combination with Bayesian neural networks we observe small acceptance rates below 10%, if the step size is kept in a reasonable range that guarantees the optimization to move forward. To get into acceptance rates of about 30% we had to reduce the trajectory length to 2 epochs and keep the step size at 10^{-4}, which was with factor 100 considerably smaller than without the AMAGOLD acceptance steps. The results from section “Ablation study” shows the problem of efficient acceptance steps was not overcome for Bayesian neural networks. The current solution even limit prediction accuracy. While AMAGOLD was evaluated on deep neural network cases, other competitors could not be evaluated on Bayesian deep learning problems in the literature. Efficient Metropolis-Hastings acceptance steps for Markov Chain Monte Carlo methods have also been approached in refs. 47–51. While the application to Bayesian deep learning is not in the focus of these studies, some problems occur evaluating them on the highly non-convex log-likelihoods in deep learning. The most prominent problem are variable batch-sizes, which slow down the inference of Bayesian neural networks considerably. These appear in refs. 47–49,51. The work of ref. 50 even requires occasional resampling on the data, which also limit the speed of each step and cannot be implemented with a fixed computational budget. Additionally, the batch-size is a relevant regularization hyperparameter. If it is not fixed, the knowledge transfer from regular neural network trainings is quite limited. Therefore, these methods have not been tested on any Bayesian deep learning problems yet. Together with the methological limits, these methods were excluded from the study in this paper, similar to previous studies on sampling algorithms for Bayesian deep learning36. We approach an algorithm that fixes the required batch-size for the gradients and acceptance step fixing the computational budget and therefore being especially compatible with deep neural networks. To do so, we rather orientate on partial momentum updates in contrast to the majority of these methods, which replace an acceptance step in an momentum resampling method. Therefore, we choose to call the method generalized SGHMC in accordance with generalized HMC.

Our proposed Metropolis-Hastings acceptance methods

Our two approaches aim to facilitate an application perspective. The first proposed method uses the similarity between classical SGHMC and Generalized HMC to construct an unbiased SGHMC algorithm. The second method is motivated from the long trajectory length parametrizations from HMC simulation experiments as in ref. 36. We ask ourselves, whether the posterior influenced by stochastic gradients might actually cause beneficial regularization, as it is well known for large ensembles64. The exploration properties of HMC, which allows to explore the posterior with the exact posterior sampling property. This section presents both algorithms, which have been experimentally analyzed.

Generalized SGHMC (GSGHMC)

The generalized SGHMC methodology uses step wise acceptance steps and is therefore composed of an batch-acceptance algorithm and an integrator, that does not result in an acceptance rate of zero. A batch-acceptance algorithm with a conservation of minima property is proposed in ref. 65. They introduce a clever notation to analyze the properties of their algorithm, which we will also use, but find the acceptance rates to be relatively small throughout their experiments. Our algorithm bases on their idea, fixing the acceptance rate by using a reversible integrator, such as OBABO. The preservation of local optima property of the algorithm is required in our case to explore the extrema of the true posterior exactly. In ref. 65 the minibatch acceptance setting is considered by an state space extension with a binary random vector \tau \in {0,1}^{|J|} such that B = {j: \tau(j)=1, j \in J}. This allows to reformulate U~ (\theta) to

U~_\tau (\theta) = - \frac{|J|}{|B|} \sum_{j \in J} \tau_j log p(y_j|x_j, \theta) - log p(\theta): (7)

uniform distribution over I_\tau = {\tau: \sum_{j \in J} \tau_j = |B|}, which we denote by Unif_B(\tau):

\pi_e(z) / e^{-c_n U~_\tau (\theta)} e^{-K(m)} Unif_B (\tau): (8)

The additional distribution over the momentum remains independent of \theta and \tau. With the Markov Chain in equation (6), we construct the augmented target distribution and observe that the distribution on \tau does not change from initial to proposal state and does not influence the relation. The acceptance probability becomes

r = min(1, \frac{\pi_e(z_p, \tau_p) t(z_p, \tau_p \to z_0, \tau_0)}{\pi_e(z_0, \tau_0) t(z_0, \tau_0 \to z_p, \tau_p)}) = min(1, \frac{e^{-c_n U~_{\tau_p} (\theta_p)} e^{-K(m_p)} t(z_p \to z_0)}{e^{-c_n U~_{\tau_0} (\theta_0)} e^{-K(m_0)} t(z_0 \to z_p)}) (9)

where each step is accepted or rejected based on a second minibatch that is selected with \tau_p. To reconstruct the posterior distribution \pi(\theta), we have to marginalize over \tau. Using the independence of the distribution on m, we directly compute the distribution on \theta,

\pi_e(\theta) / \pi(\theta) c_n \frac{1}{|B|} \sum_{\tau \in I_\tau} e^{c_n (-U~_\tau(\theta) - U(\theta))}, (10)

which slightly differs from the distribution in ref. 65 due to our different definition of U. The parameter c_n balances the tempering of the posterior distribution against resulting acceptance rates and overall bias. While in the proposing paper c_n = 20 appears as a generally good choice to trade off acceptance rate against tempering effect65, we tried to stay at c_n = 1 during our simulation. By slightly increasing the step size in comparison to the SGHMC simulation without the acceptance step, our experiments ran with acceptance rates of 50.0% \pm 3.5% both in the benchmark and the neural network case. We decided to keep c_n fixed at this value that our algorithm results in an unbiased chain, approximating the exact posterior. The independence of m guarantees the convergence results from ref. 65, which are shortly reviewed in appendix E.

The acceptance rate of a sampling algorithm often decides, whether it can be applied on large neural networks. We discuss some observations based on equation (9) that what to expect from our GSGHMC algorithm. Therefore, we take the logarithm of the acceptance probability and plug in the transition probability for OBABO as in ref. 24

log \frac{t(z_0 \to z_p)}{t(z_p \to z_0)} + log exp(-c_n U~_{\tau_p} (\theta_p)) exp(-K(m_p)) = - (K(m_0) - K(m_p) + K(m_{3/4}) - K(m_{1/4})) (11)

- c_n U~_{\tau_p} (\theta_p) + c_n U~_{\tau_0} (\theta_0) - K(m_p) + K(m_0) (12)

= c_n (U~_{\tau_0} (\theta_0) - U~_{\tau_p} (\theta_p)) - K(m_{3/4}) + K(m_{1/4}): (13)

This equation demonstrates the effect of c_n as it increases the acceptance rate, as long as the loss improves. The K terms depend on the gradients. We can write K(m_{3/4}) - K(m_{1/4}) as a function of m such that the gradient is \nabla_\theta U~_{\tau_p} (\theta_p) + \nabla_\theta U~_{\tau_0} (\theta_0), while the K terms are \epsilon > 0 = \nabla_\theta U~_{\tau_p} (\theta_p) - \nabla_\theta U~_{\tau_0} (\theta_0): (14)

The first term in this equation shows as long the series of gradients agree on the direction as they are captured in m_{1/4}, the first term is large and dominates the second term, which is scaled by the square of the small step size. This behavior is valuable as a method to prevent U-turns in the parameter trajectory was proposed earlier in ref. 13. In case of unstable gradients and small changes in the values of U, the acceptance behavior will mostly depend on the random choice.

Algorithm 1. Hamiltonian Trajectory Ensemble (HTE)

HTE(N, U, K, \theta_0, t_L, \epsilon, M, integrator)

S = { } \triangleright Posterior samples

\theta = \theta_0 \triangleright Initialize position

while |S| < N do

m ~ N(0, M) \triangleright Resample momentum

z = {\theta, m} \triangleright Store initial value

for t = 0 to t_L do

integrator.step(\theta, m, \epsilon, M, U~, T=0)

endfor

r = min{1, exp(U(z(\theta)) - U(\theta) - K(m) + K(z(m)))}

\nu ~ Unif(0, 1)

if r > \nu then

S = S \cup {\theta} \triangleright Accepted!

else

\theta = z(\theta) \triangleright Rejected!

endif

endwhile

Hamiltonian trajectory ensemble (HTE)

The second method we propose is motivated by the results of 36 on Hamiltonian Monte Carlo. To achieve overarching performance improvements, the Hamiltonian dynamics are simulated over a long trajectory length up to a thousand steps, using effectively gradient descent with momentum between the momentum resets, as this is what the trajectory simulation of HMC does. Although the integrators from Table 1 with stochastic gradients cannot simulate the same trajectory as the one followed by the true Hamiltonian dynamics, we know from optimization literature that the risk will reduce similarly to a gradient descent optimization (cp. 64), which is performed when simulating the Hamiltonian dynamics in between the momentum resampling steps. Therefore, we decided to try reproducing the HMC performance benefits with an OBABO stochastic gradient trajectory, still applying the same acceptance used in HMC. While this introduces an uncontrolled bias as it was discussed in ref. 66, the divergence of the trajectories is the main concern. The aim of our method is to accept only trajectories, that stayed close to the actual Hamiltonian trajetory and are therefore excepted by the Metropolis-Hastings step. The OBABO integrator is used with T=0, effectively simulating a leapfrog discretization21 with an additional momentum decay (cp. 67). A summary of the program realization is given in Alg. 1. We use a program-like notation, where the formulas of the integrator steps are implemented in a step method, which is applied to the variables in place for every trajectory step as given in Table 1. Evaluation of U refers to equation (2) and K to the quadratic kinetic energy term depending on m and M. The notation z(\cdot) refers to the choice of \theta or m from z.

Our primary objective using this naive approach was to sample from the same likelihood, which is optimized during deep neural network training. Through the usage of stochastic gradients small minima are hidden and the likelihood gets biased towards minima that facilitate generalization64. Taking the acceptance probability for the GSGHMC algorithm in formula (14), the momentum terms do differ by the gradient noise. Therefore, in contrast to Hamiltonian Monte Carlo, the Kinetic energy terms do not eliminate each other. With a quadratic kinetic energy term, there is a sum of kinetic energy terms of the noise as well as scalar product terms between the momentum and the gradient noise that we explicitly ignore and that will deformate the likelihood depending on the energy in the noise. While this was argued in ref. 66 to eventually cause acceptance rates of zero, we could not empirically validate this for the gradient noise evaluated in this paper. While this does not eliminate the bias of this method, it explores a regularized log-likelihood. Applying the algorithm on neural networks, we showed to conserve benefits from regularized neural network training.

The algorithm showed surprising robustness in the neural network setting, when the dimension of the target distribution increases, a property that is observed for regular HMC21 too. Since this method is not a Bayesian inference, but an ensemble method, we decided to use the term Hamiltonian trajectory ensemble (HTE). Advantages over regular ensemble methods from the machine learning community68 are generated through a speed up, since it is not required to retrain the complete model in different settings. By running the trajectory length after an initial burn in creates large ensembles considerably faster than regular ensemble methods, saving time proportional to the ensemble size.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Results\nStrengths\n- The proposed methods demonstrate improvements over existing approaches in terms of both accuracy and uncertainty\nestimation.\n- The authors conduct an extensive ablation study to assess the effect of hyperparameters.\nWeaknesses\n- Many experimental results lack error bars, including Tables III, IV, and V. Since the performance of deep neural networks\ncan be influenced by random seed selection, it is crucial to include error bars to assess whether the reported improvements\nare statistically significant.\n- The paper lacks a comparison to mini-batch-based Metropolis-Hastings methods such as AustereMH [1] and MHminibatch\n[1]. Especially, AustereMH has been used with stochastic gradient MCMC as shown in the original paper. Since AMAGOLD\nis asymptotically unbiased and the proposed methods appear to be asymptotically biased, these methods, which are\nasymptotically biased, still serve as reasonable baselines for comparison.\n- A cost comparison in the experiments would be valuable to better understand the trade-offs between computational\nefficiency and performance.\nMethods\nStrengths\n- The proposed methods are intuitive and easy to implement in practice.\nWeaknesses\n- It is unclear whether the proposed GSGHMC method converges to the target distribution asymptotically. Although the\nauthors claim \"bias-free sampling\" in the abstract, my understanding is that the method introduces bias due to the use of\nmini-batching in the acceptance step. If this is the case, it would be beneficial to analyze the magnitude of the asymptotic\nbias. A similar concern applies to the proposed Hamiltonian Trajectory Ensemble (HTE).\n- The novelty of the proposed methods compared to existing approaches is not clearly established. Mini-batching or\nsubsampling in the Metropolis-Hastings step is a well-known technique (e.g., [1,2,3,4]). It remains unclear what distinguishes\nthis work from prior methods.\nProviding a convergence rate analysis would strengthen the theoretical contributions of the paper, as such\nguarantees are commonly expected for MCMC-based methods.\nSummary\nThe experimental results suggest promising improvements over existing methods. However, key concerns remain regarding\nthe statistical significance of these improvements, the novelty of the proposed approaches, and their theoretical properties.\nAddressing these issues would improve the strength and clarity of the paper.\n[1] Anoop Korattikara, Yutian Chen, and Max Welling. Austerity in MCMC land: Cutting the Metropolis-Hastings budget. In\nInternational Conference on Machine Learning, 2014\n[2] Daniel Seita, Xinlei Pan, Haoyu Chen, and John Canny. An efficient minibatch acceptance test for Metropolis-Hastings.\nUncertainty in Artificial Intelligence, 2017\n[3] Zhang, Ruqi, A. Feder Cooper, and Christopher De Sa. Asymptotically Optimal Exact Minibatch Metropolis-Hastings.\nAdvances in Neural Information Processing Systems, 2020\n[4] Dougal Maclaurin and Ryan Prescott Adams. Firefly Monte Carlo: Exact MCMC with subsets of data. In Twenty-Fourth\nInternational Joint Conference on Artificial Intelligence, 2015", "reviewer_2": "# Summary\nThe paper presents two novel Metropolis-Hastings-based methods to improve sampling efficiency in Bayesian neural\nnetwork (BNN) inference using stochastic gradient Hamiltonian Monte Carlo (SGHMC). The motivation arises from the\ncomputational impracticality of exact Hamiltonian Monte Carlo (HMC) for large-scale deep learning models. The authors aim\nto address the trade-off between inference accuracy, uncertainty calibration, and computational efficiency by incorporating\napproximate Metropolis-Hastings corrections into SGHMC dynamics.\nThe first proposed method, *Generalized SGHMC (GSGHMC)*, introduces minibatch-based acceptance steps in stochastic\ngradient-driven simulations, approximating the true posterior with tempered acceptance criteria while maintaining detailed\nbalance through a reversible integrator (OBABO). This method is grounded in a theoretical formulation extending reversible\nMarkov Chain Monte Carlo (MCMC) with stochastic gradients, relying on batch-tempered acceptance schemes with\nguarantees on convergence and local mode preservation under smoothness and regularity assumptions.\nThe second method, *Hamiltonian Trajectory Ensemble (HTE)*, forgoes unbiased sampling in favor of long stochastic\ntrajectories and Metropolis-like selection steps after deterministic integration with zero-temperature OBABO dynamics. While\nit does not converge to the true posterior, it produces ensembles with strong mixing properties and high predictive accuracy.\nThe authors benchmark both methods against SGHMC, AMAGOLD, and classical Bayesian neural network baselines\n(variational inference, variational online Gauss-Newton, and deep ensembles) on synthetic Gaussian benchmarks and real\nimage classification tasks (LeNet/EMNIST and ResNet/CIFAR-10). They show that GSGHMC achieves the best calibration\n(ECE down to 0.010) while HTE achieves the highest accuracy (up to 93.7% on EMNIST), outperforming existing sampling-\nbased methods. They also demonstrate that ensemble performance plateaus with significantly reduced ensemble sizes\n(down to ~20% of original), thereby reducing inference costs without sacrificing accuracy.\nThe work demonstrates that with properly tuned Metropolis-Hastings acceptance mechanisms and careful integration of\nstochastic gradients, one can achieve scalable, high-performing BNN inference pipelines without the prohibitive\ncomputational burden of full HMC, while retaining or improving uncertainty calibration and predictive accuracy.\n# Comments\n### Calibration Trade-offs in HTE\nTable V shows that the Hamiltonian Trajectory Ensemble (HTE) method achieves the highest accuracy across both\nLeNet/EMNIST and ResNet/CIFAR-10 settings, but this comes at the cost of degraded calibration. While the Expected\nCalibration Error (ECE) for HTE (0.027 for LeNet, 0.047 for ResNet) remains comparable to deterministic baselines, it is\nconsistently worse than all other Bayesian approaches evaluated. Given that HTE does not converge to the true posterior, it\nis plausible that the increased accuracy derives from overconfident but miscalibrated predictions concentrated around a\nsubset of low-energy modes explored through long stochastic trajectories.\nDid the authors investigate whether this drop in calibration is systematic across all classes, or is it localized to particular\nregions of the input space? Additionally, were any calibration-improving post-processing techniques (e.g., temperature\nscaling) considered for HTE, and if so, what impact did they have on the observed trade-off between calibration and\naccuracy? How about out-of-distribution performance?\n### Task and Architecture Generality\nThe empirical results presented in the paper focus exclusively on image classification tasks using convolutional neural\nnetworks (CNNs), specifically LeNet on EMNIST and a small ResNet variant on CIFAR-10. While these benchmarks are\nstandard, the restriction to vision-based classification is a limiting factor, especially given the paper's strong emphasis on\nempirical validation. The architectural choices exclude evaluation on non-convolutional models, such as transformers, and\nno experiments are reported on regression or structured output tasks, where calibration and uncertainty quantification play an\nequally critical role.\nGiven that stochastic gradient-based MCMC methods are broadly applicable, it would be informative to test GSGHMC and\nHTE in alternative settings—such as time series forecasting, tabular regression, or even natural language tasks—to assess\ntheir robustness to domain shift, heteroscedastic noise, and varying data modalities. Were any such evaluations attempted\nor preliminary results obtained in non-image domains?\n### Acceptance Rate Dynamics and Interpretation\nThe paper reports an average acceptance rate of approximately 50% ± 3.5% for both the GSGHMC and HTE methods\nduring experiments (Section II-C). While this figure reflects a moderate trade-off between exploration and rejection, it looks\nsurprisingly low for methods based on SG-MCMC. A detailed temporal analysis of how the acceptance rate evolves during\nburn-in and sampling would be valuable—e.g., does the rate stabilize, decline, or exhibit oscillatory behavior? Furthermore,\nthe relationship between acceptance rate and predictive metrics (e.g., accuracy, NLL, calibration) is not explicitly addressed.\nIs there evidence that lower acceptance rates correlate with more informative exploration of wide minima, or do they instead\nreflect miscalibrated proposals due to aggressive step sizes? Were different acceptance thresholds tested, or were\ndiagnostics (e.g., autocorrelation, effective sample size) computed to assess mixing quality as a function of acceptance\nbehavior over time?\n### On the Role of Subsampling and Ensemble Compression\nThe greedy subsampling procedure proposed in Section II-D, while practical, appears conceptually orthogonal to the core\ncontributions of the paper. The main message of the work centers on improving the fidelity and efficiency of posterior\nsampling via novel Metropolis-Hastings corrections in stochastic gradient settings. Subsampling an ensemble post hoc for\ntest-time speedup is a separate compression strategy that could in principle be applied to any Bayesian ensemble—\nregardless of how it was generated—and thus may not align tightly with the specific methodological innovations of\nGSGHMC or HTE. The observation that performance remains stable even when reducing ensembles by 60–80% is\ninteresting, but does not inform the efficacy of the sampling method itself.\nMoreover, there exist more principled approaches for compressing Bayesian predictive distributions, starting with Bayesian\nDark Knowledge (Korattikara et al., 2015). Did the authors consider or compare such model distillation approaches to their\nsubsampling baseline? It would be useful to contextualize their method within the broader literature on posterior\ncompression.\n## Minor Comments\n### Paper Structure and Readability\nThe current structure—presenting results and discussion before the methods—makes the paper challenging to follow,\nespecially given its technical density and the mathematical formulation of the proposed sampling algorithms. Without a clear\nexposition of GSGHMC and HTE prior to their empirical evaluation and discussion, it is difficult to contextualize the core\ncontributions or understand how specific hyperparameter studies relate to the underlying dynamics. While the reverse\nordering may be justified for narrative reasons (e.g., to emphasize empirical gains), the technical nature of the work suggests\nthat a conventional structure placing methods before results might better support clarity.\nOf course, this may be a matter of personal preference, but I feel that the paper would benefit from a more standard structure."}, "author_response": "Technische Universität Ilmenau\nData intensive Systems and Visualization Group (dAI.SY)\nMax-Planck-Ring 14\n98684 Ilmenau\nPoint-by-point Response\nReviewer Nr. 1\nResults\nStrengths\n• The proposed methods demonstrate improvements over existing approaches in terms of both\naccuracy and uncertainty estimation.\n• The authors conduct an extensive ablation study to assess the e!ect of hyperparameters.\nWeaknesses\n• Many experimental results lack error bars, including Tables III, IV, and V. Since the performance\nof deep neural networks can be influenced by random seed selection, it is crucial to include error\nbars to assess whether the reported improvements are statistically significant.\nResponseR1.0: Thanksforraisingthisimportantconcern. Wefullyagreewiththeneedofseeded\nrepetitionsandaccordinglyextendedtable4withrepetitionson5seeds, reportingthestatisticsby\nmean and standard deviation. We also include SGHMC in table 4 to demonstrate the di!erence\nbetween SGHMC and GSGHMC. Running multiple seeds during the hyperparameter study was\nnot possible due to the enormous cost of the evaluation in general. However, to our knowledge\nthere are no reported cases on overfitting to a specific seed and the comparison remains fair.\n• The paper lacks a comparison to mini-batch-based Metropolis-Hastings methods such as Aus-\ntereMH [1] and MHminibatch [1]. Especially, AustereMH has been used with stochastic gradient\nMCMC as shown in the original paper. Since AMAGOLD is asymptotically unbiased and the\nproposed methods appear to be asymptotically biased, these methods, which are asymptotically\nbiased, still serve as reasonable baselines for comparison.\nResponse R1.1: Thanks for this valuable suggestion. While we could not include the method for\nallexperiments, weextendeditsdiscussionfromlines765-782. Inessence, thevariableacceptance\ncost with an additional while loop over the data proposes a serious drawback for Bayesian deep\nlearning. Therefore, it was not applied on Bayesian deep learning yet and remains a theoretically\nclose, but empirically less established baseline.\n• A cost comparison in the experiments would be valuable to better understand the trade-o!s\nbetween computational e\"ciency and performance.\nResponse R1.2: Thanks for this suggestion. We include this by reporting the average number of\ntotalepochsperexperimentintable4. Whilethedurationofanepochisfixedonahardwaresetup,\nthis enables to compare independently of hardware requirements. This is additionally discussed in\nlines 450-452, 497-499 and 564-566.\nMethods\nStrengths\n• The proposed methods are intuitive and easy to implement in practice.\nMatthias Schmal\n2\nmatthias.schmal@tu-ilmenau.de\nTechnische Universität Ilmenau\nData intensive Systems and Visualization Group (dAI.SY)\nMax-Planck-Ring 14\n98684 Ilmenau\nWeaknesses\n• It is unclear whether the proposed GSGHMC method converges to the target distribution asymp-\ntotically. Although the authors claim \"bias-free sampling\" in the abstract, my understanding is\nthat the method introduces bias due to the use of mini-batching in the acceptance step. If this\nis the case, it would be beneficial to analyze the magnitude of the asymptotic bias. A similar\nconcern applies to the proposed Hamiltonian Trajectory Ensemble (HTE).\nResponse R1.3: Thanks for this valuable suggestion. While it was not our intention to mix up\nterminology we stressed the sampling bias for our methods in the abstract to prevent wrong\nexpectations in lines 31, 32 and 39. Additionally, we included the bias in the bullet point list\n(introduction,line151). Weexpandedthediscussiontoequation(8)showingthebiasofGSGHMC\nwith its dependencies on the parameter c and the batch size. We extend the discussion regarding\nn\nthe bias of HTE in lines 921-924 and 931-944. In theory, we should not be able to quantify the\nbias of HTE without additional assumptions on the gradient noise (cp. [5]) since the divergence\nof the trajetories is required. A side e!ect of the regularization through gradient noise is a better\ngeneralizing model, which is in line with our observation of improved model performance and\nseems to be conveyed to the HTE method.\n• The novelty of the proposed methods compared to existing approaches is not clearly established.\nMini-batching or subsampling in the Metropolis-Hastings step is a well-known technique (e.g.,\n[1, 2, 3, 4]). It remains unclear what distinguishes this work from prior methods. Providing a\nconvergence rate analysis would strengthen the theoretical contributions of the paper, as such\nguarantees are commonly expected for MCMC-based methods.\nResponse R1.4: Thanks for raising these important concerns. We extend the discussion on the\ncurrentliteraturefromlines765-782andadditionallystressthechoiceofourmethods. Wealready\ncovered[1,2,3]inbriefbeforeandnowextendedthis. TheconvergenceofGSGHMCisaddressed\ninsection4.4andequation(10). Thenecessarytheoremsarenotedintheappendix. ForHTEthis\nissue is addressed in the changes we applied at lines 921-924 and 931-944. Furthermore, we note\nthat the main contributions of this paper are in the experimental evaluation of methods in the\nfield of deep learning rather than presenting a new baseline for any Bayesian inference problem.\nWe argue, that these tasks require a special treatment barely recognized in the literature so far.\nSummary\nThe experimental results suggest promising improvements over existing methods. However, key con-\ncerns remain regarding the statistical significance of these improvements, the novelty of the proposed\napproaches, and their theoretical properties. Addressing these issues would improve the strength and\nclarity of the paper.\nReviewer Nr. 2\nReviewer 2 (Remarks to the Author):\nSummary\nThe paper presents two novel Metropolis-Hastings-based methods to improve sampling e\"ciency in\nBayesianneuralnetwork(BNN)inferenceusingstochasticgradientHamiltonianMonteCarlo(SGHMC).\nMatthias Schmal\n3\nmatthias.schmal@tu-ilmenau.de\nTechnische Universität Ilmenau\nData intensive Systems and Visualization Group (dAI.SY)\nMax-Planck-Ring 14\n98684 Ilmenau\nThe motivation arises from the computational impracticality of exact Hamiltonian Monte Carlo (HMC)\nfor large-scale deep learning models. The authors aim to address the trade-o! between inference ac-\ncuracy, uncertainty calibration, and computational e\"ciency by incorporating approximate Metropolis-\nHastings corrections into SGHMC dynamics.\nThe first proposed method, Generalized SGHMC (GSGHMC), introduces minibatch-based accep-\ntance steps in stochastic gradient-driven simulations, approximating the true posterior with tempered\nacceptance criteria while maintaining detailed balance through a reversible integrator (OBABO). This\nmethodisgroundedinatheoreticalformulationextendingreversibleMarkovChainMonteCarlo(MCMC)\nwith stochastic gradients, relying on batch-tempered acceptance schemes with guarantees on conver-\ngence and local mode preservation under smoothness and regularity assumptions.\nThe second method, Hamiltonian Trajectory Ensemble (HTE), forgoes unbiased sampling in favor\nof long stochastic trajectories and Metropolis-like selection steps after deterministic integration with\nzero-temperature OBABO dynamics. While it does not converge to the true posterior, it produces\nensembles with strong mixing properties and high predictive accuracy.\nThe authors benchmark both methods against SGHMC, AMAGOLD, and classical Bayesian neu-\nral network baselines (variational inference, variational online Gauss-Newton, and deep ensembles) on\nsyntheticGaussianbenchmarksandrealimageclassificationtasks(LeNet/EMNISTandResNet/CIFAR-\n10). They show that GSGHMC achieves the best calibration (ECE down to 0.010) while HTE achieves\nthehighestaccuracy(upto93.7%onEMNIST),outperformingexistingsampling-basedmethods. They\nalso demonstrate that ensemble performance plateaus with significantly reduced ensemble sizes (down\nto 20% of original), thereby reducing inference costs without sacrificing accuracy.\nThe work demonstrates that with properly tuned Metropolis-Hastings acceptance mechanisms and\ncareful integration of stochastic gradients, one can achieve scalable, high-performing BNN inference\npipelines without the prohibitive computational burden of full HMC, while retaining or improving un-\ncertainty calibration and predictive accuracy.\nComments\nCalibration Trade-o!s in HTE: Table V shows that the Hamiltonian Trajectory Ensemble (HTE)\nmethod achieves the highest accuracy across both LeNet/EMNIST and ResNet/CIFAR-10 settings, but\nthis comes at the cost of degraded calibration. While the Expected Calibration Error (ECE) for HTE\n(0.027 for LeNet, 0.047 for ResNet) remains comparable to deterministic baselines, it is consistently\nworse than all other Bayesian approaches evaluated. Given that HTE does not converge to the true pos-\nterior, itisplausiblethattheincreasedaccuracyderivesfromoverconfidentbutmiscalibratedpredictions\nconcentrated around a subset of low-energy modes explored through long stochastic trajectories.\nDid the authors investigate whether this drop in calibration is systematic across all classes, or\nis it localized to particular regions of the input space? Additionally, were any calibration-improving\npost-processing techniques (e.g., temperature scaling) considered for HTE, and if so, what impact did\nthey have on the observed trade-o! between calibration and accuracy? How about out-of-distribution\nperformance?\nResponseR2.0: Thanksforthiscarefulobservation. InresponsetoR2.0,weexpandedourevaluation\nby seeded repetitions and now report aggregated results. The aggregated results do not show the clear\ntrade-o!betweenaccuracyandcalibrationanymore. Wedidruntemperaturescalingforallexperiments,\nbut found the results to change marginally keeping the non-scaled predictions. Additionally, we extend\nour results on classification tasks by an out-of-distribution evaluation in figure 4 and text changes\nwere applied 501-523. In conclusion, the ResNet-HTE configuration still performed worse in terms\nof calibration error compared to the ResNet-SGHMC configurations, even when analyzing the results\nMatthias Schmal\n4\nmatthias.schmal@tu-ilmenau.de\nTechnische Universität Ilmenau\nData intensive Systems and Visualization Group (dAI.SY)\nMax-Planck-Ring 14\n98684 Ilmenau\nwith decreasing step-sizes. We discuss these observations in lines 586-593 and refer to the log-likelihood\ngraphshownintheappendixvisualizingthedi!erencesinthetwomodels. Duetotheextremelyvariable\nlog-likelihood in the ResNet example, it is very likely to end up in di!erent minima rather than carefully\nexploring those.\nTask and Architecture Generality: The empirical results presented in the paper focus exclusively\non image classification tasks using convolutional neural networks (CNNs), specifically LeNet on EM-\nNIST and a small ResNet variant on CIFAR-10. While these benchmarks are standard, the restriction to\nvision-based classification is a limiting factor, especially given the paper’s strong emphasis on empirical\nvalidation. The architectural choices exclude evaluation on non-convolutional models, such as trans-\nformers, and no experiments are reported on regression or structured output tasks, where calibration\nand uncertainty quantification play an equally critical role.\nGiventhatstochasticgradient-basedMCMCmethodsarebroadlyapplicable, itwouldbeinformative\nto test GSGHMC and HTE in alternative settings—such as time series forecasting, tabular regression,\nor even natural language tasks—to assess their robustness to domain shift, heteroscedastic noise, and\nvarying data modalities. Were any such evaluations attempted or preliminary results obtained in non-\nimage domains?\nResponse R2.1: Thanks for this valuable suggestion. We added a time-series prediction problem to\nthe final evaluation, which is a part of table 4 and is introduced in lines 454-483 and discussed in lines\n484-491 and 594-597. While we could not run all hyperparameter evaluation steps with this problem\ntoo, we show the general applicability in a di!erent domain.\nAcceptance Rate Dynamics and Interpretation: The paper reports an average acceptance rate of\napproximately50% 3.5%forboththeGSGHMCandHTEmethodsduringexperiments(SectionII-C).\n±\nWhile this figure reflects a moderate trade-o! between exploration and rejection, it looks surprisingly\nlow for methods based on SG-MCMC. A detailed temporal analysis of how the acceptance rate evolves\nduringburn-inandsamplingwouldbevaluable—e.g.,doestheratestabilize,decline,orexhibitoscillatory\nbehavior? Furthermore, the relationship between acceptance rate and predictive metrics (e.g., accuracy,\nNLL, calibration) is not explicitly addressed.\nIs there evidence that lower acceptance rates correlate with more informative exploration of wide\nminima, or do they instead reflect miscalibrated proposals due to aggressive step sizes? Were di!erent\nacceptance thresholds tested, orwere diagnostics (e.g., autocorrelation, e!ective sample size) computed\nto assess mixing quality as a function of acceptance behavior over time?\nResponse R2.2: Thanks for pointing out these questions. We show the acceptance rate dynamics\nwith uncertainties in figure C1 as well as autocorrelation for the test cases in figure C2. The text\nof Appendix C discusses their relevance in detail. We address this issue in the main document from\nlines 401-414. Our main observation regarding the acceptance rates was, that it did not change much,\nwhile the performance metrics did change a lot. Specifically for the step-size, practical experience\nsuggest to descrease the step size for increasing the acceptance rate. However, this procedure was\nnot applicable for large-scale deep learning models increasing the burn-in phase enormously or even\npreventing convergence to compatible minima. Therefore, we did not rely on acceptance rates for the\ntuning. We were optimizing performance while making sure to converge in reasonable time not getting\nstuck at local minima. During our hyper-parameter study, the acceptance rate was only minimally\ninfluenced.\nMatthias Schmal\n5\nmatthias.schmal@tu-ilmenau.de\nTechnische Universität Ilmenau\nData intensive Systems and Visualization Group (dAI.SY)\nMax-Planck-Ring 14\n98684 Ilmenau\nOn the Role of Subsampling and Ensemble Compression: The greedy subsampling procedure\nproposed in Section II-D, while practical, appears conceptually orthogonal to the core contributions of\nthe paper. The main message of the work centers on improving the fidelity and e"}, {"round": 2, "reviewer_comments": {"reviewer_1": "- Although Table 4 provides the number of epochs, it is difficult to compare the methods and understand when each should be used. For example, in the LeNet experiment, VOGN takes 500 epochs to achieve 16.2% MCE, while GSGHMC takes 5000 epochs to achieve 14.4% MCE. It would be clearer to present the performance under the same compute budget or the budget required to achieve the same performance.\n- Regarding the bias, I suggest providing a more explicit characterization. For instance, the distance between the asymptotic distribution and the target distribution, as well as the dependency of the bound on the critical hyperparameters, would make the analysis more transparent.\n- Similarly, a more explicit characterization of the convergence rate would be helpful. For example, the distance between the estimated distribution at step t and the target distribution. It is unclear how Section 4.4 and Equation 10 provide convergence guarantees, and further clarification would strengthen the argument.\n- My concern regarding the positioning of this paper remains unresolved. It is unclear how the proposed methods differ from existing ones, both in theory and in experiments. The authors suggest that the paper's goal is to provide an empirical comparison in deep learning. If this is the case, then the paper overlooks a substantial body of literature in the area of uncertainty estimates in deep learning, which should be discussed and compared, such as [1,2,3,4].\nSummary:\nMy concerns regarding error bars and comparisons to mini-batch-based Metropolis-Hastings have been addressed.\nHowever, other issues remain.\n[1] Laplace Redux -- Effortless Bayesian Deep Learning, 2021\n[2] A Simple Baseline for Bayesian Uncertainty in Deep Learning, 2019\n[3] Cyclical Stochastic Gradient MCMC for Bayesian Deep Learning, 2020\n[4] Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning, 2016", "reviewer_3": "The paper proposes two methods for approximate Bayesian inference in in large machine learning models.\nThe proposed methodology seems sound, and the comparisons to contending methods are in my mind convincing.\nMy main objection with the paper is mainly the structure, though I understand that this is related to the journal guidelines. In particular the results before methodology is confusing as the reader lacks knowledge of what is really compared. Further, what I assume also are journal constraints limits the exposition of the proposed methods, and in particular gives the reader very few details about how the methods really functon(which would be helpful for someone implementing the methods). I.e. how does changes to tuning parameters influence the behavior of the chains, where are the potential for further work further improving the methods and so on. Going directly to large scale comparisons hides such fine detail.\nMinor details:\nLine 130: missing reference"}, "author_response": "Point-by-point Response\nReviewer Nr. 1\n• Although Table 4 provides the number of epochs, it is di!cult to compare the methods and\nunderstand when each should be used. For example, in the LeNet experiment, VOGN takes\n500 epochs to achieve 16.2% MCE, while GSGHMC takes 5000 epochs to achieve 14.4% MCE.\nIt would be clearer to present the performance under the same compute budget or the budget\nrequired to achieve the same performance.\nResponse: Thank you for this feedback. Sampling algorithms are naturally more costly than\nvariational inference or even maximum a posteriori optimization. However, they come with the\nadvantageoflessmodelassumptionsandthepossibilitytoapproximatemulti-modaldistributions.\nWe conducted an additional study reported as appendix G of the manuscript, which among others\n(see below) analyzes convergence properties of the proposed methods. We refer to the study\nin lines 624-627 of the manuscript as follows: “Our convergence analysis on a small benchmark\nproblem showed rapid convergence for medium and large step sizes and multi modal distributions\n(cp. Appendix G). In this analysis, the proposed GSGHMC and HTE methods show improved\nmixingovertheircompetitors.” Althoughthestudyonlycomparessamplingalgorithms,theresults\nprovide more insights into the additional costs caused by sampling. We agree, that the cost is\ntherefore hard to compare, but this remains a problem of Bayesian inference with sampling in\ngeneral.\n• Regarding the bias, Isuggest providing amore explicitcharacterization. Forinstance, the distance\nbetween the asymptotic distribution and the target distribution, as well as the dependency of the\nbound on the critical hyperparameters, would make the analysis more transparent.\n• Similarly, a more explicit characterization of the convergence rate would be helpful. For example,\nthe distance between the estimated distribution at step t and the target distribution. It is unclear\nhow Section 4.4 and Equation 10 provide convergence guarantees, and further clarification would\nstrengthen the argument.\nResponse: Thanks for pointing this out. We added a further experiment to demonstrate this\nbehavior in the new appendix G, covering the concerns from both points above. Our results allow\nfor comparing the convergence in kernel stein discrepancy as a measure of distributional distance.\nFurthermore, it highlights the e\"ect of di\"erent step sizes on this convergence.\n• My concern regarding the positioning of this paper remains unresolved. It is unclear how the\nproposed methods di\"er from existing ones, both in theory and in experiments. The authors\nsuggest that the paper’s goal is to provide an empirical comparison in deep learning. If this is the\ncase, then the paper overlooks a substantial body of literature in the area of uncertainty estimates\nin deep learning, which should be discussed and compared, such as [1, 2, 3, 4].\nResponse: We agree that the suggested methods are an valuable addition to our study and now\nadditionally report results on [2, 3, 4] in Table 4. The Laplace approximation of [1] is structurally\nsimilar to the VOGN results we already report and VOGN has been demonstrated to outperform\nLaplace in [5]. Therefore, we decided to not add [1] and keep the results of the VOGN method.\nWe extended the description in lines 458 to 460 accordingly.\nSummary: My concerns regarding error bars and comparisons to mini-batch-based Metropolis-\nHastings have been addressed. However, other issues remain.\nmatthias.schmal@tu-ilmenau.de\nData intensive Systems and Visualization Group (dAI.SY)\nMax-Planck-Ring 14\n98684 Ilmenau\nReviewer Nr. 3\nThe paper proposes two methods for approximate Bayesian inference in large machine learning models.\nThe proposed methodology seems sound, and the comparisons to contending methods are in my mind\nconvincing.\nMy main objection with the paper is mainly the structure, though I understand that this is related\nto the journal guidelines. In particular the results before methodology is confusing as the reader lacks\nknowledge of what is really compared. Further, what I assume also are journal constraints limits the\nexposition of the proposed methods, and in particular gives the reader very few details about how the\nmethods really functon (which would be helpful for someone implementing the methods). I.e. how does\nchanges to tuning parameters influence the behavior of the chains, where are the potential for further\nwork further improving the methods and so on. Going directly to large scale comparisons hides such\nfine detail.\nResponse: Thanks for this feedback. We modified the introductural paragraph of the results section\nhinting towards the corresponding sections to account for this feedback (lines 170–174)\nMinor details: Line 130: missing reference\nResponse: We fixed the reference in line 130.\nReferences\n[1] Daxberger,ErikA.,AgustinusKristiadi,AlexanderImmer,RunaEschenhagen,M.BauerandPhilipp\nHennig.“LaplaceRedux-E\"ortlessBayesianDeepLearning.”NeuralInformationProcessingSystems\n(2021)\n[2] Maddox, Wesley J., T. Garipov, Pavel Izmailov, Dmitry P. Vetrov and Andrew Gordon Wilson. “A\nSimple Baseline for Bayesian Uncertainty in Deep Learning.” Neural Information Processing Systems\n(2019).\n[3] Zhang, Ruqi, Chunyuan Li, Jianyi Zhang, Changyou Chen and Andrew Gordon Wilson. “Cyclical\nStochastic Gradient MCMC for Bayesian Deep Learning.” International Conference on Learning\nRepresentations, (2020)\n[4] Gal, Yarin and Zoubin Ghahramani. “Dropout as a Bayesian Approximation: Representing Model\nUncertainty in Deep Learning.” International Conference on Machine Learning (2015).\n[5] Osawa, Swaroop et al. “Practical deep learining with Bayesian principles.” Conference on Neural\nInformation Processing Systems (2019)\nmatthias.schmal@tu-ilmenau.de\n3"}, {"round": 3, "reviewer_comments": {"reviewer_3": "I'm happy with the revision. Don't think there are further things to do given the constraints of the journal."}, "author_response": ""}]

Decision letter