ROBUST BAYESIAN NEURAL NETWORKS
Year: 2019

Abstract:
Bayesian neural networks (BNNs) hold great promise as a flexible and principled solution to deal with uncertainty when learning from finite data. Among approaches to realize probabilistic inference in deep neural networks, variational Bayes (VB) is theoretically grounded, generally applicable, and computationally efficient. With wide recognition of potential advantages, why is it that variational Bayes has seen very limited practical use for BNNs in real applications? We argue that variational inference in neural networks is fragile: successful implementations require careful initialization and tuning of prior variances, as well as controlling the variance of Monte Carlo gradient estimates. We provide two innovations that aim to turn VB into a robust inference tool for Bayesian neural networks: first, we introduce a novel deterministic method to approximate moments in neural networks, eliminating gradient variance; second, we introduce a hierarchical prior for parameters and a novel Empirical Bayes procedure for automatically selecting prior variances. Combining these two innovations, the resulting method is highly efficient and robust. On the application of heteroscedastic regression we demonstrate good predictive performance over alternative approaches.

heading: 1 INTRODUCTION
text: Bayesian approaches to neural network training marry the representational flexibility of deep neural networks with principled parameter estimation in probabilistic models. Compared to “standard” parameter estimation by maximum likelihood, the Bayesian framework promises to bring key advantages such as better uncertainty estimates on predictions and automatic model regularization (MacKay, 1992; Graves, 2011). These features are often crucial for informing downstream decision tasks and reducing overfitting, particularly on small datasets. However, despite potential advantages, such Bayesian neural networks (BNNs) are often overlooked due to two limitations: First, posterior inference in deep neural networks is analytically intractable and approximate inference with Monte Carlo (MC) techniques can suffer from crippling variance given only a reasonable computation budget (Kingma et al., 2015; Molchanov et al., 2017; Miller et al., 2017; Zhu et al., 2018). Second, performance of the Bayesian approach is sensitive to the choice of prior (Neal, 1993), and although we may have a priori knowledge concerning the function represented by a neural network, it is generally difficult to translate this into a meaningful prior on neural network weights. Sensitivity to priors and initialization makes BNNs non-robust and thus often irrelevant in practice.
In this paper, we describe a novel approach for inference in feed-forward BNNs that is simple to implement and aims to solve these two limitations. We adopt the paradigm of variational Bayes (VB) for BNNs (Hinton & van Camp, 1993; MacKay, 1995c) which is normally deployed using Monte ∗Work done during an internship at Microsoft Research, Cambridge. †Work done while at Microsoft Research, Cambridge.
Carlo variational inference (MCVI) (Graves, 2011; Blundell et al., 2015). Within this paradigm we address the two shortcomings of current practice outlined above: First, we address the issue of high variance in MCVI, by reducing this variance to zero through novel deterministic approximations to variational inference in neural networks. Second, we derive a general and robust Empirical Bayes (EB) approach to prior choice using hierarchical priors. By exploiting conjugacy we derive data-adaptive closed-form variance priors for neural network weights, which we experimentally demonstrate to be remarkably effective.
Combining these two novel ingredients gives us a performant and robust BNN inference scheme that we refer to as “deterministic variational inference” (DVI). We demonstrate robustness and improved predictive performance in the context of non-linear regression models, deriving novel closed-form results for expected log-likelihoods in homoscedastic and heteroscedastic regression (similar derivations for classification can be found in the appendix).
Experiments on standard regression datasets from the UCI repository, (Dheeru & Karra Taniskidou, 2017), show that for identical models DVI converges to local optima with better predictive loglikelihoods than existing methods based on MCVI. In direct comparisons, we show that our Empirical Bayes formulation automatically provides better or comparable test performance than manual tuning of the prior and that heteroscedastic models consistently outperform the homoscedastic models.
Concretely, our contributions are:
• Development of a deterministic procedure for propagating uncertain activations through neural networks with uncertain weights and ReLU or Heaviside activation functions. • Development of an EB method for principled tuning of weight priors during BNN training. • Experimental results showing the accuracy and efficiency of our method and applicability to
heteroscedastic and homoscedastic regression on real datasets.

heading: 2 VARIATIONAL INFERENCE IN BAYESIAN NEURAL NETWORKS
text: We start by describing the inference task that our method must solve to successfully train a BNN. Given a model M parameterized by weights w and a dataset D = (x,y), the inference task is to discover the posterior distribution p(w|x,y). A variational approach acknowledges that this posterior generally does not have an analytic form, and introduces a variational distribution q(w;θ) parameterized by θ to approximate p(w|x,y). The approximation is considered optimal within the variational family for θ∗ that minimizes the Kullback-Leibler (KL) divergence between q and the true posterior.
θ∗ = argmin θ DKL [q(w;θ)||p(w|x,y)].
Introducing a prior p(w) and applying Bayes rule allows us to rewrite this as optimization of the quantity known as the evidence lower bound (ELBO):
θ∗ = argmax θ {Ew∼q [log p(y|w,x)]−DKL [q(w;θ)||p(w)]} . (1)
Analytic results exist for the KL term in the ELBO for careful choice of prior and variational distributions (e.g. Gaussian families). However, whenM is a non-linear neural network, the first term in equation 1 (referred to as the reconstruction term) cannot be computed exactly: this is where MC approximations with finite sample size S are typically employed:
Ew∼q [log p(y|w,x)] ≈ 1
S S∑ s=1 log p(y|w(s),x), w(s) ∼ q(w;θ). (2)
Our goal in the next section is to develop an explicit and accurate approximation for this expectation, which provides a deterministic, closed-form expectation calculation, stabilizing BNN training by removing all stochasticity due to Monte Carlo sampling.

heading: 3 DETERMINISTIC VARIATIONAL APPROXIMATION
text: Figure 1 shows the architecture of the computation of Ew∼q [log p(D|w)] for a feed-forward neural network. The computation can be divided into two parts: first, propagation of activations though
ℒ
𝔼𝒘∼𝑞 log 𝑝 𝒚 𝒘, 𝒙
𝒂0 𝒉1 𝒂1
Layer 1 Layer 2 Layer 𝐿…
(a) (b)
𝒂𝐿𝒙
𝒚 𝑞 (𝑾
0 ;𝜽 )
𝑞 (𝑾
1 ;𝜽 )
𝑞 (𝒂
0 )
𝑞 (𝒂
1 ) 𝑞 (𝒂
𝐿 )
Figure 1: Architecture of a Bayesian neural network. Computation is divided into (a) propagation of activations (a) from an input x and (b) computation of a log-likelihood function L for outputs y. Weights are represented as high dimensional variational distributions (blue) that induce distributions over activations (yellow). MCVI computes using samples (dots); our method propagates a full distribution.
parameterized layers and second, evaluation of an unparameterized log-likelihood function (L). In this section, we describe how each of these stages is handled in our deterministic framework.

heading: 3.1 MOMENT PROPAGATION
text: We begin by considering activation propagation (figure 1(a)), with the aim of deriving the form of an approximation q̃(aL) to the final layer activation distribution q(aL) that will be passed to the likelihood computation. We compute aL by sequentially computing the distributions for the activations in the preceding layers. Concretely, we define the action of the lth layer that maps a(l−1) to al as follows:
hl = f(a(l−1)),
al = hlW l + bl ,
where f is a non-linearity and {W l, bl} ⊂ w are random variables representing the weights and biases of the lth layer that are assumed independent from weights in other layers. For notational clarity, in the following we will suppress the explicit layer index l, and use primed symbols to denote variables from the (l − 1)th layer, e.g. a′ = a(l−1). Note that we have made the non-conventional choice to draw the boundaries of the layers such that the linear transform is applied after the nonlinearity. This is to emphasize that al is constructed by linear combination of many distinct elements of h′, and in the limit of vanishing correlation between terms in this combination, we can appeal to the central limit theorem (CLT). Under the CLT, for a large enough hidden dimension and for variational distributions with finite first and second moments, elements ai will be normally distributed regardless of the potentially complicated distribution for hj induced by f 1. We empirically observe that this claim is approximately valid even when (weak) correlations appear between the elements of h during training (see section 3.1.1).
Having argued that a adopts a Gaussian form, it remains to compute the first and second moments. In general, these cannot be computed exactly, so we develop an approximate expression. An overview of this derivation is presented here with more details in appendix A. First, we modelW , b and h as independent random variables, allowing us to write:
〈ai〉 = 〈hj〉 〈Wji〉+ 〈bi〉 , Cov(ai, ak) = 〈hjhl〉Cov(Wji,Wlk) + 〈Wji〉Cov(hj , hl) 〈Wlk〉+ Cov(bi, bk), (3)
where we have employed the Einstein summation convention and used angle brackets to indicate expectation over q. If we choose a variational family with analytic forms for weight means and covariances (e.g. Gaussian with variational parameters 〈Wji〉 and Cov(Wji,Wlk)), then the only difficult terms are the moments of h:
〈hj〉 ∝ ∫ f(αj) exp [ − (αj−〈a ′ j〉)2
2Σ′jj
] dαj , (4)
〈hjhl〉 ∝ ∫ f(αj)f(αl) exp [ −1
2 ( αj−〈a′j〉 αl−〈a′l〉 )> ( Σ′jj Σ ′ jl Σ′lj Σ ′ ll )−1( αj−〈a′j〉 αl−〈a′l〉 )] dαj dαl , (5)
1We are also required to choose a Gaussian variational approximation for b to preserve the Gaussian distribution of a.
where we have used the Gaussian form of a′ parameterized by mean 〈a′〉 and covariance Σ′, and for brevity we have omitted the normalizing constants. Closed form solutions for the integral in equation 4 exist for Heaviside or ReLU choices of non-linearity f (see appendix A). Furthermore, for these non-linearities, the 〈 a′j 〉 → ±∞ and 〈a′l〉 → ±∞ asymptotes of the integral in equation 5 have closed form. Figure 2 shows schematically how these asymptotes can be used as a first approximation for equation 5. This approximation is improved by considering that (by definition) the residual decays to zero far from the origin in the ( 〈 a′j 〉 , 〈a′l〉) plane, and so is well modelled by a decaying function
exp[−Q( 〈 a′j 〉 , 〈a′l〉 ,Σ′)], where Q is a polynomial in 〈a′〉 with a dominant positive even term. In practice we truncate Q at the quadratic term, and calculate the polynomial coefficients by matching the moments of the resulting Gaussian with the analytic moments of the residual. Specifically, using dimensionless variables µ′i = 〈a′i〉 / √ Σ′ii and ρ ′ ij = Σ ′ ij/ √ Σ′iiΣ ′ jj , this improved approximation takes the form
〈hjhl〉 = S′jl { A(µ′j , µ ′ l, ρ ′ jl) + exp [ −Q(µ′j , µ′l, ρ′jl) ]} , (6)
where the expressions for the dimensionless asymptote A and quadratic Q are given in table table 1 and derived in appendix A.2.1 and A.2.2. The dimensionful scale factor S′jl is 1 for a Heaviside non-linearity or Σ ′ jl/ρ ′ jl for ReLU. Using equation 6 in equation 3 gives a closed form approximation for the moments of a as a function of moments of a′. Since a is approximately normally distributed by the CLT, this is sufficient information to sequentially propagate moments all the way through the network to compute the mean and covariances of q̃(aL), our explicit multivariate Gaussian approximation to q(aL). Any deep learning framework supporting special functions arcsin and Φ will immediately support backpropagation through the deterministic expressions we have presented. Below we briefly empirically verify the presented approximation, and in section 3.2 we will show how it is used to compute an approximate log-likelihood and posterior predictive distribution for regression and classification tasks.

heading: 3.1.1 EMPIRICAL VERIFICATION
text: Approximation accuracy The approximation derived above relies on three assumptions. First, that some form of CLT holds for the hidden units during training where the iid assumption of the classic CLT is not strictly enforced; second, that a quadratic truncation of Q is sufficient2; and third that there are only weak correlation between layers so that they can be represented using independent variables in the variational distribution. To provide evidence that these assumptions hold in practice, we train a small ReLU network with two hidden layers each of 128 units
2Additional Taylor expansion terms can be computed if this assumption fails.
𝑚
𝑥 𝒂0 𝒂1 𝒂𝐿=2
1 128 128 2
ℓ
2
1
−0.5
−1
𝑦
0.0 0.5
0
𝑥
Data samples Data 1-std Model 1-std
(a) (b)
0 10 20−10−20
𝑞(𝑚) 𝑞(ℓ)
𝑞(𝑚) 𝑞(ℓ)
MC ours
0 10−10
−0.6 −0.4−0.8−1.0 −2.5 −2.0−3.0−3.5 𝑚 ℓ
B ef
o re
T ra
in in g A ft er T ra in in g
(c)
Figure 3: Empirical accuracy of our approximation on toy 1-dimensional data. (a) We train a 2 layer ReLU network to perform heteroscedastic regression on the dataset shown in (b) and obtain the fit shown in blue. (c) The output distributions for the activation units m and ` evaluated at x = 0.25 are in excellent agreement with Monte Carlo (MC) integration with a large number (20k) of samples both before and after training.
to perform 1D heteroscedastic regression on a toy dataset of 500 points drawn from the distribution shown in figure 3(b). Deeper networks and skip connections are considered in appendix C. The training objective is taken from section 4, and the only detail required here is that aL is a 2-element vector where the elements are labelled as (m, `). We use a diagonal Gaussian variational family to represent the weights, but we preserve the full covariance of a during propagation. Using an input x = 0.25 (see arrow, Figure 3(b)) we compute the distributions form and ` both at the start of training (where we expect the iid assumption to hold) and at convergence (where iid does not necessarily hold). Figure 3(c) shows the comparison between aL distributions reported by our deterministic approximation and MC evaluation using 20k samples from q(w;θ). This comparison is qualitatively excellent for all cases considered.
Computational efficiency In traditional MCVI, propagation of S samples of d-dimensional activations through a layer containing a d× d-dimensional transformation requires O(Sd2) compute and O(Sd) memory. Our DVI method approximates the S →∞ limit, while only demanding O(d3) compute and O(d2) memory (the additional factor of d arises from manipulation of the quadratically large covariance matrix Cov[hj , hl]). Whereas MCVI can always trade compute and memory for accuracy by choosing a small value for S, the inherent scaling of DVI with d could potentially limit its practical use for networks with large hidden size. To avoid this limitation, we also consider the case where only the diagonal entries Cov(hj , hj) are computed and stored at each layer. We refer to this method as “diagonal-DVI” (dDVI), and in section 6 we show the surprising result that the strong test performance of DVI is largely retained by dDVI across a range of datasets. Figure 4 shows the time required to propagate activations through a single layer using the MCVI, DVI and dDVI methods
on a Tesla V100 GPU. As a rough rule of thumb (on this hardware), for layer sizes of practical relevance, we see that absolute DVI runtimes roughly equate to MCVI with S = 300 and dDVI runtime equates to S = 1.

heading: 3.2 LOG-LIKELIHOOD EVALUATION
text: To use the moment propagation procedure derived above for training BNNs, we need to build a function L that maps final layer activations aL to the expected log-likelihood term in equation 1 (see figure 1(b)). In appendix B.1 we show the intuitive result that this expected log-likelihood over q(w)
can be rewritten as an expectation over q̃(aL). Ew∼q [log p(y|x,w)] = EaL∼q(aL) [ log p(y|aL) ] . (7)
With this form we can derive closed forms for specific tasks; for brevity we focus on the regression case and refer the reader to appendices B.4 and B.5 for the classification case.
Regression Case For simplicity we consider scalar y and a Gaussian noise model parameterized by mean m(x;w) and heteroscedastic log-variance log σ2y(x) = `(x;w). The parameters of this Gaussian are read off as the elements of a 2-dimensional output layer aL = (m, `) so that p(y|aL) = N [ y|m, e` ] . Recall that these parameters themselves are uncertain and the statistics〈
aL 〉
and ΣL can be computed following section 3.1. Inserting the Gaussian forms for p(y|aL) and q(aL) into equation 7 and performing the integral (see appendix B.2) gives a closed form expression for the ELBO reconstruction term:
EaL∼q̃(aL) [ log p(y|aL) ] = − 12 [ log 2π + 〈`〉+ Σmm+(〈m〉−Σm`−y) 2 e〈`〉−Σ``/2 ] . (8)
This heteroscedastic model can be made homoscedastic by setting 〈`〉 = Σ`` = Σm` = 0. The expression in equation 8 completes the derivations required to implement the closed form approximation to the ELBO reconstruction term for training a network. In addition, we can also compute a closed form approximation to the predictive distribution that is used at test-time to produce predictions that incorporate all parameter uncertainties. By approximating the moments of the posterior predictive and assuming normality (see appendix B.3), we find:
p(y) ≈ ∫ p(y|aL) q̃(aL) daL ≈ N ( y ∣∣∣ 〈m〉 ,Σmm + e〈`〉+Σ``/2) . (9)

heading: 4 EMPIRICAL BAYES FOR VARIATIONAL BNNS
text: So far, we have described methods for deterministic approximation of the reconstruction term in the ELBO. We now turn to the KL term. For a d-dimensional Gaussian prior p(w) = N (µp,Σp), the KL divergence with the Gaussian variational distribution q = N (µq,Σq) has closed form:
DKL [q||p] = 12 [ log |Σp| |Σq| − d+ Tr ( Σ−1p Σq ) + (µp − µq)>Σ−1p (µp − µq) ] . (10)
However, this requires selection of (µp,Σp) for which there is usually little intuition beyond arguing µp = 0 by symmetry and choosing Σp to preserve the expected magnitude of the propagated activations (Glorot & Bengio, 2010; He et al., 2015). In practice, variational Bayes for neural network parameters is sensitive to the choice of prior variance parameters, and we will demonstrate this problem empirically in section 6 (figure 5).
To make variational Bayes robust we parameterize the prior hierarchically, retaining a conditional diagonal Gaussian prior and variational distribution on the weights. The hierarchical prior takes the form s ∼ p(s);w ∼ p(w|s), using an inverse gamma distribution on s as the conjugate prior to the elements of the diagonal Gaussian variance. We partition the weights into sets {λ} that typically coincide with the layer partitioning3, and assign a single element in s to each set:
sλ ∼ Inv-Gamma(α, β) , wλi ∼ N (0, sλ), (11)
for shape α and scale β, and where wλi is the i th weight in set λ.
Rather than taking the fully Bayesian approach, we adopt an empirical Bayes approach (Type-2 MAP), optimizing sλ, assuming that the integral is dominated by a contribution from this optimal value sλ = sλ∗ . We use the data to inform the optimal setting of s λ ∗ to produce the tightest ELBO:
ELBO = Ew∼q [ log p(y|hL(w)) ] − {DKL [q(w;θ)||p(w|s∗)p(s∗)]}
=⇒ sλ∗ = argmin sλ
{ DKL [ q(w;θ)||p(wλ|sλ) ] − log p(sλ) } (12)
3In general, any arbitrary partitioning can be used
Writing out the integral for the KL in equation 12, substituting in the forms of the distributions in equation 11 and differentiating to find the optimum gives
sλ∗ = Tr [ Σλq + µ λ q(µ λ q) >]+ 2β
Ωλ + 2α+ 2 , (13)
where Ωλ is the number of weights in the set λ. The influence of the data on the choice of sλ∗ is made explicit here through dependence on the learned variational parameters Σq and µq. Using sλ∗ to populate the elements of the diagonal prior variance Σp, we can evaluate the KL in equation 10 under the empirical Bayes prior. Optimization of the resulting ELBO then simultaneously tunes the variational distribution and prior.
In the experiments we will demonstrate that the proposed empirical Bayes approach works well; however, it only approximates the full Bayesian solution, and it could fail if we were to allow too many degrees of freedom. To see this, assume we were to use one prior per weight element, and we would also define a hyperprior for each prior mean. Then, adjusting both the prior variance and prior mean using empirical Bayes would always lead to a KL-divergence of zero and the ELBO objective would degenerate into maximum likelihood.

heading: 5 RELATED WORK
text: Bayesian neural networks have a rich history. In a 1992 landmark paper David MacKay demonstrated the many potential benefits of a Bayesian approach to neural network learning (MacKay, 1992); in particular, this work contained a convincing demonstration of naturally accounting for model flexibility in the form of the Bayesian Occam’s razor, facilitating comparison between different models, accurate calibration of predictive uncertainty, and to perform learning robust to overfitting. However, at the time Bayesian inference was achieved only for small and shallow neural networks using a comparatively crude Laplace approximation. Another early review article summarizing advantages and challenges in Bayesian neural network learning is (MacKay, 1995c).
This initial excitement around Bayesian neural networks led to two main methods being developed; First, Hinton & van Camp (1993) and MacKay (1995b) developed the variational Bayes (VB) approach for posterior inference. Whereas Hinton & van Camp (1993) were motivated from a minimum description length (MDL) compression perspective, MacKay (1995b) motivated his equivalent ensemble learning method from a statistical physics perspective of variational free energy minimization. Barber & Bishop (1998) extended the methodology for two-layer neural networks to use general multivariate Normal variational distributions. Second, Neal (1993) developed efficient gradient-based Monte Carlo methods in the form of “hybrid Monte Carlo”, now known as Hamiltonian Monte Carlo, and also raised the question of prior design and limiting behaviour of Bayesian neural networks.
Rebirth of Bayesian neural networks. After more than a decade of no further work on Bayesian neural networks Graves (2011) revived the field by using Monte Carlo variational inference (MCVI) to make VB practical and scalable, demonstrating gains in predictive performance on real world tasks.
Since 2015 the VB approach to Bayesian neural networks is mainstream (Blundell et al., 2015); key research drivers since then are the problems of high variance in MCVI and the search for useful variational families. One approach to reduce variance in feedforward networks is the local reparameterization trick (Kingma et al., 2015) (see appendix E). To enhance the variational families more complicated distributions such as Matrix Gaussian posteriors (Louizos & Welling, 2016), multiplicative posteriors (Kingma et al., 2015), and hierarchical posteriors (Louizos & Welling, 2017) are used. Both our methods, the deterministic moment approximation and the empirical Bayes estimation, can potentially be extended to these richer families.
Prior choice. Choosing priors in Bayesian neural networks remains an open issue. The hierarchical priors for feedforward neural networks that we use have been investigated before by Neal (1993) and MacKay (1995a), the latter proposing a “cheap and cheerful” heuristic, alternating optimization of weights and inverse variance parameters. Barber & Bishop (1998) also used a hierarchical prior and an efficient closed-form factored VB approximation; our approach can be seen as a point estimate to their approach in order to enable use of our closed-form moment approximation. Note that Barber & Bishop (1998) manipulate an expression for 〈hjhl〉 into a one-dimensional integral, whereas our
approach gives closed form approximations for this integral without need for numerical integration. Graves (2011) also used hierarchical Gaussian priors with flat hyperpriors, deriving a closed-form update for the prior mean and variance. Compared to these prior works our approach is rigorous and with sufficient data accurately approximates the Bayesian approach of integrating over the prior parameters.
Alternative inference procedures. As an alternative to variational Bayes, probabilistic backpropagation (PBP) (Hernández-Lobato & Adams, 2015) applies approximate inference in the form of assumed density filtering (ADF) to refine a Gaussian posterior approximation. Like in our work, each update to the approximate posterior requires propagating means and variances of activations through the network. (Hernández-Lobato & Adams, 2015) only consider the diagonal propagation case and homoscedastic regression. Since the original work, PBP has been generalized to classification (Ghosh et al., 2016) and richer posterior families such as the matrix variate Normal posteriors (Sun et al., 2017). Our moment approximation could be used to improve the inference accuracy of PBP, and since we handle minibatches of data rather than processing one data point at a time, our method is more computationally efficient.
Gaussianity in neural networks. Our demonstration of Gaussianity of ReLU network activations is also directly relevant to recent work on Gaussian process interpretations of deep neural networks (Matthews et al., 2018; Lee et al., 2017), validating the insight that activations in deep neural networks are closely approximated by Gaussian processes. Two recent works derived deterministic moment approximations for deep neural networks: Bibi et al. (2018), using Price’s theorem, derived exact first and second moment expressions for ReLU activations but limit themselves to the case of zero-mean Gaussian activations. Kandemir et al. (2018) also derive closed-form solutions to the ELBO for the case of diagonal Gaussian variational families. However, their approach is limited to linear layers without bias.
Markov chain Monte Carlo approaches. Another rich class of approximate inference methods for Bayesian neural networks are stochastic gradient Markov chain Monte Carlo (SG-MCMC) methods. These methods allow for approximate posterior parameter inference using unbiased log-likelihood estimates. Stochastic gradient Langevin dynamics (SGLD) was the first method in this class (Welling & Teh, 2011). SGLD is particularly simple and efficient to implement, but recent methods increase efficiency in the case of correlated posteriors by estimating the Fisher information matrix (Ahn et al., 2012) and extend Hamiltonian Monte Carlo to the stochastic gradient case (Chen et al., 2014). A complete characterization of SG-MCMC methods is given by (Ma et al., 2015; Gong et al., 2018). However, despite this progress, important theoretical questions regarding approximation guarantees for practical computational budgets remain (Nagapetyan et al., 2017). Moreover, while SG-MCMC methods work robustly in practice, they remain computationally inefficient, especially because evaluation of the posterior predictive requires evaluating an ensemble of models.
Wild approximations. The above methods are principled but often require sophisticated implementations; recently, a few methods aim to provide “cheap” approximations to the Bayes posterior. Dropout has been interpreted by Gal & Ghahramani (2016) to approximately correspond to variational inference. Likewise, Bootstrap posteriors (Lakshminarayanan et al., 2017; Fushiki et al., 2005; Harris, 1989) have been proposed as a general, robust, and accurate method for posterior inference. However, obtaining a bootstrap posterior ensemble of size k is computationally intense at k times the computation of training a single model.

heading: 6 EXPERIMENTS
text: We implement4 deterministic variational inference (DVI) as described above to train small ReLU networks on UCI regression datasets (Dheeru & Karra Taniskidou, 2017). The experiments address the claims that our methods for eliminating gradient variance and automatic tuning of the prior improve the performance of the final trained model. In Appendix D we present extended results to demonstrate that our method is competitive against a variety of models and inference schemes.
4Our implementation in TensorFlow is available at https://github.com/Microsoft/ deterministic-variational-inference
Deterministic vs. Stochastic We compare DVI with MCVI from equation 2 with S = 10 samples (we consider vanilla MCVI and discuss the local reparameterization trick in appendix E). The same model is used for each inference method: a single hidden layer of 50 units for each dataset considered, extending this to 100 units in the special case of the larger protein structure dataset, prot. Although neither DVI nor MCVI is limited to a particular choice of variational family q(w;θ), we use a factorized Gaussian family (i.e. a diagonal Cov(Wji,Wlk)). Factorization reduces the computational complexity of terms involving Cov(Wji,Wlk) in DVI5 from O(N2) to O(N), where N is the number of elements in W (see appendix A.1). Additionally, both methods use the same EB prior from equation 13 with a broad inverse Gamma hyperprior (α = 1, β = 10) and an independent sλ for each linear transformation. Each dataset is split into random training and test sets with 90% and 10% of the data respectively. This splitting process is repeated 20 times and the average test performance of each method at convergence is reported in table 2 (see also learning curves in appendix F). We see that DVI consistently outperforms MCVI, by up to 0.35 nats per data point on some datasets. The computationally efficient diagonal-DVI (dDVI) surprisingly retains much of this performance. By default we use the heteroscedastic model, and we observe that this uniformly delivers better results than a homoscedastic model (hoDVI; rightmost column in table 2) on these datasets with no overfitting issues6.
Empirical Bayes In Figure 5 we compare the performance of networks trained with manual tuning of a fixed Gaussian prior to networks trained with the automatic EB tuning. We find that the EB method consistently finds priors that produce models with competitive or significantly improved test log-likelihood relative to the best manual setting. Since this observation holds across all datasets considered, we say that our method is “robust”. Note that the EB method can outperform manual tuning because it automatically finds different prior variances for each weight matrix, whereas in the manual tuning case we search over a single hyperparameter control-
ling all prior variances. An additional ablation study showing the relative contribution of our deterministic approach and the EB prior are shown in appendix D.1.
5For MCVI with full rank covariance, Cholesky decomposition required for sampling is O(N3). 6Note that this result is non-trivial because heteroscedastic models are more complex and could result in
poorer approximate inference leading to worse test performance

heading: 7 CONCLUSION
text: We introduced two innovations to make variational inference for neural networks more robust: 1. an effective deterministic approximation to the moments of activations of a neural networks; and 2. a simple empirical Bayes hyperparameter update. We demonstrate that together these innovations make variational Bayes a competitive method for Bayesian inference in neural heteroscedastic regression models.
Bayesian neural networks have been shown to substantially improve upon standard networks in these settings where calibrated predictive uncertainty estimates, sequential decision making, or continual learning without catastrophic forgetting are required (see e.g. Oliveira et al. (2016); Gal et al. (2017); Nguyen et al. (2018)). In future work, the new innovations proposed in this paper can be applied to these areas. In the sequential decision making and continual learning applications, approximate Bayesian inference must be run as an inner loop of a larger algorithm. This requires a robust and automated version of BNN training: this is precisely where we believe the innovations in this paper will have large impact since they pave the way to automated and robust deployment of BBNs that do not involve an expert in-the-loop.

heading: APPENDIX
text: A MOMENTS OF THE ACTIVATION VARIABLES a`
Under assumption of independence of h,W and b, we can write:
〈ai〉 = 〈hjWji〉+ 〈bi〉 = 〈hj〉 〈Wji〉+ 〈bi〉 (14)
Cov(ai, ak) = Cov(hjWji, hlWlk) + Cov(bi, bk)
= 〈hjWji hlWlk〉 − 〈hjWji〉 〈hlWlk〉+ Cov(bi, bk) = 〈hjhl〉 〈WjiWlk〉 − 〈hj〉 〈hl〉 〈Wji〉 〈Wlk〉+ Cov(bi, bk) = 〈hjhl〉 [Cov(WjiWlk) + 〈Wji〉 〈Wlk〉] − 〈hj〉 〈hl〉 〈Wji〉 〈Wlk〉+ Cov(bi, bk) = 〈hjhl〉Cov(Wji,Wlk) + 〈Wji〉Cov(hj , hl) 〈Wlk〉+ Cov(bi, bk), (15)
which is seen in the main text as equation 3. For Heaviside and ReLU activation functions, closed forms exist for 〈hj〉 in equation 14:
Heaviside 〈hj〉 = 1√2πΣ′jj ∫ ∞ 0 e − 12Σ′jj (αj−〈a′j〉)2 dαj = Φ ( µ′j )
ReLU 〈hj〉 = 1√2πΣ′jj ∫ ∞ 0 αje − 12Σ′jj (αj−〈a′j〉)2 dαj = √ Σ′jjSR(µ ′ j),
where SR(x) ··= φ(x) + xΦ(x) is a “soft ReLU”, φ and Φ represent the standard Gaussian PDF and CDF, and we have introduced the dimensionless variables µ′j = 〈 a′j 〉 / √
Σ′jj . These results are is sufficient to evaluate equation 14, so in the following sections we turn to each term from equation 15.
A.1 EVALUATION OF TERM 1: 〈hjhl〉Cov(Wji,Wlk)
In the general case, we can use the results from section A.2 to evaluate off-diagonal 〈hjhl〉. However, in our experiments we always consider the the special case where Cov(Wji,Wlk) is diagonal. In this case we can write the first term in equation 15 as (reintroducing the explicit summation):∑
jl 〈hjhl〉Cov(Wji,Wlk) = ∑ jl 〈hjhl〉 δjlδikVar(Wji)
=δik ∑ j 〈zjzj〉Var(Wji)
= diag [vVar(W )]
i.e. this term is a diagonal matrix with the diagonal given by the left product of the vector vj = 〈hjhj〉 with the matrix Var(Wki). Note that 〈hjhj〉 can be evaluated analytically for Heaviside and ReLU activation functions:
Heaviside 〈hjhj〉 = 1√2πΣ′jj ∫ ∞ 0 e − 12Σ′jj (αj−〈a′j〉)2 dαj = Φ ( µ′j )
ReLU 〈hjhj〉 = 1√2πΣ′jj ∫ ∞ 0 α2je − 12Σ′jj (αj−〈a′j〉)2 dαj = Σ ′ jj [ µ′jφ(µ ′ j) + (1 + µ ′2 j )Φ(µ ′ j) ]
A.2 EVALUATION OF TERM 2: 〈Wji〉Cov(hj , hl) 〈Wlk〉
Evaluation of Cov(hj , hl) requires an expression for 〈hjhl〉. From equation 5, we write:
〈hjhl〉 ∝ ∫ f(αj)f(αl) exp [ − 12P (αj , αl;a′,Σ′) ] dαjdαl, (16)
where P is the quadratic form:
P (αj , αl;a ′,Σ′) = ( αj−〈a′j〉 αl−〈a′l〉 )> ( Σ′jj Σ ′ jl Σ′lj Σ ′ ll )−1( αj−〈a′j〉 αl−〈a′l〉 ) = ( ηj−µ′j ηl−µ′l )> ( 1 ρ′jl ρ′lj 1 )−1 ( ηj−µ′j ηl−µ′l ) .
Here we have introduced further dimensionless variables ηj = αj/ √ Σ′jj , ηl = αl/ √
Σ′ll and ρ′jl = Σ ′ jl/ √ Σ′jjΣ ′ ll. We can then rewrite equation 16 in terms of a dimensionless integral I using a scale factor S′jl that is 1 for the Heaviside non-linearity or Σ ′ jl/ρ ′ jl for ReLU:
〈hjhl〉 = S′jlI(µ′j , µ′l, ρ′jl) ; I = 1
Z
∫ f(ηj)f(ηl) exp [ − 12P (ηj , ηl;µ′,ρ′) ] dηjdηl.
The normalization constant, Z, is evaluated by integrating over e−P/2 and is explicitly written as Z = 2πρ̄′jl, where ρ̄ ′ jl = √ 1− ρ′2jl. Now, following equation 6, we have the task to write I as an asymptote A plus a decaying correction e−Q. To evaluate A and Q, we have to insert the explicit form of the non-linearity f , which we do for Heaviside and ReLU functions in the next sections.

heading: A.2.1 HEAVISIDE NON-LINEARITY
text: For the Heaviside activation, we can represent the integral I as the shaded area under the Gaussian in the upper-left quadrant shown below. In general, this integral does not have a closed form. However, for µ′j →∞, vanishing weight appears under the Gaussian in the upper-right quadrant, so we can write down the asymptote of the integral in this limit:
lim µj→∞
I = 1
Z ∫ ∞ ηj=−∞ ∫ ∞ ηl=0 exp [ − 12P (ηj , ηl;µ′, ρ′jl) ] dηjdηl = Φ(µ ′ l)
𝜇𝑗 → ∞ 𝜇𝑙
𝜂𝑙
𝜂𝑗
Here we performed the integral by noticing that the outer integral over ηj marginalizes out ηj from the bivariate Gaussian, leaving the inner integral as the definition of the Gaussian CDF. By symmetry, we also have limµ′l→∞ I = Φ(µ ′ j) and limµ′j,l→−∞ I = 0. We can then write down the following symmetrized form that satisfies all the limits required to qualify as an asymptote:
A = Φ(µ′j)Φ(µ ′ l)
To compute the correction factor we evaluate the derivatives of (I − A) at the origin up to second order to match the moments of e−Q for quadratic Q. Description of this process is found below
Zeroth derivative At the origin µj = µl = 0, we can diagonalize the quadratic form P :
P (ηj , ηl; 0, ρ ′ jl) = 1 2ρ̄′2jl
( η2j − 2ρ′jlηjηl + η2l ) = 1
4ρ̄′2jl
( ξ2+ + ξ 2 − ) ,
where ξ± = √ 1∓ ρ(η1 ± η2). Performing this change of variables in the integral gives:
I = 1
4πρ̄′2jl ∫ H exp [ 1 4ρ̄′2jl ( ξ2+ + ξ 2 − )] dξ+dξ− = ψ/π
𝜂𝑙
𝜂𝑗
𝜉−
𝜉+
ℋ ℋ
𝜓
where we integrated in polar coordinates over the region H in which the Heaviside function is non-zero. The angle ψ can be found from the coordinate transform between η and ξ as7:
ψ = arctan √ 1+ρ′jl 1−ρ′jl = π2 − 12 arccos ρ′jl.
Since A|µ′=0 = Φ(0)Φ(0) = 1/4, we can evaluate: (I −A)|µ′=0 = π2 − 12 arccos ρ′jl − 14
= 12π arcsin ρjl 7Here we use the identity cos(2 arctanx) = cos2 arctanx− sin2 arctanx = 1−x 2
1+x2
First derivative Performing a change of variables xi = ηi − µ′i, we can write I as:
I = 1
Z
∫ H(xj + µ ′ j)H(xl + µ ′ l) exp [ − 12P (xj , xl; 0, ρ′jl) ] dxjdxl
where H is the Heaviside function. Now, using ∂xH(x) = δ(x), we have:
∂ ∂µ′j ∣∣∣ µ′=0 I = 1 Z ∫ δ(xj)H(xl) exp [ − 12P (xj , xl; 0, ρ′jl) ] dxjdxl = 1 2 √ 2π . (17)
In addition, using ∂xΦ(x) = φ(x), we have:
∂ ∂µ′j ∣∣∣ µ′=0 A = 1 2 √ 2π =⇒ ∂∂µ′j ∣∣∣ µ′=0 (I −A) = 0.
By symmetry (I − A) also has zero gradient with respect to µ′l at the origin. Therefore Q has no linear term in µ′.
Second derivative Taking another derivative in equation 17 gives:
∂2
∂µ′2j ∣∣∣ µ′=0 I = − 1 Z ∫ δ(xj)H(xl) ρ′jl ρ̄′2jl xl exp [ − 12P (xj , xl; 0, ρ′jl) ] dxjdxl = − ρ ′ jl 2πρ̄′jl .
where we used the identity ∫ f(x)∂xδ(x)dx = − ∫ δ(x)∂xf(x)dx, which holds for arbitrary f . In addition, we have:
∂2
∂µ′2j ∣∣∣ µ′=0 A = 0 =⇒ ∂2 ∂µ′2j ∣∣∣ µ′=0 (I −A) = − ρ ′ jl 2πρ̄′jl .
and the same result holds for the second derivative w.r.t. µ′l. To complete the Hessian, it is a simple extension of previous results to show that:
∂2
∂µ′j∂µ ′ l ∣∣∣ µ′=0 (I −A) = 1−ρ̄ ′ jl 2πρ̄′jl .
Now that we have obtained derivatives of the residual (I − A) up to second order we propose a correction factor of the form e−Q where Q is truncated at quadratic terms:
Q = − log α2π + β ( µ′2j + µ ′2 l ) + γµ′2j µ ′2 l .
We then find the coefficients {α, β, γ} by matching ( !=) derivatives at µ = 0:
e−Q ∣∣ µ=0 = α2π ! = arcsin ρ′jl 2π =⇒ α = arcsin ρ′jl
∂ ∂µ′i e−Q ∣∣∣ µ=0 = 0 ! = 0
∂2 ∂µ′2i e−Q ∣∣∣ µ=0 = −2β α2π ! = − ρ ′ jl 2πρ̄′jl =⇒ β = ρ ′ jl 2αρ̄′jl
∂2 ∂µ′j∂µ ′ l e−Q ∣∣∣ µ=0 = γ α2π ! = 1−ρ̄′jl 2πρ̄′jl
=⇒ γ = 1−ρ̄ ′ jl
αρ̄′jl
This yields the expression seen in table 1 of the main text.

heading: A.2.2 RELU NON-LINEARITY
text: As in the Heaviside case, we begin by computing the asymptote of I by inspecting the limit as µ′j →∞:
lim µj→∞
I = 1
Z ∫ ∞ ηj=−∞ ∫ ∞ ηl=0 ηjηl exp [ − 12P (ηj , ηl;µ′, ρ′jl) ] dηjdηl
= 1√ 2π ∫ ∞ ηl=0 µ′jηle − 12 (ηl−µ ′ l) 2 dηl + ρ′jl√ 2π ∫ ∞ ηl=0 ηl(ηl − µl)e− 1 2 (ηl−µ ′ l) 2 dηl = µ′j SR(µ ′ l) + ρ ′ jlΦ(µ ′ l) (18)
Now, we construct a full 2-dimensional asymptote by symmetrizing equation 18 (using properties SR(x) → x and Φ(x) → 1 as x → ∞ to check that the correct limits are preserved after symmetrizing):
A = SR(µ′j)SR(µ ′ l) + ρ ′ jlΦ(µ ′ j)Φ(µ ′ l)
Next we compute the correction factor e−Q. The details of this procedure closely follow those for the Heaviside non-linearity of the previous section, so we omit them here (and in practice we use Mathematica to perform the intermediate calculations). The final result is presented in table 1 of the main text.

heading: B LOG-LIKELIHOOD AND POSTERIOR PREDICTIVE COMPUTATION
text: Here we give derivations of expressions quoted in section 3.2. In section B.1 we justify the intuitive result that expectation of the ELBO reconstruction term over q(w;θ) can be re-written as an expectation over q̃(aL). We then derive expected log-likelihoods and posterior predictive distributions for the cases of univariate Gaussian regression and classification. The latter sections are arranged as follows:
Regression Classification
Log-likelihood section B.2 section B.4 Posterior predictive section B.3 section B.5
B.1 LOG-LIKELIHOODS: FROM Ew TO EaL
We begin by rewriting the reconstruction term for data point (x, y) in terms of aL: Ew∼q [log p(y|w)] = ∫ q(w) log p(y|w) dw = ∫ q(aL)q(w|aL) log p(y|w) dw daL
where we have suppressed explicit conditioning on x for brevity. Our goal now is to perform the integral over w, leaving the expectation in terms of aL only, thus allowing it to be evaluated using the approximation q̃(aL) from section 3.1.
To eliminate w, consider the case where the output of the model is a distribution p(y|aL) that is a parameter-free transformation of aL (e.g. aL are logits of a softmax distribution for classification or the moments of a Gaussian for regression). Since the model output is conditioned only on aL, we must have p(y|w) = p(y|aL) for all configurationsw that satisfy the deterministic transformation aL =M(x;w), whereM is the neural network (i.e p(y|w) = p(y|aL) for allw where q(w|aL) is non-zero). This allows us to write:∫
q(w|aL) log p(y|x,w) dw = log p(y|aL) ∫ q(w|aL) dw = log p(y|aL),
so the reconstruction term becomes: Ew∼q [log p(y|x,w)] = ∫ q(aL) log p(y|aL)daL = EaL∼q(aL) [ log p(y|aL) ] .
This establishes the equivalence given in equation 7 in the main text. Since we are using an approximation to q, we will actually compute EaL∼q̃(aL) [ log p(y|aL) ] .

heading: B.2 UNIVARIATE REGRESSION: LOG-LIKELIHOOD
text: Here we give a derivation of equation 8 from the main text. Throughout this section we label the 2 elements of the final activation vector as aL = (m, `). We first insert the Gaussian form for p(y|aL) ∼ N [ m, e` ] into the log-likelihood expression:
EaL∼q̃(aL) [ log p(y|aL) ] = − 12EaL∼q̃(aL) [ log (2π exp(`)) + exp(−`)(y −m)2 ] = − 12 log 2π − 12 〈`〉 − EaL∼q̃(aL) [ exp(−`)(y −m)2 ] . (19)
Now we use the Gaussian form of q̃(aL): q̃(aL) ∝ exp [ − 12X>(ΣL)−1X ] ; X = ( m−〈m〉 `−〈`〉 ) . and note that∫ exp(−`) exp [ − 12X>(ΣL)−1X ] dm d`
= exp (−〈`〉) ∫ exp [ − 12X>(ΣL)−1X − (`− 〈`〉) ] dmd`
= exp (−〈`〉) ∫ exp [ − 12X>(ΣL)−1X − e>` X ] dm d` = exp (
Σ`` 2 − 〈`〉
) ∫ exp [ − 12 (X> + e>` ΣL)(ΣL)−1(X + e`ΣL) ] dm d`, (20)
where e>` = (0, 1) is the unit vector in the ` coordinate, and we completed the square to obtain the final line. Inserting equation 20 into equation 19 and marginalizing out the ` coordinate gives:
EaL∼q̃(aL) [ log p(y|aL) ] = − 12 [ log 2π + 〈`〉+ eΣ``/2−〈`〉√
2πΣmm
∫ (y −m)2 exp ( − [m−(〈m〉−Σm`)] 2
2Σmm
) dm ] .
Finally, performing the integral over m gives the result seen in equation 8.

heading: B.3 UNIVARIATE REGRESSION: POSTERIOR PREDICTIVE DISTRIBUTION
text: Here we give a derivation of equation 9 from the main text. We first calculate the first and second moments of the predictive distribution under the approximation q(aL) ≈ q̃(aL):
Ey∼p(y)[y] = ∫ yp(y|aL)q̃(aL) dy daL
= ∫ [∫ yp(y|aL) dy ] q̃(aL) daL
= ∫ mq̃(aL) daL
= 〈m〉
Var[y] = Var[Ey∼p(y|aL)(y)] + Ey∼p(y) [ Var(y|aL) ] = Var[m] + Ey∼p(y) [ e` ]
= Σmm +
∫ e`p(y|aL)q̃(aL) daL dy
= Σmm +
∫ e`q̃(aL) daL
= Σmm + exp (〈`〉+ Σ``/2) where the final integral in the variance computation is performed by inserting the Gaussian form for q̃(aL) and completing the square. Then we assume normality of the predictive distribution to obtain the result in equation 9.

heading: B.4 CLASSIFICATION: LOG-LIKELIHOOD
text: There is no exact form for the expected log-likelihood for multivariate classification with logits aL. However, using the second-order Delta method (Small, 2010), we find the expansion
EaL∼q̃(aL) [ log p(y|aL) ] = 〈 aL 〉 − EaL∼q̃(aL) [ logsumexp(aL) ] ≈ 〈 aL 〉 − logsumexp( 〈 aL 〉 )− 12 ( p> diag(ΣL)− p>ΣLp ) , (21) To derive this expansion, we first state the second order expansion for the expectation of a function g of random variable x using the Delta method as follows8:
E [g(x)] ≈ g (E[x]) + 12 ∑ ij [ Cij ∂2g ∂xi∂xj ] x=E[x] , (22)
where Cij = Cov(xi, xj). Now we note that the logsumexp function has a simple Hessian ∂2
∂xi∂xj logsumexp(x) = δijpi − pipj ,
where p = softmax(x). Putting these results together allows us to write: E [logsumexp(x)] ≈ logsumexp (E[x]) + 12 [ p> diag(C)− p>Cp ] x=E[x] ,
This result is sufficient to complete the derivation of equation 21 and enable training of a classifier using our method.
8This result is obtained by Taylor expansion inside the expectation.

heading: B.5 CLASSIFICATION: POSTERIOR PREDICTIVE DISTRIBUTION
text: Using the same second-order Delta method, we find the following expansion for the posterior predictive distribution:
p(y) = EaL∼q̃(aL) [ p(y|aL) ] ≈ p [ 1 + p>ΣLp−ΣLp + 12 diag(ΣL)− 12p> diag(ΣL) ] .
(23) where p = softmax( 〈 aL 〉 ).
For this expansion, we begin by computing the Hessian:
∇i∇j [p]k = ∂ 2 ∂xi∂xj pk = [2pipj − (δkipj + δkjpi) + δikδjk − δijpi] pk,
where p = softmax(x), and we used the intermediate result∇j [p]k = δjkpk − pjpk. Then we can form the product:
Tr [C∇∇p] = p [ 2p>Cp− 2Cp + diag(C)− p> diag(C) ] and insert this into equation 22 to obtain equation 23.
Preliminary experiments show that good results are obtained either using these approximations or a lightweight MC approximation just to perform the mapping of aL to (log)p after the deterministic heavy-lifting of computing aL. In this work we are primarily concerned with demonstrating the benefits of the moment propagation method from section 3.1, so we limit our experiments to regression examples without additional complication from approximation of the likelihood function.

heading: C DEEPER NETWORKS
text: Here we consider the applicability of our method to the regime of deep, narrow networks. This regime is challenging because for small hidden dimension the Gaussian approximation for a (reliant on the CLT) breaks down, and these errors accumulate as the net becomes deep. We empirically explore this potential problem by investigating deep networks containing 5 layers of only 5, 25 or 125 units each. Figure 6 shows results analogous to figure 3 that qualitatively illustrate how well our approximation matches the true variational distribution of output activations both at the start and end of training. We see that our CLT-based approximation is good in the 125- and 25-unit cases, but is poor in the 5-unit case. Since it is generally considered that optimization of neural networks only works well in the high dimensional setting with at least a few tens of hidden units, these empirical observations suggest that our approximation is applicable in practically relevant architectures.

heading: C.1 SKIP CONNECTIONS
text: Training deep networks is considered difficult even in the traditional maximum-likelihood setting due to the problems of exploding and vanishing gradients. A popular approach to combat these issues is to add skip connections to the architecture. Here we derive the necessary results to add skip connections to our deterministic BNN.
We consider a simple layer with skip connections of the following form:
h = f(a′),
δ = hW + b,
a = a′ + δ.
The moment propagation expressions for this layer are (using the bilinearity of Cov):
〈ai〉 = 〈a′i〉+ 〈δi〉 , Cov(ai, ak) = Cov(a ′ i, a ′ k) + Cov(δi, δk) + Cov(a ′ i, δk) + Cov(δi, a ′ k),
where 〈δi〉 and Cov(δi, δk) can be computed using analogy to equations 14 and 15. This just leaves computation of Cov(a′i, δk) and its transpose, which can be performed analytically using integral
results and methods borrowed from appendix A. Cov(a′i, δk) = Cov(a ′ i, ∑ j f(a′j)Wjk) + Cov(a ′ i, bk)
= ∑ j (〈 a′if(a ′ j) 〉 − 〈a′i〉 〈 f(a′j) 〉) 〈Wjk〉
= ∑ j 〈Wjk〉 (Σ′ij) { (Σ′jj) −1φ(µ′j) Heaviside Φ(µ′j) ReLU
Using this result, we implement a 5-layer, 25-unit network with skip connections. In figure 6(d) we qualitatively verify the validity of our approximation on this architecture by observing a good match with Monte Carlo simulations using 20k samples.
Published as a conference paperatIC L R 2019 D EXTENDED RESULTS Here we include comparison with a number of different models and inference schemes on the 9 UCI datasets considered in the main text. We report test log-likelihoods at convergence and find that our method is competitive or superior to a range of state-of-the-art techniques (reproduced from Bui et al. (2016)). Dataset bost conc ener kin8 nava powe prot wine yach |D| 506 1030 768 8192 11934 9568 45730 1588 308 dx 13 8 8 8 16 4 9 11 6
GP 50 −2.22± 0.07 −2.85± 0.02 −1.29± 0.01 1.31± 0.01 4.86± 0.04 −2.66± 0.01 −2.95± 0.05 −0.67± 0.01 −1.15± 0.03 DGP-1 50 −2.33± 0.06 −3.13± 0.03 −1.32± 0.03 0.68± 0.07 3.60± 0.33 −2.81± 0.01 −2.55± 0.03 −0.35± 0.04 −1.39± 0.14 DGP-2 50 −2.17± 0.10 −2.61± 0.02 −0.95± 0.01 1.79± 0.02 4.77± 0.32 −2.58± 0.01 −2.11± 0.04 −0.10± 0.03 −0.99± 0.07 DGP-3 50 −2.09± 0.07 −2.63± 0.03 −0.95± 0.01 1.93± 0.01 5.11± 0.23 −2.58± 0.01 −2.03± 0.07 −0.13± 0.02 −0.94± 0.05 GP 100 −2.16± 0.07 −2.65± 0.02 −1.11± 0.02 1.68± 0.01 5.51± 0.03 −2.55± 0.01 −2.52± 0.07 −0.57± 0.02 −1.26± 0.03 DGP-1 100 −2.37± 0.10 −2.92± 0.03 −1.21± 0.02 1.09± 0.04 3.75± 0.37 −2.67± 0.02 −2.18± 0.06 0.07± 0.03 −1.34± 0.10 DGP-2 100 −2.09± 0.06 −2.43± 0.02 −0.90± 0.01 2.31± 0.01 5.13± 0.27 −2.39± 0.02 −1.51± 0.09 0.37± 0.02 −0.96± 0.06 DGP-3 100 −2.13± 0.09 −2.44± 0.02 −0.91± 0.01 2.46± 0.01 5.78± 0.05 −2.37± 0.02 −1.32± 0.06 0.25± 0.03 −0.80± 0.04 VI(KW)-2 −2.64± 0.02 −3.07± 0.02 −1.89± 0.07 2.91± 0.10 6.10± 0.19 −2.28± 0.02 −0.42± 0.31 −0.85± 0.01 −1.92± 0.03 SGLD-2 −2.38± 0.06 −3.01± 0.03 −2.21± 0.01 1.68± 0.00 3.21± 0.02 −2.61± 0.01 −1.23± 0.01 0.14± 0.02 −3.23± 0.03 SGLD-1 −2.40± 0.05 −3.08± 0.03 −2.39± 0.01 1.28± 0.00 3.33± 0.01 −2.67± 0.00 −3.11± 0.02 −0.41± 0.01 −2.90± 0.02 HMC-1 −2.27± 0.03 −2.72± 0.02 −0.93± 0.01 1.35± 0.00 7.31± 0.00 −2.70± 0.00 −2.77± 0.00 −0.91± 0.02 −1.62± 0.01 PBP-1 −2.57± 0.09 −3.16± 0.02 −2.04± 0.02 0.90± 0.01 3.73± 0.01 −2.84± 0.01 −2.97± 0.00 −0.97± 0.01 −1.63± 0.02 VI(G)-1 −2.90± 0.07 −3.39± 0.02 −2.39± 0.03 0.90± 0.01 3.73± 0.12 −2.89± 0.01 −2.99± 0.01 −0.98± 0.01 −3.44± 0.16 VI(KW)-1 −2.43± 0.03 −3.04± 0.02 −2.38± 0.02 2.40± 0.05 5.87± 0.29 −2.66± 0.01 −1.84± 0.07 −0.78± 0.02 −1.68± 0.04 Dropout-1 −2.46± 0.25 −3.04± 0.09 −1.99± 0.09 0.95± 0.03 3.80± 0.05 −2.89± 0.01 −2.80± 0.05 −0.93± 0.06 −1.55± 0.12 DVI −2.41± 0.02 −3.06± 0.01 −1.01± 0.06 1.13± 0.00 6.29± 0.04 −2.80± 0.00 −2.85± 0.01 −0.90± 0.01 −0.47± 0.03 dDVI −2.42± 0.02 −3.07± 0.02 −1.06± 0.06 1.13± 0.00 6.22± 0.06 −2.80± 0.00 −2.84± 0.01 −0.91± 0.02 −0.47± 0.03 DVI-MC −2.41± 0.02 −3.05± 0.01 −1.00± 0.06 1.13± 0.00 6.25± 0.03 −2.80± 0.00 −2.85± 0.01 −0.93± 0.04 −0.55± 0.03 DVI-MC-softplus −2.42± 0.02 −3.06± 0.02 −1.03± 0.05 1.13± 0.00 6.20± 0.04 −2.80± 0.01 −2.85± 0.01 −0.89± 0.01 −0.54± 0.03 MCVI −2.46± 0.02 −3.07± 0.01 −1.03± 0.04 1.14± 0.00 5.94± 0.05 −2.80± 0.00 −2.87± 0.01 −0.92± 0.01 −0.68± 0.03 MCVI-softplus −2.47± 0.02 −3.08± 0.02 −1.02± 0.05 1.14± 0.00 5.99± 0.02 −2.80± 0.00 −2.85± 0.01 −0.92± 0.01 −0.69± 0.03
Table 3: Average test log likelihood on UCI datasets. See table 4 for model glossary and implementation references.
21

heading: D.1 ABLATION STUDY
text: Here we provide an ablation study that indicates the individual contributions of (1) the deterministic approximation and (2) the the empirical Bayes prior. We consider all combinations of DVI or MCVI with and without empirical Bayes. In the DVI-fixed and MCVI-fixed cases without empirical Bayes we use a fixed zero-mean Gaussian prior during training and we perform separate runs to tune the prior variance, reporting the best performance achieved (cf. figure 5)9. Since the EB approach requires no hyperparameter tuning between the datasets shown, these results hide the considerable computational advantaged that the EB approach brings.
E VARIANCE REDUCTION AND THE LOCAL REPARAMETERIZATION TRICK
By eliminating MC sampling and its associated variance entirely, our method directly tackles the problem of high variance gradient estimates that hinder MC approaches to training of BNNs. Alternative methods that only reduce variance have been considered, and among these, the local
9Note that the EB method can out perform manual tuning because it automatically finds different prior variances for each weight matrix, whereas manual tuning case searches over a single hyperparameter controlling all prior variances.
reparameterization trick (Kingma et al., 2015) is particularly popular. Similar to our approach, the local reparameterization trick maps the uncertainty in the weights to an uncertainty in activations, however, unlike the fully deterministic DVI, MC methods are then used to propagate this uncertainty through non-linearities. The benefits of MCVI with the reparameterization trick (rMCVI) over vanilla MCVI are two-fold:
• The variance of the gradient estimates during back propagation are reduced (see details in Kingma et al. (2015)). • Since the sampling dimension in rMCVI only appears on the activations and not on the
weights, an H ′×H linear transform can be implemented using SB×H ′ by H ′×H matrix multiplies (where S is the number of samples and B is the batch size). This contrasts with the S×B×H ′ by S×H ′×H batched matrix multiply required for MCVI. Although both of these algorithms have the same asymptotic complexityO(SBH ′H), a single large matrix multiplication is generally more efficient on GPUs than smaller batched matrix multiplies.
Figure 7 shows empirical studies of the gradient variance and runtime for rMCVI vs. MCVI applied to the model described in section 3.1.1 and figure 3. To evaluate the gradient variance, we initialize the model with partially trained weights and measure the variance of the gradient of the ELBO reconstruction term L with respect to variational parameters. Specifically, we inspect the gradient with respect to the parameters ΣLq describing the variance of the q distribution for the weight matrix in the final layer.
Gradient variance ··= mean s∈ΣLq
[ Var ( ∂L ∂s )] The plots in figure 7 serve to show that rMCVI is not fundamentally different from MCVI, and the performance of one (on either the speed or variance metric) can be transformed into the other by varying the number of samples. A comparison of DVI with rMCVI is included in table 3 using the implementation labelled as “VI(KW)-1”.

heading: F LEARNING CURVES
text: Figure 8 shows the test log-likelihood during the training of the models from table 2 using DVI and MCVI inference algorithms. Since the underlying model is identical, both methods should achieve the same test log-likelihood given infinite time and infinite MC samples (or a suitable learning rate schedule) to mitigate the increased variance of the MCVI method. However, since we use only 10 samples and do not employ a leaning rate schedule, we find that MCVI converges to a log-likelihood that is consistently worse than that achieved by DVI.

Review summaries:
The authors propose a new approach to perform deterministic variational inference for feed-forward BNN with specific nonlinear activation functions by approximating layerwise moments. Under certain conditions, the authors show that the proposed method achieves better performance than existing Monte Carlo variational inference. This paper is interesting since most of the existing works focus on Monte Carlo variational inference. The main contribution of this paper is to perform Gaussian approximation. The authors show that for specific activation functions, the Gaussian approximation is reasonable. The main concern is the cumulative error due to the Gaussian approximation. Since the authors argue that the proposed method fixes the issues of stochastic VI for BNN, the authors should also investigate/clarify the following cases. 

This work is tackling two difficulties in current VB applied to DNNs ("Bayes by backprop"). First, MC approximations of intractable expectations are replaced by deterministic approximations. While this has been done before, the solution here is new and very interesting. Second, a Gaussian prior with length scales is learned by VB empirical Bayes alongside the normal training, which is also very useful.