Variational Inference for Bayesian Neural Networks under Model and
  Parameter Uncertainty
Bayesian neural networks (BNNs) have recently regained a significant amount
of attention in the deep learning community due to the development of scalable
approximate Bayesian inference techniques. There are several advantages of
using a Bayesian approach: Parameter and prediction uncertainties become easily
available, facilitating rigorous statistical analysis. Furthermore, prior
knowledge can be incorporated. However, so far, there have been no scalable
techniques capable of combining both structural and parameter uncertainty. In
this paper, we apply the concept of model uncertainty as a framework for
structural learning in BNNs and hence make inference in the joint space of
structures/models and parameters. Moreover, we suggest an adaptation of a
scalable variational inference approach with reparametrization of marginal
inclusion probabilities to incorporate the model space constraints.
Experimental results on a range of benchmark datasets show that we obtain
comparable accuracy results with the competing models, but based on methods
that are much more sparse than ordinary BNNs.
Aliaksandr Hubin, Geir Storvik
# Variational Inference for Bayesian Neural Networks under Model and Parameter Uncertainty

###### Abstract

Bayesian neural networks (BNNs) have recently regained a significant amount of attention in the deep learning community due to the development of scalable approximate Bayesian inference techniques. There are several advantages of using a Bayesian approach: Parameter and prediction uncertainties become easily available, facilitating rigorous statistical analysis. Furthermore, prior knowledge can be incorporated. However, so far, there have been no scalable techniques capable of combining both structural and parameter uncertainty. In this paper, we apply the concept of model uncertainty as a framework for structural learning in BNNs and hence make inference in the joint space of structures/models and parameters. Moreover, we suggest an adaptation of a scalable variational inference approach with reparametrization of marginal inclusion probabilities to incorporate the model space constraints. Experimental results on a range of benchmark datasets show that we obtain comparable accuracy results with the competing models, but based on methods that are much more sparse than ordinary BNNs.

Keywords:Bayesian neural networks; Structural learning; Model selection; Model averaging; Approximate Bayesian inference; Predictive uncertainty.

## 1 Introduction

In recent years, frequentist deep learning procedures have become extremely popular and highly successful in a wide variety of real-world applications ranging from natural language to image analyses (Goodfellow et al., 2016). These algorithms iteratively apply some nonlinear transformations aiming at optimal prediction of response variables from the outer layer features. This yields high flexibility in modelling complex conditional distributions of the responses. Each transformation yields another hidden layer of features which are also called neurons. The architecture/structure of a deep neural network includes the specification of the nonlinear intra-layer transformations (_activation functions_), the number of layers (_depth_), the number of features at each layer (_width_) and the connections between the neurons (_weights_). In the standard (frequentist) settings, the resulting model is trained using some optimization procedure (e.g. stochastic gradient descent) with respect to its parameters in order to fit a particular objective (like minimization of the root mean squared error or negative log-likelihood). Very often deep learning procedures outperform traditional statistical models, even when the latter are carefully designed and reflect expert knowledge (Refenes et al., 1994; Razi and Athapplily, 2005; Adya and Colopy, 1998; Sargent, 2001; Kanter and Veeramachaneni, 2015). However, typically one has to use huge datasets to be able to produce generalizable neural networks and avoid overfitting issues. Even though several regularization techniques (\(L_{1}\) and \(L_{2}\) penalties on the weights, dropout, batch normalization, etc.) have been developed for deep learning procedures to avoid overfitting to training datasets, the success of such approaches is not obvious. Unstuetured pruning the network, either by putting some weights to zero or by removing some nodes has been shown to be possible.
As an alternative to frequentist deep learning approaches, Bayesian neural networks represent a very flexible class of models, which are quite robust to overfitting (Neklyudov et al., 2018). However, they often remain heavily over-parametrized. There are several implicit approaches for the sparsification of BNNs by shrinkage of weights through priors (Jylanki et al., 2014; Blundell et al., 2015; Molchanov et al., 2017; Ghosh et al., 2018; Neklyudov et al., 2017). For example, Blundell et al. (2015) suggest a mixture of two Gaussian densities and then perform a fully factorizable mean-field variational approximation. Ghosh et al. (2019); 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). Some algorithmic procedures can also be seen to correspond to specific Bayesian priors, e.g. Molchanov et al. (2017) show that Gaussian dropout corresponds to BNNs with log uniform priors on the weight parameters.

In this paper, we consider a formal Bayesian approach for jointly taking into account _structural uncertainty_ and _parameter uncertainty_ in BNNs as a generalization of the methods developed for Bayesian model selection within linear regression models. The approach is based on introducing latent binary variables corresponding to the inclusion-exclusion of particular weights within a given architecture. This is done by means of introducing spike-and-slab priors. Such priors for the BNN setting were suggested in Polson and Rockova (2018) and Hubin (2018). A computational procedure for inference in such settings was proposed in an early version of this paper (Hubin and Storvik, 2019) and was further used without any changes in Bai et al. (2020). An asymptotic theoretical result for the choice of prior inclusion probabilities is additionally presented in Bai et al. (2020), however, it only is valid when the number of parameters goes to infinity and when the prior variance of the slab components is fixed. Here, we go further, in that we introduce hyperpriors on the prior inclusion probabilities and variance of the slab components making them stochastic, we also allow for more flexible variational approximations based on the multivariate Gaussian structures for inclusion indicators. Additionally, we consider several alternative prediction procedures for fully Bayesian model averaging, including posterior mean-based models, and the median probability model, and perform a comprehensive experimental study comparing the suggested approach with several competing algorithms and several data sets.

Using a Bayesian formalization in the space of models allows adapting the whole machinery of Bayesian inference in the joint model-parameter 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, 2021) and _posterior mean_ model-based inference for BNNs. Sparsifying properties of BMS (in particular the median probability model) are also addressed within the experiments. 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 inference algorithm is based on scalable stochastic variational inference.

The suggested approach has similarities to binary dropout that has become very popular (Srivastava et al., 2014). However, while standard binary dropout can be seen as a Bayesian approximation to a Gaussian process model where only parameter estimation is taken into account (Gal, 2016), our approach 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 of having the same dropout probabilities for the whole layer. Furthermore, the variational approximation procedure applied in Gal et al. (2017) has not been studied in the model uncertainty context.

At the same time, it is important to state explicitly that our approach does not currently aim at interpretable inference on Bayesian neural networks in most of the cases with an exception of having a special case of model selection in GLM models. Also, interpretable models could in principle be feasible if direct connections from all the layers to the responses are allowed. However, such cases are not addressed in this paper.
The rest of the paper is organized as follows: Section 2 gives some background and discuss related work. The class of BNNs and the corresponding model space are mathematically defined in Section 3. In Section 4, we describe the algorithm for training the suggested class of models using the reparametrization of marginal inclusion probabilities. Section 4.3 discusses several predictive inference possibilities. In Section 5, the suggested approach is applied to the two classical benchmark datasets MNIST, FMNIST (for image classifications) as well as PHONEME (for sound classification). We also compare the results with some of the existing approaches for inference on BNNs. Finally, in Section 6 some conclusions and suggestions for further research are given. Additional results are provided in the supplementary materials to the paper.

## 2 Background and related work

Bayesian neural networks (BNNs) were already introduced a few decades ago by Neal (1992); MacKay (1995); Bishop (1997). BNNs take advantage of the rigorous Bayesian approach and are able to properly handle parameter and prediction uncertainty and can in principle also incorporate prior knowledge. In many cases, this leads to more robust solutions with less overfitting. However, this comes at a price of extremely high computational costs. Until recently, inference on BNNs could not scale to large and high-dimensional data due to the limitations of standard MCMC approaches, the main numerical procedure in use. Several attempts based on subsampling techniques for MCMC, which are either approximate (Bardenet et al., 2014, 2017; Korattikara et al., 2014; Quiroz et al., 2014; Welling and Teh, 2011) or exact (Quiroz et al., 2016; Maclaurin and Adams, 2014; Liu et al., 2015; Welling and Teh, 2011) have been proposed, but none of them is able to explore the parameter spaces efficiently in ultrahigh-dimensional settings.

An alternative to the MCMC technique is to perform approximate Bayesian inference through variational Bayes, also known as variational inference (Jordan et al., 1999). Due to the fast convergence properties of the variational methods, variational inference algorithms are typically orders of magnitude faster than MCMC algorithms in high-dimensional problems (Ahmed et al., 2012). The variational inference has various applications in latent variable models, such as mixture models (Humphreys and Titterington, 2000), hidden Markov models (MacKay, 1997) and graphical models (Attias, 2000) in general. Graves (2011) suggested the methodology for scalable variational inference to Bayesian neural networks. This methodology was further improved by incorporating various variance reduction techniques, which are discussed in Gal (2016).

As mentioned in the introduction, it has been shown that the majority of the weight parameters in neural netwooks can be pruned out from the model without a significant loss of predictive accuracy. However, pruning is typically done implicitly by deleting the weights via ad-hoc thresholding. Yet, learning which parameters to include in a model, can also be framed as a structure learning or a model selection problem.

As discussed in Claeskens and Hjort (2008), Steel (2020) or Hansen and Yu (2001) (among other venues), model selection and model averaging in statistics generally assumes a discrete and countable (finite or infinite) set of models living on a corresponding model space. Models within a model space can differ in terms of the likelihood used, the link functions addressed, or which parameters are included in a linear or non-linear predictor. The purpose of model selection is to choose a single (best in some sense) model from a model space. This choice can lead to more interpretable models in some use cases (like variable selection in linear regression, Kuo and Mallick, 1998) or simply best models for some purpose (like prediction) in others (Geisser and Eddy, 1979), sometimes both coincide (Hubin et al., 2021), but that is not always the case (Breiman, 2001). Sparsity in model selection may or may not be of interest dependent on the context, although parsimony in some sense is typically desired. At the same time, the model selection often leads to problems with uncertainty handling resulting in too narrow confidence/credible intervals of the parameters (Heinze et al., 2018) and thus often may result in similar problems for predictions, i.e. overfitting. Model averaging (if model uncertainty is
properly addressed) can resolve these issues (Bornkamp et al., 2017) and has other advantages (Steel, 2020). Also, model uncertainty aware model selection, e.g. using the median probability model is more robust to overfitting (Ghosh, 2015).

There have been numerous works showing the efficiency and accuracy of model selection/averaging related to parameter selection through introducing latent variables corresponding to different discrete model configurations. In the Bayesian context, the posterior distribution can then be used to both select the best sparse configuration and address the joint model-and-parameters-uncertainty explicitly (George and McCulloch, 1993; Clyde et al., 2011; Frommlet et al., 2012; Hubin and Storvik, 2018; Hubin et al., 2020, 2021). Spike-and slab priors (Mitchell and Beauchamp, 1988) are typically used in this setting. All of these approaches have demonstrated both good predictive performance of the obtained sparse models and the ability to recover meaningful complex nonlinearities. They are however based on adaptations of Markov chain Monte Carlo (MCMC) and do not scale well to large high-dimensional 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 high-dimensional problems. Polson and Rockova (2018) study the use of the spike-and-slab approach in BNNs from a theoretical standpoint.

Logsdon et al. (2010); Carbonetto et al. (2012) suggest a fully-factorized variational distribution capable of efficiently and precisely "linearizing" the computational burden 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 of his Ph.D. 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.

## 3 The model

A neural network model links (possibly multidimensional) observations \(\mathbf{y}_{i}\in\mathcal{R}^{r}\) and explanatory variables \(\mathbf{x}_{i}\in\mathcal{R}^{p}\) via a probabilistic functional mapping with a mean parameter vector \(\mathbf{\mu}_{i}=\mathbf{\mu}_{i}(\mathbf{x}_{i})\in\mathcal{R}^{r}\):

\[\mathbf{y}_{i}\sim\mathfrak{f}\left(\mathbf{\mu}_{i}(\mathbf{x}_{i}),\phi \right),\quad i\in\{1,...,n\}, \tag{3.1}\]

where \(\mathfrak{f}\) is some observation distribution, typically from the exponential family, while \(\phi\) is a dispersion parameter. To construct the vector of mean parameters \(\mathbf{\mu}_{i}\), one builds a sequence of building blocks of hidden layers through semi-affine transformations:

\[z_{ij}^{(l+1)}= g_{j}^{(l)}\left(\beta_{0j}^{(l)}+\sum_{k=1}^{p^{(l)}}\beta_{ kj}^{(l)}z_{ik}^{(l)}\right),l=1,...,L-1,j=1,...,p^{(l+1)}, \tag{3.2}\]

with \(\mu_{ij}=z_{ij}^{(L)}\). Here, \(L\) is the number of layers, \(p^{(l)}\) is the number of nodes within the corresponding layer while \(g_{j}^{(l)}\) is a univariate function (further referred to as the _activation function_). Further, \(\beta_{kj}^{(l)}\in\mathcal{R},k>0\) are the weights (slope coefficients) for the inputs \(z_{ik}^{(l)}\) of the \(l\)-th layer (note that \(z_{ik}^{(1)}=x_{ik}\) and \(p^{(1)}=p\)). For \(k=0\), we obtain the intercept/bias terms. Finally, we introduce latent binary indicators \(\gamma_{kj}^{(l)}\in\{0,1\}\) switching the corresponding weights on and off such that \(\beta_{kj}^{(l)}=0\) if \(\gamma_{kj}^{(l)}=0\).

In our notation, we explicitly differentiate between discrete structural/model configurations defined by the vectors \(\mathbf{\gamma}=\{\gamma_{kj}^{(l)},j=1,..,p^{(l+1)},k=0,...,p^{(l)},l=1,...,L-1\}\) (further referred to as models) constituting the model space \(\Gamma\) and
parameters of the models, conditional on these configurations \(\mathbf{\theta}|\mathbf{\gamma}=\{\mathbf{\beta},\phi|\mathbf{\gamma}\}\), where only those \(\beta^{(l)}_{kj}\) for which \(\gamma^{(l)}_{kj}=1\) are included. This approach is (in statistical science literature) a rather standard way to explicitly specify the model uncertainty in a given class of models and is used in e.g. Clyde et al. (2011); Frommlet et al. (2012); Hubin et al. (2021).

A Bayesian approach is completed by specification of model priors \(p(\mathbf{\gamma})\) and parameter priors for each model \(p(\mathbf{\beta}|\mathbf{\gamma},\phi)\). If the dispersion parameter is present in the distribution of the outcomes, one also has to define a prior \(p(\phi|\mathbf{\gamma})\). Many kinds of priors on \(p(\mathbf{\beta}|\mathbf{\gamma},\phi)\) can be considered, including the mixture of Gaussians prior (Blundell et al., 2015), the Horseshoe prior (Ghosh et al., 2019; Louizos et al., 2017), or mixtures of g-priors (Li and Clyde, 2018), which could give further penalties to the weight parameters. We first following our early preprint (Hubin and Storvik, 2019) as well as even earlier ideas from Hubin (2018); Polson and Rockova (2018) and consider the independent Gaussian spike-and-slab weight priors combined with independent Bernoulli priors for the latent inclusion indicators being equal to \(1\). This choice of the priors corresponds to marginal spike-and-slab priors for the weights (Clyde et al., 2011):

\[p(\beta^{(l)}_{kj}|\sigma^{2}_{\beta,l},\gamma^{(l)}_{kj})= \gamma^{(l)}_{kj}\mathcal{N}(0,\sigma^{2}_{\beta,l})+(1-\gamma^{( l)}_{kj})\delta_{0}(\beta^{(l)}_{kj}), \tag{3.3a}\] \[p(\gamma^{(l)}_{kj})= \text{Bernoulli}(\psi^{(l)}), \tag{3.3b}\]

Here, \(\delta_{0}(\cdot)\) is the delta mass or "spike" at zero, \(\sigma^{2}_{\beta,l}\) is the prior variance of \(\beta^{(l)}_{kj}\), whilst \(\psi^{(l)}\in(0,1)\) is the prior probability for including the weight \(\beta^{(l)}_{kj}\) into the model.

To automatically infer the prior variance and the prior probability for including a weight, we assume for \(\sigma^{2}_{\beta,l}\) a standard inverse Gamma hyperprior with hyperparameters \(a^{(l)}_{\beta},b^{(l)}_{\beta}\), and for \(\psi^{(l)}\) a Beta\((a^{(l)}_{\psi},b^{(l)}_{\psi})\) prior:

\[p(\beta^{(l)}_{kj}|\sigma^{2}_{\beta,l},\gamma^{(l)}_{kj})= \gamma^{(l)}_{kj}\mathcal{N}(0,\sigma^{2}_{\beta,l})+(1-\gamma^{ (l)}_{kj})\delta_{0}(\beta^{(l)}_{kj}), \tag{3.4a}\] \[p(\sigma^{2}_{\beta,l})= \text{Inv-Gamma}(a^{(l)}_{\beta},b^{(l)}_{\beta}),\] (3.4b) \[p(\gamma^{(l)}_{kj})= \text{Bernoulli}(\psi^{(l)}),\] (3.4c) \[p(\psi^{(l)})= \text{Beta}(a^{(l)}_{\psi},b^{(l)}_{\psi}). \tag{3.4d}\]

Here, \(\delta_{0}(\cdot)\) is the delta mass or "spike" at zero, \(\sigma^{2}_{\beta,l}\) is the prior variance of \(\beta^{(l)}_{kj}\) with a standard inverse Gamma hyperprior with hyperparameters \(a^{(l)}_{\beta},b^{(l)}_{\beta}\), whilst \(\psi^{(l)}\in[0,1]\) is the prior probability for including the weight \(\beta^{(l)}_{kj}\) into the model. Further, \(\psi^{(l)}\) is assumed Beta\((a^{(l)}_{\psi},b^{(l)}_{\psi})\) distributed. With the presence of a dispersion parameter, an additional prior is needed, see Dey et al. (2000). We will refer to our model as the Latent Binary Bayesian Neural Network (LBBNN) model.

## 4 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 \(\Delta\) (for example the distribution of a new observation \(y^{*}\) conditional on new covariates \(\mathbf{x}^{*}\)) based on data \(\mathcal{D}\):

\[p(\Delta|\mathcal{D})=\sum_{\mathbf{\gamma}\in\Gamma}\int_{\mathbf{\theta}\in\Theta_{ \gamma}}p(\Delta|\mathbf{\theta},\mathbf{\gamma},\mathcal{D})p(\mathbf{\theta},\mathbf{\gamma} |\mathcal{D})\mathrm{d}\mathbf{\theta}, \tag{4.1}\]
where \(\Theta_{\gamma}\) is the parameter space defined through \(\mathbf{\gamma}\). Standard procedures for dealing with complex posteriors is to apply Monte Carlo methods which involve simulations from \(p(\mathbf{\theta},\mathbf{\gamma}|\mathcal{D})\). For the model defined by (3.1)-(3.4) with many hidden variables, such simulations become problematic.

The idea behind variational inference (Graves, 2011; Blei et al., 2017) is to apply the approximation

\[\tilde{p}(\Delta|\mathcal{D})=\sum_{\mathbf{\gamma}\in\Gamma}\int_{\mathbf{\theta}\in \Theta_{\gamma}}p(\Delta|\mathbf{\theta},\mathbf{\gamma},\mathcal{D})q_{\mathbf{\eta}}( \mathbf{\theta},\mathbf{\gamma})\mathrm{d}\mathbf{\theta}. \tag{4.2}\]

for some suitable (parametric) distribution \(q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{\gamma})\) which, with appropriate choices of the parameters \(\mathbf{\eta}\), approximates the posterior well and is _simple_ to sample from. The specification of \(\mathbf{\eta}\) is typically obtained through the minimization of the Kullback-Leibler divergence from the variational family distribution to the posterior distribution:

\[\mathrm{KL}(q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{\gamma})||p(\mathbf{\theta},\mathbf{\gamma} |\mathcal{D}))=\sum_{\mathbf{\gamma}\in\Gamma}\int_{\Theta_{\gamma}}q_{\mathbf{\eta}} (\mathbf{\theta},\mathbf{\gamma})\log\tfrac{q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{\gamma})}{ p(\mathbf{\theta},\mathbf{\gamma}|\mathcal{D})}\mathrm{d}\mathbf{\theta}, \tag{4.3}\]

with respect to the variational parameters \(\mathbf{\eta}\). Compared to standard variational inference approaches, the setting is extended to include the discrete model identifiers \(\mathbf{\gamma}\). For an optimal choice \(\hat{\mathbf{\eta}}\) of \(\mathbf{\eta}\), inference on \(\Delta\) is performed through Monte Carlo estimation of (4.2) inserting \(\hat{\mathbf{\eta}}\) for \(\mathbf{\eta}\). The main challenge then becomes choosing a suitable variational family and a computational procedure for minimizing (4.3). Note that although this minimization is still a computational challenge, it will typically be much easier than directly obtaining samples from the true posterior. The final Monte Carlo estimation will be simple, provided the variational distribution \(q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{\gamma})\) is selected such that it is simple to sample from.

As in standard settings of variational inference, minimization of the divergence (4.3) is equivalent to maximization of the evidence lower bound (ELBO)

\[\mathcal{L}_{VI}(\mathbf{\eta})= \sum_{\mathbf{\gamma}\in\Gamma}\int_{\Theta_{\gamma}}q_{\mathbf{\eta}}( \mathbf{\theta},\mathbf{\gamma})\log p(\mathcal{D}|\mathbf{\theta},\mathbf{\gamma})\mathrm{d} \mathbf{\theta}-\mathrm{KL}\left(q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{\gamma})||p(\mathbf{ \theta},\mathbf{\gamma})\right) \tag{4.4}\]

through the equality

\[\mathcal{L}_{VI}(\mathbf{\eta})= p(\mathcal{D})-\mathrm{KL}\left(q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{ \gamma})||p(\mathbf{\theta},\mathbf{\gamma}|\mathcal{D})\right),\]

which also shows that \(\mathcal{L}_{VI}(\mathbf{\eta})\) is a lower bound of the marginal likelihood \(p(\mathcal{D})\).

### Variational distributions

We will consider a variational family previously proposed for linear regression (Logsdon et al., 2010; Carbonetto et al., 2012) which we extend to the LBBNN setting. Assume

\[q_{\mathbf{\eta}}(\mathbf{\theta},\mathbf{\gamma})=q_{\mathbf{\eta}_{0}}(\phi)\prod_{l=1}^{L-1} \prod_{j=1}^{p^{(l+1)}}\prod_{k=0}^{p^{(l)}}q_{\kappa_{kj},\tau_{kj}}(\beta_{kj }^{(l)}|\gamma_{kj}^{(l)})q_{\alpha_{kj}^{(l)}}(\gamma_{kj}^{(l)}), \tag{4.5}\]

where \(q_{\mathbf{\eta}_{0}}(\phi)\) is some appropriate distribution for the dispersion parameter,

\[q_{\kappa_{kj}^{(l)},\tau_{kj}^{(l)}}\left(\beta_{kj}^{(l)}|\gamma_{kj}^{(l)} \right)= \gamma_{kj}^{(l)}\mathcal{N}(\kappa_{kj}^{(l)},\tau^{2}{}_{kj}^{( l)})+(1-\gamma_{kj}^{(l)})\delta_{0}(\beta_{kj}^{(l)}), \tag{4.6}\]
and

\[q_{\alpha_{kj}^{(l)}}(\gamma_{kj}^{(l)})= \text{Bernoulli}(\alpha_{kj}^{(l)}). \tag{4.7}\]

With probability, \(\alpha_{kj}^{(l)}\in[0,1]\), the posterior of parameters of weight \(\beta_{kj}^{(l)}\) will be approximated by a normal distribution with some mean and variance ("slab"), and otherwise, the weight is put to zero. Thus, \(\alpha_{kj}^{(l)}\) will approximate the marginal posterior inclusion probability of the weight \(\beta_{kj}^{(l)}\). Here, \(\mathbf{\eta}=\{\mathbf{\eta}_{0},(\kappa_{kj}^{(l)},\tau_{kj}^{2(l)},\alpha_{kj}^{(l )}),l=1,...,L-1,k=1,...,p^{(l+1)},j=1,...,p^{(l)}\}\).

A similar variational distribution has also been considered within BNN through the dropout approach (Srivastava et al., 2014). For dropout, however, the final network is dense but trained through a Monte Carlo average of sparse networks. In our approach, the target distribution is different in the sense of including the binary variables \(\{\gamma_{kj}^{(l)}\}\) as part of the _model_. 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.

The variational distribution (4.5)-(4.7), corresponding to the commonly applied mean-field approximation, can be seen as a rather crude approximation, which completely ignores all posterior dependence between the model structures or parameters. Consequently, the resulting conclusions can be misleading or inaccurate as the posterior probability of one weight might be highly affected by the inclusion of others. Such a dependence structure can be built into the variational approximation either through the \(\gamma\)'s or through the \(\beta\)'s (or both). Here, we only consider dependence structures in the inclusion variables. We still assume independence between layers, but within layers, we introduce a dependence structure by defining \(\mathbf{\alpha}^{l}=\{\alpha_{kj}^{(l)}\}\) now to be a stochastic vector, which on logit-scale follows a multivariate normal distribution:

\[\text{logit}(\mathbf{\alpha}^{(l)})\sim\text{MVN}(\mathbf{\xi}^{(l)},\mathbf{\Sigma}^{(l) }). \tag{4.8}\]

Here, either a full covariance matrix \(\mathbf{\Sigma}^{(l)}\) or a low-rank parametrization for the covariance is possible. For the latter, \(\mathbf{\Sigma}^{(l)}=\mathbf{F}^{(l)}\mathbf{F}^{(l)}{}^{T}+\mathbf{D}^{(l)}\) with \(\mathbf{F}^{(l)}\) being the factor part of low-rank form of covariance matrix and \(D^{(l)}\) is the diagonal part of low-rank form of covariance matrix. This drastically reduces the number of parameters and allows for efficient computations of the determinant and inverse matrix. A particularly interesting case is when \(\mathbf{F}^{(l)}\) has rank zero in which case we retain independence between the components but some penalization in the variability of the \(\alpha_{kj}^{(l)}\)'s. Under the parametrization (4.6)-(4.8), the parameters \(\{\mathbf{\xi}^{(l)},l=1,...,L-1\}\) and \(\{\mathbf{\Sigma}^{(l)},l=1,...,L-1\}\) are added to the parameter vector \(\mathbf{\eta}\). Then the reparametrization trick is also performed for these parameters using the default representations for MVN and LFMVN, available out of the box in PyTorch probabilities.

### Optimization by stochastic gradient

For simplicity, we assume here that there is no dispersion parameters, so the target distribution is \(p(\mathbf{\beta},\mathbf{\gamma}|\mathcal{D})\). We can rewrite the ELBO (4.4) as

\[\mathcal{L}_{VI}(\mathbf{\eta})= \sum_{\mathbf{\gamma}\in\Gamma}\int_{\beta_{\gamma}}q_{\mathbf{\eta}}( \mathbf{\beta},\mathbf{\gamma})[\log p(\mathcal{D}|\mathbf{\beta},\mathbf{\gamma})-\log\tfrac {q_{\mathbf{\eta}}(\mathbf{\beta},\mathbf{\gamma})}{p(\mathbf{\beta},\mathbf{\gamma})}]\text{d} \mathbf{\beta}. \tag{4.9}\]

Due to the huge computational cost in the computation of gradients when \(\Gamma\) and \(\mathcal{D}\) are large, stochastic gradient methods using Monte Carlo estimates for obtaining unbiased estimates of the gradients have become the standard approach for
variational inference in such situations. Both the reparametrization trick and minibatching (Kingma et al., 2015; Blundell et al., 2015) are further applied.

Another complication in our setting is the discrete nature of \(\mathbf{\gamma}\). Following Gal et al. (2017), we relax the Bernoulli distribution (4.7) with the _Concrete distribution_:

\[\tilde{\gamma}=\gamma_{tr}(\nu,\delta;\alpha)=\text{sigmoid}(( \text{logit}(\alpha)-\text{logit}(\nu))/\delta),\quad\nu\sim\text{Unif}[0,1], \tag{4.10}\]

where \(\delta\) is a tuning parameter, which is selected to take some small value. In the zero limit, \(\tilde{\gamma}\) reduces to a Bernoulli(\(\alpha\)) variable. Combined with the reparametrization of the \(\beta\)'s,

\[\beta=\beta_{tr}(\varepsilon;\kappa,\tau)=\kappa+\tau\varepsilon, \quad\varepsilon\sim N(0,1) \tag{4.11}\]

we define the following approximation to the ELBO:

\[\mathcal{L}_{VI}^{\delta}(\mathbf{\eta}):=\int_{\mathbf{\nu}}\int_{\mathbf{ \varepsilon}}q_{\mathbf{\nu},\mathbf{\varepsilon}}(\mathbf{\nu},\mathbf{\varepsilon})[ \log p(\mathcal{D}|\beta_{tr}(\mathbf{\varepsilon},\mathbf{\kappa},\mathbf{\tau}),\gamma _{tr}(\mathbf{\nu},\mathbf{\alpha},\delta))-\] \[\log\frac{q_{\mathbf{\eta}}(\beta_{tr}(\mathbf{\varepsilon},\mathbf{ \kappa},\mathbf{\tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta))}{p(\beta_{tr}( \mathbf{\varepsilon},\mathbf{\kappa},\mathbf{\tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha}, \delta))}]\text{d}\mathbf{\varepsilon}\text{d}\mathbf{\nu} \tag{4.12}\]

where the transformations on vectors are performed elementwise. Further, due to that \(\text{d}\mathbf{\varepsilon}\text{d}\mathbf{\nu}\) does not depend on \(\mathbf{\eta}\), we can change the order of integration and differentiation when taking the gradient of \(\mathcal{L}_{VI}^{\delta}(\mathbf{\eta})\):

\[\nabla_{\mathbf{\eta}}\mathcal{L}_{VI}^{\delta}(\mathbf{\eta})=\int_{ \mathbf{\nu}}\int_{\mathbf{\varepsilon}}q_{\mathbf{\nu},\mathbf{\varepsilon}}(\mathbf{\nu},\mathbf{ \varepsilon})\nabla_{\mathbf{\eta}}[ \log p(\mathcal{D}|\beta_{tr}(\mathbf{\varepsilon},\mathbf{\kappa},\mathbf{ \tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta)-\] \[\log\frac{q_{\mathbf{\eta}}(\beta_{tr}(\mathbf{\varepsilon},\mathbf{\kappa}, \mathbf{\tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta))}{p(\beta_{tr}(\mathbf{ \varepsilon},\mathbf{\kappa},\mathbf{\tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta)) }]\text{d}\mathbf{\varepsilon}\text{d}\mathbf{\nu}. \tag{4.13}\]

An unbiased estimator of \(\nabla_{\mathbf{\eta}}\mathcal{L}_{VI}^{\delta}(\mathbf{\eta})\) is then given in Proposition 1.

**Proposition 1**.: _Assume for all \(m=1,...,M\)\(\left(\mathbf{\nu}^{(m)},\mathbf{\eta}^{(m)}\right)\sim q_{\mathbf{\nu},\mathbf{ \varepsilon}}(\mathbf{\nu},\mathbf{\eta})\)and \(S\) is a random subset of indices \(\{1,...,n\}\) of size N. Also, assume the observations to be conditionally independent. Then, for any \(\delta>0\), an unbiased estimator for the gradient of \(\mathcal{L}_{VI}^{\delta}(\mathbf{\eta})\) is given by_

\[\widetilde{\nabla}_{\mathbf{\eta}}\mathcal{L}_{VI}^{\delta}(\mathbf{\eta })=\frac{1}{M}\sum_{m=1}^{M}\Big{[}\frac{n}{N}\sum_{i\in S}\nabla_{\mathbf{\eta}} \log p(\mathbf{y}_{i}|\mathbf{x}_{i},\beta_{tr}(\mathbf{\varepsilon}^{(m)},\mathbf{\kappa}, \mathbf{\tau}),\gamma_{tr}(\mathbf{\nu}^{(m)},\mathbf{\alpha},\delta))-\] \[\nabla_{\mathbf{\eta}}\log\frac{q_{\mathbf{\eta}}(\beta_{tr}(\mathbf{ \varepsilon}^{(m)},\mathbf{\kappa},\mathbf{\tau}),\gamma_{tr}(\mathbf{\nu}^{(m)},\mathbf{ \alpha},\delta))}{p(\beta_{tr}(\mathbf{\varepsilon}^{(m)},\mathbf{\kappa},\mathbf{\tau}), \gamma_{tr}(\mathbf{\nu}^{(m)},\mathbf{\alpha},\delta))}\Big{]}. \tag{4.14}\]

Proof.: From (4.13) we have that

\[\frac{1}{M}\sum_{m=1}^{M}\nabla_{\mathbf{\eta}}[\log p(\mathcal{D} |\beta_{tr}(\mathbf{\varepsilon}^{(m)},\mathbf{\kappa},\mathbf{\tau}),\gamma_{tr}(\mathbf{\nu} ^{(m)},\mathbf{\alpha},\delta)-\log\frac{q_{\mathbf{\eta}}(\beta_{tr}(\mathbf{ \varepsilon}^{(m)},\mathbf{\kappa},\mathbf{\tau}),\gamma_{tr}(\mathbf{\nu}^{(m)},\mathbf{ \alpha},\delta))}{p(\beta_{tr}(\mathbf{\varepsilon}^{(m)},\mathbf{\kappa},\mathbf{\tau}), \gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta))}]\]

is an unbiased estimate of the gradient. Further, since we assume the observations to be conditionally independent, we have

\[\nabla_{\mathbf{\eta}}\log p(\mathcal{D}|\beta_{tr}(\mathbf{\varepsilon},\mathbf{\kappa}, \mathbf{\tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta))=\sum_{i=1}^{n}\nabla_{\mathbf{ \eta}}\log p(\mathbf{y}_{i}|\mathbf{x}_{i};\beta_{tr}(\mathbf{\varepsilon},\mathbf{\kappa}, \mathbf{\tau}),\gamma_{tr}(\mathbf{\nu},\mathbf{\alpha},\delta)),\]

for which an unbiased estimator can be constructed through a random subset, showing the result.
```
sample\(N\) indices uniformly from \(\{1,...,n\}\) defining \(S\); for\(m\)in\(\{1,...,M\}\)do for\((k,j,l)\in\mathcal{B}\)do sample\(\nu_{kl}^{(l)}\sim\text{Unif}[0,1]\) and \(\varepsilon_{kj}^{(l)}\sim N(0,1)\) endfor endfor calculate\(\widetilde{\nabla}_{\eta}\mathcal{L}_{VI}^{\delta}(\mathbf{\eta})\) according to (4.14) update\(\mathbf{\eta}\leftarrow\mathbf{\eta}+\mathbf{A}\widetilde{\nabla}_{\mathbf{\eta}}\mathcal{L}_{VI}^{ \delta}(\mathbf{\eta})\)
```

**Algorithm 1** Doubly stochastic variational inference step

Algorithm 1 describes one iteration of a doubly stochastic variational inference approach where updating is performed on the parameters for the case of mean field assumption (for simplicity). The set \(\mathcal{B}\) is the collection of all combinations \(j,k,l\) in the network. The matrix of learning rates \(\mathbf{A}\) will always be diagonal, allowing for different step sizes on the parameters involved. Following Blundell et al. (2015), constraints of \(\tau_{kj}^{(l)}\) are incorporated by means of the reparametrization \(\tau_{kj}^{(l)}=\log(1+\exp(\rho_{kj}^{(l)}))\) where \(\rho_{kj}^{(l)}\in\mathcal{R}\). Typically, updating is performed over a full _epoch_, in which case the observations are divided into \(n/N\) subsets and updating is performed sequentially over all subsets.

In case of the dependence structure (4.8), \(\mathbf{\alpha}\) is sampled instead while \(\mathbf{\xi}^{(l)}\) and the components of \(\mathbf{\Sigma}^{(l)}\) go into \(\mathbf{\eta}\). Constraints on \(\alpha_{kj}^{(l)}\) are incorporated by means of the reparametrization \(\alpha_{kj}^{(l)}=(1+\exp(-\omega_{kj}^{(l)}))^{-1}\) with \(\omega_{kj}^{(l)}\in\mathcal{R}\)

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 can be calculated by the usual backpropagation algorithm on a neural network. This algorithm assumes a known dispersion parameter \(\phi\), but can be easily generalized to include learning about \(\phi\) as well.

### Prediction

Once the estimates \(\widehat{\mathbf{\eta}}\) of the parameters \(\mathbf{\eta}\) of the variational approximating distribution are obtained, we go back to the original discrete model for \(\mathbf{\gamma}\) (setting \(\delta=0\)). Then, there are several ways to proceed with predictive inference. We list these below.

```
sample\(N\) indices uniformly from \(\{1,...,n\}\) defining \(S\); for\(m\)in\(\{1,...,M\}\)do for\((k,j,l)\in\mathcal{B}\)do sample\(\nu_{kl}^{(l)}\sim\text{Unif}[0,1]\) and \(\varepsilon_{kj}^{(l)}\sim N(0,1)\) endfor endfor calculate\(\widetilde{\nabla}_{\eta}\mathcal{L}_{VI}^{\delta}(\mathbf{\eta})\) according to (4.14) update\(\mathbf{\eta}\leftarrow\mathbf{\eta}+\mathbf{A}\widetilde{\nabla}_{\mathbf{\eta}}\mathcal{L}_{VI}^{ \delta}(\mathbf{\eta})\)
```

**Algorithm 2** Doubly stochastic variational inference step

#### Fully Bayesian model averaging

In this case, define

\[\hat{p}(\Delta|\mathcal{D})=\frac{1}{R}\sum_{r=1}^{R}p(\Delta|\mathbf{\beta}^{r}, \mathbf{\gamma}^{r}) \tag{4.15}\]

where \((\mathbf{\beta}^{r},\mathbf{\gamma}^{r})\sim q_{\widehat{\mathbf{\eta}}}(\mathbf{\beta},\mathbf{ \gamma})\). This procedure takes uncertainty in both the model structure \(\mathbf{\gamma}\) and the parameters \(\mathbf{\beta}\) into account in a formal Bayesian setting. A bottleneck of this approach is that we have to both sample from a huge approximate posterior distribution of parameters and models _and_ keep all of the components of \(\widehat{\mathbf{\eta}}\) stored during the prediction phase, which might be computationally and memory inefficient.
The posterior mean based model (Wasserman, 2000)In this case we put \(\beta^{(l)}_{kj}=\hat{E}\{\beta^{(l)}_{kj}|\mathcal{D}\}\) where

\[E\{\beta^{(l)}_{kj}|\mathcal{D}\}= p(\gamma^{(l)}_{kj}=1|\mathcal{D})E\{\beta^{(l)}_{kj}|\gamma^{(l)}_{ kj}=1,\mathcal{D}\}\approx\hat{\alpha}^{(l)}_{kj}\hat{\kappa}^{(l)}_{kj}.\]

Here \(\hat{\alpha}^{(l)}_{kj}\) is either the estimate of \(\alpha^{(l)}_{kj}\) obtained through the variational inference procedure or, in case of dependence structure \(\boldsymbol{\alpha}\)\(E\{\alpha^{(l)}_{kj}|\mathcal{D}\}\). In the latter case, one formally would want to integrate out \(\boldsymbol{\alpha}\) instead but this is not quite feasible in practice and extra sampling is avoided. This approach specifies one dense model \(\hat{\boldsymbol{\gamma}}\) with no sparsification. At the same time, no sampling is needed.

The median probability model (Barbieri et al., 2004)This approach is based on the notion of a median probability model, which has been shown to be optimal in terms of predictions in the context of simple linear models. Here, we set \(\gamma^{(l)}_{kj}=\mathrm{I}(\hat{\alpha}^{(l)}_{kj}>0.5)\) while \(\beta^{(l)}_{kj}\sim\gamma^{(l)}_{kj}N(\hat{\kappa}^{(l)}_{kj},\hat{\tau}^{2( l)}_{kj})\). A model averaging approach similar to (4.15) is then applied. 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.

Median probability model-based inference combined with parameter posterior meanHere, again, we set \(\gamma^{(l)}_{kj}=\mathrm{I}(\hat{\alpha}^{(l)}_{kj}>0.5)\) but now we use \(\beta^{(l)}_{kj}=\gamma^{(l)}_{kj}\hat{\kappa}^{(l)}_{kj}\). Similarly to the posterior mean-based model, no sampling is needed but in addition, we only need to store the variational parameters of \(\widehat{\boldsymbol{\eta}}\) corresponding to marginal inclusion probabilities above 0.5. Hence, we significantly sparsify the BNN of interest and reduce the computational cost of the predictions drastically.

Post-trainingOnce it is decided to make inference based on a selected model, one might take several additional iterations of the training algorithm concerning the parameters of the models, having the architecture-related parameters fixed. This might 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 \(\boldsymbol{\alpha}\) any longer. Moreover, the number of weights \(\beta^{(l)}_{jk}\)'s corresponding to \(\gamma^{(l)}_{jk}=1\) to make inference on is typically significantly reduced due to the sparsity induced by using the selected median probability model. It is also possible to keep the \(\alpha^{(l)}_{jk}\)'s fixed but still allow the \(\gamma^{(l)}_{jk}\)'s to be random.

Other model selecting criteria and alternative thresholdingThe median probability model is (at least in theory) 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 (the 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 post-train the distributions of the parameters. Other model selection criteria, including DIC and WAIC, can be used in the same way as the most probable model. Another heuristic way to tackle the issue is to replace conditioning on \(\mathrm{I}(\gamma^{(l)}_{kj}>0.5)\) with \(\mathrm{I}(\gamma^{(l)}_{kj}>\lambda)\), where \(\lambda\) is a tuning parameter. The latter might also improve predictive performance in case too conservative priors on the model configurations are used. At the same time, we are not addressing the methods described in this paragraph in our experiments and rather leave them for further research.
## 5 Applications

In-depth studies of the suggested variational approximations in the context of _linear_ regression models have been performed in earlier studies, including multiple synthetic and real data examples with the aims of both recovering meaningful relations and predictions (Carbonetto et al., 2012; Hernandez-Lobato et al., 2015). The results from these studies show that the approximations based on the suggested variational family distributions are reasonably precise and indeed scalable, but can be biased. We will not address toy examples and simulation-based examples in this article and rather refer the curious readers to the very detailed and comprehensive studies in the references mentioned above, whilst we will address some more complex examples here. In particular, we will address the classification of MNIST (LeCun et al., 1998) and fashion-MNIST (FMNIST Xiao et al., 2017) images as well as the PHONEME data (Hastie et al., 1995). Both MNIST and FMNIST datasets comprise of 70 000 grayscale images (size 28x28) from 10 categories (handwritten digits from 0 to 9, and "Top", "Trouser", "Pullover", "Dress", "Coat", "Sandal", "Shirt", "Sneaker", "Bag", and "Ankle Boot" Zalando's fashion items respectively), with 7 000 images per category. The training sets consist of 60 000 images, and the test sets have 10 000 images. For the PHONEME dataset, we have 256 covariates and 5 classes in the responses. In this dataset, we have 3 500 observations in the training set and 1 000 - in the test set. The PHONEME data are extracted from the TIMIT database (TIMIT Acoustic-Phonetic Continuous Speech Corpus, NTIS, US Dept of Commerce), which is a widely used resource for research in speech recognition. This dataset was formed by selecting five phonemes for classification based on a digitized speech from this database. The phonemes are transcribed as follows: "sh" as in "she", "dcl" as in "dark", "iy" as the vowel in "she", "aa" as the vowel in "dark", and "ao" as the first vowel in "water".

Experimental designFor all the datasets, we address a dense neural network with the ReLU activation function, and multinomially distributed observations. For the two first examples, we have 10 classes and 784 input explanatory variables (pixels) while for the third one, we have 256 input variables and 5 classes. In all three cases, the network has 2 hidden layers with 400, and 600 neurons correspondingly. Priors for the parameters and model indicators were chosen according to (3.4) with parameters in the priors specified through an empirical Bayes approach. The inference was performed using the suggested doubly stochastic variational inference approach (Algorithm 1) on 250 epochs with a batch size of 100. \(M\) was set to 1 to reduce computational costs and due to the fact that this choice of \(M\) is argued to be sufficient in combination with the reparametrization trick (Gal, 2016). Up to 20 first epochs were used for pre-training of the models and parameters as well as empirically (aka Empirical Bayes) estimating the hyperparameters of the priors \((a_{\psi},b_{\psi},a_{\beta},b_{\beta})\) through adding them into the computational graph. After that, the main training cycle began (with fixed hyperparameters on the priors). We used the ADAM stochastic gradient ascent optimization (Kingma and Ba, 2014) with the diagonal matrix \(\mathbf{A}\) in Algorithm 1 and the diagonal elements specified in Tables 1 and S-1 for pre-training, and the main training stage. Typically, one would maximize the marginal likelihood for the Empirical Bayes methods, but since we do not have it available for the addressed models, its lower bound (ELBO) was used in the pre-training stage. After 250 training epochs, post-training was performed. When post-training the parameters, either with fixed marginal inclusion probabilities or with the median probability model, we ran additional 50 epochs of the optimization routine with \(\mathbf{A}\) specified in the bottom rows of Tables 1 and S-1. For the fully Bayesian model averaging approach, we used both \(R=1\) and \(R=10\). Even though \(R=1\) can give a poor Monte Carlo estimate of the prediction distribution, it can be of interest due to high sparsification. All the PyTorch implementations used in the experiments are available in our GitHub repository.

We report results for _our model_**LBBNN** applied with the Gaussian priors (**GP**) for the slab components of \(\beta\)'s combined with variational inference based on Mean-Field (**MF**), MVN (**MVN**) and Low Factor MVN (for **LFMVN**, the predictions' results are reported in the supplemental material) dependence structures between the latent indicators. We use the combined names **LBBNN-GP-MF**, **LBBNN-GP-MVN**, and **LBBNN-GP-LFMVN** respectively to denote combination of model, prior and variational distribution. In Section 4 of the supplementary materials to the paper, results
from our early preprint (Hubin and Storvik, 2019) for MNIST and FMNIST datasets where the hyperparameters are fixed (corresponding to the setting also considered by Bai et al., 2020) are reported.

\begin{table}
\begin{tabular}{l|l} \hline
**Model** & **Meaning** \\ \hline BNN & Bayesian neural network \\ LBBNN & Latent binary Bayesian neural network \\ \hline
**Parameters prior** & **Meaning** \\ \hline GP & Independent Gaussian priors for weights \\ MGP & Independent mixture of Gaussians prior for weights \\ HP & Independent horseshoe priors for weights \\ \hline
**Inference** & **Meaning** \\ \hline MF & Mean-field variational inference \\ MVN & Multivariate Gaussian structure for the inclusion probabilities \\ LFMVN & Low factor for the covariance of MVN structure for the inclusion probabilities \\ \hline
**\(\mathbf{\gamma}\)** & **Meaning** \\ \hline SIM & Inclusion of the weights is drawn from the posterior of inclusion indicators \\ ALL & All weights are used \\ MED & Weights corresponding to the median probability model are used \\ PRN & Not pruned using a threshold-based rule weights are used \\ \hline
**\(\mathbf{\beta}\)** & **Meaning** \\ \hline SIM & The included weights are drawn from their posterior \\ MEA & Posterior means of the weights are used \\ \hline \(R\) & **Meaning** \\ \hline
10 & 10 samples are drawn \\
1 & 1 sample is drawn or posterior means are used \\ \hline
**Evaluation metric** & **Meaning** \\ \hline All cl Acc & Accuracy computed for all samples in the test set \\
0.95 threshold Acc & Accuracy computed for those samples in the test set where the maximum \\  & (across classes) model averaged predictive posterior exceeds 0.95 \\
0.95 threshold & Number of samples in the test set where the maximum \\ Num.cl & (across classes) model averaged predictive posterior exceeds 0.95 \\ Dens. level & Fraction of weights that are used to make predictions \\ Epo. time & Average time elapsed per epoch of training \\ \hline \end{tabular}
\end{table}
Table 1: Specifications of diagonal elements of \(\mathbf{A}\) matrices for the step sizes of optimization routines for LBBNN-GP-MF and LBBNN-GP-MVN, see Table 2 for explanation of the abbreviations. Note that \(A_{\omega}\) is only used in LBBNN-GP-MF, while \(A_{\xi}\) and \(A_{\Sigma}\) are only used in LBBNN-GP-MVN. For tuning parameters of LBBNN-GP-LFMVN, see Table S-1 in the supplementary materials to the paper.

\begin{table}
\begin{tabular}{l c c c c c c} \hline  & \(A_{\beta}\), \(A_{\rho}\) & \(A_{\xi}\) & \(A_{\omega}\) & \(A_{\Sigma}\) & \(A_{a_{\omega}}\),\(A_{b_{\omega}}\) & \(A_{a_{\beta}}\),\(A_{b_{\beta}}\) \\ \hline Pre-training & 0.00010 & 0.10000 & 0.10000 & 0.10000 & 0.00100 & 0.00001 \\ Training & 0.00010 & 0.01000 & 0.00010 & 0.00010 & 0.00000 & 0.00000 \\ Post-training & 0.00010 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 \\ \hline \end{tabular}
\end{table}
Table 1: Specifications of diagonal elements of \(\mathbf{A}\) matrices for the step sizes of optimization routines for LBBNN-GP-MF and LBBNN-GP-MVN, see Table 2 for explanation of the abbreviations. Note that \(A_{\omega}\) is only used in LBBNN-GP-MF, while \(A_{\xi}\) and \(A_{\Sigma}\) are only used in LBBNN-GP-MVN. For tuning parameters of LBBNN-GP-LFMVN, see Table S-1 in the supplementary materials to the paper.
In addition, we also used several relevant _baselines_. In particular, we addressed a standard Dense BNN with Gaussian priors and mean-field variational inference (Graves, 2011), denoted as **BNN-GP-MF**, which can be seen as a special case of our original model with all \(\gamma_{kj}^{(l)}\) being fixed and equal to 1 and no prior on the variance components of weights. This model is important in measuring how predictive power is changed due to introducing sparsity. Furthermore, we report the results for a Dense BNN with mixture priors (**BNN-MGP-MF**) with two Gaussian components of the mixtures (Blundell et al., 2015) with probabilities of 0.5 for each and variances equal to 1 and \(e^{-6}\) correspondingly. Additionally, we have addressed two popular sparsity-inducing approaches, in particular, a dense network with Concrete dropout (**BNN-GP-CMF**) (Gal et al., 2017) and a dense network with Horseshoe priors (**BNN-HP-MF**) (Louizos et al., 2017). Finally, a frequentist fully connected neural network (**FNN**) (with posthoc weight pruning) was used as a more basic baseline. We only report the results for FNN in the supplementary materials to make the experimental design cleaner. All of the baseline methods (including the FNN) also have 2 hidden layers with 400, and 600 neurons correspondingly. They were trained for 250 epochs with an Adam optimizer (with a learning rate \(a=0.0001\) for all involved parameters) and a batch size equal to 100. For the BNN with Horseshoe priors, we are reporting statistics separately before and after ad-hoc pruning (PRN) of the weights. Post-training (when necessary) was performed for additional 50 epochs. For FNN, for all three experiments, we performed weight and neuron pruning (Blalock et al., 2020) to have the same sparsity levels as those obtained by the Bayesian approaches to make them directly comparable. Pruning of FNN was based on removing the corresponding share of weights/neurons having the smallest magnitude (absolute value). No uncertainty was taken into consideration and neither was structure learning considered for FNNs.

For prediction, several methods were described in Section 4.3. All essentially boils down to choices on how to treat the model parameters \(\mathbf{\gamma}\) and the weights \(\mathbf{\beta}\). For \(\mathbf{\gamma}\), we can either simulate (**SIM**) from the (approximate) posterior or use the median probability model (**MED**). An alternative for **BNN-HP-MF** here is the pruning method (**PRN**) applied in Louizos et al. (2017). We also consider the choice of including all weights for some of the baseline methods (**ALL**). For \(\mathbf{\beta}\), we consider either sampling from the (approximate) posterior (**SIM**) or using the posterior mean (**MEA**). Under this notation, the fully Bayesian model averaging from Section 4.3 is denoted as **SIM SIM**, whilst the posterior mean based model as **ALL MEA**, the median probability model as **MED SIM**, the and the median probability model combined with parameter posterior mean as **MED MEA**.

We then evaluated accuracies (**Acc** - the proportion of the correctly classified images). Accuracies based on the median probability model (through either \(R=1\) or \(R=10\)) and the posterior mean models were also obtained. Finally, accuracies based on post-training of the parameters with fixed marginal inclusion probabilities and post-training of the median probability model were evaluated. For the cases when model averaging is addressed (\(R=10\)), we are additionally reporting accuracies when classification is only performed if 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, we are reporting the overall density level (the fraction of \(\gamma_{kj}^{(l)}\)'s equal to one within at least one of the simulations), for different approaches. To guarantee reproducibility, summaries (medians, minimums, maximums) across 10 independent runs of the described experiment \(s\in\{1,...,10\}\) were computed for all of these statistics. Estimates of the marginal inclusion probabilities \(\hat{p}(\gamma_{kj}^{(l)}=1|\mathcal{D})\) 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 each layer \(l\) as \(\rho(\gamma^{(l)}|\mathcal{D}):=\frac{1}{p^{(l+1)}p^{(l)}}\sum_{kj}\hat{p}( \gamma_{kj}^{(l)}=1|\mathcal{D})\), summarized in Table 6. To make the abbreviations used in the reported results more clear, we provide Table 2 with their short summaries.

MnistThe results reported in Tables 3 and 6 (with some additional results on LBBNN-GP-LFMVN and post-training reported in Tables S-2 and S-8 in the supplementary material) show that within our LBBNN approach: a) model averaging
across different BNNs (\(R=10\)) gives significantly higher accuracy than the accuracy of a random individual BNN from the model space (\(R=1\)); b) the median probability model and posterior mean based model also perform significantly better than a randomly sampled model. The performance of the median probability model and posterior mean-based model is in fact on par with full model averaging; c) according to Table 6 and Figure 1, for the mean-field variational distribution, the majority of the weights of the models have very low marginal inclusion probabilities for the weights at layers 1 and 2, while more weights have high marginal inclusion probabilities at layer 3 (although a significant reduction also at this layer). This resembles 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; d) for the MVN with full rank structure within variational approximation, the input layer is the most dense, followed by extreme sparsification in the second layer and a moderate sparsification at layer 3; e) the MVN approach with a low factor parametrization of the covariance matrix (results in the supplementary) only provides very moderate sparsification not exceeding 50% of the weight parameters; f) variations of all of the performance metrics across simulations are low, showing stable behavior across the repeated experiments; g) 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 h) both using the posterior mean-based model and using sample averaging improves accuracy compared to a single sample from the parameter space; i) variability in the estimates of the target parameters is 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 j) dense approaches outperform sparse approaches in terms of the accuracy in general; k) Concrete dropout marginally outperforms other approaches in terms of median accuracy, however, it exhibits large variance, whilst our full BNN and the compressed BNN with horseshoe priors yield stable performance across experiments; l) neither our approach nor baselines managed to reach state of the art results in terms of hard classification accuracy of predictions (Palvanov and Im Cho, 2018); m) 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 a 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; n) this might mean that the thresholds need to be calibrated towards the specific methods; o) our approach under the mean-field variational approximation and the full rank MVN structure of variational approximation yields the highest sparsity of weights when using the median probability model. Also, q) post-training (results in the supplementary) does not seem to significantly improve either the predictive quality of the models or uncertainty handling; o) all BNN for all considered sparsity levels on a given configuration of the network depth and widths are significantly outperforming the frequentist counterpart (with the corresponding same sparsity levels) in terms of the generalization error. Finally, in terms of computational time, r) as expected FNNs were the fastest in terms of time per epoch, whilst for the Bayesian approaches we see a strong positive correlation between the number of parameters and computational time, where BNN-GP-CMF is the fastest method and LBBNN-GP-MVN is the slowest. All times were obtained whilst training our models on a GeForce RTX 2080 Ti GPU card. Having said that, it is important to notice that the speed difference between the fastest and slowest Bayesian approach is less than 3 times, which given the fact that the time is also influenced by the implementation of different methods and a potentially different load of the server when running the experiments might be considered quite a tolerable difference in practice.
\begin{table}
\begin{tabular}{c c c c|c c c c c} \hline \multicolumn{2}{c|}{Prediction} & \multicolumn{2}{c|}{Model-Prior-} & \multicolumn{2}{c}{All cl} & \multicolumn{2}{c}{0.95 threshold} & Dens. & Epo. \\ \(\boldsymbol{\gamma}\) & \(\boldsymbol{\beta}\) & Method & \(R\) & Acc & Acc & Num.cl & level & time \\ \hline SIM & SIM & LBBNN-GP-MF & 1 & 0.968 (0.966,0.970) & - & - & 0.090 & 8.363 \\ SIM & SIM & LBBNN-GP-MF & 10 & 0.981 (0.979,0.982) & 0.999 & 8322 & 1.000 & 8.363 \\ ALL & MEA & LBBNN-GP-MF & 1 & 0.981 (0.980,0.983) & - & - & 1.000 & 8.363 \\ MED & SIM & LBBNN-GP-MF & 1 & 0.969 (0.968,0.974) & - & - & 0.079 & 8.363 \\ MED & SIM & LBBNN-GP-MF & 10 & 0.980 (0.979,0.982) & 0.999 & 8444 & 0.079 & 8.363 \\ MED & MEA & LBBNN-GP-MF & 1 & 0.981 (0.980,0.983) & - & - & 0.079 & 8.363 \\ \hline SIM & SIM & LBBNN-GP-MVN & 1 & 0.965 (0.964,0.966) & - & - & 0.180 & 9.651 \\ SIM & SIM & LBBNN-GP-MVN & 10 & 0.978 (0.976,0.979) & 1.000 & 7818 & 1.000 & 9.651 \\ ALL & MEA & LBBNN-GP-MVN & 1 & 0.978 (0.976,0.980) & - & - & 1.000 & 9.651 \\ MED & SIM & LBBNN-GP-MVN & 1 & 0.968 (0.966,0.969) & - & - & 0.163 & 9.651 \\ MED & SIM & LBBNN-GP-MVN & 10 & 0.977 (0.975,0.979) & 1.000 & 7928 & 0.163 & 9.651 \\ MED & MEA & LBBNN-GP-MVN & 1 & 0.974 (0.972,0.976) & - & - & 0.163 & 9.651 \\ \hline ALL & SIM & BNN-GP-MF & 1 & 0.965 (0.965,0.966) & - & - & 1.000 & 5.094 \\ ALL & SIM & BNN-GP-MF & 10 & 0.984 (0.982,0.985) & 0.999 & 8477 & 1.000 & 5.094 \\ ALL & MEA & BNN-GP-MF & 1 & 0.984 (0.982,0.985) & - & - & 1.000 & 5.094 \\ \hline ALL & SIM & BNN-MGP-MF & 1 & 0.965 (0.964,0.967) & - & - & 1.000 & 5.422 \\ ALL & SIM & BNN-MGP-MF & 10 & 0.982 (0.981,0.983) & 0.999 & 8329 & 1.000 & 5.422 \\ ALL & MEA & BNN-MGP-MF & 1 & 0.983 (0.981,0.984) & - & - & 1.000 & 5.422 \\ \hline SIM & SIM & BNN-GP-CMF & 1 & 0.982 (0.894,0.984) & - & - & 0.226 & 3.477 \\ SIM & SIM & BNN-GP-CMF & 10 & 0.984 (0.896,0.986) & 0.995 & 9581 & 1.000 & 3.477 \\ ALL & MEA & BNN-GP-CMF & 1 & 0.984 (0.893,0.986) & - & - & 1.000 & 3.477 \\ \hline SIM & SIM & BNN-HP-MF & 1 & 0.964 (0.962,0.967) & - & - & 1.000 & 4.254 \\ SIM & SIM & BNN-HP-MF & 10 & 0.982 (0.981,0.983) & 1.000 & 0003 & 1.000 & 4.254 \\ ALL & MEA & BNN-HP-MF & 1 & 0.966 (0.963,0.968) & - & - & 1.000 & 4.254 \\ PRN & SIM & BNN-HP-MF & 1 & 0.965 (0.962,0.969) & - & - & 0.194 & 4.254 \\ PRN & SIM & BNN-HP-MF & 10 & 0.982 (0.981,0.983) & 1.000 & 0002 & 0.194 & 4.254 \\ PRN & MEA & BNN-HP-MF & 1 & 0.965 (0.963,0.968) & - & - & 0.194 & 4.254 \\ \hline \end{tabular}
\end{table}
Table 3: Performance metrics for the MNIST data for the compared approaches. All results are medians across 10 repeated experiments (with min and max included in parentheses). No post-training is used. For further details see Table 2.
this sense compared to the previous example. For FNN, the same conclusions as those obtained for the MNIST data set are valid.

PhonemeFinally, the same set of approaches, model specifications (except for having 256 input covariates and 5 classes of the responses), and tuning parameters of the algorithms as in the MNIST and FMNIST examples were used for the classification of PHONEME data. The results a)- r) for the PHONEME data, based on Tables 5, 6, and Tables S-4 and S-10 in the supplementary are also overall consistent with the results from the MNIST and FMNIST experiments, however, predictive performances for all of the approaches are better than on FMNIST yet poorer than on MNIST. All of the methods, where sparsifications are possible, gave a lower sparsity level for this example. Yet, rather considerable sparsification is still shown to be feasible. For FNN, the same conclusions as those obtained for MNIST and FMNIST data sets are valid, though the deterioration of performance of FNN here was less drastic. Also, as we demonstrate in Figures S-5 - S-8 in the supplementary materials, the conclusions are consistent across various width configurations of Bayesian

\begin{table}
\begin{tabular}{l l l l|c c c c c} \hline \multicolumn{2}{c|}{Prediction} & \multicolumn{2}{c|}{Model-Prior-} & \multicolumn{2}{c}{All cl} & 0.95 threshold & Dens. & Epo. \\ \(\boldsymbol{\gamma}\) & \(\boldsymbol{\beta}\) & Method & \(R\) & Acc & Acc & Num.cl & level & time \\ \hline SIM & SIM & LBBNN-GP-MF & 1 & 0.864 (0.861,0.866) & - & - & 0.120 & 7.969 \\ SIM & SIM & LBBNN-GP-MF & 10 & 0.883 (0.881,0.886) & 0.995 & 4946 & 1.000 & 7.969 \\ ALL & MEA & LBBNN-GP-MF & 1 & 0.882 (0.879,0.887) & - & - & 1.000 & 7.969 \\ MED & SIM & LBBNN-GP-MF & 1 & 0.867 (0.864,0.871) & - & - & 0.108 & 7.969 \\ MED & SIM & LBBNN-GP-MF & 10 & 0.883 (0.880,0.886) & 0.995 & 5025 & 0.108 & 7.969 \\ MED & MEA & LBBNN-GP-MF & 1 & 0.880 (0.877,0.886) & - & - & 0.108 & 7.969 \\ \hline SIM & SIM & LBBNN-GP-MVN & 1 & 0.858 (0.854,0.859) & - & - & 0.156 & 9.504 \\ SIM & SIM & LBBNN-GP-MVN & 10 & 0.879 (0.874,0.880) & 0.995 & 4503 & 1.000 & 9.504 \\ ALL & MEA & LBBNN-GP-MVN & 1 & 0.875 (0.873,0.876) & - & - & 1.000 & 9.504 \\ MED & SIM & LBBNN-GP-MVN & 1 & 0.865 (0.860,0.866) & - & - & 0.129 & 9.504 \\ MED & SIM & LBBNN-GP-MVN & 10 & 0.877 (0.875,0.879) & 0.995 & 4694 & 0.129 & 9.504 \\ MED & MEA & LBBNN-GP-MVN & 1 & 0.871 (0.868,0.875) & - & - & 0.129 & 9.504 \\ \hline ALL & SIM & BNN-GP-MF & 1 & 0.864 (0.863,0.866) & - & - & 1.000 & 5.368 \\ ALL & SIM & BNN-GP-MF & 10 & 0.893 (0.890,0.894) & 0.997 & 5089 & 1.000 & 5.368 \\ ALL & MEA & BNN-GP-MF & 1 & 0.886 (0.882,0.888) & - & - & 1.000 & 5.368 \\ \hline ALL & SIM & BNN-MGP-MF & 1 & 0.867 (0.866,0.868) & - & - & 1.000 & 4.803 \\ ALL & SIM & BNN-MGP-MF & 10 & 0.893 (0.892,0.897) & 0.996 & 5151 & 1.000 & 4.803 \\ ALL & MEA & BNN-MGP-MF & 1 & 0.888 (0.885,0.890) & - & - & 1.000 & 4.803 \\ \hline SIM & SIM & BNN-GP-CMF & 1 & 0.896 (0.820,0.902) & - & - & 0.094 & 3.369 \\ SIM & SIM & BNN-GP-CMF & 10 & 0.897 (0.823,0.901) & 0.942 & 8825 & 1.000 & 3.369 \\ ALL & MEA & BNN-GP-CMF & 1 & 0.896 (0.821,0.901) & - & - & 1.000 & 3.369 \\ \hline SIM & SIM & BNN-HP-MF & 1 & 0.864 (0.863,0.869) & - & - & 1.000 & 4.613 \\ SIM & SIM & BNN-HP-MF & 10 & 0.887 (0.886,0.889) & 1.000 & 0181 & 1.000 & 4.613 \\ ALL & MEA & BNN-HP-MF & 1 & 0.867 (0.861,0.868) & - & - & 1.000 & 4.613 \\ PRN & SIM & BNN-HP-MF & 1 & 0.865 (0.860,0.868) & - & - & 0.302 & 4.613 \\ PRN & SIM & BNN-HP-MF & 10 & 0.887 (0.884,0.888) & 1.000 & 0179 & 0.302 & 4.613 \\ PRN & MEA & BNN-HP-MF & 1 & 0.865 (0.862,0.869) & - & - & 0.302 & 4.613 \\ \hline \end{tabular}
\end{table}
Table 4: Performance metrics for the FMNIST data for the suggested in the article Bayesian approaches to BNN. For further details, see Table 2 and the caption of Table 3.
neural networks. Also, sparsification increases with increased width for all the methods, yet the growth in sparsity is not proportional to the growth of width.

Out-of-domain experimentsFollowing the example of measuring the in and out-of-domain uncertainty suggested in Nitarshan (2018), we will first look at the ability of the LBBNN-GP-MF 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 dataset 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 \(s=10\). As can be seen in Figure 2, the samples from LBBNN-GP-MF give highly confident predictions for the MNIST dataset with almost no variance in the samples from the posterior predictive distribution. At the same time, the out-of-domain 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

\begin{table}
\begin{tabular}{l l l l|c c c c c} \hline \hline \multicolumn{2}{c|}{Prediction} & \multicolumn{2}{c|}{Model-Prior-} & \multicolumn{2}{c}{All cl} & 0.95 threshold & Dens. & Epo. \\ \(\boldsymbol{\gamma}\) & \(\boldsymbol{\beta}\) & Method & \(R\) & Acc & Acc & Num.cl & level & time \\ \hline SIM & SIM & LBBNN-GP-MF & 1 & 0.913 (0.898,0.929) & - & - & 0.371 & 0.433 \\ SIM & SIM & LBBNN-GP-MF & 10 & 0.927 (0.923,0.933) & 0.992 & 690 & 1.000 & 0.433 \\ ALL & MEA & LBBNN-GP-MF & 1 & 0.925 (0.921,0.933) & - & - & 1.000 & 0.433 \\ MED & SIM & LBBNN-GP-MF & 1 & 0.923 (0.910,0.928) & - & - & 0.307 & 0.433 \\ MED & SIM & LBBNN-GP-MF & 10 & 0.925 (0.912,0.934) & 0.984 & 757 & 0.307 & 0.433 \\ MED & MEA & LBBNN-GP-MF & 1 & 0.925 (0.913,0.932) & - & - & 0.307 & 0.433 \\ \hline SIM & SIM & LBBNN-GP-MVN & 1 & 0.919 (0.911,0.927) & - & - & 0.255 & 0.505 \\ SIM & SIM & LBBNN-GP-MVN & 10 & 0.929 (0.927,0.935) & 0.995 & 649 & 1.000 & 0.505 \\ ALL & MEA & LBBNN-GP-MVN & 1 & 0.926 (0.918,0.931) & - & - & 1.000 & 0.505 \\ MED & SIM & LBBNN-GP-MVN & 1 & 0.925 (0.916,0.929) & - & - & 0.225 & 0.505 \\ MED & SIM & LBBNN-GP-MVN & 10 & 0.929 (0.925,0.933) & 0.995 & 668 & 0.225 & 0.505 \\ MED & MEA & LBBNN-GP-MVN & 1 & 0.924 (0.921,0.928) & - & - & 0.225 & 0.505 \\ \hline ALL & SIM & BNN-GP-MF & 1 & 0.915 (0.907,0.919) & - & - & 1.000 & 0.203 \\ ALL & SIM & BNN-GP-MF & 10 & 0.919 (0.900,0.929) & 0.966 & 834 & 1.000 & 0.203 \\ ALL & MEA & BNN-GP-MF & 1 & 0.917 (0.901,0.922) & - & - & 1.000 & 0.203 \\ \hline ALL & SIM & BNN-MGP-MF & 1 & 0.913 (0.910,0.925) & - & - & 1.000 & 0.208 \\ ALL & SIM & BNN-MGP-MF & 10 & 0.916 (0.912,0.926) & 0.969 & 833 & 1.000 & 0.208 \\ ALL & MEA & BNN-MGP-MF & 1 & 0.921 (0.914,0.926) & - & - & 1.000 & 0.208 \\ \hline SIM & SIM & BNN-GP-CMF & 1 & 0.879 (0.706,0.906) & - & - & 0.509 & 0.103 \\ SIM & SIM & BNN-GP-CMF & 10 & 0.922 (0.918,0.930) & 0.965 & 187 & 1.000 & 0.103 \\ ALL & MEA & BNN-GP-CMF & 1 & 0.873 (0.712,0.904) & - & - & 1.000 & 0.103 \\ \hline SIM & SIM & BNN-HP-MF & 1 & 0.921 (0.915,0.929) & - & - & 1.000 & 0.136 \\ SIM & SIM & BNN-HP-MF & 10 & 0.921 (0.915,0.926) & 0.895 & 019 & 1.000 & 0.136 \\ ALL & MEA & BNN-HP-MF & 1 & 0.921 (0.916,0.926) & - & - & 1.000 & 0.136 \\ PRN & SIM & BNN-HP-MF & 1 & 0.919 (0.909,0.926) & - & - & 0.457 & 0.136 \\ PRN & SIM & BNN-HP-MF & 10 & 0.919 (0.916,0.927) & 0.926 & 028 & 0.457 & 0.136 \\ PRN & MEA & BNN-HP-MF & 1 & 0.920 (0.914,0.926) & - & - & 0.457 & 0.136 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Performance metrics for the PHONEME data for the suggested in the article Bayesian approaches to BNN. For further detail see Table 2 and the caption of Table 3.
Figure 1: An illustration of histograms of the marginal inclusion probabilities of the weights for the three layers (from top to bottom) of LBBNN-GP-MF from simulation \(s=10\) for MNIST (left) and FMNIST (right).
the posterior predictive distribution in this case. The reversed example of inference on FMNIST and uncertainty related to MNIST data, illustrated in Figure 3, leads to the same conclusions. More or less identical results were obtained for the LBBNN-GP-MVN and LBBNN-GP-LFMVN approaches, but they are not reported in the paper due to space constraints.

Figure 4 shows the results on more detailed out-of-domain experiments using FMNIST data for the models trained on MNIST data and vice versa. Following Louizos and 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, where the ideal is a CDF close to the lower right corner. Concrete drop-out is overconfident with a distribution of test classes being far from uniform, the horseshoe prior-based approach (both before and after pruning) is the closest to uniform (but it was

\begin{table}
\begin{tabular}{c c c} \hline \hline \multicolumn{2}{c}{**MNIST data**} & **FMNIST data** & **PHONEME data** \\ \hline
**LBBNN-GP-MF** & & \\ \hline \(\rho(\gamma^{(1)}=1|\mathcal{D})\) & 0.0844 (0.0835,0.0853) & 0.1323 (0.1291,0.1349) & 0.3806 (0.3764,0.3838) \\ \(\rho(\gamma^{(2)}=1|\mathcal{D})\) & 0.0959 (0.0942,0.0967) & 0.1005 (0.0981,0.1020) & 0.3670 (0.3641,0.3699) \\ \(\rho(\gamma^{(3)}=1|\mathcal{D})\) & 0.2945 (0.2808,0.3056) & 0.2790 (0.2709,0.2921) & 0.4236 (0.4053,0.4367) \\ \hline
**LBBNN-GP-MVN** & & \\ \hline \(\rho(\gamma^{(1)}=1|\mathcal{D})\) & 0.2975 (0.2928,0.2993) & 0.2461 (0.2410,0.2515) & 0.3201 (0.3142,0.3273) \\ \(\rho(\gamma^{(2)}=1|\mathcal{D})\) & 0.0368 (0.0363,0.0377) & 0.0392 (0.0383,0.0398) & 0.2287 (0.2235,0.2355) \\ \(\rho(\gamma^{(3)}=1|\mathcal{D})\) & 0.1394 (0.1311,0.1475) & 0.1462 (0.1368,0.1521) & 0.2763 (0.2632,0.2953) \\ \hline
**LBBNN-GP-LFMVN** & & \\ \hline \(\rho(\gamma^{(1)}=1|\mathcal{D})\) & 0.4474 (0.4448,0.4498) & 0.4589 (0.4565,0.4603) & 0.4973 (0.4965,0.4987) \\ \(\rho(\gamma^{(2)}=1|\mathcal{D})\) & 0.4525 (0.4501,0.4537) & 0.4516 (0.4493,0.4528) & 0.4972 (0.4952,0.4990) \\ \(\rho(\gamma^{(3)}=1|\mathcal{D})\) & 0.4815 (0.4685,0.4871) & 0.4805 (0.4654,0.4868) & 0.4979 (0.4925,0.5048) \\ \hline \hline \end{tabular}
\end{table}
Table 6: Medians and standard deviations of the average (per layer) marginal inclusion probability (see the text for the definition) for our model for both MNIST and FMNIST data across 10 repeated experiments.

Figure 2: Uncertainty related to the in-domain test data (MNIST, left) and out-of-domain test data (FMNIST, right) based on 100 samples from the posterior predictive distribution. Yellow lines are model-averaged posterior class probabilities (in percent). Green bars mark the correct classes, blue bars for other samples (with heights corresponding to an alternative estimate of class probabilities using hard classification within each of the replicates in the prediction procedure. The dashed black lines give the 95% threshold for making decisions with doubt possibilities). The original images are depicted to the left.
also closer to uniform for the in-domain predictions), whilst the 2 other baselines are in between, our approaches (both before pruning and after pruning with the median probability model) are on par with them showing that they handle out-of-domain uncertainty rather well.

Figure 4: Empirical CDF for the entropy of the marginal posterior predictive distributions trained on MNIST and applied to FMNIST (left) and vice versa (right) for simulation \(s=10\). Postfix S indicates the model sparsified by an appropriate method.

Figure 3: Uncertainty related to the in-domain test data (FMNIST, left) and out-of-domain test data (MNIST, right) based on 100 samples from the posterior predictive distribution. See Figure 2 for additional details.
More on misclassification uncertaintiesFigure 5 shows the misclassification uncertainties associated with posterior predictive sampling. One can see that for the majority of the cases when the LBBNN-GP-MF makes a misclassification, the class certainty of the predictions is relatively 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 misclassification illustrated here, even a human would have serious doubts about making a decision. Here, again, very similar results were obtained for LBBNN-GP-MVN and LBBNN-GP-LFMVN approaches, but they are not reported in the paper due to space constraints.

## 6 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 approximating the joint posterior of models and the parameters of these models. Approximate posterior predictive distributions, with both models and parameters marginalized out, can be easily obtained. Furthermore, marginal inclusion probabilities give proper probabilistic interpretation to Bayesian binary dropout and allow to perform model (or architecture) selection. This comes at the price of having only one additional parameter per weight included.

We provide image and sound 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. Regarding the computational costs of optimization: For the mean-field approximations, we are introducing only one additional parameter \(\alpha_{kj}^{l}\) for each weight. With underlying Gaussian structure on \(\boldsymbol{\alpha}^{(l)}\), additional parameters of the covariance matrix are further introduced. The complexity of each optimization step is proportional to the number of parameters to optimize, thus the deterioration in terms of computational time (as demonstrated in the experiments) is not at all drastic as compared to the fully connected BNN or even FNN. For the obtained predictions the complexities for different methods are proportional to the number of "active" parameters involved in predictions giving typically benefits to more sparse methods.

Regarding practical recommendations, we suggest, based on our empirical results, to use LBBNN-GP-MF if one is interested in a reasonable trade-off between sparsity, predictive accuracy, and uncertainty as well as computational costs. If sparsity is not needed standard BNN-GP-MF and BNN-MGP-MF are sufficient.

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 down-weighting 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 within simpler models (Smith and LeSage, 2004; Fahrmeir and Lang, 2001; Dobra et al., 2004). Obviously, in the context of inference in the joint parameter-model settings in BNNs, more research should be done on the choice of priors. Specifically, for image analysis, it might be of interest to develop convolution-inducing priors, whilst for recurrent models, one can think of exponentially decaying parameter priors for controlling the short-long memory.

In this work, we restrict ourselves to a subclass of BNNs, defined by the inclusion-exclusion 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
## Appendix A

Figure 5: Uncertainty based on the samples from the LBBNN-GP-MF model from the joint posterior (from simulation \(s=10\)) for 16 potentially wrongly classified (under model averaging) images for MNIST data (left) and FMNIST data (right). Yellow lines are model-averaged posterior class probabilities based on the Full BNN approach (in percent). Here, green bars are the true classes and red - the incorrectly predicted, blue bars - other samples, dashed black lines indicate the 95% threshold for making a decision when a doubt possibility is included. The original images are depicted to the left.
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), normalizing flows (Louizos and Welling, 2017), or address Jackknife to remove the bias (Nowozin, 2018). We leave these opportunities for further research.

The approach suggested in this paper can only implicitly perform width-selection of the architectures: Our approach allows to select the width under specific activation functions like ReLU as for a neuron with all weights switched off, i.e. if all \(\mathbf{\gamma}\)'s of a given neuron are put to 0, the unit is excluded from a layer thus reducing the width. The paper does not address the depth selection of neural networks, but a similar implicit approach would be possible for the depth as long as architectures with skip-connections are allowed (i.e. all layers are connected to the responses as well as the next layers of the networks). In such a case, it would be possible to have all nodes excluded in a specific layer of depth \(k\) making only the lower depth layers influence the responses. Then, the depth uncertainty could be inferred. Addressing the depth selection, thus, is an interesting possibility for follow-up research. Also, addressing depth and width selection more explicitly through assigning targeted priors as in Hubin et al. (2021) could be of interest in the future. But the ideas from Hubin et al. (2021) would impose more challenges for variational Bayes as compared to MCMC. At the same time, as shown in Hubin et al. (2021), such a procedure is much more likely to provide highly interpretable models. Also, further research on model priors will be needed in case of the explicit width-depth selection.

Last but not least, we would like to discuss some concurrent work that appeared while our paper was in the submission/review process. Firstly, two theoretical papers on posterior consistency of Bayesian variational deep learning in general and sparse contexts were published (Bhattacharya and Maiti, 2021; Cherief-Abdellatif, 2020). Also, Bai et al. (2020) justified theoretically (in an asymptotic setting) the choice of prior inclusion probabilities for the model and priors from Hubin (2018); Hubin and Storvik (2019). They used an almost identical variational inference technique as the one proposed in Hubin and Storvik (2019). The approach from (Hubin and Storvik, 2019) further recently found an application in genetic association studies (Cheng et al., 2022). Finally, Sun et al. (2022) used MCMC for inference on the model proposed in the early version of our work Hubin and Storvik (2019). These advancements show how rapidly the field develops and once again emphasize how actual and to date, the methodological developments on BNNs are in the field.

#### Acknowledgments

The authors would like to acknowledge Sean Murray (Norwegian Computing Center) for the comments on the language of the article and Dr. Pierre Lison (Norwegian Computing Center) for thoughtful discussions of the literature, potential applications, and technological tools. We also thank Dr. Petter Mostad, Department of Mathematical Sciences, The Chalmers University of Technology and the University of Gothenburg for valuable comments on Proposition 2. We also acknowledge constructive comments from the reviews and editorial comments we received at all stages of the publication of this article.