Bayesian sparsification for deep neural networks with Bayesian model
  reduction
Deep learning's immense capabilities are often constrained by the complexity
of its models, leading to an increasing demand for effective sparsification
techniques. Bayesian sparsification for deep learning emerges as a crucial
approach, facilitating the design of models that are both computationally
efficient and competitive in terms of performance across various deep learning
applications. The state-of-the-art -- in Bayesian sparsification of deep neural
networks -- combines structural shrinkage priors on model weights with an
approximate inference scheme based on stochastic variational inference.
However, model inversion of the full generative model is exceptionally
computationally demanding, especially when compared to standard deep learning
of point estimates. In this context, we advocate for the use of Bayesian model
reduction (BMR) as a more efficient alternative for pruning of model weights.
As a generalization of the Savage-Dickey ratio, BMR allows a post-hoc
elimination of redundant model weights based on the posterior estimates under a
straightforward (non-hierarchical) generative model. Our comparative study
highlights the advantages of the BMR method relative to established approaches
based on hierarchical horseshoe priors over model weights. We illustrate the
potential of BMR across various deep learning architectures, from classical
networks like LeNet to modern frameworks such as Vision Transformers and
MLP-Mixers.
Dimitrije Marković, Karl J. Friston, Stefan J. Kiebel
# Bayesian sparsification for deep neural networks with Bayesian model reduction

###### Abstract

Deep learning's immense capabilities are often constrained by the complexity of its models, leading to an increasing demand for effective sparsification techniques. Bayesian sparsification for deep learning emerges as a crucial approach, facilitating the design of models that are both computationally efficient and competitive in terms of performance across various deep learning applications. The state-of-the-art - in Bayesian sparsification of deep neural networks - combines structural shrinkage priors on model weights with an approximate inference scheme based on stochastic variational inference. However, model inversion of the full generative model is exceptionally computationally demanding, especially when compared to standard deep learning of point estimates. In this context, we advocate for the use of Bayesian model reduction (BMR) as a more efficient alternative for pruning of model weights. As a generalization of the Savage-Dickey ratio, BMR allows a post-hoc elimination of redundant model weights based on the posterior estimates under a straightforward (non-hierarchical) generative model. Our comparative study highlights the advantages of the BMR method relative to established approaches based on hierarchical horseshoe priors over model weights. We illustrate the potential of BMR across various deep learning architectures, from classical networks like LeNet to modern frameworks such as Vision Transformers and MLP-Mixers.

Bayesian model reduction, Stochastic variational inference, Deep neural networks

## 1 Introduction

Bayesian deep learning integrates the principles of Bayesian methodology with the objectives of deep learning, facilitating the training of expansive parametric models tailored for classifying and generating intricate audio-visual data, including images, text, and speech (Wang and Yeung, 2020; Wilson, 2020; Wang and Yeung, 2016). Notably, the Bayesian approach frames
the challenge of model optimization as an inference problem. This perspective is especially apt for scenarios necessitating decision-making under uncertainty (Murphy, 2022; Ghahramani, 2015). As a result, Bayesian formulations in deep learning have proven advantageous in various respects, offering enhancements in generalization (Wilson and Izmailov, 2020), accuracy, calibration (Izmailov et al., 2020; Luo and Kareem, 2020), and model compression (Louizos et al., 2017).

These functional enhancements are intrinsically tied to judiciously chosen structural priors (Fortuin, 2022). The priors, integral to the probabilistic generative model, scaffold the architecture of the network, thereby reducing the data required for the inference of optimal parametric solutions. Recent studies have highlighted the efficacy of hierarchical shrinkage priors over model weights, a specific category of structural priors, in achieving highly-sparse network representations (Nalisnick et al., 2019; Louizos et al., 2017; Seto et al., 2021; Ghosh et al., 2018). Sparse representations not only reduce redundancy but also evince additional performance benefits. However, the adoption of shrinkage priors in all deep learning models presents a conundrum: the ballooning space of latent parameters and the diminishing scalability of prevailing approximate inference schemes (Snoek et al., 2015; Krishnan et al., 2019; Izmailov et al., 2020; Daxberger et al., 2021).

In line with ongoing research on scalable Bayesian inference, we introduce an approximate inference scheme rooted in Bayesian model reduction (BMR). In essence, BMR extends the foundational principles of the Savage-Dickey Density Ratio method (Cameron, 2013). BMR is typically conceptualized as a combinatorial model comparison framework, enabling swift estimations of model evidence across an extensive array of models, that differ in their prior assumptions, to identify the most probable one. Originally conceived for model comparison within the dynamical causal modeling framework (Rosa et al., 2012; Friston and Penny, 2011), the scope of BMR has since broadened. Subsequent works expanded its methodology (Friston et al., 2016, 2017, 2018) and adapted it for structure learning (Smith et al., 2020). More recently, BMR has found applications in Bayesian nonlinear regression and classification tasks using Bayesian neural networks with variance backpropagation (Beckers et al., 2022; Haussmann et al., 2020).

The BMR method is intimately connected with the spike-and-slab prior, a type of shrinkage prior (Mitchell and Beauchamp, 1988). Intriguingly, this specific structured shrinkage prior has parallels with Dropout regularization (Nalisnick et al., 2019). Such an association spurred researchers in Bayesian deep learning to formulate sparsification methods based on a different type of shrinkage prior--the hierarchical horseshoe prior (Piironen and Vehtari, 2017)--as a tool for automated depth determination. Subsequent studies suggested that merging horseshoe priors with structured variational approximations yields robust, highly sparse representations (Ghosh et al., 2018). The allure of continuous shrinkage priors (e.g., horseshoe priors) stems from the computational challenges associated with model inversion reliant on spike-and-slab priors (Nalisnick et al., 2019; Piironen and Vehtari, 2017). However, continuous shrinkage priors necessitate a considerably more expansive parameter space, to represent the approximate posterior, compared to optimizing neural networks using the traditional point estimate method.

In this work, we reexamine the spike-and-slab prior within the framework of BMR-based sparsification, highlighting its efficiency. Notably, this approach circumvents the need to expand the approximate posterior beyond the conventional fully factorised mean-field ap
proximation, making it more scalable than structured variational approximations (Ghosh et al., 2018). In this light, BMR can be seen as a layered stochastic and black-box variational inference technique, which we term _stochastic BMR_. We subject the stochastic BMR to rigorous validation across various image classification tasks and network architectures, including LeNet-5 (LeCun et al., 1989), Vision Transformers (Dosovitskiy et al., 2020), and MLP-Mixers (Tolstikhin et al., 2021).

Central to our study is an empirical comparison of stochastic BMR with methods anchored in hierarchical horseshoe priors. Through multiple metrics - from Top-1 accuracy to expected calibration error and negative log-likelihood - we establish the competitive performance of stochastic BMR. We argue its computational efficiency, and remarkable sparsification rate, position BMR as an appealing choice, enhancing the scalability and proficiency of contemporary deep learning networks across diverse machine learning challenges, extending well beyond computer vision. We conclude with a discussion on potential avenues of future research that could further facilitate of BMR based pruning of deep neural networks.

## 2 Methods

In this section, we first describe the methods and techniques used in our research to address the problem of efficient Bayesian sparsification of deep neural networks. We provide a detailed overview of our approach, starting with variational inference methods, followed by the formulation of the Bayesian model reduction (BMR), Bayesian neural networks with shrinkage priors, and the description of corresponding approximate posterior.

### Variational inference

Given a joint density of latent variables, represented as \(\mathbf{z}=\left(z_{1},\ldots,z_{k}\right)\), and a dataset of \(n\) observations \(\mathbf{\mathcal{D}}=\left(y_{1},\ldots,y_{n}\right)\) we can express the joint density, that is, the generative model, as

\[p\left(\mathbf{\mathcal{D}},\mathbf{z}\right)=p\left(\mathbf{z}\right)p\left(\mathbf{\mathcal{ D}}|\mathbf{z}\right).\]

The posterior density is then obtained, following the Bayes rule, as

\[p\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\propto p\left(\mathbf{z}\right)p\left(\mathbf{ \mathcal{D}}|\mathbf{z}\right). \tag{1}\]

For complex generative models, direct inference as described above becomes computationally prohibitive. To circumvent this, we approximate the exact posterior \(p\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\), constraining it to a distribution \(q\left(z\right)\) that belongs to a named distribution family \(\mathcal{Q}\). We then seek \(q^{*}\left(z\right)\in\mathcal{Q}\), an approximate solution that minimizes the following Kullback-Leibler divergence (Blei et al., 2017)

\[q^{*}\left(z\right)=\underset{q\in\mathcal{Q}}{\text{argmin}}D_{KL}\left(q \left(\mathbf{z}\right)||p\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\right)=\underset{ q\in\mathcal{Q}}{\text{argmin}}F\left[q\right],\]

where \(F\left[q\right]\) stands for the variational free energy (VFE), defined as

\[F\left[q\right]=E_{q\left(\mathbf{z}\right)}\left[\ln q\left(\mathbf{z}\right)-\ln p \left(\mathbf{\mathcal{D}},\mathbf{z}\right)\right]\]

VFE serves as an upper bound on the marginal log-likelihood

\[F\left[q\right]=D_{KL}\left(q\left(\mathbf{z}\right)||p\left(\mathbf{z}|\mathbf{\mathcal{ D}}\right)\right)-\ln p\left(\mathbf{\mathcal{D}}\right)\geq-\ln p\left(\mathbf{ \mathcal{D}}\right)\]
As KL-divergence is always greater or equal to zero, minimizing VFE brings the approximate solution as close as possible to the true posterior, without having to compute the exact posterior.

The most straightforward way to obtain the approximate posterior \(q^{*}\left(\mathbf{z}\right)\), is to minimize the VFE along its negative gradient:

\[\dot{\mathbf{\phi}}=-\nabla_{\mathbf{\phi}}F\left[q\right]\]

where \(\mathbf{\phi}\) signifies the parameters of the approximate posterior \(q_{\mathbf{\phi}}\left(\mathbf{z}\right)=q\left(\mathbf{z}|\mathbf{\phi}\right)\). Thus, variational inference reframes the inference problem highlighted in eq. (1) as an optimization problem Beal (2003).

### Stochastic and black-box variational inference

_Stochastic variational inference_ (SVI) improves the computational efficiency of gradient descent by approximating the variational free energy using a subset--\(\mathbf{\mathcal{K}}_{i}=\left(y_{s_{1}^{i}},\ldots,y_{s_{k}^{i}}\right);\,k\ll n\)--of the entire data set \(\mathbf{\mathcal{D}}\). This approach fosters a stochastic gradient descent (SGD) mechanism, capable of managing large datasets (Hoffman et al., 2013). Crucially, at every iteration step \(i\) of the SGD process, the subset \(\mathbf{\mathcal{K}}_{i}\) undergoes re-sampling.

_Black-box Variational Inference_ (BBVI) facilitates the optimization of any (named or unnamed) posterior density \(q_{\mathbf{\phi}}\left(\mathbf{z}\right)\), through the integration of Monte Carlo estimates for variational gradients (Ranganath et al., 2014). This can be formulated as the following relation

\[\nabla_{\mathbf{\phi}}F\left[q\right]\approx\nabla_{\mathbf{\phi}}\tilde{F}\left[q \right]=\frac{1}{S}\sum_{s=1}^{S}\nabla_{\mathbf{\phi}}\ln q_{\mathbf{\phi}}\left(\bm {z}\right)\left[\ln\frac{q_{\mathbf{\phi}}\left(\mathbf{z}\right)}{p\left(\mathbf{ \mathcal{D}},\mathbf{z}\right)}+1\right];\quad\mathbf{z}_{s}\sim q\left(\mathbf{z}|\mathbf{ \phi}\right) \tag{2}\]

which is known as the REINFORCE estimator (Williams, 1992). To mitigate the variance inherent to Monte Carlo gradient estimations, we employ Rao-Blackwellization (Schulman et al., 2015), with an implementation sourced from NumPyro (Bingham et al., 2019). For optimizing the variational objective stochastically, we leverage the AdaBelief optimizer (Zhuang et al., 2020). As an adaptive algorithm, AdaBelief ensures swift convergence, robust generalization, and steady optimization. Notably, we use AdaBelief's implementation from the Optax package within the JAX ecosystem (Babuschkin et al., 2020).

### Bayesian model reduction

Let us consider two generative processes for the data: a full model

\[p\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\propto p\left(\mathbf{\mathcal{D}}|\mathbf{z} \right)p\left(\mathbf{z}\right)\]

and a reduced model 1 in which the original prior \(p\left(\mathbf{z}\right)\) is replaced with a more informative prior \(\tilde{p}\left(\mathbf{z}\right)=p\left(\mathbf{z}|\mathbf{\theta}\right)\) that depends on hyper-parameters \(\mathbf{\theta}\). This change leads to a different posterior

Footnote 1: the reduction here implies applying constraints of any form to the prior to obtain a posterior with reduced entropy.

\[\tilde{p}\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\propto p\left(\mathbf{\mathcal{D}}| \mathbf{z}\right)\tilde{p}\left(\mathbf{z}\right)\]
Noting that as the following relation holds:

\[1=\int\mathrm{d}\mathbf{z}\tilde{p}\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)=\int\mathrm{d }\mathbf{z}p\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\frac{\tilde{p}\left(\mathbf{z}\right)p \left(\mathbf{\mathcal{D}}\right)}{p\left(\mathbf{z}\right)\tilde{p}\left(\mathbf{ \mathcal{D}}\right)},\]

we can express the link between the models as:

\[\begin{split}-\ln\tilde{p}\left(\mathbf{\mathcal{D}}\right)& =-\ln p\left(\mathbf{\mathcal{D}}\right)-\ln\int d\mathbf{z}p\left(\mathbf{z}| \mathbf{\mathcal{D}}\right)\frac{\tilde{p}\left(\mathbf{z}\right)}{p\left(\mathbf{z} \right)}\\ &\approx F\left(\mathbf{\phi}^{*}\right)-\ln\int d\mathbf{z}q_{\mathbf{\phi}^ {*}}\left(\mathbf{z}\right)\frac{\tilde{p}\left(\mathbf{z}\right)}{p\left(\mathbf{z} \right)}\end{split} \tag{3}\]

where we assumed the approximate posterior for the full model corresponds to \(p\left(\mathbf{z}|\mathbf{\mathcal{D}}\right)\approx q_{\mathbf{\phi}^{*}}\left(\mathbf{z}\right)\), and that \(-\ln p\left(\mathbf{\mathcal{D}}\right)\approx F\left(\mathbf{\phi}^{*}\right)\).

From eq. (3) we obtain the free energy of the reduced model as

\[-\ln\tilde{p}\left(\mathbf{\mathcal{D}}\right)\approx-\ln E_{q}\left[\frac{\tilde {p}\left(\mathbf{z}\right)}{p\left(\mathbf{z}\right)}\right]+F\left(\mathbf{\phi}^{*} \right)=-\Delta F\left(\mathbf{\theta}\right). \tag{4}\]

where \(\Delta F\left(\mathbf{\theta}\right)\) denotes the change in the free energy of going from the full model to the reduced model, given hyper-parameters \(\mathbf{z}_{H}\). Note that for \(\Delta F\left(\mathbf{\theta}\right)>0\) the reduced model has a better variational free energy compared to the flat model. Consequently, the reduced model offers a model with a greater marginal likelihood; i.e., a better explanation for the data and improved generalization capabilities. Heuristically, this can be understood as minimising model complexity, without sacrificing accuracy (because log evidence can be expressed as accuracy minus complexity, where complexity is the KL divergence between posterior and prior beliefs). This relationship is pivotal in formulating efficient pruning criteria, especially for extensive parametric models commonly employed in deep learning.

### Bayesian neural networks

In a general (nonlinear) regression problem, we model the relationship between predictors \(\mathbf{X}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\right)\) and target variables \(\mathbf{Y}=\left(\mathbf{y}_{1},\ldots,\mathbf{y}_{n}\right)\) using a likelihood distribution from an exponential family as

\[\mathbf{y}_{i}\sim p\left(\mathbf{y}|\mathbf{\mathcal{W}},\mathbf{x}_{i}\right)=h(y)\exp \left[\mathbf{\eta}\left(\mathbf{f}\left(\mathbf{\mathcal{W}},\mathbf{x}_{i}\right)\right) \cdot\mathbf{T}\left(\mathbf{y}\right)-A\left(\mathbf{f}\left(\mathbf{\mathcal{W}},\mathbf{x}_{i} \right)\right)\right]. \tag{5}\]

Functions \(h(\cdot),\mathbf{\eta}(\cdot),\mathbf{T}(\cdot),A(\cdot)\) are known and selected depending on the task. For example in a regression problem the likelihood will correspond to a multivariate normal distribution and in a classification problem to a categorical distribution. In this work, we will only consider a categorical likelihood, as it is the most suitable for image classification tasks.

The mapping \(\mathbf{f}\left(\mathbf{\mathcal{W}},\mathbf{x}_{i}\right)\) represents a generic deep neural network of depth \(L\) defined as

\[\mathbf{\mathcal{W}} =\left(\mathbf{W}_{1},\ldots,\mathbf{W}_{L}\right)\] \[\mathbf{h}_{i}^{0} =\mathbf{x}_{i}\] \[\mathbf{h}_{i}^{l} =\mathbf{g}\left(\mathbf{W}_{l}\cdot\left[\mathbf{h}_{i}^{l-1};1\right]\right)\] \[\mathbf{f}\left(\mathbf{\mathcal{W}},\mathbf{x}_{i}\right) =\mathbf{W}_{L}\cdot\left[\mathbf{h}_{i}^{L-1};1\right]\]
A probabilistic formulation of the deep learning task, that is, inferring model weights, introduces implicit bias to the parameters \(\boldsymbol{\mathcal{W}}\) of an artificial neural network in the form of a prior distribution \(p\left(\boldsymbol{\mathcal{W}}\right)\). Hence, parameter estimation is cast as an inference problem where

\[p\left(\boldsymbol{\mathcal{W}}|\boldsymbol{\mathcal{D}}\right)\propto p\left( \boldsymbol{\mathcal{W}}\right)\prod_{i=1}^{n}p\left(\boldsymbol{y}_{i}| \boldsymbol{\mathcal{W}},\boldsymbol{x}_{i}\right)\]

The choice of the prior distribution is crucial for optimal task performance, and a prior assumption of structural sparsity is essential for inferring sparse representations of over-parameterised models, such as deep neural networks.

### Bayesian neural networks with shrinkage priors

Shrinkage priors instantiate a prior belief about the sparse structure of model parameters. Here, we will investigate two well-established forms of shrinkage priors for network weight parameters, a canonical spike-and-slab prior (George and McCulloch, 1993; Mitchell and Beauchamp, 1988) defined as

\[w_{ijl} \sim\mathcal{N}\left(0,\lambda_{ijl}^{2}\gamma_{0}^{2}\right)\] \[\lambda_{ijl} \sim\text{Bernoulli}\left(\pi_{l}\right)\] \[\pi_{l} \sim\mathcal{B}e\left(\alpha_{0},\beta_{0}\right)\]

and a regularised-horseshoe prior (Piironen and Vehtari, 2017)

\[w_{ijl} \sim\mathcal{N}\left(0,\gamma_{il}^{2}\right)\] \[\gamma_{il}^{2} =\frac{c_{l}^{2}v_{l}^{2}\tau_{il}^{2}}{c_{l}^{2}+\tau_{il}^{2}v_ {l}^{2}} \tag{6}\] \[c_{l}^{-2} \sim\Gamma\left(2,6\right)\] \[\tau_{il} \sim\mathcal{C}^{+}(0,1)\] \[v_{l} \sim\mathcal{C}^{+}(0,\tau_{0})\]

where \(i\in\{1,\ldots,K_{l}\}\), \(j\in[1,\ldots,K_{l-1}+1]\), and where \(w_{ijl}\) denotes \(ij\)th element of the weight matrix at depth \(l\). The symbols \(\mathcal{B}e\), and \(\mathcal{C}^{+}\) denote a Beta distribution and a half-Cauchy distribution, respectively.

Importantly, the spike-and-slab prior relates to dropout regularisation, which is commonly introduced as a sparsification method in deep learning (Nalisnick et al., 2019; Mobiny et al., 2021). This type of prior is considered the gold standard in shrinkage priors and has been used in many recent applications of Bayesian sparsification on neuronal networks (Bai et al., 2020; Hubin and Storvik, 2023; Jantre et al., 2021; Sun et al., 2022; Ke and Fan, 2022) showing excellent sparsification rates. However, the inversion of the resulting hierarchical model is challenging and requires carefully constructed posterior approximations. Moreover, their dependence on discrete random variables renders them unsuitable for Markov-Chain Monte Carlo-based sampling schemes. As a result, researchers often use continuous formulations of the shrinkage-prior, with the horseshoe prior being a notable example.
In contexts that involve sparse learning with scant data, the regularised horseshoe prior has emerged as one of the preferred choices within shrinkage prior families (Ghosh et al., 2019). A distinct advantage of this prior is its ability to define both the magnitude of regularisation for prominent coefficients and convey information about sparsity. It is worth noting a dependency highlighted in Ghosh et al. (2018): for \(v_{l}\tau_{il}\ll 1\) the equation simplifies to \(\gamma_{il}\approx v_{l}\tau_{il}\) recovering the original horseshoe prior. In contrast, for \(v_{l}\tau_{il}\gg 1\), the equation becomes \(\gamma_{il}^{2}\approx c_{l}^{2}\). In this latter scenario, the prior over the weights is defined as \(w_{ijl}\sim\mathcal{N}\left(0,c_{l}^{2}\right)_{i}\), with \(c_{l}\) serving as a weight decay hyper-parameter for layer \(l\).

### Approximate posterior for Bayesian neural networks

To benchmark stochastic BMR, we explore two forms of prior distribution \(p\left(\boldsymbol{\mathcal{W}}\right)\)--a flat and a hierarchical structure--in conjunction with a fully factorised mean-field approximation.

Firstly, let us consider the flat prior over model weights, represented in a non-centered parameterization:

\[\begin{split} c_{l}^{-2}&\sim\Gamma(2,2)\\ \hat{w}_{ijl}&\sim\mathcal{N}\left(0,1\right)\\ w_{ijl}&=\gamma_{0}c_{l}\hat{w}_{ijl}\end{split} \tag{7}\]

where we set \(\gamma_{0}=0.1\). Note that in the flat prior we incorporate a layer specific scale parameter, which we found to stabilise variational inference. Based on this, we describe a fully factorised approximate posterior as a composite of Normal and Log-Normal distributed random variables. Hence,

\[\begin{split} q\left(\boldsymbol{\hat{\mathcal{W}}},\boldsymbol {c}\right)&=\prod_{l}q\left(c_{l}^{-2}\right)\prod_{i}\prod_{j}q \left(\hat{w}_{ijl}\right)\\ q\left(\hat{w}_{ijl}\right)&=\mathcal{N}\left( \mu_{ijl},\sigma_{ijl}^{2}\right)\\ q\left(c_{l}^{-2}\right)&=\mathcal{L}\mathcal{N} \left(\mu_{c,l},\sigma_{c,l}^{2}\right).\end{split} \tag{8}\]

When inverting a hierarchical generative model over weights of artificial neural network, we exclusively apply stochastic black-box variational inference to the model variant with the regularised horseshoe prior. This choice is motivated by its documented superiority over the spike-and-slab prior, as established in Ghosh et al. (2018). We express the hierarchical prior
in the non-centered parameterization as:

\[a_{il},b_{il} \sim\Gamma\left(\frac{1}{2},1\right)\] \[\hat{a}_{l},\hat{b}_{l} \sim\Gamma\left(\frac{1}{2},1\right)\] \[\tau_{il} =\sqrt{\frac{a_{il}}{b_{il}}}\] \[v_{l} =\tau_{0}\sqrt{\frac{\hat{a}_{l}}{\hat{b}_{l}}}\] \[\hat{w}_{ijl} \sim\mathcal{N}\left(0,1\right)\] \[w_{ijl} =\gamma_{il}\hat{w}_{ijl}\]

Note that the expressions above involve a reparameterization of Half-Cauchy distributed random variables as the square-root of the quotient of two Gamma distributed random variables, a strategy drawn from Wand et al. (2011) (see Appendix B for additional details). Such a reparameterization of the Half-Cauchy ensures capturing of fat-tails in the posterior, even when leveraging a fully-factorised mean-field posterior approximation, as referenced in Ghosh et al. (2018).

For the fully-factorised mean-field approximation, the approximate posterior is portrayed as a composite of Normal and Log-Normal distributed random variables, expressed as:

\[q\left(\hat{\mathbf{\mathcal{W}}},\mathbf{a},\mathbf{b},\mathbf{\hat{a}},\mathbf{ \hat{b}},\mathbf{c}\right) =\prod_{l}q\left(c_{l}^{-2}\right)q\left(\hat{a}_{l}\right)q \left(\hat{b}_{l}\right)\prod_{i}q\left(a_{il}\right)q\left(b_{il}\right)\prod _{j}q\left(\hat{w}_{ijl}\right)\] \[q\left(c_{l}\right) =\mathcal{LN}\left(\mu_{c,l},\sigma_{c,l}^{2}\right)\] \[q\left(\hat{a}_{l}\right) =\mathcal{LN}\left(\hat{\mu}_{a,l},\hat{\sigma}_{a,l}^{2}\right)\] \[q\left(\hat{b}_{l}\right) =\mathcal{LN}\left(\hat{\mu}_{b,l},\hat{\sigma}_{b,l}^{2}\right)\] \[q\left(a_{il}\right) =\mathcal{LN}\left(\mu_{a,il},\sigma_{a,il}^{2}\right)\] \[q\left(b_{il}\right) =\mathcal{LN}\left(\mu_{b,il},\sigma_{b,il}^{2}\right)\] \[q\left(\hat{w}_{ijl}\right) =\mathcal{N}\left(\mu_{w,ijl},\sigma_{w,ijl}^{2}\right)\]

### Application of stochastic BMR to Bayesian neural networks

To apply BMR to Bayesian neural networks, we commence by estimating an approximate posterior for the flat model, as detailed in eq. (7). To retain high computational efficiency, we pair BMR solely with the fully factorised approximate posterior, as presented in eq. (8). While it is feasible to use this method alongside the structured posterior (Ghosh et al., 2018), it requires considerably more computationally intensive estimations of the reduced free energy. As shown below, we obtain satisfactory results with a fully factorised posterior. Therefore, we defer the exploration of BMR with a structured posterior to future endeavours.
Given a fully factorised approximate posterior, we can determine the change in variational free energy, \(\Delta F\)--after substituting the prior \(\mathcal{N}\left(0,1\right)\) with \(\mathcal{N}\left(0,\theta_{ijl}^{2}\right)\) for the weight \(\hat{w}_{ijl}\)--as:

\[\Delta F\left(\theta_{ijl}\right) =-\frac{1}{2}\ln\rho_{ijl}^{2}-\frac{1}{2}\frac{\mu_{ijl}^{2}}{ \sigma_{ijl}^{2}}\left(1-\frac{\theta_{ijl}^{2}}{\rho_{ijl}^{2}}\right)\] \[\rho_{ijl}^{2} =\theta_{ijl}^{2}+\sigma_{ijl}^{2}-\theta_{ijl}^{2}\sigma_{ijl}^{2}\]

For the second hierarchical level of the approximate posterior, we aim to minimize the following form for the variational free energy:

\[F=\sum_{l=1}^{L}E_{q\left(\boldsymbol{\theta}_{l}\right)}\left[-\sum_{i,j} \Delta F(\theta_{ijl})+\ln\frac{q\left(\boldsymbol{\theta}_{l}\right)}{p\left( \boldsymbol{\theta}_{l}\right)}\right] \tag{9}\]

This minimization is done with respect to \(q\left(\boldsymbol{\Theta}\right)=\prod_{l}q\left(\boldsymbol{\theta}_{l}\right)\), the approximate posterior over hyper-parameters. Note the application of eq.4 in substituting the marginal log-likelihood with the change in the variational free energy.

For the spike-and-slab prior we can write the following relation:

\[\boldsymbol{\theta}_{l} =\left[\pi_{l},\lambda_{ijl}\right]\text{ for }i\in\left\{1,\ldots,K_{l}\right\}, \text{ and }j\in\left\{1,\ldots,K_{l-1}+1\right\}\right]\] \[\theta_{ijl} =\lambda_{ijl}\]

Consequently, the approximate posterior at the second level of the hierarchy can be approximated as:

\[q\left(\boldsymbol{\Theta}\right) =\prod_{l}q\left(\pi_{l}\right)\prod_{ij}q\left(\lambda_{ijl}\right)\] \[q\left(\lambda_{ijl}\right) =q_{ijl}^{\lambda_{ijl}}\left(1-q_{ijl}\right)^{1-\lambda_{ijl}}\] \[q\left(\pi_{l}\right) =\mathcal{B}\left(\alpha_{l},\beta_{l}\right)\]

The iterative update to obtain the minimum of the simplified variational free energy (eq.9) is then:

\[q_{ijl}^{k+1} =\frac{1}{1+e^{-\left[\zeta_{l}^{k}-\Delta F\left(\lambda_{ijl}= 0\right)\right]}}\] \[\zeta_{l}^{k} =\psi(\alpha_{l}^{k})-\psi(\beta_{l}^{k})\] \[\alpha_{l}^{k+1} =\sum_{i,j}q_{ijl}^{k+1}+\alpha_{0}\] \[\beta_{l}^{k+1} =\sum_{i,j}\left(1-q_{ijl}^{k+1}\right)+\beta_{0}\]

Here, \(\alpha_{l}^{0}=\alpha_{0}\), \(\beta_{l}^{0}=\beta_{0}\), \(\Delta F\left(\lambda_{ijl}=0\right)=-\frac{1}{2}\left[\ln\sigma_{ijl}^{2}+ \frac{\mu_{ijl}^{2}}{\sigma_{ijl}^{2}}\right]\), and \(\psi\left(\cdot\right)\) refers to the digamma function. The efficiency of this inference scheme is remarkable, typically achieving convergence after a few iterations. In practice, we cap the maximum number of iterations at \(k_{max}=4\).
Finally, we use the following pruning heuristics to eliminate model weights and sparsify network structure

\[\text{if }q_{ijl}^{k_{max}}<\frac{1}{2}\text{, set }\hat{w}_{ijl}=0.\]

To achieve the high sparsification rate presented in the next section, we adopt an iterative optimisation and pruning approach proposed in Beckers et al. (2022). We perform weight pruning at the beginning of each epoch (except the first one), and further optimisation for 500 iterations, completing one epoch. In total, we apply iterative pruning and optimisation for fifty epochs in all examples below.

The complete implementation of stochastic BMR is available at an online repository github.com/dimarkov/bmr4pml with notebooks and scripts necessary to recreate all result figures.

## 3 Results

In this section, we present the outcomes of our experiments and analyses conducted to evaluate the performance and efficiency of the stochastic Bayesian model reduction in the context of Bayesian sparsification of deep neural networks. Our results are structured to provide insights into the capabilities and advantages of our approach.

### Performance Comparison

The training regimen used a batch size of \(N_{B}=128\) and the AdaBelief algorithm with learning rate set to \(\alpha=10^{-3}\) in the case of the MAP estimate, \(\alpha=5\cdot 10^{-3}\) in the case of the mean-field methods, and \(\alpha=10^{-2}\) in the case of stochastic BMR (the exponential decay rates were kept at default values \(\beta_{1}=0.9\), and \(\beta_{2}=0.999\)). Figure 1 charts the

Figure 1: Classification performance comparison on FashionMNIST dataset for different neuronal architectures and approximate inference schemes.
epoch-wise evolution of ACC, ECE, and NLL for each architecture, under five distinct approximate inference strategies: (i) Maximum a posteriori (MAP) estimate for the flat generative model, akin to traditional deep learning point estimates coupled with weight decay. (ii) A fully factorised posterior approximation for the flat generative model (Flat-FF). (iii) A fully factorised posterior approximation of the hierarchical generative model with a regularised horseshoe prior (Tiered-FF). (iv) The stochastic BMR algorithm augmented with a spike-and-slab prior (BMR-S&S). Each epoch is defined by 500 stochastic gradient steps, with each step randomly drawing \(N_{B}\) data instances from the training pool.

Interestingly, all approximate inference methods demonstrate comparable top-1 accuracy scores. However, the stochastic BMR method followed by the Tiered-FF approximatio (with a single exception), consistently resulted in the lowest ECE and NLL scores across the majority of DNN architectures and datasets (see Figure S1 for CIFAR10 dataset and Figure S2 for CIFAR100 dataset).The implicit reduction in model complexity suggests that--as anticipated--Stochastic BMR furnishes a model of the data that has the greatest evidence or marginal likelihood (not shown). In this setting, the NLL of the test data can be regarded as a proxy for (negative log) marginal likelihood.

### Learning of sparse representations

Figure 2 depicts the fraction of pruned model parameters for different DNN architectures and datasets. It is noteworthy to observe the substantive sparsity achieved by the stochastic BMR algorithm. This sparsity is consistent across datasets and architectures, with the exception of the LeNet-5 structure when used for the FashionMNIST dataset, because by default LeNet-5 architecture is already sparse and contains relatively low-number of model weights (for other data sets we substantially increased the dimensionality of hidden layers as detailed in Appendix A).

To delve deeper into the pruning behavior across varying network depths, Figure 3 presents a per-layer cumulative distribution function (CDF) for model parameters, highlighting the proportion of parameters whose absolute mean posterior estimate falls below a given threshold. When juxtaposing the BMR CDF trajectories with those obtained from the Tiered-FF method (sparsification is induced by the regularised half-cauchy prior), it is evident that BMR furnishes more pronounced sparsification. This distinction is crucial,

Figure 2: Total fraction of pruned model parameters obtained with the stochastic BMR algorithm across different DNN architectures and datasets.
as the stochastic BMR not only matches or surpasses the performance of the Tiered-FF algorithm but also averages a 30% faster stochastic gradient descent.

To illustrate the structural learning variations among algorithms, Figure 4 presents heatmaps of posterior expectations obtained using the four different methods. The Figure 4 reveals subtle differences between inferred representations of the MLP and LeNet-5

Figure 4: Posterior expectations (color coded) over model parameters obtained using different approximate inference schemes at the first layer of (a) MLP architecture, and (b) LeNet architectures.

Figure 3: Cumulative Distribution Function (CDF) of absolute posterior parameter expectations at different layers of MLP (top row), and LeNet architectures (bottom row). The y-axis represents the fraction of parameters with values less than or equal to the value on the x-axis.
architecture's input layers trained on the Fashion MNIST dataset. Divergent compression rates among the algorithms indicate inherent trade-offs between efficiency and performance. It is evident that the stochastic BMR strikes a balance between compression advantages and performance, as it is less prone to over-pruning as compared to the Tiered-FF method (two featured of the LeNet-5 input layer are effective removed - see Figure 4(b)).

## 4 Discussion

In this study, we presented a novel algorithm--stochastic Bayesian model reduction--designed for an efficient Bayesian sparsification of deep neural networks. Our proposed method seamlessly integrates stochastic and black-box variational inference with Bayesian model reduction (BMR), a generalisation of the Savage-Dickey ratio. Through the stochastic BMR strategy, we enable iterative pruning of model parameters, relying on posterior estimates acquired from a straightforward variational mean-field approximation to the generative model. This model is characterized by Gaussian priors over individual parameters and layer-specific scale parameters. The result is an efficient pruning algorithm for which the computational demand of the pruning step is negligible compared to the direct stochastic black-box optimization of the full hierarchical model.

Over recent years, the Bayesian sparsification of neural networks has gained momentum, primarily driven by the spike-and-slab prior Bai et al. (2020); Hubin and Storvik (2023); Jantre et al. (2021); Sun et al. (2022); Ke and Fan (2022). These works have showcased the remarkable sparsification capabilities inherent to such shrinkage priors. Nevertheless, when juxtaposed with the stochastic BMR algorithm, they often necessitate supplementary assumptions related to the approximate posterior. These assumptions, in turn, lead to a more computation-intensive model inversion. Moreover, in contrast to related approaches, the versatility of stochastic BMR allows its integration with more efficient optimization techniques, like variational Laplace Daxberger et al. (2021) and proximal-gradient methods Khan et al. (2018), provided the resulting approximate posterior in the form of a normal distribution is apt for the application at hand.

The insights obtained here pave the way for a deeper exploration of the potential applications of Bayesian model reduction across a wider array of architectures and tasks in probabilistic machine learning, such as audiovisual and natural language processing tasks. A more detailed fine tuning of the core dynamics of these algorithms, in terms of iterations steps, learning rates, and other free-parameters, might be the key to unveiling even more proficient Bayesian deep learning methodologies in the near future.

We thank Conor Heins, Magnus Koudahl, and Beren Millidge for valuable discussions during the initial stages of this work. SK acknowledges support by DFG TRR 265/1 (Project ID 402170461, B09) and Germany's Excellence Strategy--EXC 2050/1 (Project ID 390696704)--Cluster of Excellence "Centre for Tactile Internet with Human-in-the-Loop" (CeTI) of Technische Universitat Dresden.
## Appendix A

For the simple multi-layer perceptron, we configure the architecture with five hidden layers, each comprising 400 neurons. The chosen activation function is the Swish activation function (Ramachandran et al., 2017).

For the LeNet-5 architecture, we adhere to the original design, which includes three convolutional layers, average pooling following the initial two convolutional layers, and two linear layers. The activation function used is the hyperbolic tangent. The convolutional layers employ a kernel size of \(5\times 5\), while the average pooling uses a window of shape \(2\times 2\). For the FashionMNIST dataset, the feature counts of the convolutional layers are designated as (6, 16, 120), and the two linear layers have neuron counts of (84, 10). However, for the CIFAR10 and CIFAR100 datasets, we elevate the feature counts of the convolutional layers to (18, 48, 360), with linear layer neuron counts set to (256, 10) for CIFAR10 and (256, 100) for CIFAR100.

For the MlpMixer architecture we employ six layers and a patch resolution of \(4\times 4\). Across all datasets, we maintain constant values for hidden size (\(C\)), sequence length (\(S\)), MLP channel dimension (\(D_{C}\)), and MLP token dimension (\(D_{S}\)); specifically \(C=256\), \(S=64\), \(D_{C}=512\) and \(D_{S}=512\) for all datasets.

For the VisionTransformer architecture, we adopt a slightly modified version of the ViT-Tiny setup: we use six layers, eight heads for each attention block, an embedding dimension of 256, and a hidden dimension of 512. The patch resolution of \(4\times 4\) is consistent with the MlpMixer. In both MlpMixer and VisionTransformer architectures, the GeLU activation function is used (Hendrycks and Gimpel, 2016).

For training using the maximum a posteriori estimate (Flat-MAP), dropout regularization, with dropout probability set to 0.2, is applied to all linear layers across all architectures, with the exception of the MlpMixer.

## Appendix B

In the centered parameterization of a generative model, Stochastic Variational Inference (SVI) with a fully factorized posterior yields a non-sparse solution, undermining the objective of employing shrinkage priors (Ghosh et al., 2019). Typically, this limitation is addressed by adopting the non-centered parameterization of the prior.

Consider the unique property of the half-Cauchy distribution: given \(x\sim C^{+}(0,1)\), and \(z=bx\) the resulting probability distribution for \(z\) is \(z\sim C^{+}(0,b)\). Therefore, the non-centred parameterization is formulated as

\[\hat{\tau}_{i}^{l} \sim\mathcal{C}^{+}(0,1)\] \[\hat{\lambda}_{ij}^{l} \sim\mathcal{C}^{+}(0,1)\] \[\hat{w}_{ij}^{l} \sim\mathcal{N}\left(0,1\right)\] \[\left[\gamma_{ij}^{l}\right]^{2} =\frac{\left[c^{l}\tau_{0}^{l}\tau_{i}^{l}\hat{\lambda}_{ij}^{l} \right]^{2}}{\left[c^{l}\right]^{2}+\left[\tau_{0}^{l}\tau_{i}^{l}\hat{ \lambda}_{ij}^{l}\right]^{2}}\] \[w_{ij}^{l} =\gamma_{ij}^{l}\hat{w}_{ij}^{l}\]
However, while the half-Cauchy distribution is frequently chosen for sampling-based inference, it poses challenges in variational inference (Piironen and Vehtari, 2017). Firstly, exponential family-based approximate posteriors (e.g., Gamma or log-Normal distributions) inadequately capture the half-Cauchy distribution's fat tails. Secondly, using a Cauchy approximating family for the posterior results in high variance gradients during stochastic variational inference (Ghosh et al., 2019). Hence, in the context of stochastic variational inference, the half-Cauchy distribution undergoes a reparameterization, as described in (Ghosh et al., 2018):

\[x\sim\mathcal{C}^{+}(0,b)\equiv x=\sqrt{\frac{1}{u}},u\sim\Gamma\left(\frac{1} {2},\frac{1}{v}\right),v\sim\Gamma\left(\frac{1}{2},b^{2}\right)\]

or, when represented in the non-centered parameterization:

\[x=b\sqrt{\frac{v}{u}},u\sim\Gamma\left(\frac{1}{2},1\right),v\sim\Gamma\left( \frac{1}{2},1\right) \tag{10}\]

## Appendix C
