{'fraction_non_alphanumeric': 0.03802797772647019, 'fraction_numerical': 0.012178007152881525, 'mean_word_length': 4.839280993920169, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gradient Regularization as Approximate Variational Inference

We developed Variational Laplace for Bayesian neural networks (BNNs), which exploits a local approximation of the curvature of the likelihood to estimate the ELBO without the need for stochastic sampling of the neural-network weights. The Variational Laplace objective is simple to evaluate, as it is the log-likelihood plus weight-decay, plus a squared-gradient regularizer. Variational Laplace gave better test performance and expected calibration errors than maximum a posteriori inference and standard sampling-based variational inference, despite using the same variational approximate posterior. Finally, we emphasize the care needed in benchmarking standard VI, as there is a risk of stopping before the variance parameters have converged. We show that early-stopping can be avoided by increasing the learning rate for the variance parameters.


Introduction
Neural networks are increasingly being used in safety-critical settings, such as selfdriving cars [1] and medical diagnosis [2]. In these settings, it is critical to be able to reason about uncertainty in the parameters of the network-for instance, so that the system is able to call for additional human input when necessary [3]. Several approaches to Bayesian inference in neural networks are available, including stochastic gradient Langevin dynamics [4] Laplace's method [5][6][7] and variational inference [8,9].
Here, we focus on combining the advantages of Laplace's method [5][6][7] and variational inference (VI; [10]). In particular, Laplace's method is very fast, as it begins by finding a mode using a standard gradient descent procedure, and then computes a local Gaussian approximate of the mode by performing a second-order Taylor expansion. However, as the mode is discovered by standard gradient descent, it may be a narrow mode that generalizes poorly [11]. In contrast, variational inference (VI; [8]) is slower, as it requires stochastic sampling of the weights, but that stochastic sampling forces it to find a broad, flat mode that presumably generalizes better.
We developed a new Variational Laplace (VL) method that combines the best of both worlds, giving a method that finds broad, flat modes even in the absence of the stochastic sampling. The resulting objective is composed of the log-likelihood, standard weight-decay regularization and a squared-gradient regularizer, which is weighted by the variance of the approximate posterior. VL displayed improved performance over VI and MAP in standard benchmark tasks.
Our squared gradient regularizer relates strongly to the effectiveness of stochastic gradient descent. In particular, recent work has shown that gradient descent implicitly uses a squared gradient regularizer [12], and that full-batch gradient descent with a squared gradient regularizer can recover much of the benefits of implicit regularization from minibatched stochastic gradient descent [13]. Our work implies that these regularizers can be interpreted as a form of approximate inference over the neural-network weights.

Variational Inference (VI) for Bayesian Neural Networks
To perform variational inference for neural networks, we follow the usual approach [8,14], by using independent Gaussian priors, P and approximate posteriors Q for all parameters, w: where µ λ and σ 2 λ are learned parameters of the approximate posterior, and where Σ is a diagonal matrix, with Σ λλ = σ 2 λ . We fit the approximate posterior by optimizing the evidence lower bound objective (ELBO) with respect to parameters of the variational posterior, µ λ and σ 2 λ : Here, x is all training inputs, y is all training outputs, and β is the tempering parameter which is 1 for a close approximation to Bayesian inference, but is often set to smaller values to "temper" the posterior, which often improves empirical performance [15,16] and has theoretical justification as accounting for the data-curation process [17]. We need to optimize the expectation in Equation (3) with respect to the parameters of Q(w), the distribution over which the expectation is taken. To perform this optimization efficiently, the usual approach is to use the reparameterization trick [8,18,19]-we write w in terms of : where λ ∼ N (0, 1). Thus, the ELBO can be written as an expectation over : where the distribution over is now fixed. Critically, now the expected gradient of the term inside the expectation is equal to the gradient of L VI , so we can use samples of to estimate the expectation.

Laplace's Method
Laplace's method [5][6][7] first finds a mode by doing gradient ascent on the log-joint: and uses a Gaussian approximate posterior around that mode, where H(w * ) is Hessian of the log-joint at w * .

Related Work
There is past work on Variational Laplace [20][21][22], which learns the mean parameters, µ, of a Gaussian approximate posterior, Q µ (w) = N w; µ, −H −1 (µ) (8) and obtains the covariance matrix as a function of the mean parameters using the Hessian, as in Laplace's method. However, instead of taking the approximation to be centered around a MAP solution, w * , they take the approximate posterior to be centered on learned mean parameters, µ. Importantly, they simplify the ELBO by substituting this approximate posterior into Equation (3), and approximating the log-joint using its Taylor series expansion. Ultimately they obtain L VI ≈ log P(y|w=µ, x) + log P(w=µ) − 1 2 log|H(µ)| + const.
However, there are two problems with this approach when applied to neural networks. First, the algebraic manipulations required to derive Equation (9) require the full N × N Hessian, H(µ), for all N parameters, and neural networks have too many parameters for this to be feasible. Second, the log|H(µ)| term in Equation (9) cannot be minibatched, as we need the full sum over minibatches inside the log to compute the Hessian: where H j (µ) is the contribution to the Hessian from an individual minibatch. Due to these issues, past Variational Laplace methods did not scale to large neural networks. An alternative deterministic approach to variational inference in Bayesian neural networks, approximates the distribution over activities induced by stochasticity in the weights [23]. Unfortunately, it is important to capture the covariance over features induced by stochasticity in the weights. In fully connected networks, this is feasible, as we usually have a small number of features at each layer. However, in convolutional networks, we have a large number of features, channels × height × width. In the lower layers of a ResNet, we may have 64 channels and a 32 × 32 feature map, resulting in 64 × 32 2 = 65,536 features and a 65,536 × 65,536 covariance matrix. These scalability issues prevented them from applying their approach to convolutional networks. In contrast, our approach is highly scalable and readily applicable to the convolutional setting. Subsequent work such as Haußmann et al. [24] introduced other deterministic approximations, based on decomposing the relu into a linear function and a Heaviside step. However, their approach had errors of ∼30% on CIFAR-10.
Ritter et al. [7] and MacKay [25] used Laplace's method in Bayesian neural networks, by first finding the mode by doing gradient ascent on the log-joint probability, and expanding around that mode. As usual for Laplace's method, they risk finding a narrow mode that generalizes poorly. In contrast, we find a mode using an approximation to the ELBO that takes the curvature into account and hence is biased towards broad, flat modes that presumably generalise better.
Our approach gives a squared-gradient regularizer that is similar to those discovered in past work [12,26]. They showed that squared-gradient regularizers connect to gradient descent, in that approximation errors due to finite-step sizes in gradient-descent imply an effective squared gradient regularization. The similarity of our objectives raises profound questions about the extent to which gradient descent can be said to perform Bayesian inference. That said there are three key differences. First, their approach connects fullbatch gradient descent to squared-gradient regularizers. Of course, most neural network training is stochastic gradient descent based on minibatches. Given that the stationary distribution of SGD is isotropic Gaussian (in a quadratic loss-function; Appendix A), we are able to connect stochastic gradient descent to squared gradient regularizers, and hence to approximate variational inference. This is especially important in light of recent work showing full-batch gradient descent gives poor performance, but that performance can be improved by including explicit squared gradient regularization. Our work indicates that explicit squared gradient regularization is mimicking the implicit regularization from stochastic gradient descent. First, our method uses the Fisher, (i.e., the gradients for data sampled from the model) whereas their approach uses the empirical Fisher, (i.e., gradients for the observed data) to form the squared gradient regularizer [27]. Second, our approach gives a principled method to learn a separate weighting for the squared-gradient for each parameter, whereas the connection to SGD forces Barrett and Dherin [12] to use a uniform weighting across all parameters.
Our work differs from, e.g., Khan et al. [28] by explicitly providing a simply implemented loss-function in terms of a squared-gradient regularizer, instead of working with NTK-inspired approximations to the Hessian.
Other approaches include "Broad Bayesian Learning" [29], which optimizes the architecture of a Bayesian neural network, exploiting information from previously trained but different networks. Of course, quantification of uncertainty for Bayesian neural networks is always fraught [30]. As such, we followed standard practice in the literature of reporting OOD detection performance and a measure of calibration accuracy [31].

Methods
To combine the best of VI and Laplace's method, we begin by noting that the ELBO can be rewritten in terms of the KL divergence between the prior and approximate posterior: where the KL-divergence can be evaluated analytically: As such, the only term we need to approximate is the expected log-likelihood.
To approximate the expectation, we begin by taking a second-order Taylor series expansion of the log-likelihood around the current settings of the mean parameters, µ: (13) where B is the number of minibatches, g j is the gradient for minibatch j and H is the Hessian for the full dataset: Here, x and y are the the inputs and outputs for the full dataset, whereas x j and y j are the inputs and outputs for minibatch j. Now we consider the expectation of each of these terms under the approximate posterior, Q(w). The first term is constant and independent of w.
The second (linear) term is zero, because the expectation of (w − µ) under the approximate posterior is zero The third (quadratic) term might at first appear difficult to evaluate because it involves H, the N × N matrix of second derivatives, where N is the number of parameters in the model. However, using properties of the trace, and noting that the expectation of (w − µ)(w − µ) T is the covariance of the approximate posterior, we obtain writing the trace in index notation, and substituting for the (diagonal) posterior covariance, Σ: Thus, our first approximation of the expected log-likelihood is and substituting this into Equation (11) gives . (20) This resolves most of the issues with the original Variational Laplace method: it requires only the diagonal of the Hessian, it can be minibatched and it does not blow up if H λλ is zero.

Pathological Optima When Using the Hessian
However, a new issue arises: H λλ is usually negative, in which case the approximation in Equation (20) can be expected to work well. However there is nothing to stop H λλ from becoming positive. Usually if we, e.g., took the log-determinant of the negative Hessian, this would immediately break the optimization process (as we would be taking the logarithm of a negative number). However, in our context, there is no immediate issue as Equation (20) takes on a well-defined value even when one or more H λλ 's are positive. That said, we rapidly encounter similar issues as we get pathological optimal values of σ 2 λ . In particular, picking out the terms in the objective that depend on σ 2 λ , absorbing the other terms into the constant, and taking β = 1 for simplicity, we have Thus, the gradient wrt a single variance parameter is In the typical case, H λλ is negative so 1 s 2 λ − H λλ is positive, and we can find the optimum by solving for the value of σ 2 λ where the gradient is zero:  (23) undefined, but the gradient is always positive:

Avoiding Pathologies with the Fisher
To avoid pathologies arising from the fact that the Hessian is not necessarily negative definite, a common approach is to approximate the Hessian using the Fisher information matrix: Importantly,g is the gradient of the log-likelihood for data sampled from the model,ỹ j , not for the true data:g This gives us the Fisher, which is a commonly used and well-understood approximation to the Hessian [27]. Importantly, this contrasts with the empirical Fisher [27], which uses the gradient conditioned on the actual data (and not data sampled from the model): which is problematic, because there is a large rank-1 component in the direction of the mean gradient, which disrupts the estimated matrix specifically in the direction of interest for problems such as optimization [27].
Using the Fisher information (Equation (25)) in Equation (19), we obtain an approximate expected log-likelihood: Substituting this into Equation (11) gives us the final VL objective, L VL , which is an approximation of the ELBO: In practice, we typically take the objective for a minibatch, divided by the number of datapoints in a minibatch, S: where 1 Sg j;λ are the gradients of the log-likelihood for the minibatch averaged across datapoints, i.e., the gradient of 1 S log P ỹ j |x j , w=µ . Remember B is the number of minibatches so SB is the total number of training datapoints.

Constraints on the Network Architecture
Importantly, here the regularizer is the squared gradient of the loss with respect to the parameters. As such, computing the loss implicitly involves a second-derivative of the log-likelihood, and we therefore cannot use piecewise linear activation functions such as ReLU, which have pathological second derivatives. In particular, the second derivative has a delta-function "spike" at zero: where φ is the relu nonlinearity, Θ(x) is the Heaviside step function which is zero for x < 0 and one for 0 < x, and δ(x) is the Dirac delta function. As the function is almost never evaluated at exactly zero, it is not possible to sensibly take into account the contribution of the infinitely high spike in the second derivative at zero. Interestingly, this issue is very similar to the one that turns up when differentiating step (i.e., Θ(x)) activationsthe derivative is well-defined and zero almost everywhere. The issue is that there are delta-function spikes in the gradient at zero that gradient descent cannot reasonably work with. Instead, we used a softplus activation function, but any activation with well-behaved second derivatives is admissible.

Results
We compared MAP, VI and our method (VL) on four different datasets (CIFAR-10, CIFAR-100 [32], SVHN [33] and fashion-MNIST [34] MIT Licensed) using a PreactResNets-18 [35] with an initial learning rate of 1E-4, which decreased by a factor of 10 after 100 and 150 epochs and a batch size of 128 with all the other optimizer hyperparameters set to their default values. We tried two variants of variational inference: evaluating test-performance using the mean network, VI (mean), and evaluating test performance by drawing 10 samples from the approximate posterior, VI (sampled). We swept across different degrees of posterior tempering, β. Using β < 1 is normatively justified in the Bayesian framework as accounting for the effect of data curation [17]. For many values of β VL gave better test accuracies, test log-likelihoods and expected calibration errors [31,36] than VI or MAP inference (Figure 1). Importantly though, for the optimal value of β, VL almost always gave better performance on these metrics (Table 1). These experiments took ∼480 GPU hours, and were run on a mixture of nVidia 1080 and 2080 GPUs in an internal cluster.   [31,36] and OOD detection metric (AUROC) for CIFAR-10, CIFAR-100, SVHN and fashion MNIST. Downsampled Imagenet [37] was used as OOD data. See Appendix B.
The runtimes of the methods are listed in Table 2. VL or gradient regularization was around a factor of three slower than either VI or MAP due to the need to compute secondderivatives. It is still eminently feasible, especially in comparison to past methods for deterministic variational inference that have fundamental difficulties in scaling to convolutional networks [23]. Furthermore, we did not find that increasing the number of epochs improved performance either for VI or MAP, as we were already training for convergence.

Early-Stopping and Poor Performance in VI
Before performing comparisons where we learn the approximate posterior variance, it is important to understand the pitfalls when optimizing variational Bayesian neural networks using adaptive optimizers such as Adam. In particular, there is a strong danger of stopping the optimization before the variances have converged. To illustrate this risk, note that Adam [38] updates take the form where η is the learning rate, m is an unbiased estimator of the mean gradient, g , v is an unbiased estimator of the squared gradient, g 2 and is a small positive constant to avoid divide-by-zero. The magnitude of the updates, |∆θ|, is maximized by having exactly the same gradient on each step, in which case, neglecting , we have |∆θ| = η. As such, with a learning rate of η = 10 −4 , a training set of 50,000 and a batch size of 128 parameters can move at most 50,000/128 ×10 −4 ≈ 0.04 per epoch. Doing 100 epochs at this learning rate, a parameter can change by at most 4 over the 100 epochs before the first learning rate step. This is fine for the weights, which typically have very small values. However, the underlying parameters used for the variances typically take on larger values. In our case, we will use log σ λ as the parameter, and initialize it to three less than the prior standard deviation, log s λ − 3. To ensure reasonable convergence, log σ λ should be able to revert back to the prior, implying that it must be able to change by at least three during the course of training. Unfortunately, 3 is very close to the maximum possible change of 4, raising the possibility that the variance parameters will not actually converge. To check whether early-stopping was indeed an issue, we plotted the (tempered) ELBO for VI ( Figure 2A) and VL ( Figure 2B). For VI (Figure 2A) with the standard setup (lightest line with a learning rate multiplier of 1), the ELBO clearly has not converged at 100 epochs, indicating earlystopping. Notably, this was still an issue with VL ( Figure 2B), especially if we were to train for fewer epochs. However, the effect was smaller for VL, which may have been because the gradients were more consistent, as it did not sample the weights. These issues can be rectified by increasing the learning rate specifically for the log σ λ parameters (darker lines).
We then plotted the test log-likelihood ( Figure 2C), test accuracy ( Figure 2D) and ELBO ( Figure 2E) against the learning rate multiplier. Again, the performance for VL (orange) was reasonably robust to changes in the learning rate multiplier. However, the performance of VI (blue) was very sensitive to the multiplier: as the multiplier increased, test performance fell but the ELBO rose. As we ultimately care about test performance, these results would suggest that we should use the lowest multiplier (1), and accept the possibility of early-stopping. That may be a perfectly good choice in many cases. However, VI is supposed to be an approximate Bayesian method, and using an alternative form for the ELBO, L VI = log P(y|x) − D KL (Q(w)|| P(w|y, x)), (33) we can see that the ELBO measures KL-divergence between the true and approximate posterior, and hence the quality of our approximate Bayesian inference. As such, very poor ELBOs imply that the KL-divergence between the true and approximate posterior is very large, and hence the "approximate posterior" is no longer actually approximating the true posterior. As such, if we are to retain a Bayesian interpretation of VI, we need to use larger learning rate multipliers which give better values for the ELBO ( Figure 2E). However, in doing that, we get worse test performance ( Figure 2C,D). This conflict between approximate posterior quality and test performance is very problematic: the Bayesian framework would suggest that as Bayesian inference becomes more accurate, performance should improve, whereas for VI, performance gets worse. Concretely, by initializing log σ λ to a small value and then early-stopping, we leave log σ λ at a small value through training, in which case VI becomes equivalent to MAP inference with a negligibly small amount of noise added to the weights. We would therefore expect early-stopped VI to behave (and be) very similar to MAP inference.
In subsequent experiments, we chose to use a learning rate multiplier of 10, as this largely eliminated early-stopping (though see VI with β = 0.1; Figure 2E). Overall, this indicates that we have to be very careful to avoid early stopping when running standard, sampling-based variational inference.

Conclusions
We gave a novel Variational Laplace approach to inference in Bayesian neural networks which combines the best of previous approaches based on variational inference and Laplace's Method. This method gave excellent empirical performance compared to VI.