1 Introduction
Bayesian neural networks represent a very flexible class of models, which are quite robust to overfitting. However they often remain heavily overparametrized. There are several implicit approaches for BNNs sparsification through shrinkage of weights (Blundell et al., 2015; Molchanov et al., 2017; Ghosh et al., 2018; Neklyudov et al., 2017). For example, Blundell et al. (2015) suggest a mixtureofGaussians spike and slab prior on the weights and then perform fully factorisable meanfield variational approximation. Ghosh & DoshiVelez (2017); Louizos et al. (2017) independently generalize this approach by means of suggesting Horseshoe priors (Carvalho et al., 2009)
for the weights, providing even stronger shrinkage and automatic specification of the mixture component variances required in
Blundell et al. (2015). Molchanov et al. (2017) suggest that the interpretation of Gaussian dropout as performing variational approximations for a Bayesian neural network with log uniform prior over the weight parameters leads to improved sparsity in the latter. Some of these approaches show that the majority of the parameters can be pruned out from the network without significant loss of predictive accuracy. At the same time pruning is done in an implicit manner by deleting the weights via adhoc thresholding.In Bayesian model selection problems there have been numerous works showing efficiency and accuracy of model selection by means of introducing latent variables corresponding to different discrete model configurations, and then conditioning on their marginal posterior to both select the best sparse configuration and address the joint modelandparametersuncertainty explicitly (George & McCulloch, 1993; Clyde et al., 2011; Frommlet et al., 2012; Hubin & Storvik, 2018; Hubin et al., 2018b, a). For instance, Hubin et al. (2018a) address inference in the class of deep Bayesian regression models (DBRM), which generalizes the class of Bayesian neural networks. They show both good predictive performance of the obtained sparse models and the ability to recover meaningful complex nonlinearities. The approach suggested in Hubin et al. (2018a)
is based on adaptations of Markov chain Monte Carlo (MCMC) and does not scale well to large highdimensional data samples.
Louizos et al. (2017) also warn about the complexity of explicit discretization of model configuration within BNNs, as it causes an exponential explosion with respect to the total number of parameters, and hence infeasibility of inference for highdimensional problems. At the same time, Polson & Rockova (2018) encourage the use of spike and slab approach in BNNs from a theoretical standpoint.Logsdon et al. (2010); Carbonetto et al. (2012) suggested a fullyfactorized variational distribution capable of efficiently and precisely "linearizing" the problem of Bayesian model selection in the context of linear models with an ultrahigh number of potential covariates, typical for genome wide association studies (GWAS). In the discussion to his PhD thesis Hubin (2018) proposed combining the approaches of Logsdon et al. (2010); Carbonetto et al. (2012) and Graves (2011) for scalable approximate Bayesian inference on the joint space of models and parameters in deep Bayesian regression models. We develop this idea further in this article.
More specifically, we introduce a formal Bayesian approach for jointly taking into account model (structural) uncertainty and parameter uncertainty
in BNNs. The approach is based on introducing latent binary variables corresponding to inclusionexclusion of particular weights within a given architecture. Using Bayesian formalization in the space of models allows to adapt the whole machinery of Bayesian inference in the joint modelparameter settings, including
Bayesian model averaging (BMA) (across all models) or Bayesian model selection (BMS) of one “best” model with respect to some model selection criterion (Claeskens et al., 2008). In this paper we study BMA as well as the median probability model (Barbieri et al., 2004, 2018) and posterior mean model based inference for BNNs. Spasifying properties of BMA and the median probability model are also addressed. Finally, following Hubin (2018) we will link the obtained marginal inclusion probabilities to binary dropout rates, which gives proper probabilistic reasoning for the latter. The suggested inference approach is based on scalable stochastic variational inference.The approach has similarities to binary dropout that has become very popular (Srivastava et al., 2014)
. However, while standard binary dropout can only be seen as a Bayesian approximation to a Gaussian process model where only parameter estimation is taken into account
(Gal, 2016), our approach also explicitly models structural uncertainty. In this sense it is closely related to Concrete dropout (Gal et al., 2017). However, the model proposed by Gal et al. (2017) does not allow for BMS: The median probability model will either select all weights or nothing due to a strong assumption on having the same dropout probabilities for the whole layer. Furthermore, it uses variational approximations, which were not studied in the model uncertainty context.The rest of the paper is organized as follows: The class of BNNs and the corresponding model space are mathematically defined in Section 2. In Section 3 we describe the algorithm for making inference on the suggested class of models using reparametrization of marginal inclusion probabilities. In Section 4 the suggested approach is applied to a classical benchmark data set MNIST and also to a data set FMNIST, there we also compare the results with some of the existing approaches for inference on BNNs. Finally, in Section 5 some conclusions and suggestions for further research are given. Additional results, discussions, and proofs can be found in the web supplement to the paper.
2 The model
A neural network model links observations and explanatory variables via a probabilistic functional mapping of the form (1), with a possibly multidimensional mean parameter , which can be written as:
(1) 
where
is a dispersion parameter. To construct the vector of mean parameters
of the distribution of interest, one builds a sequence of building blocks of hidden layers through semiaffine transformations:(2) 
with . Here is the number of layers, is the number of nodes within the layer while are univariate functions (further referred to as activation functions). Further, are the weights (slope coefficients) for the inputs of the th layer (note that ; for we obtain the intercepts or the bias terms). Finally, are latent binary indicators switching the corresponding weights on and off.
In our notation we explicitly differentiate between discrete model configurations defined by the vectors (further referred to as models) constituting the model space and parameters of the models, conditional on these configurations , where only those for which are included. This approach is a rather standard (in statistical science literature) way to explicitly specify the model uncertainty in a given class of models and is used in Clyde et al. (2011); Frommlet et al. (2012); Hubin & Storvik (2018); Hubin et al. (2018b, a). A Bayesian approach is obtained by specification of model priors and parameter priors for each model . If the dispersion parameter is present in the distribution of the outcomes, one also has to define . Many kinds of priors on can be considered, including the mixture of Gaussians prior (Blundell et al., 2015), Horseshoe prior (Ghosh & DoshiVelez, 2017; Louizos et al., 2017), or mixtures of gpriors (Li & Clyde, 2018), but for simplicity of representation we consider the multivariate Gaussian with the diagonal covariance matrix weight prior combined with independent Bernoulli priors for the latent inclusion indicators (Clyde et al., 2011):
(3) 
The parameter is the penalty for including the weight into the model. To induce AICtype model complexity penalization, one typically uses , whilst for BICtype penalization one uses , where is the sample size.
3 Bayesian inference
The main goal of inference with uncertainty in both models and parameters is to infer the posterior marginal distribution of some parameter of interest (for example probability of a new observation conditional on new covariates ) based on data :
(4) 
The posterior joint distribution of models and parameters
will be approximated by means of scalable variational inference (Graves, 2011) for the joint parametermodel setting. Consider the following variational families, parameterized with , with independence across weight components:(5) 
where the factorized components follow the form suggested (in a simpler setting) by Logsdon et al. (2010); Carbonetto et al. (2012):
(6) 
Here is the delta mass or "spike" at zero and . Thus, with probability , the posterior of parameters of a weight
will be approximated by a normal distribution with some mean and variance ("slab"), and otherwise the edge is considered to have no effect on the observations
. Hence will approximate the marginal inclusion probability of the weight . Note that while the variational distribution is similar to the dropout approach, the target distribution, which we aim at approximating, is different in the sense of including the binary variables as well. Hence our marginal inclusion probabilities can serve as a particular case of dropout rates with a proper probabilistic interpretation in terms of structural model uncertainty.To perform variational approximation of the posterior we aim at minimizing the KullbackLeibler divergence between the variational family distribution and the posterior distribution.
(7) 
with respect to the variational parameters . Constraints on are incorporated by means of the reparametrization
For numerical stability we also reparameterize the standard deviations of the weights
. In Proposition 1 we show minimization of the divergence (7) to be equivalent to maximization of the evidence lower bound .Proposition 1.
Minimization of is equivalent to maximization of the evidence lower bound (ELBO)
Proof of Proposition 1 is given in Section A of the web supplement.
Following Blundell et al. (2015) we approximate the objective
with an unbiased estimate based on minibatching. Note that the observations are assumed conditionally independent, which means
. We will hence sample minibatches of size from the full data, giving the unbiased estimateHowever, the cardinality of is , where is the total number of potential weights in the neural network, hence iterating through all of the models in
for every iteration of stochastic gradient descent optimization is infeasible for reasonably complex Bayesian neural networks. It is also infeasible to integrate over all of the weights for a given configuration. To resolve this, we adopt a yet another unbiased MonteCarlo based approximation:
(8) 
where for .
Stochastic gradient based methods for optimization require calculation of an unbiased estimate of the gradient (Bottou et al., 2018). One can further use sampling to obtain such an estimate (Graves, 2011; Blundell et al., 2015). Some extra care is needed in this case, due to the discrete variables. In Proposition 2 we suggest a generalization of Blundell et al. (2015) and Mnih & Gregor (2014), allowing to make inference in the joint space of models and parameters. Remark 2.1 allows to use control variates (Weaver & Tao, 2001) to reduce the variance of (Mnih & Gregor, 2014).
Proposition 2.
Assume for and is a random subset of of size N. Then an unbiased estimator for the gradient of is given by:
(9) 
Remark 2.1.
Without inducing bias, a term can be added to .
Algorithm 1 describes one iteration of a doubly stochastic variational inference approach where updating is performed on the parameters in the transformed space. The set is the collection of all combinations
in the network. Note that in the suggested algorithm partial derivatives with respect to marginal inclusion probabilities, as well as mean and standard deviation terms of the weights are shared and coincide with the gradients found by the usual backpropagation algorithm on a neural network. This algorithm assumes no dispersion parameter
, but can be easily generalized to include it. Once the estimates of the parameters of the variational approximating distribution are obtained, there are several ways to proceed with the inference. Algorithm 2 describes a procedure for inference onbased on the marginal posterior predictive distribution (
4), taking uncertainty in both the model structure and the parameters into account. sample indices uniformly from defining ; for in do for do set , ; sample ; sample ; end for end for set with from (9); foreach prediction of interest do for in do for do sample ; sample ; end for calculate ; end for set . end foreachA bottleneck of Algorithm 2 is that we have to both sample from a huge approximate posterior distribution of parameters and models and keep all of the components of stored during the inference, which might be computationally inefficient. The posterior mean based model (Wasserman, 2000), which sets and or the median probability model (Barbieri et al., 2004), which sets , can be used as simplifying alternatives. Both approaches specify one model but the former model results in a dense model while the latter results in a sparse one.
Given a specific model there are several options with respect to : Simulation conditional on as in Algorithm 2 or inserting the posterior means conditional on . A main benefit with the latter is that a single set is predefined so that calculation of only one term needs to be done in the prediction step. In both of these cases, one can either use directly the learned variational approximation or perform posttraining of the parameters. By posttraining we mean that only distributions of the parameters of the models (means and variances of the weights and dispersion parameter of the distribution) are optimized, whilst the distributions of the model configurations are fixed. This might improve inferential properties of the posterior predictive distribution and reduce its variance. We provide a comparison of several alternatives in Section 4 and give a further discussion of the simplifying inference techniques in Section B of the web supplement.
4 Applications
In this section we will address the classification of MNIST (LeCun et al., 1998) images and classification of fashionMNIST (FMNIST Xiao et al., 2017) images. Both of the datasets comprise of 28x28 grayscale images of items from 10 categories, with images per category. The training sets consist of images and the test sets have images.
Experimental design
For both of the data sets we addressed a dense neural network with the ReLU activation function, multinomially distributed observations with 10 classes and 784 input explanatory variables (pixels). The network has 3 hidden layers with 400, 600 and 600 neurons correspondingly. Priors for the parameters and model indicators were chosen according to (
3) with and. The inference was performed by means of the suggested doubly stochastic variational inference approach on 250 epochs with batch size 100. We used the ADAM stochastic gradient descent optimization
(Kingma & Ba, 2014) with (allowing the parameter in the last step of Algorithm 1 to have separate values for the different sets of parameters). When posttraining the parameters, either with fixed marginal inclusion probabilities or with the median probability model, we ran additional epochs of the optimization routine with , and . In Algorithm 2 we used both and .We then evaluated accuracies (Acc
 the proportion of the correctly classified images) for each vector of predictions and recorded their median. Accuracies based on a single sample (
) from the median probability model and the model obtained by inserting the posterior mean of the parameters were also obtained. Finally, accuracies based on posttraining of the parameters with fixed marginal inclusion probabilities and posttraining of the median probability model were evaluated. For the cases when model averaging is addressed () we are additionally reporting accuracies when classification is only performed when the maximum model averaged class probability exceeds 95% as suggested by Posch et al. (2019). Otherwise, a doubt decision is made (Ripley, 2007, sec 2.1). In this case we both report the accuracy within the classified images as well as the number of classified images. Finally, in column Density of Table 1 we are reporting the median of the overall density level (proportion of used at least once weights in the prediction stage to the total number of weights in a BNN), for different approaches. To guarantee reproducibility, summaries across 10 independent runs of the described experiment were computed for all of these statistics. Summaries of these statistics are reported in Table 1. Estimates of the marginal inclusion probabilities based on the suggested variational approximations were also computed for all of the weights. In order to compress the presentation of the results we only present the mean marginal inclusion probabilities for all of the layers as , summarized in Table 2.In addition to our approach (denoted as Full BNN, Gaussian priors in Table 1) we also used several baselines. In particular, we addressed a standard Dense BNN with Gaussian priors (Graves, 2011), which is a particular case of our original model with all being fixed and equal to 1. This model is important in measuring how much of the predictive power we might loose due to introducing sparsity. Furthermore, we report the results for a Dense BNN with mixture priors with two Gaussian components of the mixtures (Blundell et al., 2015) with probabilities 0.5 for each and variances equal to 1 and correspondingly. Additionally we have addressed two popular sparsity inducing approaches, in particular, BNN with Concrete dropout (Gal et al., 2017) and BNN with Horseshoe priors (Louizos et al., 2017). All of the baseline methods have, similarly to the full BNN, 3 hidden layers with 400, 600 and 600 neurons correspondingly. They were trained for 250 epochs with Adam optimizer () and batch size equal to 100. For the BNN with Horseshoe priors we are reporting statistics separately before and after adhoc pruning of the weights. Posttraining (when necessary) was performed for additional 50 epochs.
MNIST  FMNIST  
Method  All cl  0.95 threshold  Density  All cl  0.95 threshold  Density  
PT  Acc  Acc  Num.cl  level  Acc  Acc  Num.cl  level  
Full BNN, Gaussian priors  
SIM  SIM  N  1  0.958 (0.954,0.960)      0.056  0.854 (0.850,0.858)      0.066 
SIM  SIM  Y  1  0.971 (0.969,0.973)      0.056  0.868 (0.863,0.872)      0.066 
SIM  SIM  N  10  0.967 (0.966,0.971)  0.999  7064  0.084  0.867 (0.863,0.870)  0.996  4097  0.083 
SIM  SIM  Y  10  0.978 (0.976,0.980)  0.999  8366  0.084  0.880 (0.875,0.882)  0.994  4933  0.083 
ALL  EXP  N  1  0.969 (0.967,0.970)      1.000  0.866 (0.864,0.874)      1.000 
ALL  EXP  Y  1  0.979 (0.978,0.980)      1.000  0.880 (0.877,0.884)      1.000 
MED  SIM  N  1  0.961 (0.957,0.964)      0.051  0.858 (0.854,0.865)      0.065 
MED  SIM  Y  1  0.973 (0.971,0.977)      0.051  0.872 (0.870,0.875)      0.065 
MED  SIM  N  10  0.964 (0.962,0.967)  0.998  7441  0.051  0.863 (0.859,0.869)  0.993  4347  0.065 
MED  SIM  Y  10  0.977 (0.976,0.979)  0.999  8645  0.051  0.878 (0.876,0.881)  0.992  5223  0.065 
MED  EXP  N  1  0.965 (0.963,0.968)      0.051  0.863 (0.859,0.870)      0.065 
MED  EXP  Y  1  0.978 (0.976,0.979)      0.051  0.879 (0.876,0.882)      0.065 
Dense BNN, Gaussian priors  
ALL  SIM  N  1  0.965 (0.965,0.966)      1.000  0.864 (0.863,0.866)      1.000 
ALL  SIM  N  10  0.984 (0.982,0.985)  0.999  8477  1.000  0.893 (0.890,0.894)  0.997  5089  1.000 
ALL  EXP  N  1  0.984 (0.982,0.985)      1.000  0.886 (0.882,0.888)      1.000 
Dense BNN, mixture priors  
ALL  SIM  N  1  0.965 (0.964,0.967)      1.000  0.867 (0.866,0.868)      1.000 
ALL  SIM  N  10  0.982 (0.981,0.983)  0.999  8329  1.000  0.893 (0.892,0.897)  0.996  5151  1.000 
ALL  EXP  N  1  0.983 (0.981,0.984)      1.000  0.888 (0.885,0.890)      1.000 
BNN, Concrete dropout  
SIM  SIM  N  1  0.982 (0.894,0.984)      0.226  0.896 (0.820,0.902)      0.094 
SIM  SIM  N  10  0.984 (0.896,0.986)  0.995  9581  0.820  0.897 (0.823,0.901)  0.942  8825  0.447 
SIM  SIM  Y  1  0.982 (0.894,0.984)      0.226  0.897 (0.820,0.899)      0.094 
SIM  SIM  Y  10  0.984 (0.896,0.986)  0.995  9586  0.820  0.897 (0.823,0.902)  0.943  8826  0.447 
ALL  EXP  N  1  0.983 (0.896,0.984)      1.000  0.896 (0.821,0.901)      1.000 
ALL  EXP  Y  1  0.983 (0.894,0.984)      1.000  0.896 (0.820,0.901)      1.000 
BNN, horseshoe priors  
SIM  SIM  N  1  0.964 (0.962,0.967)      1.000  0.864 (0.863,0.869)      1.000 
SIM  SIM  N  10  0.982 (0.981,0.983)  1.000  0003  1.000  0.887 (0.886,0.889)  1.000  0181  1.000 
ALL  EXP  N  1  0.966 (0.963,0.968)      1.000  0.867 (0.861,0.868)      1.000 
PRN  SIM  N  1  0.965 (0.962,0.969)      0.194  0.865 (0.860,0.868)      0.302 
PRN  SIM  N  10  0.982 (0.981,0.983)  1.000  0002  0.194  0.887 (0.884,0.888)  1.000  0179  0.302 
PRN  EXP  N  1  0.965 (0.963,0.968)      0.194  0.865 (0.862,0.869)      0.302 
PRN  SIM  Y  1  0.967 (0.965,0.968)      0.194  0.867 (0.864,0.871)      0.302 
PRN  SIM  Y  10  0.982 (0.981,0.983)  1.000  0007  0.194  0.888 (0.887,0.890)  1.000  0147  0.302 
PRN  EXP  Y  1  0.966 (0.964,0.969)      0.194  0.868 (0.864,0.869)      0.302 
MNIST data  FMNIST data  

Layer  Med.  SD.  Med.  SD. 
0.0520  0.0005  0.0665  0.0004  
0.0598  0.0003  0.0613  0.0005  
0.2217  0.0064  0.2013  0.0051 
Mnist
The results reported in Tables 1 and 2 show that within the suggested fully Bayesian approach: a) model averaging across different BNNs () gives significantly higher accuracy than the accuracy of a random individual BNN from the model space (), additional accuracy improvement can be achieved by using more samples from the joint posterior for model averaging as shown in Figure 1
; b) the median probability model and posterior mean based model also perform significantly better than a randomly sampled model, they perform somewhat worse than model averaging, but equally well to model averaging when posttraining is performed; c) the majority of the weights of the models have very low marginal inclusion probabilities for the weights at layers 1 and 2, and significantly more weights have higher marginal inclusion probabilities at layer 3, resembling the structure of convolutional neural networks (CNN) where typically one first has a set of sparse convolutional layers, followed by a few fully connected layers, unlike CNNs the structure of sparsification is learned automatically within our approach (see also Section
C.1 in the web supplement for further details on this issue); d) variations of all of the performance metrics across simulations are low, showing stable behavior across the repeated experiments; e) inference with a doubt option gives almost perfect accuracy, however this comes at a price of rejecting to classify some of the items.For other approaches it is also the case that f) both posterior mean based model and using sample averaging improves accuracy compared to a single sample from the parameter space; g) variations of the target parameters are low for the dense BNNs with Gaussian/mixture of Gaussians priors and BNN with horseshoe priors and rather high for the Concrete dropout approach. When it comes to comparing our approach to baselines we notice that h) dense approaches outperform sparse approaches in terms of the accuracy in general; i) Concrete dropout marginally outperforms other sparse approaches in terms of median accuracy, however it exhibits large variance, whilst our full BNN and the compressed BNN with horseshoe priors yield equivalent performance across experiments; j) neither our approach nor baselines managed to reach state of the art results in terms of hard classification accuracy of predictions (Palvanov & Im Cho, 2018); k) including a 95% threshold for making a classification results in a very low number of classified cases for the horseshoe priors (it is extremely underconfident), the Concrete dropout approach seems to be overconfident when doing inference with the doubt option (resulting in lower accuracy but larger number of decisions), the full BNN, and BNN with Gaussian and mixture of Gaussian priors give less classified cases than the Concrete dropout approach but reach significantly higher accuracy; l) this might mean that the thresholds need to be calibrated towards the specific methods; m) our approach yields the highest sparsity of weights, both for model averaging and for the median probability model.
Fmnist
The same set of approaches, model specifications and tuning parameters of the algorithms as in the MNIST example were used for this application. The results a) m) for FMNIST data, based on Tables 1, 2 and Figure 1 are completely consistent with the results from the MNIST experiment, however predictive performance for all of the approaches is poorer on FMNIST. Also whilst full BNN and BNN with horseshoe priors on FMNIST obtain lower sparsity levels than on MNIST, Concrete dropout here improves in this sense compared to the previous example.
Outofsample experiments
Figure 2 shows the results on outofsample experiments using FMNIST data on a BNN trained by MNIST data and vice versa. Following Louizos & Welling (2017)
, the goal now is to obtain as inconclusive results as possible (reaching ideally a uniform distribution across classes), corresponding to a large entropy. The plot shows the empirical cumulative distribution function (CDF) of the entropies over the classified samples. Concrete dropout is over confident with a distribution of test classes being far from uniform, the horseshoe prior based approach is the closest to uniform (but it was also closer to uniform for the indomain predictions), whilst the 2 other baselines are in between, our approach is on par with them for the first case (left graph) and even becomes close to the horseshoe prior based approach for the second case (right graph) showing that it handles outofsample uncertainty rather well (just as it does handle the indomain uncertainty). See also Section
C.3 in the web supplement for further details on this issue.5 Discussion
In this paper we have introduced the concept of Bayesian model (or structural) uncertainty in BNNs and suggested a scalable variational inference technique for fitting the approximation to the joint posterior of models and parameters of these models. Posterior predictive distributions, with both models and parameters marginalized out, can be easily obtained. Furthermore, marginal inclusion probabilities suggested in our approach give proper probabilistic interpretation to Bayesian binary dropout and allow to perform model selection. This comes at the price of having just one additional parameter per weight included. We further provide image classification applications of the suggested technique showing that it both allows to significantly sparsify neural networks without noticeable loss of predictive power and accurately handle the predictive uncertainty.
Currently, fairly simple prior distributions for both models and parameters are used. These prior distributions are assumed independent across the parameters of the neural network, which might not always be reasonable. Alternatively, both parameter and model priors can incorporate joint dependent structures, which can further improve the sparsification of the configurations of neural networks. When it comes to the model priors with local structures and dependencies between the variables (neurons), one can mention the so called dilution priors (George et al., 2010). These priors take care of the similarities between models by means of downweighting the probabilities of the models with highly correlated variables. There are also numerous approaches to incorporate interdependencies between the model parameters via priors in different settings (Smith & LeSage, 2004; Fahrmeir & Lang, 2001; Dobra et al., 2004). Obviously, in the context of inference in the joint parametermodel settings in BNNs, more research should be done on the choice of priors. Specifically, for image analysis, it might be of interest to develop convolutioninducing priors, whilst for recurrent models one can think of exponentially decaying parameter priors for controlling the shortlong memory.
In this work we restrict ourselves to a subclass of BNNs, defined by inclusionexclusion of particular weights within a given architecture. In the future it can be of particular interest to extend the approach to the choice of the activation functions as well as the maximal depth and width of each layer of the BNN. A more detailed discussion of these possibilities and ways to proceed is given in Hubin (2018). Finally, studies of the accuracy of variational inference within these complex nonlinear models should be performed. Even within linear models Carbonetto et al. (2012) have shown that the results can be strongly biased. Various approaches for reducing the bias in variational inference are developed. One can either use more flexible families of variational distributions by for example introducing auxiliary variables (Ranganath et al., 2016; Salimans et al., 2015) or address Jackknife to remove the bias (Nowozin, 2018).
References
 (1)
 Barbieri et al. (2018) Barbieri, M., Berger, J. O., George, E. I. & Rockova, V. (2018), ‘The median probability model and correlated variables’, arXiv preprint arXiv:1807.08336 .
 Barbieri et al. (2004) Barbieri, M. M., Berger, J. O. et al. (2004), ‘Optimal predictive model selection’, The annals of statistics 32(3), 870–897.
 Blundell et al. (2015) Blundell, C., Cornebise, J., Kavukcuoglu, K. & Wierstra, D. (2015), ‘Weight uncertainty in neural networks’, arXiv preprint arXiv:1505.05424 .

Bottou et al. (2018)
Bottou, L., Curtis, F. E. & Nocedal, J. (2018), ‘Optimization methods for largescale machine learning’,
SIAM Review 60(2), 223–311.  Carbonetto et al. (2012) Carbonetto, P., Stephens, M. et al. (2012), ‘Scalable variational inference for Bayesian variable selection in regression, and its accuracy in genetic association studies’, Bayesian analysis 7(1), 73–108.

Carvalho et al. (2009)
Carvalho, C. M., Polson, N. G. & Scott, J. G. (2009), Handling sparsity via the horseshoe, in
‘Artificial Intelligence and Statistics’, pp. 73–80.
 Claeskens et al. (2008) Claeskens, G., Hjort, N. L. et al. (2008), ‘Model selection and model averaging’, Cambridge Books .
 Clyde et al. (2011) Clyde, M. A., Ghosh, J. & Littman, M. L. (2011), ‘Bayesian adaptive sampling for variable selection and model averaging’, Journal of Computational and Graphical Statistics 20(1), 80–101.

Dobra et al. (2004)
Dobra, A., Hans, C., Jones, B., Nevins, J. R., Yao, G. & West, M.
(2004), ‘Sparse graphical models for
exploring gene expression data’,
Journal of Multivariate Analysis
90(1), 196–212.  Fahrmeir & Lang (2001) Fahrmeir, L. & Lang, S. (2001), ‘Bayesian inference for generalized additive mixed models based on Markov random field priors’, Journal of the Royal Statistical Society: Series C (Applied Statistics) 50(2), 201–220.
 Frommlet et al. (2012) Frommlet, F., Ljubic, I., Arnardóttir Helga, B. & Bogdan, M. (2012), ‘QTL Mapping Using a Memetic Algorithm with Modifications of BIC as Fitness Function’, Statistical Applications in Genetics and Molecular Biology 11(4), 1–26.
 Gal (2016) Gal, Y. (2016), Uncertainty in Deep Learning, PhD thesis, University of Cambridge.
 Gal et al. (2017) Gal, Y., Hron, J. & Kendall, A. (2017), Concrete dropout, in ‘Advances in Neural Information Processing Systems’, pp. 3581–3590.
 George & McCulloch (1993) George, E. I. & McCulloch, R. E. (1993), ‘Variable selection via Gibbs sampling’, Journal of the American Statistical Association 88(423), 881–889.
 George et al. (2010) George, E. I. et al. (2010), Dilution priors: Compensating for model space redundancy, in ‘Borrowing Strength: Theory Powering Applications–A Festschrift for Lawrence D. Brown’, Institute of Mathematical Statistics, pp. 158–165.
 Ghosh & DoshiVelez (2017) Ghosh, S. & DoshiVelez, F. (2017), ‘Model selection in Bayesian neural networks via horseshoe priors’, arXiv preprint arXiv:1705.10388 .
 Ghosh et al. (2018) Ghosh, S., Yao, J. & DoshiVelez, F. (2018), ‘Structured Variational Learning of Bayesian Neural Networks with Horseshoe Priors’, arXiv preprint arXiv:1806.05975 .
 Graves (2011) Graves, A. (2011), Practical variational inference for neural networks, in ‘Advances in Neural Information Processing Systems’, pp. 2348–2356.
 Hubin (2018) Hubin, A. (2018), Bayesian model configuration, selection and averaging in complex regression contexts, PhD thesis, University of Oslo.
 Hubin & Storvik (2018) Hubin, A. & Storvik, G. (2018), ‘Mode jumping MCMC for Bayesian variable selection in GLMM’, Computational Statistics & Data Analysis .
 Hubin et al. (2018a) Hubin, A., Storvik, G. & Frommlet, F. (2018a), ‘Deep Bayesian regression models’, arXiv preprint arXiv:1806.02160 .

Hubin et al. (2018b)
Hubin, A., Storvik, G. & Frommlet, F. (2018b), ‘A novel algorithmic approach to Bayesian logic
regression’, Bayesian Anal. .
Advance publication.
https://doi.org/10.1214/18BA1141  Kingma & Ba (2014) Kingma, D. P. & Ba, J. (2014), ‘Adam: A method for stochastic optimization’, arXiv preprint arXiv:1412.6980 .
 LeCun et al. (1998) LeCun, Y., Cortes, C. & Burges, C. (1998), ‘MNIST handwritten digit database, 1998’, URL http://www. research. att. com/~ yann/ocr/mnist .
 Li & Clyde (2018) Li, Y. & Clyde, M. A. (2018), ‘Mixtures of gpriors in generalized linear models’, Journal of the American Statistical Association 113(524), 1828–1845.
 Logsdon et al. (2010) Logsdon, B. A., Hoffman, G. E. & Mezey, J. G. (2010), ‘A variational Bayes algorithm for fast and accurate multiple locus genomewide association analysis’, BMC bioinformatics 11(1), 58.
 Louizos et al. (2017) Louizos, C., Ullrich, K. & Welling, M. (2017), Bayesian compression for deep learning, in ‘Advances in Neural Information Processing Systems’, pp. 3288–3298.
 Louizos & Welling (2017) Louizos, C. & Welling, M. (2017), Multiplicative normalizing flows for variational bayesian neural networks, in ‘Proceedings of the 34th International Conference on Machine LearningVolume 70’, JMLR. org, pp. 2218–2227.
 Mnih & Gregor (2014) Mnih, A. & Gregor, K. (2014), ‘Neural variational inference and learning in belief networks’, arXiv preprint arXiv:1402.0030 .
 Molchanov et al. (2017) Molchanov, D., Ashukha, A. & Vetrov, D. (2017), Variational dropout sparsifies deep neural networks, in ‘Proceedings of the 34th International Conference on Machine LearningVolume 70’, JMLR. org, pp. 2498–2507.
 Neklyudov et al. (2017) Neklyudov, K., Molchanov, D., Ashukha, A. & Vetrov, D. P. (2017), Structured bayesian pruning via lognormal multiplicative noise, in ‘Advances in Neural Information Processing Systems’, pp. 6775–6784.
 Nitarshan (2018) Nitarshan, R. (2018), ‘Weight Uncertainty in Neural Networks’, https://www.nitarshan.com/bayesbybackprop/.

Nowozin (2018)
Nowozin, S. (2018), ‘Debiasing evidence approximations: On importanceweighted autoencoders and jackknife variational inference’.
 Palvanov & Im Cho (2018) Palvanov, A. & Im Cho, Y. (2018), ‘Comparisons of deep learning algorithms for MNIST in realtime environment’, International Journal of Fuzzy Logic and Intelligent Systems 18(2), 126–134.
 Polson & Rockova (2018) Polson, N. & Rockova, V. (2018), ‘Posterior concentration for sparse deep learning’, arXiv preprint arXiv:1803.09138 .
 Posch et al. (2019) Posch, K., Steinbrener, J. & Pilz, J. (2019), ‘Variational inference to measure model uncertainty in deep neural networks’, arXiv preprint arXiv:1902.10189 .
 Ranganath et al. (2016) Ranganath, R., Tran, D. & Blei, D. (2016), Hierarchical variational models, in ‘International Conference on Machine Learning’, pp. 324–333.
 Ripley (2007) Ripley, B. D. (2007), Pattern recognition and neural networks, Cambridge university press.
 Salimans et al. (2015) Salimans, T., Kingma, D. & Welling, M. (2015), Markov chain Monte Carlo and variational inference: Bridging the gap, in ‘International Conference on Machine Learning’, pp. 1218–1226.
 Smith & LeSage (2004) Smith, T. E. & LeSage, J. P. (2004), A Bayesian probit model with spatial dependencies, in ‘Spatial and spatiotemporal econometrics’, Emerald Group Publishing Limited, pp. 127–160.
 Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. & Salakhutdinov, R. (2014), ‘Dropout: a simple way to prevent neural networks from overfitting’, The Journal of Machine Learning Research 15(1), 1929–1958.
 Wasserman (2000) Wasserman, L. (2000), ‘Bayesian model selection and model averaging’, Journal of mathematical psychology 44(1), 92–107.

Weaver & Tao (2001)
Weaver, L. & Tao, N. (2001), The optimal reward baseline for gradientbased reinforcement learning,
in ‘Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence’, Morgan Kaufmann Publishers Inc., pp. 538–545.  Xiao et al. (2017) Xiao, H., Rasul, K. & Vollgraf, R. (2017), ‘Fashionmnist: a novel image dataset for benchmarking machine learning algorithms’, arXiv preprint arXiv:1708.07747 .
Web Supplement
Appendix A Proofs of propositions
Proof of Proposition 1.
We have
from which the result follows. ∎
Proof of Proposition 2.
Consider the continuous relaxation and reparametrization of the model
and similarly an approximative prior distribution given by  
with and . For this model converges to the original one. Note however that in this model we have a continuous dependence on the ’s. Direct application of the pathwise derivative estimation (Blundell et al., 2015; Gal, 2016; Ghosh & DoshiVelez, 2017) can then be performed for any fixed value of , similarly to Gal et al. (2017). However we are interested in the limiting case.
Define to be the expectation with respect to , then for this parametrization the ELBO becomes:
and its gradient can be written as:
To approximate the expectation part and the sum over all observations we use MonteCarlo sampling. Draw for and as a random subsample of size from . The unbiased estimator of the gradient of ELBO then becomes:
where . When , estimator converges to . ^{1}^{1}1Note that the gradient of the last term in actually could have been derived directly, leading to the same results. ∎
Proof of Remark 2.1.
The identity
leads to the result. ∎
Appendix B Other inference possibilities
Modelparameter posterior mean based inference
Approximate as , where , which within our variational Bayes approach simplifies to . Here no sampling (as in Algorithm 2) is needed, but no sparsification is achieved either.
Median probability model based inference combined with sampling
This approach is based on the notion of a median probability model, which was shown to be optimal in terms of predictions in the context of simple linear models (Barbieri et al., 2004). This corresponds to applying Algorithm 2 with the sampling of replaced by putting . Within this approach we significantly sparsify the network and only sample from the distributions of those weights that have marginal inclusion probabilities above 0.5. The rest of the weights are simply replaced with deterministic values of zero.
Median probability model based inference combined with parameter posterior mean
Approximate as , where , which within our variational Bayes approach simplifies to . Here again no sampling is needed and we only need to store the variational parameters of corresponding to marginal inclusion probabilities above 0.5. Hence we are significantly sparsifying the BNN of interest and reduce the computational cost of inference drastically.
Median probability model: Posttraining
Once it is decided to make inference based on the median probability model, one might take a number of additional iterations of the training algorithm with respect to the parameters of the models, having the architecture fixed. This might often give additional improvements in terms of the quality of inference as well as make the training steps much easier, since the number of parameters is reduced dramatically. This is so, since one does not have to estimate marginal inclusion probabilities any longer. Moreover, the number of weights ’s corresponding to to make inference on is significantly reduced due to the sparsity induced by using the median probability model.
Infeasibility remark: Other model selecting criteria and alternative thresholding
The median probability model is not always feasible in the sense that one needs at least one connected path across all of the layers with all of the weights linking the neurons having marginal inclusion above 0.5. One way to resolve the issue is to use the most probable model (model with the largest marginal posterior probability) instead of the median probability model. Then conditionally on its configuration one can sample from the distribution of the parameters, select the mean (mode) of the parameters or posttrain the distributions of the parameters. Other model selection criteria including DIC, WAIC, and FIC
(Claeskens et al., 2008)can be used in the same way as the most probable model. Another heuristic way to tackle the issue is to replace conditioning on
with , where is an arbitrary threshold. The latter will also improve predictive performance in case too conservative priors on the model configurations are used.Appendix C Extensions of the applications
c.1 Sparsity in the posterior
In the top row of Figure 3 we are reporting histograms of the marginal inclusion probabilities for the weights at all of the three hidden layers for MNIST data. Just like in Table 2 the histograms show that our approach yields extremely high sparsification levels for layers 1 and 2 and a more moderate sparsification at layer 3. A similar structure is seen in the bottom row for the FMNIST data.
Figures 4 and 5 show the sparse structure of typically sampled models from the model space and the corresponding weights. The final sample from the joint posterior is obtained by multiplying the model mask and weights elementwisely. There seems to be some pattern in the structure of the obtained weights for layer 1. This requires additional research and might allow to further compress the BNN.
c.2 Missclassification uncertainties
Figure 6
shows the missclassification uncertainties associated with posterior predictive sampling. One can clearly see that for almost all of the cases, when the BNN makes a missclassification, class certainty of the predictions is very low, indicating that the network is unsure. Moreover, even in these cases the truth is typically within the 95% credible interval of the predictions, which following
Posch et al. (2019) can be read from whether less than 95 out of 100 samples belong to a wrong class and at least 6 out of 100 samples belong to the right one. Also notice that in many of the cases of missclassification illustrated here, even a human would have serious doubts in making a decision.c.3 Inoutofdomain results
Following the example of measuring the in and outofdomain uncertainty suggested in Nitarshan (2018) we will test the ability of the approach to give confidence in its predictions by means of trying to classify a sample from FMNIST images with samples from the posterior predictive distribution based on the joint posterior of models and parameters trained on MNIST data set and compare this to the results for a sample of images from the test set of MNIST data. The results are reported for the joint posterior (of models and parameters) obtained in experiment run . As can be seen in Figure 7, the samples from BNN give highly confident predictions for the MNIST data set with almost no variance in the samples from the posterior predictive distribution. At the same time, the outofdomain uncertainty, related to the samples from the posterior predictive distribution based on FMNIST data, is typically high (with some exceptions) showing low confidence of the samples from the posterior predictive distribution in this case. The reversed example of inference on FMNIST and uncertainty related to MNIST data, illustrated in Figure 8, leads to exactly the same conclusions.
Comments
There are no comments yet.