Title: DETERMINISTIC VARIATIONAL INFERENCE FOR ROBUST BAYESIAN NEURAL NETWORKS

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.

1. INTRODUCTION:

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

2. VARIATIONAL INFERENCE IN BAYESIAN NEURAL NETWORKS:

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 θ D KL [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):

EQUATION

Analytic results exist for the KL term in the ELBO for careful choice of prior and variational distributions (e.g. Gaussian families). However, when M 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:

E w∼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. 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.

3. DETERMINISTIC VARIATIONAL APPROXIMATION:



3.1. MOMENT PROPAGATION:

We begin by considering activation propagation (figure 1 (a)), with the aim of deriving the form of an approximation q(a L ) to the final layer activation distribution q(a L ) that will be passed to the likelihood computation. We compute a L by sequentially computing the distributions for the activations in the preceding layers. Concretely, we define the action of the l th layer that maps a (l-1) to a l as follows:

h l = f (a (l-1) ), a l = h l W l + b l ,

where f is a non-linearity and {W l , b l } ⊂ w are random variables representing the weights and biases of the l th 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 a l 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 a i will be normally distributed regardless of the potentially complicated distribution for h j 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 model W , b and h as independent random variables, allowing us to write:

EQUATION

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 W ji and Cov(W ji , W lk )), then the only difficult terms are the moments of h:

EQUATION

h j h l ∝ 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) A(µ 1 , µ 2 , ρ) Q(µ 1 , µ 2 , ρ) Heaviside Φ(µ 1 )Φ(µ 2 ) -log( g h ρ 2π ) + ρ 2g h ρ µ 2 1 + µ 2 2 -2ρ 1+ ρ µ 1 µ 2 + O(µ 4 ) ReLU SR(µ 1 )SR(µ 2 ) + ρΦ(µ 1 )Φ(µ 2 ) -log( gr 2π ) + ρ 2gr(1+ ρ) µ 2 1 + µ 2 2 -arcsin ρ-ρ grρ µ 1 µ 2 + O(µ 4 )

Table 1 : Forms for the components of the approximation in equation 6 for Heaviside and ReLU non-linearities. Φ is the CDF of a standard Gaussian, SR is a "soft ReLU" that we define as SR(x) = φ(x) + xΦ(x) where φ is a standard Gaussian, ρ = 1ρ2 , g h = arcsin ρ and g r = g h + ρ

1+ ρ

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 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(a L ), our explicit multivariate Gaussian approximation to q(a L ). 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.

EQUATION

3.1.1. EMPIRICAL VERIFICATION:

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 sufficient 2 ; 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 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 a L 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 for m 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 a L distributions reported by our deterministic approximation and MC evaluation using 20k samples from q(w; θ). This comparison is qualitatively excellent for all cases considered. We show the time to propagate a batch of 10 activation vectors through a single d × d layer. For MCVI we label curves with the number of samples used, and we show quadratic and cubic scaling guides-to-the-eye (black). Black dots indicate where our implementation runs out of memory (16GB).

Computational efficiency In traditional MCVI, propagation of S samples of d-dimensional activations through a layer containing a d × d-dimensional transformation requires O(Sd 2 ) compute and O(Sd) memory. Our DVI method approximates the S → ∞ limit, while only demanding O(d 3 ) compute and O(d 2 ) memory (the additional factor of d arises from manipulation of the quadratically large covariance matrix Cov[h j , h l ]). 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(h j , h j ) 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.

3.2. LOG-LIKELIHOOD EVALUATION:

To use the moment propagation procedure derived above for training BNNs, we need to build a function L that maps final layer activations a L 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(a L ).

EQUATION

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 σ 2 y (x) = (x; w). The parameters of this Gaussian are read off as the elements of a 2-dimensional output layer a L = (m, ) so that p(y|a L ) = N y|m, e . Recall that these parameters themselves are uncertain and the statistics a L and Σ L can be computed following section 3.1. Inserting the Gaussian forms for p(y|a L ) and q(a L ) into equation 7 and performing the integral (see appendix B.2) gives a closed form expression for the ELBO reconstruction term:

EQUATION

)

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|a L ) q(a L ) da L ≈ N y m , Σ mm + e +Σ /2 . (9)

4. EMPIRICAL BAYES FOR VARIATIONAL BNNS:

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:

EQUATION

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:

EQUATION

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 = E w∼q log p(y|h L (w)) -{D KL [q(w; θ)||p(w|s * )p(s * )]} =⇒ s λ * = argmin s λ D KL q(w; θ)||p(w λ |s λ ) -log p(s λ ) (12)

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

EQUATION

)

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.

5. RELATED WORK:

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 h j h l 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) (Herná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. (Herná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.

6. EXPERIMENTS:

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

7. CONCLUSION:

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.

APPENDIX A MOMENTS OF THE ACTIVATION VARIABLES a:

Under assumption of independence of h, W and b, we can write:

EQUATION

)

which is seen in the main text as equation 3. For Heaviside and ReLU activation functions, closed forms exist for h j in equation 14:

Heaviside h j = 1 √ 2πΣ jj ∞ 0 e - 1 2Σ jj (αj-a j ) 2 dα j = Φ µ j ReLU h j = 1 √ 2πΣ jj ∞ 0 α j e - 1 2Σ jj (αj-a j ) 2 dα j = Σ jj SR(µ 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: h j h l Cov(W ji , W lk )

In the general case, we can use the results from section A.2 to evaluate off-diagonal h j h l . However, in our experiments we always consider the the special case where Cov(W ji , W lk ) is diagonal. In this case we can write the first term in equation 15 as (reintroducing the explicit summation):

jl h j h l Cov(W ji , W lk ) = jl h j h l δ jl δ ik Var(W ji ) =δ ik j z j z j Var(W ji ) = diag [vVar(W )]

i.e. this term is a diagonal matrix with the diagonal given by the left product of the vector v j = h j h j with the matrix Var(W ki ). Note that h j h j can be evaluated analytically for Heaviside and ReLU activation functions:

Heaviside h j h j = 1 √ 2πΣ jj ∞ 0 e - 1 2Σ jj (αj-a j ) 2 dα j = Φ µ j ReLU h j h j = 1 √ 2πΣ jj ∞ 0 α 2 j e - 1 2Σ jj (αj-a j ) 2 dα j = Σ jj µ j φ(µ j ) + (1 + µ 2 j )Φ(µ j ) A.2 EVALUATION OF TERM 2: W ji Cov(h j , h l ) W lk

Evaluation of Cov(h j , h l ) requires an expression for h j h l . From equation 5, we write:

EQUATION

)

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:

h j h l = S jl I(µ j , µ l , ρ jl ) ; I = 1 Z f (η j )f (η l ) exp -1 2 P (η j , η l ; µ , ρ ) dη j dη l .

The normalization constant, Z, is evaluated by integrating over e -P/2 and is explicitly written as Z = 2π ρ jl , where ρ jl = 1ρ 2 jl . 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.

A.2.1 HEAVISIDE NON-LINEARITY:

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 -1 2 P (η j , η l ; µ , ρ jl ) dη j dη 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 ρ 2 jl η 2 j -2ρ jl η j η l + η 2 l = 1 4 ρ 2 jl ξ 2 + + ξ 2 -,

where

ξ ± = √ 1 ∓ ρ(η 1 ± η 2 ).

Performing this change of variables in the integral gives:

I = 1 4π ρ 2 jl H exp 1 4 ρ 2 jl ξ 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 ξ as 7 :

ψ = arctan 1+ρ jl 1-ρ jl = π 2 -1 2 arccos ρ jl .

Since A| µ =0 = Φ(0)Φ(0) = 1/4, we can evaluate:

(I -A)| µ =0 = π 2 -1 2 arccos ρ jl -1 4 = 1 2π arcsin ρ jl 7

Here we use the identity cos(2 arctan x) = cos 2 arctan x -

sin 2 arctan x = 1-x 2 1+x 2

First derivative Performing a change of variables x i = η iµ i , we can write I as:

I = 1 Z H(x j + µ j )H(x l + µ l ) exp -1 2 P (x j , x l ; 0, ρ jl ) dx j dx l

where H is the Heaviside function. Now, using ∂ x H(x) = δ(x), we have:

EQUATION

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 ∂µ 2 j µ =0 I = - 1 Z δ(x j )H(x l ) ρ jl ρ 2 jl x l exp -1 2 P (x j , x l ; 0, ρ jl ) dx j dx l = - ρ jl 2π ρ jl .

where we used the identity f (x)∂ x δ(x)dx =δ(x)∂ x f (x)dx, which holds for arbitrary f . In addition, we have:

∂ 2 ∂µ 2 j µ =0 A = 0 =⇒ ∂ 2 ∂µ 2 j µ =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π + β µ 2 j + µ 2 l + γµ 2 j µ 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 ∂µ 2 i 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.

A.2.2 RELU NON-LINEARITY:

As in the Heaviside case, we begin by computing the asymptote of I by inspecting the limit as µ j → ∞:

EQUATION

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.

B LOG-LIKELIHOOD AND POSTERIOR PREDICTIVE COMPUTATION:

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(a L ). 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: We begin by rewriting the reconstruction term for data point (x, y) in terms of a L :

E w∼q [log p(y|w)] = q(w) log p(y|w) dw = q(a L )q(w|a L ) log p(y|w) dw da L 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 a L only, thus allowing it to be evaluated using the approximation q(a L ) from section 3.1.

To eliminate w, consider the case where the output of the model is a distribution p(y|a L ) that is a parameter-free transformation of a L (e.g. a L are logits of a softmax distribution for classification or the moments of a Gaussian for regression). Since the model output is conditioned only on a L , we must have p(y|w) = p(y|a L ) for all configurations w that satisfy the deterministic transformation a L = M(x; w), where M is the neural network (i.e p(y|w) = p(y|a L ) for all w where q(w|a L ) is non-zero). This allows us to write:

q(w|a L ) log p(y|x, w) dw = log p(y|a L ) q(w|a L ) dw = log p(y|a L ), so the reconstruction term becomes:

E w∼q [log p(y|x, w)] = q(a L ) log p(y|a L )da L = E a L ∼q(a L ) log p(y|a L ) .

This establishes the equivalence given in equation 7 in the main text. Since we are using an approximation to q, we will actually compute E a L ∼q(a L ) log p(y|a L ) .

B.2 UNIVARIATE REGRESSION: LOG-LIKELIHOOD:

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 a L = (m, ). We first insert the Gaussian form for p(y|a L ) ∼ N m, e into the log-likelihood expression:

EQUATION

)

Now we use the Gaussian form of q(a L ):

q(a L ) ∝ exp -1 2 X (Σ L ) -1 X ; X = m-m - .

and note that

EQUATION

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:

E a L ∼q(a L ) log p(y|a L ) = -1 2 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.

B.3 UNIVARIATE REGRESSION: POSTERIOR PREDICTIVE DISTRIBUTION:

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(a L ) ≈ q(a L ):

E y∼p(y) [y] = yp(y|a L )q(a L ) dy da L = yp(y|a L ) dy q(a L ) da L = mq(a L ) da L = m Var[y] = Var[E y∼p(y|a L ) (y)] + E y∼p(y) Var(y|a L ) = Var[m] + E y∼p(y) e

= Σ mm + e p(y|a L )q(a L ) da L dy = Σ mm + e q(a L ) da L = Σ mm + exp ( + Σ /2) where the final integral in the variance computation is performed by inserting the Gaussian form for q(a L ) and completing the square. Then we assume normality of the predictive distribution to obtain the result in equation 9.

B.4 CLASSIFICATION: LOG-LIKELIHOOD:

There is no exact form for the expected log-likelihood for multivariate classification with logits a L . However, using the second-order Delta method (Small, 2010), we find the expansion

EQUATION

) 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 :

EQUATION

where C ij = Cov(x i , x j ). Now we note that the logsumexp function has a simple Hessian ∂ 2 ∂xi∂xj logsumexp(x) = δ ij p ip i p j , where p = softmax(x). Putting these results together allows us to write:

E [logsumexp(x)] ≈ logsumexp (E[x]) + 1 2 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.

B.5 CLASSIFICATION: POSTERIOR PREDICTIVE DISTRIBUTION:

Using the same second-order Delta method, we find the following expansion for the posterior predictive distribution:

EQUATION

)

where p = softmax( a L ).

For this expansion, we begin by computing the Hessian:

∇ i ∇ j [p] k = ∂ 2 ∂xi∂xj p k = [2p i p j -(δ ki p j + δ kj p i ) + δ ik δ jk -δ ij p i ] p k ,

where p = softmax(x), and we used the intermediate result ∇ j [p] k = δ jk p kp j p k . 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 a L to (log)p after the deterministic heavy-lifting of computing a L . 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.

C DEEPER NETWORKS:

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.

C.1 SKIP CONNECTIONS:

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):

a i = a i + δ i , Cov(a i , a k ) = 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 Published as a conference paper at ICLR 2019 results and methods borrowed from appendix A.

Cov(a:

i , δ k ) = Cov(a i , j f (a j )W jk ) + Cov(a i , b k ) = j a i f (a j ) -a i f (a j ) W jk = j W jk (Σ 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. 

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

F LEARNING CURVES:

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.

:

We are also required to choose a Gaussian variational approximation for b to preserve the Gaussian distribution of a.

Additional Taylor expansion terms can be computed if this assumption fails.

In general, any arbitrary partitioning can be used

Our implementation in TensorFlow is available at https://github.com/Microsoft/ deterministic-variational-inference

For MCVI with full rank covariance, Cholesky decomposition required for sampling is O(N 3 ).

Note that this result is non-trivial because heteroscedastic models are more complex and could result in poorer approximate inference leading to worse test performance

This result is obtained by Taylor expansion inside the expectation.

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

Figures and Tables:

Figure1shows the architecture of the computation of E w∼q [log p(D|w)] for a feed-forward neural network. The computation can be divided into two parts: first, propagation of activations though

Figure 2: Approximation of hjh l using an asymptote and Gaussian correction for (a) Heaviside and (b) ReLU non-linearities. Yellow functions have closed-forms, and blue indicates residuals. The examples are plotted for -6 < µ < 6 and ρ jl = 0.5, and the relative magnitude of each correction term is indicated on the vertical axis.

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.

Figure 4: Runtime performance of VI methods.

Figure 5: Comparison of converged test log-likelihood with a manually tuned prior variance (orange) or empirical Bayes (blue).

Figure 6: Empirical accuracy of our approximation for 5-layer networks trained analogously to figure 3. Progressively narrower networks of (a) 125 unit (b) 25 unit and (c) 5 unit are trained and our CLT-based approximation is only seen to significantly break down in the 5-unit case. (d) Qualitative verification of our approximation applied to an architecture with skip connections (orange).

Figure 8: Learning trajectories for the models from table 2.

Average test log-likelihood on UCI datasets. |D| is the dataset size, and dx is the input dimension.Deterministic vs. StochasticWe 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

Table content: <table><tr><td>Dataset</td><td>|D|</td><td>dx</td><td>DVI</td><td>dDVI</td><td>MCVI</td><td>hoDVI</td></tr><tr><td>bost</td><td>506</td><td>13</td><td>-2.41 ± 0.02</td><td>-2.42 ± 0.02</td><td>-2.46 ± 0.02</td><td>-2.58 ± 0.04</td></tr><tr><td>conc</td><td>1030</td><td>8</td><td>-3.06 ± 0.01</td><td>-3.07 ± 0.02</td><td>-3.07 ± 0.01</td><td>-3.23 ± 0.01</td></tr><tr><td>ener</td><td>768</td><td>8</td><td>-1.01 ± 0.06</td><td>-1.06 ± 0.06</td><td>-1.03 ± 0.04</td><td>-2.09 ± 0.06</td></tr><tr><td>kin8</td><td>8192</td><td>8</td><td>1.13 ± 0.00</td><td>1.13 ± 0.00</td><td>1.14 ± 0.00</td><td>1.01 ± 0.01</td></tr><tr><td>nava</td><td>11934</td><td>16</td><td>6.29 ± 0.04</td><td>6.22 ± 0.06</td><td>5.94 ± 0.05</td><td>5.84 ± 0.06</td></tr><tr><td>powe</td><td>9568</td><td>4</td><td>-2.80 ± 0.00</td><td>-2.80 ± 0.00</td><td>-2.80 ± 0.00</td><td>-2.82 ± 0.00</td></tr><tr><td>prot</td><td>45730</td><td>9</td><td>-2.85 ± 0.01</td><td>-2.84 ± 0.01</td><td>-2.87 ± 0.01</td><td>-2.94 ± 0.00</td></tr><tr><td>wine</td><td>1588</td><td>11</td><td>-0.90 ± 0.01</td><td>-0.91 ± 0.02</td><td>-0.92 ± 0.01</td><td>-0.96 ± 0.01</td></tr><tr><td>yach</td><td>308</td><td>6</td><td>-0.47 ± 0.03</td><td>-0.47 ± 0.03</td><td>-0.68 ± 0.03</td><td>-1.41 ± 0.03</td></tr></table>

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

Table content: <table><tr><td>Dataset</td><td>bost</td><td>conc</td><td>ener</td><td>kin8</td><td>nava</td><td>powe</td><td>prot</td><td>wine</td><td>yach</td></tr><tr><td>|D|</td><td>506</td><td>1030</td><td>768</td><td>8192</td><td>11934</td><td>9568</td><td>45730</td><td>1588</td><td>308</td></tr><tr><td>dx</td><td>13</td><td>8</td><td>8</td><td>8</td><td>16</td><td>4</td><td>9</td><td>11</td><td>6</td></tr></table>

Average test log likelihood on UCI datasets. See table 4 for model glossary and implementation references.

Table content: <table/>

Peer review (iclr2019):
 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. (1) A deep BNN to show that the cumulative error is negligible as the number of the hidden layers increases (2) Small latent dimension since CLT may not hold (3) A heavy-tailed variational distribution since the second moment may not be finite (4) Other nonlinear activations since the Gaussian approximation may not be accurate due to (generalized) Berry-Esseen theorem (5) A BNN with skip connections since a Bayesian multiplayer perceptron with skip connections is also a feed-forward BNN Among these cases, I am eager to see some results on a deep thin BNN. For example, a BNN with 5 hidden layers, where the latent dimension at each layer is less than 32. Furthermore, I would like to see some empirical comparison on real-world datasets between DVI and MCVI under a *fixed* prior since such comparison demonstrates the approximation accuracy of DVI and rule out the confounding factor introduced by the empirical Bayes approach.