Towards Efficient MCMC Sampling in Bayesian Neural Networks by
  Exploiting Symmetry
Bayesian inference in deep neural networks is challenging due to the
high-dimensional, strongly multi-modal parameter posterior density landscape.
Markov chain Monte Carlo approaches asymptotically recover the true posterior
but are considered prohibitively expensive for large modern architectures.
Local methods, which have emerged as a popular alternative, focus on specific
parameter regions that can be approximated by functions with tractable
integrals. While these often yield satisfactory empirical results, they fail,
by definition, to account for the multi-modality of the parameter posterior. In
this work, we argue that the dilemma between exact-but-unaffordable and
cheap-but-inexact approaches can be mitigated by exploiting symmetries in the
posterior landscape. Such symmetries, induced by neuron interchangeability and
certain activation functions, manifest in different parameter values leading to
the same functional output value. We show theoretically that the posterior
predictive density in Bayesian neural networks can be restricted to a
symmetry-free parameter reference set. By further deriving an upper bound on
the number of Monte Carlo chains required to capture the functional diversity,
we propose a straightforward approach for feasible Bayesian inference. Our
experiments suggest that efficient sampling is indeed possible, opening up a
promising path to accurate uncertainty quantification in deep learning.
Jonas Gregor Wiese, Lisa Wimmer, Theodore Papamarkou, Bernd Bischl, Stephan Günnemann, David Rügamer
# Towards Efficient MCMC Sampling in Bayesian Neural Networks by Exploiting Symmetry

###### Abstract

Bayesian inference in deep neural networks is challenging due to the high-dimensional, strongly multi-modal parameter posterior density landscape. Markov chain Monte Carlo approaches asymptotically recover the true posterior but are considered prohibitively expensive for large modern architectures. Local methods, which have emerged as a popular alternative, focus on specific parameter regions that can be approximated by functions with tractable integrals. While these often yield satisfactory empirical results, they fail, by definition, to account for the multi-modality of the parameter posterior. In this work, we argue that the dilemma between exact-but-unaffordable and cheap-but-inexact approaches can be mitigated by exploiting symmetries in the posterior landscape. Such symmetries, induced by neuron interchangeability and certain activation functions, manifest in different parameter values leading to the same functional output value. We show theoretically that the posterior predictive density in Bayesian neural networks can be restricted to a symmetry-free parameter reference set. By further deriving an upper bound on the number of Monte Carlo chains required to capture the functional diversity, we propose a straightforward approach for feasible Bayesian inference. Our experiments suggest that efficient sampling is indeed possible, opening up a promising path to accurate uncertainty quantification in deep learning.

Keywords:Uncertainty quantification Predictive uncertainty Bayesian inference Monte Carlo sampling Posterior symmetry

## 1 Introduction

Despite big data being the dominant paradigm in deep learning, the lack of infinitely many observations makes uncertainty quantification (UQ) an important problem in the field. A key component of UQ in Bayesian learning is the parameter posterior density that assigns a posterior probability to each parameter value5[23]. Between the extreme cases of all posterior probability mass
concentrating on a single value, indicating complete certainty about the model parameters, and being distributed uniformly over all possible values in a reflection of total ignorance, the shape of the parameter posterior density is central to the quantification of predictive uncertainty. However, the parameter posterior for Bayesian neural networks (BNNs) is typically highly multi-modal and rarely available in closed form. The classical Markov chain Monte Carlo (MCMC) approach asymptotically recovers the true posterior but is considered prohibitively expensive for BNNs, as the large number of posterior modes prevents a reasonable mixing of chains [25]. Popular approximation techniques, such as Laplace approximation (LA; [9, 29]) or deep ensembles (DE; [28]), therefore focus on local regions of the posterior landscape. While these methods are faster than traditional MCMC and perform well in many applications, they systematically omit regions of the parameter space that might be decisive for meaningful UQ [25] (also shown in Section 5.3).

In this work, we challenge the presumed infeasibility of MCMC for NNs and propose to exploit the - in this context, rarely considered - unidentifiability property of NNs, i.e., the existence of two or more equivalent parameter values that describe the same input-output mapping. We refer to these equivalent values as _equivoutput parameter states_. Equiooutput parameter states emerge from certain activation functions [6, 27, 37], as well as the free permutability of neuron parameters in hidden layers [20], and can be transformed into one another.

The functional redundancy arising from this phenomenon grows rapidly with the depth and width of a network (cf. Figure 1) and induces symmetries in the posterior density. For exact inference (up to a Monte Carlo error), we need to incorporate all non-equivo output parameter states that lead to distinct input-output mappings. Considering only these _functionally diverse_ mappings means, in turn, that our effective parameter space makes up a comparatively small fraction of the network's original parameter space. Since their numerous equiooutput counterparts do not contribute any new information to predictive uncertainty, we need much fewer MCMC samples when approximating the posterior predictive density (PPD) via Monte Carlo integration. By explicitly removing symmetries from samples _post-hoc_, we can even expose the functionally relevant part of the posterior and provide an opportunity for interpretation and analytical approximation in the reduced effective parameter space.

#### 1.0.1 Our Contributions.

We analyze the role of posterior space redundancies in quantifying BNN uncertainty, making the following contributions: 1) We show that the full PPD can be obtained from a substantially smaller reference set containing uniquely identified parameter states in function space. 2) We propose an estimation procedure for the number of Monte Carlo chains required to discover functionally diverse modes, providing a practical guideline for sampling from the parameter space of multi-layer perceptrons (MLPs). 3) We supply experimental evidence that our approach yields superior predictive performance compared to standard MCMC and local approximation methods. 4) Lastly, we demonstrate the posterior interpretability and analytic approximation that can
be obtained from explicitly removing symmetries _post-hoc_, for which we propose an algorithmic proof-of-concept.

## 2 Related Work

In this section, we review the literature on the existence, removal, and utilization of symmetries in MLPs, as well as the relation between such symmetries and statistical inference.

#### 2.0.1 Existence of Parameter State Symmetries.

Non-unique network parameter states have been considered in the literature before. [20] were among the first to note that equioutput states induce symmetries in the parameter space of MLPs. Focusing, within the general linear group of the parameter space, on the subgroup of transformations that leave the input-output mapping unchanged, they derived equivalence classes of equioutput parameter states and showed that, for every MLP, there exists a minimal and complete set of representatives covering all functionally different parameter states. [42] and [27] continued along this line of work to study single-hidden-layer MLPs with specific activation functions, advancing from \(\tanh\) to more general self-affine activations. An extension to MLPs of arbitrary depth was studied by [6] in the context of \(\tanh\) activations. More recently, [37] characterized equioutput parameter states for ReLU activations, again focusing on the case of a single hidden layer, and [1] classified all \(\mathbb{G}\)-invariant single-hidden-layer MLPs with ReLU activation for any finite orthogonal group \(\mathbb{G}\). Lastly, [45] generalized much of the above in a framework addressing the identifiability of affine symmetries in arbitrary architectures.

#### 2.0.2 Symmetry Removal.

Symmetries in the parameter posterior density of Bayesian models can produce adverse effects that have been addressed in several research areas of statistics and machine learning. A prominent example is _label switching_ in finite mixture models, where the permutability of label assignments to the mixture components induces symmetries similar to those in BNNs. To make mixture models identifiable, [3] introduced an adaptive Metropolis algorithm with online relabeling, effectively removing permutation symmetries by optimizing over discrete sets of permutation matrices. Such exhaustive-search approaches, however, scale poorly to modern NNs with many parameters, as the amount of equioutput states rises exponentially with the number of parameters.

In BNNs, symmetries have been known to slow down MCMC convergence to the stationary parameter posterior density due to budget wasted on visiting symmetric modes [31, 34]. [25], reporting results from extensive and large-scale experiments, indeed find that MCMC chains tend to mix better in function space than in parameter space. Consequently, reducing the effect of symmetries by imposing parameter constraints and defining anchoring points for subsets of the latent variables has been shown to improve mixing [41]. A proposal for constrained sampling can be found, for example, in [39], with application to ReLU-activated MLPs.
Utilizing Symmetries.Symmetries in the parameter posterior density are not, however, necessarily a nuisance. Quite on the contrary, they can be useful to enhance generalization and make inference affordable. An increasing body of work has been exploring the use of symmetry removal in the context of _mode connectivity_, an approach to find more robust NN solutions by retrieving connected areas of near-constant loss rather than isolated local optima [10, 17]. Focusing on equioutput permutations of hidden-layer neurons, [43] and [2], among others, propose to align the layer-wise embeddings of multiple networks and thus improve upon the performance of individual models. Following a similar idea, [38] apply a _post-hoc_ standardization on parameter vectors in ReLU-activated NNs that draws from the notion of equioutput equivalence classes.

In the field of Bayesian deep learning, the idea of utilizing - exact or approximate - parameter symmetries, represented by permutation groups, has led to the development of _lifted MCMC_[33, 5]. Orbital Markov chains have been introduced to leverage parameter symmetries in order to reduce mixing times [32, 5]. Lifted MCMC has been considered mainly in the context of probabilistic graphical models. There is scope to harness lifted MCMC in the context of MLPs since these can be cast as graphical models [26, 45].

We believe that equioutput symmetries have the potential to facilitate MCMC inference, despite the apparent complexity they introduce in the parameter posterior density. From the insight that a vast part of the sampling space is made up of symmetric copies of some minimal search set, we conclude that running multiple short MCMC chains in parallel, each of which can sample a functionally different mode, represents a more efficient use of the available budget than collecting a large number of samples from a single chain. In Section 4, we propose an upper bound for MCMC chains necessary to observe all functionally diverse posterior modes, which is a key criterion for successful inference. The perspective of parameter posterior symmetries thus lends a new theoretical justification to previous efforts in multi-chain MCMC that are motivated mainly by the exploitation of parallel computing resources [30, 40]. Our experiments in Section 5 suggest that this approach is indeed more effective than single-chain MCMC in BNNs with many symmetries in their parameter posterior density. In agreement with [25], our findings advocate to focus on function-space instead of parameter-space mixing during MCMC sampling.

We thus view the existence of equioutput parameter states as a benign phenomenon. That said, there are still benefits to be gained from removing the symmetries: with a parameter posterior density reduced to the minimal parameter set sufficient to represent its full functional diversity, we get an opportunity for better interpretation, and possibly even analytical approximation. We demonstrate the potential of symmetry removal in Section 5.3 by means of a custom algorithm (Supplementary Material B). In the following section, we provide the mathematical background and introduce the characterization and formal notation of equioutput transformations.
## 3 Background and Notation

#### 3.0.1 MLP Architectures.

In this work, we consider NNs of the following form. Let \(f:\mathcal{X}\rightarrow\mathcal{Y}\) represent an MLP with \(K\) layers, where layer \(l\in\{1,\ldots,K\}\) consists of \(M_{l}\) neurons, mapping a feature vector \(\mathbf{x}=(x_{1},\ldots,x_{n})^{\top}\in\mathcal{X}\subseteq\mathbb{R}^{n},n\in \mathbb{N}\), to an outcome vector

\[f(\mathbf{x})=:\hat{\mathbf{y}}=(\hat{y}_{1},\ldots,\hat{y}_{m})^{\top}\in\mathcal{Y} \subseteq\mathbb{R}^{m},\ \ m\in\mathbb{N},\]

to estimate \(\mathbf{y}=(y_{1},\ldots,y_{m})^{\top}\in\mathcal{Y}\). The \(i\)-th neuron in the \(l\)-th layer of the MLP is associated with the weights \(w_{lij},\ j=1,\ldots,M_{l-1}\), and the bias \(b_{li}\). We summarize all the MLP parameters in the vector

\[\mathbf{\theta}:=(w_{211},\ldots,w_{KM_{K}M_{K-1}},b_{21},\ldots,b_{KM_{K}})^{\top }\in\Theta\subseteq\mathbb{R}^{d}\]

and write \(f_{\mathbf{\theta}}\) to make clear that the MLP is parameterized by \(\mathbf{\theta}\). For each hidden layer \(l\in\{2,\ldots,K-1\}\), the inputs are linearly transformed and then activated by a function \(a\). More specifically, we define the pre-activations of the \(i\)-th neuron in the \(l\)-th hidden layer as

\[o_{li}=\sum_{j=1}^{M_{l-1}}w_{lij}z_{(l-1)j}+b_{li}\]

with post-activations \(z_{(l-1)i}=a(o_{(l-1)i})\) from the preceding layer. For the input layer, we have \(z_{1i}=x_{i},i=1,\ldots,n\), and for the output layer, \(z_{Ki}=\hat{y}_{i},i=1,\ldots,M_{K}\).

#### 3.0.2 Predictive Uncertainty.

In the Bayesian paradigm, a prior density \(p(\mathbf{\theta})\) is imposed on the parameters, typically as part of a Bayesian model of the data. Using Bayes' rule, the parameter posterior density

\[p(\mathbf{\theta}|\mathcal{D})=\frac{p(\mathcal{D}|\mathbf{\theta})p(\mathbf{\theta})}{p( \mathcal{D})}\]

updates this prior belief based on the information provided by the data \(\mathcal{D}\) and encoded in the likelihood \(p(\mathcal{D}|\mathbf{\theta})\). In supervised learning, the data are typically given by a set of \(N\) feature vectors \(\mathbf{x}\in\mathcal{X}\) and outcome vectors \(\mathbf{y}\in\mathcal{Y}\), forming the dataset \(\mathcal{D}=\{(\mathbf{x}^{(1)},\mathbf{y}^{(1)}),\ldots,(\mathbf{x}^{(N)},\mathbf{y}^{(N)})\}\). The PPD \(p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathcal{D})\) quantifies the predictive or functional uncertainty of the model for a new observation \((\mathbf{x}^{*},\mathbf{y}^{*})\in\mathcal{X}\times\mathcal{Y}\). Since

\[p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathcal{D})=\int_{\Theta}p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{ \theta})p(\mathbf{\theta}|\mathcal{D})\,\mathrm{d}\mathbf{\theta},\]

deriving this uncertainty requires access to the posterior density \(p(\mathbf{\theta}|\mathcal{D})\), which can be estimated from MCMC sampling.
Equioutput Transformations.Let us now characterize the notion of equioutput parameter states, and the transformations to convert between them, more formally. Two parameter states \(\mathbf{\theta},\mathbf{\theta}^{\prime}\) are considered _equioutput_ if the maps \(f_{\mathbf{\theta}},f_{\mathbf{\theta}^{\prime}}\) yield the same outputs for all possible inputs from \(\mathcal{X}\). We denote this equivalence relation (see proof in Supplementary Material A.2) by \(\sim\) and write:

\[\mathbf{\theta}\sim\mathbf{\theta}^{\prime}\Longleftrightarrow f_{\mathbf{\theta}}(\mathbf{x})= f_{\mathbf{\theta}^{\prime}}(\mathbf{x})\,\forall\,\mathbf{x}\in\mathcal{X},\ \ \mathbf{\theta},\mathbf{\theta}^{\prime}\in\varTheta.\]

The equioutput relation is always defined with respect to a particular MLP \(f\), which we omit in our notation when it is clear from the context.

All MLPs with more than one neuron in at least one hidden layer exhibit such equioutput parameter states that arise from permutation invariances of the input-output mapping [20, 27]. Since the operations in the pre-activation of the \(i\)-th neuron in the \(l\)-th layer commute6, the \(M_{l}>1\) neurons of a hidden layer \(l\) can be freely interchanged by permuting their associated parameters. In addition, equioutput transformations can arise from the use of certain activation functions with inherent symmetry properties. For example, in the case of \(\tanh\), the signs of corresponding parameters can be flipped using \(\tanh(x)=-\tanh(-x)\). For ReLU activations, a scaling transformation can be applied such that the mapping of the network remains unchanged, i.e., \(\text{ReLU}(x)=c^{-1}\cdot\text{ReLU}(c\cdot x)\) for \(|c|>0\) (see also Supplementary Material A.1).

Footnote 6: Recall that the pre-activation of neuron \(i\) in layer \(l\) is \(o_{li}=\sum_{j=1}^{M_{l-1}}w_{lij}z_{(l-1)j}+b_{li}\). By the commutative property of sums, any permutation \(\pi:J\to J\) of elements from the set \(J=\{1,\ldots,M_{l-1}\}\) will lead to the same pre-activation:

\[\mathcal{F}_{\mathbf{T}}:\varTheta\to\varTheta,\mathbf{\theta}\mapsto\mathbf{T}\mathbf{\theta },\ \ \mathbf{T}\in\mathbb{R}^{d\times d},\]

be an activation-related transformation of a parameter vector that might, for instance, encode an output-preserving sign flip. \(\mathcal{F}_{\mathbf{T}}\) constitutes an _equioutput_ transformation if \(f_{\mathbf{\theta}}(\cdot)=f_{\mathcal{F}_{\mathbf{T}}(\mathbf{\theta})}(\cdot)\). We collect all output-preserving transformation matrices \(\mathbf{T}\) in the set \(\mathcal{T}\), i.e.,

\[\mathcal{T}=\left\{\mathbf{T}\in\mathbb{R}^{d\times d}\mid f_{\mathbf{\theta}}(\cdot) =f_{\mathcal{F}_{\mathbf{T}}(\mathbf{\theta})}(\cdot)\right\}.\]

Similarly, let

\[\mathcal{F}_{\mathbf{P}}:\varTheta\to\varTheta,\mathbf{\theta}\mapsto\mathbf{P}\mathbf{ \theta},\ \ \mathbf{P}\in\{0,1\}^{d\times d},\]

be a transformation that permutes elements in the parameter vector. We define the set of permutation matrices that yield equioutput parameter states as

\[\mathcal{P}=\left\{\mathbf{P}\in\mathbb{R}^{d\times d}\mid f_{\mathbf{\theta}}(\cdot) =f_{\mathcal{F}_{\mathbf{P}}(\mathbf{\theta})}(\cdot)\right\}.\]
The cardinality of \(\mathcal{P}\) is at least \(\prod_{l=2}^{K-1}M_{l}!\)[20] when traversing through the NN from the first layer in a sequential manner, applying to each layer permutations that compensate for permutations in its predecessor.

Since activation functions operate neuron-wise, activation- and permutation-related equioutput transformations do not interact (for instance, we could permute the associated weights of two neurons and later flip their sign). We can, therefore, define arbitrary combinations of activation and permutation transformations as

\[\mathcal{E}=\left\{\boldsymbol{E}=\boldsymbol{TP}\in\mathbb{R}^{d\times d}, \boldsymbol{T}\in\mathcal{T},\boldsymbol{P}\in\mathcal{P}\ |\ f_{\boldsymbol{\theta}}(\cdot)=f_{\mathcal{F}_{\boldsymbol{E}}( \boldsymbol{\theta})}(\cdot)\right\}.\]

The transformation matrices in \(\mathcal{E}\) will exhibit a block-diagonal structure with blocks corresponding to network layers. This is due to the permutations \(\boldsymbol{P}\) affecting both incoming and outgoing weights, but only in the sense that two incoming and two outgoing weights swap places, never changing layers. The activation-related sign flips or re-scalings occur neuron-wise, making \(\boldsymbol{T}\) a diagonal matrix that does not alter the block-diagonal structure of \(\boldsymbol{P}\).

For the cardinality of the set \(\mathcal{E}\) of equioutput transformations, we can establish a lower bound that builds upon the minimum cardinality of \(\mathcal{P}\):

\[|\mathcal{E}|\geq\prod_{l=2}^{K-1}M_{l}!\cdot|\mathcal{T}_{l}|,\]

where \(|\mathcal{T}_{l}|\) denotes the number of activation-related transformations applicable to neurons in layer \(l\). From this, it becomes immediately clear that the amount of functional redundancy increases rapidly with the network size (see also Figure 1). As a result of equioutput parameter states, the MLP parameter posterior density exhibits functional redundancy in the form of symmetries for commonly used priors (see Section 4):

\[p(\boldsymbol{\theta}|\mathcal{D})=\frac{p(\mathcal{D}|\boldsymbol{\theta})p( \boldsymbol{\theta})}{p(\mathcal{D})}\stackrel{{ 10,13}}{{=}}\frac{p(\mathcal{D}|\boldsymbol{E \theta})p(\boldsymbol{E\theta})}{p(\mathcal{D})}=p(\boldsymbol{E\theta}| \mathcal{D}),\ \ \boldsymbol{\theta}\in\Theta,\boldsymbol{E}\in\mathcal{E}. \tag{1}\]

## 4 Efficient Sampling

These symmetric structures in the parameter posterior density suggest that sampling can be made more efficient. In the following section, we show that the PPD can theoretically be obtained from a small reference set of non-equivalent parameter states and propose an upper bound on Markov chains that suffice to sample all non-symmetric posterior modes.

### Posterior Reference Set

As introduced in Section 3, for each parameter state \(\boldsymbol{\theta}\) of an NN, there are functionally redundant counterparts \(\boldsymbol{\theta}^{\prime}\) related to \(\boldsymbol{\theta}\) by an equioutput transformation, such that \(f_{\boldsymbol{\theta}}(\cdot)=f_{\boldsymbol{\theta}^{\prime}}(\cdot)\). We can use this equivalence relation to dissect the parameter space \(\Theta\) into disjoint equivalence classes. For this, let the _reference set_
\(\mathcal{S}_{1}\) be a minimal set of representatives of each equivalence class (cf. _open minimal sufficient search sets_ in [6]). All parameter states in \(\mathcal{S}_{1}\) are functionally diverse, i.e., \(\boldsymbol{\theta},\tilde{\boldsymbol{\theta}}\in\mathcal{S}_{1}\Rightarrow \boldsymbol{\theta}\not\sim\tilde{\boldsymbol{\theta}}\), and each element in \(\Theta\) is equivalent to exactly one element in \(\mathcal{S}_{1}\). For a finite amount of equioutput transformations, as in the case of \(\tanh\)-activated MLPs (finite possibilities of sign-flip combinations of hidden neurons), the NN parameter space can then be dissected into \(|\mathcal{E}|\) disjoint _representative sets_, which contain equioutput transformations of the elements of the reference set, in the following way.

Proposition 1 (Parameter space dissection): _Let \(\mathcal{S}_{1}\) be the reference set of uniquely identified network parameter states. Then, for a finite number of equioutput transformations, it holds that the parameter space can be dissected into \(|\mathcal{E}|\) disjoint, non-empty representative sets up to a set \(\mathcal{S}^{0}\subset\Theta\), i.e.,_

\[\Theta=\left(\ \dot{\bigcup}_{j=1}^{|\mathcal{E}|}\mathcal{S}_{j}\right) \dot{\cup}\mathcal{S}^{0}\text{, where }\mathcal{S}_{j}\cong\{\boldsymbol{\theta}\ |\ \boldsymbol{\theta}=\boldsymbol{E}_{j}\boldsymbol{\theta}^{\prime}\quad \forall\boldsymbol{\theta}^{\prime}\in\mathcal{S}_{1},\boldsymbol{E}_{j}\in \mathcal{E}\}\,, \tag{2}\]

_where \(\dot{\cup}\) denotes the union over disjoint sets. We use \(\mathcal{S}^{0}\) as a residual quantity to account for cases that cannot be assigned unambiguously to one of the sets \(\mathcal{S}_{j}\) because they remain unchanged even under a transformation with non-identity matrices \(\boldsymbol{E}_{j}\in\mathcal{E}\)._

The edge cases that make up \(\mathcal{S}^{0}\) exist, e.g., on the boundary of two classes [6] or in degenerated cases such as the zero vector [42, 45]. For a characterization of the involved sets, as well as a proof sketch, see Supplementary MaterialA.4.

Equioutput parameter states have the same posterior probabilities \(p(\boldsymbol{\theta}|\mathcal{D})=p(\boldsymbol{E}\boldsymbol{\theta}| \mathcal{D})\) if the prior is transformation-invariant; see Supplementary Material A.3, Equations (10)-(13). Moreover, equioutput parameter states produce by definition the same predictions \(p(\boldsymbol{y}^{*}|\boldsymbol{x}^{*},\boldsymbol{\theta})=p(\boldsymbol{y} ^{*}|\boldsymbol{x}^{*},\boldsymbol{E}\boldsymbol{\theta})\) for any \(\boldsymbol{E}\in\mathcal{E}\). Thus, the following corollary holds.

Figure 1: Example of \(\tanh\)-activated MLPs. _Left_: Cardinality lower bound of the equioutput transformation set for a single hidden layer with 1 to 128 neurons; the functional redundancy factor for a network with 128 neurons is \(1.31\cdot 10^{254}\). _Right_: A ten-dimensional MLP parameter posterior (top-right corner, depicted as bivariate marginal density) exhibits symmetries, such that all red sample clusters are equioutput-related to the green cluster. The associated function spaces are identical, i.e., many posterior modes are redundant.
Corollary 1 (Reformulated posterior predictive density): _Let \(\mathcal{E}\) be finite. As in Proposition 1, consider the disjoint non-empty sets \(\mathcal{S}_{j},j\in\{1,\ldots,|\mathcal{E}|\}\), and residual space \(\mathcal{S}^{0}\). If the prior density \(p(\boldsymbol{\theta})\) is transformation-invariant, then the posterior predictive density expresses as_

\[p(\boldsymbol{y}^{*}|\boldsymbol{x}^{*},\mathcal{D}) =\int_{\Theta}p(\boldsymbol{y}^{*}|\boldsymbol{x}^{*},\boldsymbol {\theta})p(\boldsymbol{\theta}|\mathcal{D})\,\mathrm{d}\boldsymbol{\theta} \tag{4}\] \[=|\mathcal{E}|\int_{\mathcal{S}_{j}}p(\boldsymbol{y}^{*}| \boldsymbol{x}^{*},\boldsymbol{\theta})p(\boldsymbol{\theta}|\mathcal{D})\, \mathrm{d}\boldsymbol{\theta}+\int_{\mathcal{S}^{0}}p(\boldsymbol{y}^{*}| \boldsymbol{x}^{*},\boldsymbol{\theta})p(\boldsymbol{\theta}|\mathcal{D})\, \mathrm{d}\boldsymbol{\theta}\] \[\approx|\mathcal{E}|\int_{\mathcal{S}_{j}}p(\boldsymbol{y}^{*}| \boldsymbol{x}^{*},\boldsymbol{\theta})p(\boldsymbol{\theta}|\mathcal{D})\, \mathrm{d}\boldsymbol{\theta}.\]

The proof of Corollary 1 is given in Supplementary Material A.3. It follows from Proposition 1 and the assumption of transformation-invariant prior densities, which is often satisfied in practice (e.g., for widely-applied isotropic Gaussian priors). We can further approximate (4) by (4) as the set \(\mathcal{S}^{0}\subset\mathbb{R}^{d}\) is of negligible size (depending on \(\Theta\), potentially even with zero Lebesgue measure).

As a consequence of Corollary 1, the PPD can be obtained up to the residual set by only integrating over uniquely identified parameter states from one of the sets \(\mathcal{S}_{j}\), with a multiplicative factor \(|\mathcal{E}|\) that corrects the probability values by the amount of redundancy in the posterior. In other words, only a fraction \(1/|\mathcal{E}|\) of the posterior must be sampled in order to infer a set of uniquely identified parameter states of the NN, and thus, to obtain the full PPD. This reduces the target sampling space drastically, as illustrated in Figure 1. For example, it allows the posterior space of a single-layer, tanh-activated network with 128 neurons to be effectively reduced to a \(10^{254}\)-th of its original size.

In the case of an infinite amount of equioutput transformations, such as in ReLU-activated MLPs (the scaling factor \(|c|>0\) can be chosen arbitrarily), we can use similar reasoning. Only one representative set of the posterior density needs to be observed in order to capture the full functional diversity of a network because the integrals over two representative sets are identical. For a more in-depth discussion of ReLU symmetries, see, for example, [4].

#### 4.2.2 How to Obtain a Representative Set?

In practice, when using Monte Carlo to approximate Equation (4), it is not necessary to actually constrain the sampling procedure to a specific set \(\mathcal{S}_{j}\), which might indeed not be straightforward7. Since any equioutput transformation is known _a priori_, we just need to be aware of the fact that each sample can theoretically be mapped to different representative sets after running the sampling procedure. Hence, for the calculation of the PPD integral, the samples can remain scattered across the various representative sets as long as they cover all functionally diverse parameter states. For
the purpose of providing better interpretability and analytic approximation of the posterior, it may still be worthwhile to explicitly remove the symmetries. In Section 5.3, we demonstrate such symmetry removal using a custom algorithm for \(\tanh\)-activated networks (Supplementary Material B).

### An Upper Bound for Markov Chains

The question remains how many samples are needed to approximate a set of uniquely identified parameter states sufficiently well. Even in a symmetry-free setting, BNN posteriors can exhibit multiple functionally diverse modes representing structurally different hypotheses, depending on the network architecture and the underlying data-generating process. For example, in Section 5.3, we discuss the case of an under-parameterized network that preserves three distinctive modes caused by its restricted capacity.

In the following, we assume \(\nu\in\mathbb{N}\) functionally diverse modes with the goal of visiting every mode or its local proximity at least once when running MCMC. As the ability to switch from one mode to another within a chain depends on various factors, such as the acceptance probability and the current state of other parameters, increasing the number of samples per chain does not necessarily correlate with the number of visited modes. We, therefore, propose to focus on the number of independent chains, rather than the number of samples per chain, to effectively control the number of visited modes.

This further allows us to derive an upper bound for the number of independent chains that are required to visit every mode at least once. The number of samples from each chain will then ultimately determine the approximation quality. In the computation of the PPD, we formulate the Monte Carlo integration over all samples from all chains simultaneously [30]. In practice, given a user-defined number of maximal resources \(\rho\) (e.g., CPU cores), the following proposition provides a lower bound on the probability that the number of chains \(\mathcal{G}\) necessary to visit every mode remains below the resource limit of the user (i.e., \(\mathcal{G}<\rho\)).

Proposition 2 (Probabilistic bound for sufficient number of Markov chains): _Let \(\pi_{1},\ldots,\pi_{\nu}\) be the respective probabilities of the \(\nu\) functionally diverse modes to be visited by an independently started Markov chain and \(\Pi_{J}:=\sum_{j\in J}\pi_{j}\). Then, given \(\rho\) chains,_

\[\mathbb{P}(\mathcal{G}<\rho)\geq 1-\rho^{-1}\left\{\sum_{q=0}^{\nu-1}(-1)^{ \nu-1-q}\sum_{J:|J|=q}(1-\Pi_{J})^{-1}\right\}. \tag{5}\]

The proof can be found in Supplementary Material A.6. Note that this bound is independent of the NN architecture and only depends on the assumptions about the number and probabilities of functionally diverse modes \(\nu\), disregarding symmetric copies. Proposition 2 can be used to calculate the number of MCMC chains given certain assumptions - for example, from domain knowledge, or in a worst-case scenario calculation - and thus provides practical guidance for MCMC sampling of MLPs. Judging by the comparably high predictive performance of
local approximations such as LA and DE [28, 29], we conclude that a small amount of functional modes is reasonable to assume in practice. Our qualitative experiments in Section 5 support this supposition.

As an example of applying Proposition 2, assume \(\nu=3\) functionally diverse modes in a reference set with \(\pi_{1}=0.57,\pi_{2}=0.35,\pi_{3}=0.08\) (chosen to represent a rather diverse functional mode set). An upper bound of \(\rho=1274\) chains ensures that we observe all functionally diverse modes with probability \(\mathbb{P}(\mathcal{G}<\rho)\geq 0.99\).

## 5 Experiments

We now investigate our theoretical findings and compare the resulting approach to single-chain MCMC and DE. In all experiments8, we employ a Bayesian regression model with a normal likelihood function, standard normal prior for parameters \(\mathbf{\theta}\), and a truncated standard normal prior restricted to the positive real line for the variance of the normal likelihood, which we treat as a nuisance parameter. Depending on the task, we either use a No-U-Turn sampler [22] with \(2^{10}\) warmup steps to collect a single sample from the posterior or derive the maximum-a-posteriori estimator using a gradient-based method (details are given in Supplementary Material C.2, C.3).

Footnote 8: [https://anonymous.4open.science/r/efficient_mcmc_sampling_in_bnns_symmetry-6FAF/](https://anonymous.4open.science/r/efficient_mcmc_sampling_in_bnns_symmetry-6FAF/)

### Performance Comparison

In our first experiment, we demonstrate the predictive performance of BNNs, where the PPD is calculated based on MCMC sampling, using the derived upper bound for the number of chains (ours). In this case, we collect one sample per chain for \(G\) chains, and thus \(G\) samples in total. This is compared to MCMC sampling collecting \(G\) samples from a single chain (s.c.), and DE with ten ensemble members on three synthetic datasets (\(\mathcal{D}_{S}\), \(\mathcal{D}_{I}\), and \(\mathcal{D}_{R}\)) as well as benchmark data from [11] (for dataset details and additional results on LA,

\begin{table}
\begin{tabular}{l|c c c c|c c c}  & \multicolumn{3}{c}{**Smaller network \(f_{1}\)**} & \multicolumn{3}{c}{**Larger network \(f_{2}\)**} \\  & MCMC (ours) & MCMC (s.c.) & DE & MCMC (ours) & MCMC (s.c.) & DE \\ \hline \(\mathcal{D}_{S}\) & **-0.53** (\(\pm\) 0.09) & -0.56 (\(\pm\) 0.11) & -0.58 (\(\pm\) 0.11) & **-0.59** (\(\pm\) 0.12) & **-0.59** (\(\pm\) 0.12) & -2.13 (\(\pm\) 0.03) \\ \(\mathcal{D}_{I}\) & **0.79** (\(\pm\) 0.06) & 0.65 (\(\pm\) 0.07) & 0.56 (\(\pm\) 0.06) & **0.91** (\(\pm\) 0.09) & **0.91** (\(\pm\) 0.09) & -2.02 (\(\pm\) 0.02) \\ \(\mathcal{D}_{R}\) & 0.64 (\(\pm\) 0.10) & **0.75** (\(\pm\) 0.11) & -1.46 (\(\pm\) 0.06) & **0.95** (\(\pm\) 0.08) & **0.95** (\(\pm\) 0.08) & -2.20 (\(\pm\) 0.02) \\ Airfoil & **-0.74** (\(\pm\) 0.04) & -0.80 (\(\pm\) 0.05) & -1.62 (\(\pm\) 0.03) & **0.92** (\(\pm\) 0.05) & 0.72 (\(\pm\) 0.10) & -2.17 (\(\pm\) 0.01) \\ Concrete & **-0.41** (\(\pm\) 0.05) & -0.44 (\(\pm\) 0.06) & -1.59 (\(\pm\) 0.03) & **0.26** (\(\pm\) 0.07) & 0.25 (\(\pm\) 0.07) & -2.03 (\(\pm\) 0.01) \\ Diabetes & **-1.20** (\(\pm\) 0.07) & **-1.20** (\(\pm\) 0.07) & -1.47 (\(\pm\) 0.07) & **-1.18** (\(\pm\) 0.08) & -1.22 (\(\pm\) 0.09) & -2.09 (\(\pm\) 0.04) \\ Energy & **0.92** (\(\pm\) 0.04) & 0.69 (\(\pm\) 0.12) & -1.76 (\(\pm\) 0.02) & 2.07 (\(\pm\) 0.46) & **2.38** (\(\pm\) 0.11) & -1.99 (\(\pm\) 0.02) \\ ForestF & **-1.37** (\(\pm\) 0.07) & **-1.37** (\(\pm\) 0.07) & -1.60 (\(\pm\) 0.06) & **-1.43** (\(\pm\) 0.45) & -1.69 (\(\pm\) 0.49) & -2.20 (\(\pm\) 0.02) \\ Yacht & **1.90** (\(\pm\) 0.16) & 1.29 (\(\pm\) 0.56) & -1.14 (\(\pm\) 0.14) & **3.31** (\(\pm\) 0.21) & 0.15 (\(\pm\) 0.09) & -2.18 (\(\pm\) 0.03) \\ \end{tabular}
\end{table}
Table 1: Mean log pointwise predictive density (L PPD) values on test sets (larger is better; one standard error in parentheses). The highest performance per dataset and network is highlighted in bold.
see Supplementary Material C.1 and E.1, respectively). We use a smaller NN \(f_{1}\) with a single hidden layer containing three neurons and a larger network \(f_{2}\) with three hidden layers having 16 neurons each, both with \(\tanh\) activation. As in Section 4.2, we assume three functionally diverse modes \(\nu=3\) and mode probabilities \(\pi_{1}=0.57,\pi_{2}=0.35,\pi_{3}=0.08\) as in the given example. To demonstrate the performance of our MCMC-based PPD approximation, we measure the goodness-of-fit on the test data using the log point-wise predictive density (LPPD; [18])

\[\text{LPPD}=\log\int_{\Theta}p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{\theta})p(\mathbf{\theta} |\mathcal{D})\,\mathrm{d}\mathbf{\theta}\approx\log\Bigg{(}\tfrac{1}{G}\sum_{g=1} ^{G}p(\mathbf{y}^{*}|\mathbf{x}^{*},\mathbf{\theta}^{(g)})\Bigg{)}, \tag{6}\]

where \(\mathbf{\theta}^{(1)},\ldots,\mathbf{\theta}^{(G)}\) are \(G\) samples obtained across all chains via MCMC sampling from the parameter posterior density \(p(\mathbf{\theta}|\mathcal{D})\). Equation (6) is evaluated at each test point \((\mathbf{x}^{*},\mathbf{y}^{*})\). Table 1 reports the mean LPPD across \(N^{*}\) independent test points for each combination of dataset and sampling scheme (see Appendix C.1 for details). Our results clearly indicate that using only a moderate amount of Markov chains, following our approach, yields equal or even better performance than single-chain MCMC and DE in all but two experiments.

### Practical Evaluation of Corollary 1

Next, we investigate the property derived in Corollary 1 using our proposed upper bound of chains, again with the assumption from the example in Section 4.2. To this end, we analyze the PPD for dataset \(\mathcal{D}_{I}\), using network \(f_{2}\). For every newly collected sample in the MCMC run, the updated PPD is computed approximately on a two-dimensional (input/output) grid. Then, the Kullback-Leibler (KL) divergence between consecutive densities is averaged over the grid of input values of \(f_{2}\) (details in Supplementary Material D.1). As shown in Figure 2, despite the size of the network \(f_{2}\) and the high amount of equioutput parameter states \(|\mathcal{E}|=\left(16!\cdot 2^{16}\right)^{3}\approx 2.58\cdot 10^{54}\), the PPD converges after notably fewer than \(|\mathcal{E}|\) samples and plots of the function space indicate the saturation of functional diversity already after 1274 samples from as many chains.

### Posterior Symmetry Removal

So far we were mainly concerned with the predictive performance of MCMC sampling. Yet, mapping all samples to a joint representative set, as characterized in Section 4.1, has the potential to reduce the effective weight space enormously, facilitating interpretability and possibly even analytical approximation. For this, we propose a custom algorithm for \(\tanh\)-activated MLPs as a proof-of-concept. Our algorithm removes symmetries in a data-dependent manner and thus minimizes the number of remaining modes in the representative set (details in Supplementary Material B).
We demonstrate the efficacy of the approach for \(f_{1}\) in two experiments A and B for datasets \(\mathcal{D}_{S}\) and \(\mathcal{D}_{I}\), respectively. For experiment A (Figure 3, top) we visualize the neuron parameter space along the steps of the proposed algorithm (for details, see Algorithms 1-4 in Supplementary Material B). Different colors encode the current neuron index (i.e., one of \(\{1,2,3\}\)) in the hidden layer of the respective neuron parameter vector. Initially, symmetries of the posterior densities are clearly noticeable (Figure 3a), and the neuron parameter vectors

Figure 3: For experiment A (top row), the parameter subspaces of hidden neurons are visualized in their initial state as obtained when running MCMC (a), and after their transformation and reassignment during the steps of Algorithm 4 (b, c). Different colors encode the current neuron index in the hidden layer. By reassigning neurons, Algorithm 4 effectively finds an optimal reference set and allows to separate the multi-modal complex univariate marginal density (d) into three univariate densities. For experiment B (bottom row), the MLP parameter posterior density after symmetry removal results in a tri-modal system, here illustrated as a bivariate plot (f). Investigating these modes, we find that all are functionally diverse, i.e., represent a different hypothesis of the dataset (g-i). Combined, they form the full function space (j).

Figure 2: Convergence of MCMC depicted as the change in KL-divergence on original (black) and log-scale (blue) when consecutively adding another sample from a new and independent chain and re-estimating the posterior density. Small overlaying plots: approximated PPD of the network after \(2^{0}\), \(2^{4}\), \(2^{8}\), and \(G=1274\) samples; darker colors correspond to higher probabilities.
are distributed identically (Figure 2(d)). Upon the first step of the algorithm (Figure 2(b)), the parameter space is effectively halved as a consequence of removing the \(\tanh\)-induced sign-flip symmetry, and three clusters remain. The second (clustering) step removes all permutation symmetries from the full posterior density by assigning areas of this parameter subspace to the neurons of the model. This is depicted in Figure 2(c) and clearly shows the separation of states by different cluster colors. In Figures 2(d) and 2(e), the univariate marginal density of each neuron's incoming weight \(w_{2i1}\) reveals their reassignment to the parameter space. After running our algorithm, each neuron exhibits a distinct unimodal density, resulting in a unimodal density in the full parameter space (visualization in Supplementary Material E.3). We can conclude that only one functionally diverse mode exists in this case, which should also be recovered by local approximations like LA.

For experiment B (Figure 3, bottom), we focus on the visualization of the full parameter posterior density obtained after the application of the symmetry removal algorithm (Figure 2(f)). Three functionally diverse modes, represented by different colors, remain in the BNN posterior. We can now interpret these modes in function space by further clustering the transformed samples using a spectral clustering approach (details in Supplementary Material D.2). Figures 2(g)-j visualize the network parameter states in the function space, revealing three functionally diverse hypotheses the network can potentially learn from the given dataset. Such knowledge allows for a better-suited approximation by, e.g., a mixture of Laplace approximations (MoLA; [13]), as shown in Supplementary Material E.2. Note that approaches focusing on a single mode, such as standard LA, would have captured only a third of the functional diversity. For meaningful UQ, it is thus imperative that all functionally diverse modes are accounted for.

## 6 Discussion

We showed that the PPD for Bayesian MLPs can be obtained from just a fraction of the parameter space due to the existence of equioutput parameter states. Together with an upper bound on the number of MCMC chains to guarantee the recovery of every functionally diverse mode, our approach paves the way towards exact uncertainty quantification (up to a Monte Carlo error) in deep learning. Furthermore, we demonstrate the use of symmetry removal and present a proof-of-concept approach to map samples of \(\tanh\)-activated MLPs to a representative set. This _post-hoc_ procedure improves the interpretability of the symmetry-free posterior density drastically and facilitates analytical approximations. As a future research direction, we plan to investigate whether our MCMC sampling approach can be improved by initializing the sampling states in an informative way via ensemble training, building upon insights in [19, 36, 46].

#### Acknowledgments.

LW is supported by the DAAD programme Konrad Zuse Schools of Excellence in Artificial Intelligence, sponsored by the German Federal Ministry of Education and Research.