Adversarial Robustness Certification for Bayesian Neural Networks
We study the problem of certifying the robustness of Bayesian neural networks
(BNNs) to adversarial input perturbations. Given a compact set of input points
$T \subseteq \mathbb{R}^m$ and a set of output points $S \subseteq
\mathbb{R}^n$, we define two notions of robustness for BNNs in an adversarial
setting: probabilistic robustness and decision robustness. Probabilistic
robustness is the probability that for all points in $T$ the output of a BNN
sampled from the posterior is in $S$. On the other hand, decision robustness
considers the optimal decision of a BNN and checks if for all points in $T$ the
optimal decision of the BNN for a given loss function lies within the output
set $S$. Although exact computation of these robustness properties is
challenging due to the probabilistic and non-convex nature of BNNs, we present
a unified computational framework for efficiently and formally bounding them.
Our approach is based on weight interval sampling, integration, and bound
propagation techniques, and can be applied to BNNs with a large number of
parameters, and independently of the (approximate) inference method employed to
train the BNN. We evaluate the effectiveness of our methods on various
regression and classification tasks, including an industrial regression
benchmark, MNIST, traffic sign recognition, and airborne collision avoidance,
and demonstrate that our approach enables certification of robustness and
uncertainty of BNN predictions.
Matthew Wicker, Andrea Patane, Luca Laurenti, Marta Kwiatkowska
# Adversarial Robustness Certification

###### Abstract

We study the problem of certifying the robustness of Bayesian neural networks (BNNs) to adversarial input perturbations. Given a compact set of input points \(T\in\mathbb{R}^{n}\) and a set of output points \(S\subseteq\mathbb{R}^{n}\), we define two notions of robustness for BNNs in an adversarial setting; probabilistic robustness and decision robustness. Probabilistic robustness is the probability that for all points in \(T\) the output of a BNN sampled from the posterior is in \(S\). On the other hand, decision robustness considers the optimal decision of a BNN and checks if for all points in \(T\) the optimal decision of the BNN for a given loss function lies within the output set \(S\). Although exact computation of these robustness properties is challenging due to the probabilistic and non-convex nature of BNNs, we present a unified computational framework for efficiently and formally bounding them. Our approach is based on weight interval sampling, integration, and bound propagation techniques, and can be applied to BNNs with a large number of parameters, and independently of the (approximate) inference method employed to train the BNN. We evaluate the effectiveness of our methods on various regression and classification tasks, including an industrial regression benchmark, MNIST, traffic sign recognition, and airborne collision avoidance, and demonstrate that our approach enables certification of robustness and uncertainty of BNN predictions.

 Certification, Bayesian Neural Networks, Adversarial Robustness, Classification, Regression, Uncertainty

## I Introduction

While neural networks (NNs) regularly obtain state-of-the-art performance in many supervised machine learning problems [1, 2], they have been found to be vulnerable to adversarial attacks, i.e., imperceptible modifications of their inputs that trick the model into making an incorrect prediction [3]. Along with several other vulnerabilities [4], the discovery of adversarial examples has made the deployment of NNs in real-world, safety-critical applications - such as autonomous driving or healthcare - increasingly challenging. The design and analysis of methods that can mitigate such vulnerabilities of NNs, or provide guarantees for their worst-case behaviour in adversarial conditions, has thus become of critical importance [5, 6].

While retaining the advantages intrinsic to deep learning, Bayesian neural networks (BNNs), i.e., NNs with a probability distribution placed over their weights and biases [7], enable probabilistically principled evaluation of model uncertainty. Since adversarial examples are intuitively related to uncertainty [8], the application of BNNs is particularly appealing in safety-critical scenarios. In fact, model uncertainty of a BNN can, in theory, be taken into account at prediction time to enable safe decision-making [9, 10, 11, 12]. Various techniques have been proposed for the evaluation of their robustness, including generalisation of gradient-based adversarial attacks (i.e., non

Fig. 1: Certifications for a traffic sign recognition benchmark with two classes: speed limit (spd. lmt.) and warning sign (wam.). We plot original images, the upper and lower-bound class probabilities as red and blue horizontal lines, respectively, and a description of the result. **Top Row:** A 50 km/hr sign from the GTSRB dataset. As the lower bound class probability is 0.81, we certify that all images in the ball are classified correctly as speed limit signs and therefore no adversarial examples exist. **Bottom Row:** A nonsense traffic sign. As the upper bound probability for all classes is less than a threshold (0.75), we certify that the BNN is uncertain.
Bayesian) [13], statistical verification techniques [14], as well as approaches based on pointwise (i.e., for a specific test point \(x^{*}\)) uncertainty evaluation [15]. However, to the best of our knowledge, a systematic approach for computing formal (i.e., with certified bounds) guarantees on the behaviour of BNNs and their decisions against adversarial input perturbations is still missing.

In this work, we develop a novel algorithmic framework to quantify the adversarial robustness of BNNs. In particular, following existing approaches for quantifying the robustness of deterministic neural networks [16, 17, 18], we model adversarial robustness as an _input-output specification_ defined by a given compact set of input points \(T\subseteq\mathbb{R}^{m}\) and a given convex polytope output set \(S\subseteq\mathbb{R}^{n}\). A neural network satisfies this specification if all points in \(T\) are mapped into \(S\), called a safe set. Modelling specifications in this way encompasses many other practical properties such as classifier monotonicity [19] and individual fairness [20]. For a particular specification, we focus on two main properties of a BNN of interest for adversarial prediction settings: _probabilistic robustness_[21, 14] and _decision robustness_[18, 22]. The former, probabilistic robustness, is defined as the probability that a network sampled from the posterior distribution is robust (e.g., satisfies a robustness specification defined by a given \(T\) and \(S\)). Probabilistic robustness attempts to provide a general measure of robustness of a BNN; in contrast, _decision robustness_ focuses on the decision step, and evaluates the robustness of the optimal decision of a BNN. That is, a BNN satisfies decision robustness for a property if, for all points in \(T\), the expectation of the output of the BNN in the case of regression, or the argmax of the expectation of the softmax w.r.t. the posterior distribution for classification, are contained in \(S\).

Unfortunately, evaluating probabilistic and decision robustness for a BNN is not trivial, as it involves computing distributions and expectations of high-dimensional random variables passed through a non-convex function (the neural network architecture). Nevertheless, we derive a unified algorithmic framework based on computations over the BNN weight space that yields _certified lower_ and _upper bounds_ for both properties. Specifically, we show that probabilistic robustness is equivalent to the measure, w.r.t. the BNN posterior, of the set of weights for which the resulting deterministic NN is _robust_, i.e., it maps all points of \(T\) to a subset of \(S\). Computing upper and lower bounds for the probability involves sampling compact sets of weights according to the BNN posterior, and propagating each of these weight sets, \(H\), through the neural network architecture, jointly with the input region \(T\), to check whether all the networks instantiated by weights in \(H\) are safe. To do so, we generalise bound propagation techniques developed for deterministic neural networks to the Bayesian settings and instantiate explicit schemes for Interval Bound Propagation (IBP) and Linear Bound Propagation (LBP) [23]. Similarly, in the case of decision robustness, we show that formal bounds can be obtained by partitioning the weight space into different weight sets, and for each weight set \(J\) of the partition we again employ bound propagation techniques to compute the maximum and minimum of the decision of the NN for all input points in \(T\) and different weight configurations in \(J\). The resulting extrema are then averaged according to the posterior measure of the respective weight sets to obtain sound lower and upper bounds on decision robustness.

We perform a systematic experimental investigation of our framework on a variety of tasks. We first showcase the behaviour of our methodology on a classification problem from an airborne collision avoidance benchmark [24] and on two safety-critical industrial regression benchmarks [25]. We then consider image recognition tasks and illustrate how our method can scale to verifying BNNs on medium-sized computer vision problems, including MNIST and a two-class subset of the German Traffic Sign Recognition Benchmark (GTSRB) dataset [26]. On small networks, such as those used for airborne collision avoidance (\(\sim\)\(5000\) parameters), our method is able to verify key properties in under a second, thus enabling comprehensive certification over a fine partition of the entire state space. Moreover, when employed in conjunction with adversarial training [27], we are able to obtain non-trivial certificates for convolutional NNs (471,000 parameters) on full-colour GTSRB images (2,352 dimensions).1 As an example, we demonstrate the bounds on decision robustness in Figure 1, where we plot the upper and lower bound class probabilities (in red and blue respectively) for a BNN trained on a two-class traffic sign recognition benchmark. The bounds are computed for all images within a \(\ell_{\infty}\) ball with radius 2/255 of the two images in the left column of the figure. For the top image of a speed limit sign, our lower bound allows us to verify that the all images within the 2/255 are correctly classified by the BNN as a 50 km/hr sign. For the bottom image of a nonsense traffic sign, our upper bound allows us to verify that the BNN is uncertain for this image and all images in the ball.

Footnote 1: An implementation to reproduce all the experiments can be found at: [https://github.com/matthewwicker/AdversarialRobustnessCertificationForBNNs](https://github.com/matthewwicker/AdversarialRobustnessCertificationForBNNs).

In summary, this paper makes the following contributions.

* We present an algorithmic framework based on convex relaxation techniques for the robustness analysis of BNNs in adversarial settings.
* We derive explicit lower- and upper-bounding procedures based on IBP and LBP for the propagation of input and weight intervals through the BNN posterior function.
* We empirically show that our method can be used to certify BNNs consisting of multiple hidden layers and with hundreds of neurons per layer.

A preliminary version of this paper appeared as [21]. This work extends [21] in several aspects. In contrast to [21], which focused only on probabilistic robustness, here we also tackle decision robustness and embed the calculations for the two properties in a common computational framework. Furthermore, while the method in [21] only computes lower bounds, in this paper we also develop a technique for upper bounds computation. Finally, we substantially extend the empirical evaluation to include additional datasets, evaluation of convolutional architectures and scalability analysis, as well as certification of out-of-distribution uncertainty.
## II Related Work

Bayesian uncertainty estimates have been shown to empirically flag adversarial examples, often with remarkable success [15, 28]. These techniques are, however, empirical and can be circumvented by specially-tailored attacks that also target the uncertainty estimation [29]. Despite these attacks, it has been shown that BNN posteriors inferred by Hamiltonian Monte Carlo tend to be more robust to attacks than their deterministic counterparts [10]. Further, under idealised conditions of infinite data, infinitely-wide neural networks and perfect training, BNNs are provably robust to gradient-based adversarial attacks [11]. However, while showing that BNNs are promising models for defending against adversarial attacks, the arguments in [10] and [11] do not provide concrete bounds or provable guarantees for when an adversarial example does not exist for a given BNN posterior.

In [9, 14], the authors tackle similar properties of BNNs to those discussed in this paper. Yet these methods only provide bounds on probabilistic robustness and the bounds are _statistical_, i.e., only valid up to a user-defined, finite probability \(1-\delta\), with \(\delta>0\). In contrast, the method in this paper covers both probabilistic and decision robustness and computes bounds that are sound for the whole BNN posterior (i.e., hold with probability \(1\)). In [27], the authors incorporate worst-case information via bound propagation into the likelihood in order to train BNNs that are more adversarially robust; while that work develops a principled defense for BNNs against attack, it does not develop or study methods for analyzing or guaranteeing their robustness.

Since the publication of our preliminary work [21], the study of [30] has further investigated certifying the posterior predictive of BNNs. The definition in [30] corresponds to a subset of what we refer to as _decision_ robustness, but their method only applies to BNNs whose posterior support has been clipped to be in a finite range. Here, we pose a more general problem of certifying decision and probabilistic robustness of BNNs, and can handle posteriors on continuous, unbounded support, which entails the overwhelming majority of those commonly employed for BNNs. Furthermore, following the preliminary version of this paper [21, 31] introduced a technique for probabilistic robustness certification implemented via a recursive algorithm the that operates over the state-space of a model-based control scenario. [32] uses similar methods to those presented in [21] to study infinite-time horizon robustness properties of BNN control policies by checking for safe weight sets and modifying the posterior so that only safe weights have non-zero posterior support.

Most existing certification methods in the literature are designed for deterministic NNs. Approaches studied include abstract interpretation [23], mixed integer linear programming [33, 34, 35, 36], game-based approaches [37, 38], and SAT/SMT [24, 39]. In particular, [40, 18, 41] employ relaxation techniques from non-convex optimisation to compute guarantees over deterministic neural network behaviours, specifically using Interval Bound Propagation (IBP) and Linear Bound Propagation (LBP) approaches. However, these methods cannot be used for BNNs because they all assume that the weights of the networks are deterministic, i.e., fixed to a given value, while in the Bayesian setting we need to certify the BNN for a continuous range of values for weights that are not fixed, but distributed according to the BNN posterior.

In the context of Bayesian learning, methods to compute adversarial robustness measures have been explored for Gaussian processes (GPs), both for regression [42] and classification tasks [43, 44]. However, because of the non-linearity in NN architectures, GP-based approaches cannot be directly employed for BNNs. Furthermore, the vast majority of approximate Bayesian inference methods for BNNs do not employ Gaussian approximations over latent space [45]. In contrast, our method is specifically tailored to take into account the non-linear nature of BNNs and can be directly applied to a range of approximate Bayesian inference techniques used in the literature.

## III Background

In this section, we overview the necessary background and introduce the notation we use throughout the paper. We focus on neural networks (NNs) employed in a supervised learning scenario, where we are given a dataset of \(n_{\mathcal{D}}\) pairs of inputs and labels, \(\mathcal{D}=\{(x_{i},y_{i})\}_{i=1}^{n_{\mathcal{D}}}\), with \(x_{i}\in\mathbb{R}^{m}\), and where each target output \(y\in\mathbb{R}^{n}\) is either a one-hot class vector for classification or a real-valued vector for regression.

### _Bayesian Deep Learning_

Consider a feed forward neural network \(f^{w}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}\), parametrised by a vector \(w\in\mathbb{R}^{n_{w}}\) containing all its weights and biases. We denote with \(f^{w,1},...,f^{w,K}\) the \(K\) layers of \(f^{w}\) and take the activation function of the \(i\)th layer to be \(\sigma^{(i)}\), abbreviated to just \(\sigma\) in the case of the output activation.2 Throughout this paper, we will use \(f^{w}(x)\) to represent pre-activation of the last layer.

Footnote 2: We assume that the activation functions have a finite number of inflection points, which holds for activation functions commonly used in practice [46].

Bayesian learning of deep neural network starts with a prior distribution, \(p_{\mathbf{w}}(w)\), over the vector of random variables associated to the weights, \(\mathbf{w}\). Placing a distribution over the weights defines a stochastic process indexed by the input space, which we denote as \(f^{\mathbf{w}}\). After the data set \(\mathcal{D}\) has been observed, the BNN prior distribution is updated according to the likelihood, \(p(\mathcal{D}|w)=\prod_{i=1}^{n_{\mathcal{D}}}p(y_{i}|x_{i},w)\), which models how likely (probabilistically speaking) we are to observe an output under the stochasticity of our model parameters and observational noise [47]. The likelihood function, \(p(y_{i}|x_{i},w)\), generally takes the shape of a softmax for multiclass classification and a multivariate Gaussian for regression. The posterior distribution over the weights given the dataset is then computed by means of the Bayes formula, i.e., \(p(w|\mathcal{D})\propto p(\mathcal{D}|w)p(w)\). The cumulative distribution of \(p(w|\mathcal{D})\) we denote as \(P(\cdot)\), so that for \(R\subseteq\mathbb{R}^{n_{w}}\) we have:

\[\int_{R}p(w|\mathcal{D})dw=P(R). \tag{1}\]
The posterior \(p(w|\mathcal{D})\) is in turn used to calculate the output of a BNN on an unseen point, \(x^{*}\). The distribution over outputs is called the posterior predictive distribution and is defined as:

\[p(y^{*}|x^{*},\mathcal{D})=\int p(y^{*}|x^{*},w)p(w|\mathcal{D})dw. \tag{2}\]

Equation (2) defines a distribution over the BNN output. When employing a Bayesian model, the overall final prediction is taken to be a single value, \(\hat{y}\), that minimizes the Bayesian risk of an incorrect prediction according to the posterior predictive distribution and a loss function \(\mathcal{L}\). Formally, the final decision of a BNN is computed as

\[\hat{y}=\operatorname*{arg\,min}_{y}\int_{\mathbb{R}^{n}}\mathcal{L}(y,y^{*})p (y^{*}|x^{*},\mathcal{D})dy^{*}.\]

This minimization is the subject of Bayesian decision theory [48], and the final form of \(\hat{y}\) clearly depends on the specific loss function \(\mathcal{L}\) employed in practice. In this paper, we focus on two standard loss functions widely employed for classification and regression problems.3

Footnote 3: In Appendix B we discuss how our method can be generalised to other losses commonly employed in practice.

Classification DecisionsThe 0-1 loss, \(\ell_{0-1}\), assigns a penalty of 0 to the correct prediction, and 1 otherwise. It can be shown that the optimal decision in this case is given by the class for which the predictive distribution obtains its maximum, i.e.:

\[\hat{y}=\operatorname*{arg\,max}_{i=1,\ldots,n}p_{i}(y^{*}|x^{*},\mathcal{D}) =\operatorname*{arg\,max}_{i=1,\ldots,n}\mathbb{E}_{w\sim p(w|\mathcal{D})} \left[\sigma_{i}(f^{w}(x))\right],\]

where \(\sigma_{i}\) represents the \(i\)th output component of the softmax function.

Regression DecisionsThe \(\ell_{2}\) loss assigns a penalty to a prediction according to its \(\ell_{2}\) distance from the ground truth. It can be shown that the optimal decision in this case is given by the expected value of the BNN output over the posterior distribution, i.e.:

\[\hat{y}=\mathbb{E}_{w\sim p(w|\mathcal{D})}\left[f^{w}(x)\right].\]

Unfortunately, because of the non-linearity of neural network architectures, the computation of the posterior distribution over weights, \(p(w|\mathcal{D})\), is generally intractable [7]. Hence, various approximation methods have been studied to perform inference with BNNs in practice. Among these methods, we will consider Hamiltonian Monte Carlo (HMC) [7] and Variational Inference (VI) [45], which we now briefly describe.

#### Iii-B1 Hamiltonian Monte Carlo (HMC)

HMC proceeds by defining a Markov chain whose invariant distribution is \(p_{\mathbf{w}}(w|\mathcal{D}),\) and relies on Hamiltonian dynamics to speed up the exploration of the space. Differently from VI discussed below, HMC does not make parametric assumptions on the form of the posterior distribution, and is asymptotically correct [7]. The result of HMC is a set of samples that approximates \(p_{\mathbf{w}}(w|\mathcal{D})\).

#### Iii-B2 Variational Inference (VI)

VI proceeds by finding a Gaussian approximating distribution over the weight space \(q(w)\sim p_{\mathbf{w}}(w|\mathcal{D})\) in a trade-off between approximation accuracy and scalability. The core idea is that \(q(w)\) depends on some hyperparameters that are then iteratively optimized by minimizing a divergence measure between \(q(w)\) and \(p_{\mathbf{w}}(w|\mathcal{D})\). Samples can then be efficiently extracted from \(q(w)\).

For simplicity of notation, in the rest of the paper we will indicate with \(p(w|\mathcal{D})\) the posterior distribution estimated by either of the two methods, and clarify the methodological differences when they arise.

## IV Problem Statements

We focus on local specifications defined over an input compact set \(T\subseteq\mathbb{R}^{m}\) and output set \(S\subseteq\mathbb{R}^{n}\) in the form of a convex polytope:

\[S=\{y\in\mathbb{R}^{n}\,|\,C_{S}y+d_{S}\geq 0\}, \tag{3}\]

where \(C_{S}\in\mathbb{R}^{n_{S}xn}\) and \(d_{S}\in\mathbb{R}^{n_{S}}\) are the matrix and vector encoding the polytope constraints, and \(n_{S}\) is the number of output constraints considered. For simplicity of presentation, we assume that \(T\) is defined as a box (axis-aligned linear constraints).4 However, we stress that all the methods in this paper can be extended to the more general case where \(T\) is a convex polytope. Our formulation of input-output specifications can be used to capture important properties such as classifier monotonicity [49] and individual fairness [20], but in this work we focus exclusively on adversarial robustness. Targeted adversarial robustness, where one aims to force the neural network into a particular wrong classification, is captured in this framework by setting \(T\) to be an over-approximation of an \(\ell_{p}\) ball around a given test input, and setting \(C_{S}\) to an \(n_{S}\times n\) matrix of all zeros with a \(-1\) entry in the diagonal entry corresponding to the true class and a \(1\) on the diagonal entry corresponding to the target class or classes. For regression, one uses \(C_{S}\) to encode the absolute deviation from the target value and \(d_{S}\) to encode the maximum tolerable deviation. Throughout the paper we will refer to an input-output set pair, \(T\) and \(S\), as defined above as a _robustness specification_.

Footnote 4: Note that, where a specification is not in this form already, one can first compute a bounding box \(R=[x^{L},x^{U}]\) (or a finite sequence of them) such that \(T\subseteq R\), and then proving that the output specification holds for \(R\) also proves that it holds for \(T\). In the case that we do not prove that an output specification holds, then we cannot guarantee it is violated by nature of our method being sound but not complete.

### _Probabilistic Robustness_

Probabilistic robustness accounts for the probabilistic behaviour of a BNN in adversarial settings.

**Definition 1** (Probabilistic robustness).: _Given a Bayesian neural network \(f^{\mathbf{w}}\), an input set \(T\subseteq\mathbb{R}^{m}\) and a set \(S\subseteq\mathbb{R}^{n}\) of safe outputs, then probabilistic robustness is defined as_

\[P_{\text{safe}}(T,S):=Prob_{w\sim p(w|\mathcal{D})}\big{(}\forall x\in T,f^{w} (x)\in S\big{)}. \tag{4}\]
Given \(\eta\in[0,1]\), we then say that \(f^{\mathbf{w}}\) is probabilistically robust, or safe, for specifications \(T\) and \(S\), with probability at least \(\eta\) iff:_

\[P_{\text{safe}}(T,S)\geq\eta\]

Probabilistic robustness considers the adversarial behaviour of the model while accounting for the uncertainty arising from the posterior distribution. \(P_{\text{safe}}(T,S)\) represents the (weighted) proportion of neural networks sampled from \(f^{\mathbf{w}}\) that satisfy a given input-output specification (captured by \(T\) and \(S\)) and can be used directly as a measure of compliance for Bayesian neural networks. As such, probabilistic robustness is particularly suited to quantification of the robustness of a BNN to adversarial perturbations [9, 22, 50]. Exact computation of \(P_{\text{safe}}(T,S)\) is hindered by both the size and non-linearity of neural networks. As \(P_{\text{safe}}(T,S)\) cannot be computed exactly for general BNNs, in this work we tackle the problem of computing provable bounds on probabilistic robustness.

**Problem 1** (Bounding probabilistic robustness).: _Given a Bayesian neural network \(f^{\mathbf{w}}\), an input set \(T\subseteq\mathbb{R}^{m}\) and a set \(S\subseteq\mathbb{R}^{n}\) of safe outputs, compute (non-trivial) \(P_{\text{safe}}^{L}\) and \(P_{\text{safe}}^{U}\) such that_

\[P_{\text{safe}}^{L}\leq P_{\text{safe}}(T,S)\leq P_{\text{safe}}^{U}. \tag{5}\]

We highlight the difference between this problem definition and those discussed in prior works [9, 14]. In particular, prior works compute upper and lower bounds that hold with probability \(1-\delta\) for some pre-specified \(\delta\). While such statistical bounds can provide an estimation for \(P_{\text{safe}}(T,S)\), these may not be sufficient in safety-critical contexts where strong, worst-case guarantees over the full behaviour of the BNN are necessary. The problem statement above holds with probability \(1\) and represents sound guarantees on \(P_{\text{safe}}(T,S)\).

### _Decision Robustness_

While \(P_{\text{safe}}\) attempts to measure the compliance of all functions in the support of a BNN posterior, we are often interested in evaluating robustness w.r.t. a specific decision. In order to do so, we consider _decision robustness_, which is computed over the final decision of the BNN. In particular, given a loss function and a decision \(\hat{y}\) we have the following.

**Definition 2** (Decision robustness).: _Consider a Bayesian neural network \(f^{\mathbf{w}}\), an input set \(T\subseteq\mathbb{R}^{m}\) and a set \(S\subseteq\mathbb{R}^{n}\) of safe outputs. Assume that the decision for a loss \(\mathcal{L}\) for \(x\in\mathbb{R}^{m}\) is given by \(\hat{y}(x)\). Then, the Bayesian decision is considered to be robust if:_

\[\forall x\in T\quad\hat{y}(x)\in S. \tag{6}\]

We notice that, since the specific form of the decision depends on the loss employed in practice, the definition of decision robustness takes different form depending on whether the BNN is used for classification or for regression. In particular, we instantiate the definition for the two cases of standard loss discussed in Section III.

In the regression case, using the mean square loss we have that \(\hat{y}(x)=\mathbb{E}[f^{w}(x)]\), so that if we find upper and lower bounds on \(\mathbb{E}[f^{w}(x)]\) for all \(x\in T\), i.e., for \(i=1,...,m\):

\[D_{\text{safe},i}^{L}\leq\min_{x\in T}\mathbb{E}\left[f_{i}^{w}(x)\right],\;D_ {\text{safe},i}^{U}\geq\max_{x\in T}\mathbb{E}\left[f_{i}^{w}(x)\right],\]

we can then simply check whether these are within \(S\).

For the classification case, where the decision is given by the \(\arg\max\) of the predictive posterior, note that, in order to check the condition in Definition 2, it suffices to compute lower and upper bounds on the posterior predictive in \(T\), i.e.:

\[D_{\text{safe},i}^{L}\leq\min_{x\in T}\mathbb{E}\left[\sigma_{i}(f^{w}(x)) \right],\;D_{\text{safe},i}^{U}\geq\max_{x\in T}\mathbb{E}\left[\sigma_{i}(f ^{w}(x))\right],\]

for \(i=1,\ldots,m\). It is easy to see that the knowledge of \(D_{\text{safe},i}^{L}\) and \(D_{\text{safe},i}^{U}\) for all \(i=1,\ldots,m\) can be used to provide guarantees of the model decision belonging to \(S\), as defined in Equation (3), by simply propagating these bounds through the equations. Therefore, for both classification and regression we have to bound an expectation of the BNN output over the posterior distribution, with the additional softmax computations for classification. We thus arrive at the following problem for bounding decision robustness.

**Problem 2** (Bounding decision robustness).: _Let \(f^{\mathbf{w}}\) be a BNN with posterior distribution \(p(w|\mathcal{D})\). Consider an input-output

Fig. 2: A diagram illustrating a single iteration of the computational flow for the certification process of a BNN w.r.t. decision robustness (green) and probabilistic robustness (purple). This process is summarised in Algorithm 1.
specification (\(T\), \(S\)) and assume \(\mathcal{L}=\ell_{0-1}\) for classification or \(\mathcal{L}=\ell_{2}\) for regression. We aim at computing (non-trivial) lower and upper bounds \(D^{L}_{\text{safe}}\) and \(D^{U}_{\text{safe}}\) such that:_

\[D^{L}_{\text{safe}}\lesssim\mathbb{E}[s(f^{w}(x))]\leq D^{U}_{\text{safe}} \quad\forall x\in T,\]

_where \(s=\sigma\) for classification and \(s=\mathbb{I}\) for regression._

Note that, while for regression we bound the decision directly, for classification we compute the bounds on the predictive posterior and use these to compute bounds on the final decision. This is similar to what is done for deterministic neural networks, where in the case of classification the bounds are often computed over the logit, and then used to provide guarantees for the final decision [18]. As with probabilistic robustness, our bounds on decision robustness are sound guarantees and do not have a probability of error as with statistical bounds.

### _Outline of our Approach:_

We design an algorithmic framework for worst-case and best-case bounds on local robustness properties in Bayesian neural networks, taking account of both the posterior distribution (\(P^{L}_{\text{safe}}\) and \(P^{U}_{\text{safe}}\)) and the overall model decision (\(D^{L}_{\text{safe}}\) and \(D^{L}_{\text{safe}}\)). First, we show how the two robustness properties of Definitions 1 and 2 can be reformulated in terms of computation over weight intervals. This allows us to derive a unified approach to the bounding of the robustness of the BNN posterior (i.e., probabilistic robustness) and of the robustness of the overall model decision (i.e., decision robustness) that is based on _bound propagation_ and _posterior integral_ computation over hyper-rectangles. A visual outline for our framework is presented in Figure 2. We organise the presentation of our framework by first introducing a general theoretical framework for bounding the robustness quantities of interest (Section V). We will then show how the required integral computations can be achieved for Bayesian posterior inference techniques commonly used in practice (Section VI-A). Hence, we will show how to extend bound propagation techniques to deal with both input variable intervals and intervals over the weight space, and will instantiate approaches based on Interval and Linear Bound Propagation techniques (Section VI-B). Finally (Section VII), we will present an overall algorithm that produces the desired bounds.

## V Formulating BNN Adversarial Robustness via Weight Sets

In this section, we show how a single computational framework can be leveraged to compute bounds on both definitions of BNN robustness. We start by converting the computation of robustness into weight space and then define a family of weight intervals that we leverage to bound the integrations required by both definitions. Interestingly, we find that the resulting theoretical bounds in both cases depend on the same quantities. Proofs for the main results in this section are presented in Appendix C.

### _Bounding Probabilistic Robustness_

We first show that the computation of \(P_{\text{safe}}(T,S)\) is equivalent to computing a maximal set of safe weights \(H\) such that each network associated to weights in \(H\) is safe w.r.t. the robustness specification at hand.

**Definition 3** (Maximal safe and unsafe sets).: _We say that \(H\subseteq\mathbb{R}^{n_{w}}\) is the maximal safe set of weights from \(T\) to \(S\), or simply the maximal safe set of weights, iff \(H=\{w\in\mathbb{R}^{n_{w}}\,|\,\forall x\in T,f^{w}(x)\in S\}\). Similarly, we say that \(K\subseteq\mathbb{R}^{n_{w}}\) is the maximal unsafe set of weights from \(T\) to \(S\), or simply the maximal unsafe set of weights, iff \(K=\{w\in\mathbb{R}^{n_{w}}\,|\,\exists x\in T,f^{w}(x)\not\in S\}\)._

Intuitively, \(H\) and \(K\) simply encode the input-output specifications \(S\) and \(T\) in the BNN weight space. The following lemma, which trivially follows from Equation 4, allows us to directly relate the maximal sets of weights to the probability of robustness.

**Lemma 1**.: _Let \(H\) and \(K\) be the maximal safe and unsafe sets of weights from \(T\) to \(S\). Assume that \(w\sim p(w|\mathcal{D})\). Then, it holds that_

\[P(H)=\int_{H}p(w|\mathcal{D})dw=P_{\text{safe}}(T,S)= \tag{7}\] \[1-\int_{K}p(w|\mathcal{D})dw=1-P(K).\]

Lemma 1 simply translates the robustness specification from being concerned with the input-output behaviour of the BNN to an integration on the weight space.

An exact computation of sets \(H\) and \(K\) is infeasible in general. However, we can easily compute subsets of \(H\) and \(K\). Such subsets can then be used to compute upper and lower bounds on the value of \(P_{\text{safe}}(T,S)\) by considering subsets of the maximal safe and unsafe weights.

**Definition 4** (Safe and unsafe sets).: _Given a maximal safe set \(H\) or a maximal unsafe set \(K\) of weights, we say that \(\hat{H}\) and \(\hat{K}\) are a safe and unsafe set of weights from \(T\) to \(S\) iff \(\hat{H}\in H\) and \(\hat{K}\subseteq K\), respectively._

\(\hat{H}\) and \(\hat{K}\) can include _any_ safe and unsafe weights, respectively, without requiring they are _maximal_. Without maximality, we no longer have strict equality in Lemma 1, but instead we arrive at bounds on the value of probabilistic robustness.

We proceed by defining \(\hat{H}\) and \(\hat{K}\) as the union of a family of disjoint weight intervals, as these can provide flexible approximations of \(H\) and \(K\). That is, we consider \(\mathcal{H}=\{H_{i}\}_{i=1}^{n_{H}}\), with \(H_{i}=[w_{i}^{L,H},w_{i}^{U,H}]\) and \(\mathcal{K}=\{K_{i}\}_{i=1}^{n_{K}}\), with \(K_{i}=[w_{i}^{L,K},w_{i}^{U,K}]\) such that \(H_{i}\in H\) and \(K_{i}\subset K\), \(\hat{H}=\bigcup_{i=1}^{n_{H}}H_{i}\), \(K=\bigcup_{i=1}^{n_{K}}K_{i}\), and \(H_{i}\cap H_{j}=\emptyset\) and \(K_{i}\cap K_{j}=\emptyset\), for any \(i\neq j\). Hence, as a consequence of Lemma 1, and by the fact that \(\hat{H}=\bigcup_{i=1}^{n_{H}}H_{i}\subset H\) and \(\hat{K}=\bigcup_{i=1}^{n_{K}}K_{i}\subset K\), we obtain the following.

**Proposition 1** (Bounds on probabilistic robustness).: _Let \(H\) and \(K\) be the maximal safe and unsafe sets of weights from \(T\) to \(S\). Consider two families of pairwise disjoint weight intervals \(\mathcal{H}=\{H_{i}\}_{i=1}^{n_{H}}\), \(\mathcal{K}=\{K_{i}\}_{i=1}^{n_{K}}\), where for all \(i\):_

\[H_{i}\subseteq H,\quad K_{i}\subseteq K. \tag{8}\]
_Let \(\hat{H}\subseteq H\) and \(\hat{K}\subseteq K\) be non-maximal safe and unsafe sets of weights, with \(\hat{H}=\bigcup_{i=1}^{n_{H}}H_{i}\) and \(\hat{K}=\bigcup_{i=1}^{n_{K}}K_{i}\). Assume that \(w\sim p(w|\mathcal{D})\). Then, it holds that_

\[P_{\text{safe}}^{L}:=\sum_{i=1}^{n_{H}}P(H_{i})\leq P_{\text{safe}}(T,S,\mathbf{ w})\leq 1-\sum_{i=1}^{n_{K}}P(K_{i})=:P_{\text{safe}}^{U}, \tag{9}\]

_that is, \(P_{\text{safe}}^{L}\) and \(P_{\text{safe}}^{U}\) are, respectively, lower and upper bounds on probabilistic robustness._

Through the use of Proposition 1 we can thus bound probabilistic robustness by performing computation over sets of safe and unsafe intervals. Note that the bounds are given in the case where \(\mathcal{H}\) and \(\mathcal{K}\) are families of pairwise disjoint weight sets. The general case can be tackled by using the Bonferroni bound, which is discussed in Appendix D-D for hyper-rectangular weight sets.

Before explaining in detail how such bounds can be explicitly computed, first, in the next section, we show how a similar derivation leads us to analogous bounds and computations for decision robustness.

### _Bounding Decision Robustness_

The key difference between our formulation of probabilistic robustness and that of decision robustness is that, for the former, we are only interested in the behaviour of neural networks extracted from the BNN posterior that satisfy the robustness requirements (hence the distinction between \(H\)- and \(K\)-weight intervals), whereas for the computation of bounds on decision robustness we need to take into account the overall worst-case behaviour of an expected value computed for the BNN predictive distribution in order to compute sound bounds. As such, rather than computing safe and unsafe sets, we only need a family of weight sets, \(\mathcal{J}=\{J_{i}\}_{i=1}^{n_{J}}\), and rely on that for bounding \(D_{\text{safe}}(T,S)\). We explicitly show how this can be done for classification with likelihood \(\sigma\). The bound for regression follows similarly by using the identity function as \(\sigma\).

**Proposition 2** (Bounding decision robustness).: _Let \(\mathcal{J}=\{J_{i}\}_{i=1}^{n_{J}}\), with \(J_{i}\subset\mathbb{R}^{n_{\omega}}\) be a family of disjoint weight intervals. Let \(\sigma^{L}\) and \(\sigma^{U}\) be vectors that lower and upper bound the domain of the final activation function, and \(c\in\{1,\dots,m\}\) an index spanning the BNN output dimension. Define:_

\[D_{\text{safe},c}^{L} :=\sum_{i=1}^{n_{J}}P(J_{i})\min_{\begin{subarray}{c}x\in T\\ w\in J_{i}\end{subarray}}\sigma_{c}(f^{w}(x))+\sigma^{L}\Bigg{(}1-\sum_{i=1}^ {n_{J}}P(J_{i})\Bigg{)} \tag{10}\] \[D_{\text{safe},c}^{U} :=\sum_{i=1}^{n_{J}}P(J_{i})\max_{\begin{subarray}{c}x\in T\\ w\in J_{i}\end{subarray}}\sigma_{c}(f^{w}(x))+\sigma^{U}\Bigg{(}1-\sum_{i=1}^ {n_{J}}P(J_{i})\Bigg{)}. \tag{11}\]

_Consider the vectors \(D_{\text{safe}}^{L}=[\,D_{\text{safe},1}^{L},\dots,D_{\text{safe},m}^{L}]\) and \(D_{\text{safe}}^{U}=[\,D_{\text{safe},1}^{U},\dots,D_{\text{safe},m}^{U}]\), then it holds that:_

\[D_{\text{safe}}^{L}\leq\mathbb{E}_{p(w|\mathcal{D})}[\sigma(f^{w}(x))]\leq D_ {\text{safe}}^{U}\quad\forall x\in T,\]

_that is, \(D_{\text{safe}}^{L}\) and \(D_{\text{safe}}^{U}\) are lower and upper bounds on the predictive distribution in \(T\)._

Intuitively, the first terms in the bounds of Equations (10) and (11) consider the worst-case output for the input set \(T\) and each interval \(J_{i}\), while the second term accounts for the worst-case value of the posterior mass not captured by the family of intervals \(\mathcal{J}\) by taking a coarse, overall bound on that region. The provided bound is valid for any family of intervals \(\mathcal{J}\). Ideally, however, the partition should be finer around regions of high probability mass of the posterior distribution, as these make up the dominant term in the computation of the posterior predictive. We will discuss in Section VI how we select these intervals in practice so as to empirically obtain non-vacuous bounds.

### _Computation of the Bounds_

We now propose a unified approach to computing these lower and upper bounds. We first observe that the bounds in Equations (9), (10) and (11) require the integration of the posterior distribution over weight intervals, i.e., \(P(H_{i})\), \(P(K_{i})\) and \(P(J_{i})\). While this is in general intractable, we have built the bound so that \(H_{i}\), \(K_{i}\) and \(J_{i}\) are axis-aligned hyper-rectangles, and so the computation can be done exactly for approximate Bayesian inference methods used in practice. This will be the topic of Section VI-A, where, given a rectangle in weight space of the form \(R=[w^{L},w^{U}]\), we will show how to compute \(P(R)=\int_{R}p(w|\mathcal{D})dw\).

For the explicit computation of decision robustness, the only missing ingredient is then the computation of the minimum and maximum \(\sigma(f^{w}(x))\) for \(x\in T\) and \(w\in J_{i}\). We do this by bounding the BNN output for any given rectangle in the weight space \(R\). That is, we will compute upper and lower bounds \(y^{L}\) and \(y^{U}\) such that:

\[y^{L}\leq\min_{\begin{subarray}{c}x\in T\\ w\in R\end{subarray}}f^{w}(x)\quad y^{U}\geq\max_{\begin{subarray}{c}x\in T \\ w\in R\end{subarray}}f^{w}(x), \tag{12}\]

which can then be used to bound \(\sigma(f^{w}(x))\) by simple propagation over the softmax (if needed). The derivation of such bounds will be the subject of Section VI-B.

Finally, observe that, whereas for decision robustness we can simply select any weight interval \(J_{i}\), for probabilistic robustness one needs to make a distinction between safe sets (\(H_{i}\)) and unsafe sets (\(K_{i}\)). It turns out that this can be done by bounding the output of the BNN in each of these intervals. For example, in the case of the safe sets, by definition we have that \(\forall w\in H_{i},\forall x^{\prime}\in T\) it follows that \(f^{w}(x^{\prime})\in S\). By defining \(y^{L}\) and \(y^{U}\) as in Equation (12), we can see that it suffices to check whether \([y^{L},y^{U}]\subseteq S\). Hence, the computation of probabilistic robustness also depends on the computation of such bounds (again, discussed in Section VI-B).

Therefore, once we have shown how the computation of \(P(R)\) for any weight interval and \(y^{L}\) and \(y^{U}\) can be done, the bounds in Proposition 1 and Proposition 2 can be computed explicitly, and we can thus bound probabilistic and decision robustness. Section VII will assemble these results together into an overall computational flow of our methodology.

## VI Explicit Bound Computation

In this section, we provide details on the specific computations needed to calculate the theoretical bound presented in Section V for probabilistic and decision robustness. We start
by discussing how a weight intervals family can be generated in practice, and how to integrate over them in Section VI-A. In Section VI-B, we then derive a scheme based on convex-relaxation techniques for bounding the output of BNNs.

### _Integral Computation over Weight Intervals_

Key to the computation of the bounds derived in Section V is the ability to compute the integral of the posterior distribution over a combined set of weight intervals. Crucially, the shape of the weight sets \(\mathcal{H}=\{H_{i}\}_{i=1}^{n_{H}}\), \(\mathcal{K}=\{K_{i}\}_{i=1}^{n_{K}}\) and \(\mathcal{J}=\{J_{i}\}_{i=1}^{n_{J}}\) is a parameter of the method, so that it can be chosen to simplify the integral computation depending on the particular form of the approximate posterior distribution used. We build each weight interval as an axis-aligned hyper-rectangle of the form \(R=[w^{L},w^{U}]\) for \(w^{L}\) and \(w^{U}\in\mathbb{R}^{n_{w}}\).

Weight Intervals for Decision RobustnessIn the case of decision robustness it suffices to sample any weight interval \(J_{i}\) to compute the bounds we derived in Proposition 2. Clearly, the bound is tighter if the \(\mathcal{J}\) family is finer around the area of high probability mass for \(p(w|\mathcal{D})\). In order to obtain such a family we proceed as follows. First, we define a _weight margin_\(\gamma>0\) that has the role of parameterising the radius of the weight intervals we define. We then iteratively sample weight vectors \(w_{i}\) from \(p(w|\mathcal{D})\), for \(i=1,\ldots,n_{J}\), and finally define \(J_{i}=[w_{i}^{L},w_{i}^{U}]=[w_{i}-\gamma,w_{i}+\gamma]\). As such, thus defined weight intervals naturally hover around the area of greater density for \(p(w|\mathcal{D})\), while asymptotically covering the whole support of the distribution.

Weight Intervals for Probabilistic RobustnessOn the other hand, for the computation of probabilistic robustness one has to make a distinction between safe weight intervals \(H_{i}\) and unsafe ones \(K_{i}\). As explained in Section V-C, this can be done by bounding the output of the BNN in each of these intervals. For example, in the case of the safe sets, by definition, \(H_{i}\) is safe if and only if \(\forall w\in H_{i},\forall x^{\prime}\in T\) we have that \(f^{w}(x^{\prime})\in S\). Thus, in order to build a family of safe (respectively unsafe) weight intervals \(H_{i}\) (\(K_{i}\)), we proceed as follows. As for decision robustness, we iteratively sample weights \(w_{i}\) from the posterior used to build hyper-rectangles of the form \(R_{i}=[w_{i}-\gamma,w_{i}+\gamma]\). We then propagate the BNN through \(R\) and check whether the output of the BNN in \(R\) is (is not) a subset of \(S\). The derivation of such bounds on propagation will be the subject of Section VI-B.

Once the family of weights is computed, there remains the computation of the cumulative distribution over such sets. The explicit computations depend on the particular form of Bayesian approximate inference that is employed. We discuss explicitly the case of Gaussian variational approaches, and of sample-based posterior approximation (e.g., HMC), which entails the majority of the approximation methods used in practice [51].

Variational InferenceFor variational approximations, \(p(w|\mathcal{D})\) takes the form of a multi-variate Gaussian distribution over the weight space. The resulting computations reduce to the integral of a multi-variate Gaussian distribution over a finite-sized axis-aligned rectangle, which can be computed using standard methods from statistics [52]. In particular, under the common assumption of variational inference with a Gaussian distribution with diagonal covariance matrix [53], i.e., \(p(w|\mathcal{D})=\mathcal{N}(\mu,\Sigma)\), with \(\Sigma=\text{diag}(\Sigma_{1},\ldots,\Sigma_{n_{w}})\), we obtain the following result for the posterior integration:

\[P(R)=\int_{R} \!\!\!\!P(w|\mathcal{D})dw= \tag{13}\] \[\prod_{j=1}^{n_{w}}\frac{1}{2}\left(\text{erf}\left(\frac{\mu_{ j}-w_{i}^{L}}{\sqrt{2\Sigma_{j}}}\right)-\text{erf}\left(\frac{\mu_{j}-w_{i}^{u }}{\sqrt{2\Sigma_{j}}}\right)\right).\]

By plugging this into the bounds of Equation (9) for \(P(H_{i})\) and \(P(K_{i})\) for probabilistic robustness and in Equations (10) and (11) for decision robustness, one obtains a closed-form formula for the bounds given weight set interval families \(\mathcal{H}\), \(\mathcal{K}\) and \(\mathcal{J}\).

Sample-based approximationsIn the case of sample-based posterior approximation (e.g., HMC), we have that \(p(w|\mathcal{D})\) defines a distribution over a finite set of weights. In this case we can simplify the computations by selecting the weight margin \(\gamma=0\), so that each sampled interval will be of the form \(R=[w_{i},w_{i}]\) and its probability under the discrete posterior will trivially be:

\[P(R_{i})=p(w_{i}|\mathcal{D}). \tag{14}\]

### _Bounding Bayesian Neural Networks' Output_

Given an input specification, \(T\), and a weight interval, \(R=[w^{L},w^{U}]\), the second key step in computing probabilistic and decision robustness is the bounding of the output of the BNN over \(R\) given \(T\). That is, we need to derive methods to compute \([y^{L},y^{U}]\) such that, by construction, \(\forall w\in[w^{L},w^{U}]\), \(\forall x^{\prime}\in T\) it follows that \(f^{w}(x^{\prime})\in[y^{L},y^{U}]\).

In this section, we discuss interval bound propagation (IBP) and linear bound propagation (LBP) as methods for computing the desired output set over-approximations. Before discussing IBP and LBP in detail, we first introduce common notation for the rest of the section. We consider feed-forward neural networks of the form:

\[z^{(0)} =x \tag{15}\] \[\zeta_{i}^{(k+1)} =\sum_{j=1}^{n_{k}}W_{ij}^{(k)}z_{j}^{(k)}+b_{i}^{(k)} i=0,\ldots,n_{k+1}\] (16) \[z_{i}^{(k)} =\sigma(\zeta_{i}^{(k)}) i=0,\ldots,n_{k} \tag{17}\]

for \(k=1,\ldots,K\), where \(K\) is the number of hidden layers, \(\sigma(\cdot)\) is a pointwise activation function, \(W^{(k)}\in\mathbb{R}^{n_{k}\times n_{k-1}}\) and \(b^{(k)}\in\mathbb{R}^{n_{k}}\) are the matrix of weights and vector of biases that correspond to the \(k\)th layer of the network and \(n_{k}\) is the number of neurons in the \(k\)th hidden layer. Note that, while Equations (15)-(17) are written explicitly for fully-connected layers, convolutional layers can be accounted for by embedding them in fully-connected form [41].

We write \(W_{:}^{(k)}\) for the vector comprising the elements from the \(i\)th row of \(W^{(k)}\), and similarly \(W_{:j}^{(k)}\) for that comprising the elements from the \(j\)th column. \(\zeta^{(K+1)}\) represents the final output of the network (or the logit in the case of classification networks), that is, \(\zeta^{(K+1)}=f^{w}(x)\). We write \(W^{(k),L}\) and \(W^{(k),U}\) for the lower and upper bound induced
by \(R\) for \(W^{(k)}\) and \(b^{(k),L}\) and \(b^{(k),U}\) for those of \(b^{(k)}\), for \(k=0,\ldots,K\). Observe that \(z^{(0)}\), \(\zeta_{i}^{(k+1)}\) and \(z_{i}^{(k)}\) are all functions of the input point \(x\) and of the combined vector of weights \(w=[W^{(0)},b^{(0)},\ldots,W^{(K)},b^{(K)}]\). We omit the explicit dependency for simplicity of notation. Finally, we remark that, as both the weights and the input vary in a given set, Equation (16) defines a quadratic form.

Interval Bound Propagation (IBP)IBP has already been employed for fast certification of deterministic neural networks [18]. For a deterministic network, the idea is to propagate the input box around \(x\), i.e., \(T=[x^{L},x^{U}]\), through the first layer, so as to find values \(z^{(1),L}\) and \(z^{(1),U}\) such that \(z^{(1)}\in\left[z^{(1),L},z^{(1),U}\right]\), and then iteratively propagate the bound through each consecutive layer for \(k=1,\ldots,K\). The final box constraint in the output layer can then be used to check for the specification of interest [18]. The only adjustment needed in our setting is that at each layer we also need to propagate the interval of the weight matrix \([W^{(k),L},W^{(k),U}]\) and that of the bias vector \([b^{(k),L},b^{(k),U}]\). This can be done by noticing that the minimum and maximum of each term of the bi-linear form of Equation (16), that is, of each monomial \(W^{(k)}_{ij}z^{(k)}_{j}\), lies in one of the four corners of the interval \([W^{(k),L}_{ij},W^{(k),U}_{ij}]\times[z^{(k),L}_{j},z^{(k),U}_{j}]\), and by adding the minimum and maximum values respectively attained by \(b^{(k)}_{i}\). As in the deterministic case, interval propagation through the activation function proceeds by observing that generally employed activation functions are monotonic, which permits the application of Equation (17) to the bounding interval. Where monotonicity does not hold, we can bound any activation function that has finitely many inflection points by splitting the function into piecewise monotonic functions. This is summarised in the following proposition.

**Proposition 3**.: _Let \(f^{w}(x)\) be the network defined by the set of Equations (15)-(17), let for \(k=0,\ldots,K\):_

\[t^{(k),L}_{ij}=\min\{W^{(k),L}_{ij}z^{(k),L}_{j},W^{(k),U}_{ij}z^ {(k),L}_{j},\] \[W^{(k),L}_{ij}z^{(k),U}_{j},W^{(k),U}_{ij}z^{(k),U}_{j}\}\] \[t^{(k),U}_{ij}=\max\{W^{(k),L}_{ij}z^{(k),L}_{j},W^{(k),U}_{ij}z^ {(k),L}_{j},\] \[W^{(k),L}_{ij}z^{(k),U}_{j},W^{(k),U}_{ij}z^{(k),U}_{j}\}\]

_where \(i=1,\ldots,n_{k+1}\), \(j=1,\ldots,n_{k}\), and \(z^{(k),L}=\sigma(\zeta^{(k),L})\), \(z^{(k),U}=\sigma(\zeta^{(k),U})\) and:_

\[\zeta^{(k+1),L}=\sum_{j}t^{(k),L}_{ij}+b^{(k),L}\] \[\zeta^{(k+1),U}=\sum_{j}t^{(k),U}_{j}+b^{(k),U}.\]

_Then we have that \(\forall x\in T\) and \(\forall w\in R\):_

\[f^{w}(x)=\zeta^{(K+1)}\in\left[\zeta^{(K+1),L},\zeta^{(K+1),U}\right].\]

The proposition above, whose proof is Appendix C-C (Appendix C-C), yields a bounding box for the output of the neural network in \(T\) and \(R\).

Linear Bound Propagation (LBP)We now discuss how LBP can be used to lower-bound the BNN output over \(T\) and \(R\) as an alternative to IBP. In LBP, instead of propagating bounding boxes, one finds lower and upper Linear Bounding Functions (LBFs) for each layer and then propagates them through the network. As the bounding function has an extra degree of freedom w.r.t. the bounding boxes obtained through IBP, LBP usually yields tighter bounds, though at an increased computational cost. Since in deterministic networks non-linearity comes only from the activation functions, LBFs in the deterministic case are computed by bounding the activation functions, and propagating the bounds through the affine function that defines each layer.

Similarly, in our setting, given \(T\) in the input space and \(R\) for the first layer in the weight space, we start with the observation that LBFs can be obtained and propagated through commonly employed activation functions for Equation (17), as discussed in [41].

**Lemma 2**.: _Let \(f^{w}(x)\) be defined by Equations (15)-(17). For each hidden layer \(k=1,\ldots,K\), consider a bounding box in the pre-activation function, i.e. such that \(\zeta_{i}^{(k)}\in[\zeta_{i}^{(k),L},\zeta_{i}^{(k),U}]\) for \(i=1,\ldots,n_{k}\). Then there exist coefficients \(\alpha_{i}^{(k),L}\), \(\beta_{i}^{(k),L}\), \(\alpha_{i}^{(k),U}\) and \(\beta_{i}^{(k),U}\) of lower and upper LBFs on the activation function such that for all \(\zeta_{i}^{(k)}\in[\zeta_{i}^{(k),L},\zeta_{i}^{(k),U}]\) it holds that:_

\[\alpha_{i}^{(k),L}\zeta_{i}^{(k)}+\beta_{i}^{(k),L}\leq\sigma(\zeta_{i}^{(k)}) \leq\alpha_{i}^{(k),U}\zeta_{i}^{(k)}+\beta_{i}^{(k),U}.\]

The lower and upper LBFs can thus be minimised and maximised to propagate the bounds of \(\zeta^{(k)}\) in order to compute a bounding interval \([z^{(k),L},z^{(k),U}]\) for \(z^{(k)}=\sigma(\zeta^{(k)})\). Then, LBFs for the monomials of the bi-linear form of Equation (16) can be derived using McCormick's inequalities [54]:

\[W^{(k)}_{ij}z^{(k)}_{j}\geq W^{(k),L}_{ij}z^{(k)}_{j}+W^{(k)}_{ ij}z^{(k),L}_{j}-W^{(k),L}_{ij}z^{(k),L}_{j} \tag{18}\] \[W^{(k)}_{ij}z^{(k)}_{j}\leq W^{(k),U}_{ij}z^{(k)}_{j}+W^{(k)}_{ ij}z^{(k),L}_{j}-W^{(k),U}_{ij}z^{(k),L}_{j} \tag{19}\]

for every \(i=1,\ldots,n_{k}\), \(j=1,\ldots,n_{k-1}\) and \(k=1,\ldots,K\). The bounds of Equations (18)-(19) can thus be used in Equation (16) to obtain LBFs on the pre-activation function of the following layer, i.e. \(\zeta^{(k+1)}\). The final linear bound can be obtained by iterating the application of Lemma 2 and Equations (18)-(19) through every layer. This is summarised in the following proposition, which is proved in Appendix C along with an explicit construction of the LBFs.

**Proposition 4**.: _Let \(f^{w}(x)\) be the network defined by the set of Equations (15)-(17). Then for every \(k=0,\ldots,K\) there exists lower and upper LBFs on the pre-activation function of the form:_

\[\zeta_{i}^{(k+1)}\geq\mu_{i}^{(k+1),L}\cdot x+\sum_{l=0}^{k-1}(\nu_ {i}^{(l,k+1),L},W^{(l)})+\] \[\nu_{i}^{(k,k+1),L}\cdot W_{i}^{(k)}+\lambda_{i}^{(k+1),L}\quad \text{for }i=1,\ldots,n_{k+1}\] \[\zeta_{i}^{(k+1)}\leq\mu_{i}^{(k+1),U}\cdot x+\sum_{l=0}^{k-1}(\nu_ {i}^{(l,k+1),U},W^{(l)})+\] \[\nu_{i}^{(k-1,k+1),U}\cdot W_{i}^{(k)}+\lambda_{i}^{(k+1),U} \quad\text{for }i=1,\ldots,n_{k+1}\]

_where \(\langle\cdot,\cdot\rangle\) is the Frobenius product between matrices, \(\cdot\) is the dot product between vectors, and the explicit formulas for the LBF coefficients, i.e., \(\mu_{i}^{(k+1),L}\), \(\nu_{i}^{(l,k+1),L}\), \(\lambda_{i}^{(k+1),L}\), \(\mu_{i}^{(k+1),U}\), \(\nu_{i}^{(l,k+1),U}\), are given in Appendix C-D._
_Now let \(\zeta_{i}^{(k),L}\) and \(\zeta_{i}^{(k),U}\), respectively, be the minimum and the maximum of the right-hand side of the two equations above; then we have that \(\forall x\in T\) and \(\forall w\in R\):_

\[f^{w}(x)=\zeta^{(K+1)}\in\left[\zeta^{(K+1),L},\zeta^{(K+1),U}\right].\]

```
1:\(\#\)\(\mathcal{H}\) is a set of known safe weight intervals
2:\(\mathcal{H}\leftarrow\varnothing\)
3:\(\#\) Elementwise product to obtain width of weight margin
4:\(v\leftarrow\gamma\cdot I\cdot\Sigma\)
5:for\(i\gets 0\) to \(N\)do
6:\(w^{(i)}\gets p(w[\mathcal{D}])\)
7:\(\#\) Assume weight intervals are built to be disjoint
8:\([w^{(i),L},w^{(i),U}]\leftarrow[w_{i}-v,w_{i}+v]\)
9:\(f\) Interval/Linear Bound Propagation, Section VI-B
10:\(y^{i},y^{U}\leftarrow\) Propagate(\(f,T,[w^{(i),L},w^{(i),U}]\))
11:if\([y^{L},y^{U}]\subset S\)then
12:\(\mathcal{H}\leftarrow\mathcal{H}\bigcup\{[w^{(i),L},w^{(i),U}]\}\)
13:endif
14:endfor
15:\(P_{\text{safe}}^{L}\gets 0.0\)
16:for\([w^{(i),L},w^{(i),U}]\in\mathcal{H}\)do
17:\(\#\) Compute safe weight probs, Section VI-A
18:\(P_{\text{safe}}^{L}=P_{\text{safe}}^{L}+P([w^{(i),L},w^{(i),U}])\)
19:endfor
20:return\(P_{\text{safe}}^{L}\)
```

**Algorithm 1** Lower Bounds for BNN Probabilistic Robustness

## VII Complete Bounding Algorithm

Using the results presented in Section VI, it is possible to explicitly compute the bounds on probabilistic and decision robustness derived in Section V. In this section, we bring together all the results discussed so far, and assemble complete algorithms for the computation of bounds on \(P_{\text{safe}}(T,S)\) and \(D_{\text{safe}}(T,S)\). We discuss the procedure to lower bound \(P_{\text{safe}}(T,S)\) in Algorithm 1. We then discuss the details of upper bounds and bounds on \(D_{\text{safe}}(T,S)\), leaving the algorithms and their description for these bounds to Appendix D.

### _Lower Bounding Algorithm_

We provide a step-by-step outline for how to compute lower bounds on \(P_{\text{safe}}(T,S)\) in Algorithm 1. We start on line 2 by initializing the family of safe weight sets \(\mathcal{H}\) to be the empty set and by scaling the weight margin with the posterior weight scale (line 4). We then iteratively (line 5) proceed by sampling weights from the posterior distribution (line 6), building candidate weight boxes (line 8), and propagate the input and weight box through the BNN (line 10). We next check whether the propagated output set is inside the safe output region \(S\), and if so update the family of weights \(\mathcal{H}\) to include the weight box currently under consideration (lines 11 and 12). Finally, we rely on the results in Section VI-A to compute the overall probabilities over all the weight sets in \(\mathcal{H}\), yielding a valid lower bound for \(P_{\text{safe}}(T,S)\). For clarity of presentation, we assume that all the weight boxes that we sample in lines 6-8 are pairwise disjoint, as this simplifies the probability computation. The general case with overlapping weight boxes relies on the Bonferroni bound and is given in Appendix D-D. The algorithm for the computation of the lower bound on \(D_{\text{safe}}(T,S)\) (listed in the Appendix D as Algorithm 2) proceeds in an analogous way, but without the need to perform the check in line 11, and by adjusting line 18 to the formula from Proposition 2.

### _Upper Bounding Algorithm_

Upper bounding \(P_{\text{safe}}(T,S)\) and \(D_{\text{safe}}(T,S)\) follows the same computational flow as Algorithm 1. The pseudocode outlines computation of probabilistic and decision robustness are listed respectively in Algorithm 3 and 4 in Appendix D subsection A and Appendix D subsection B. We again proceed by sampling a rectangle around weights, propagate bounds through the NN, and compute the probabilities of weight intervals. The key change to the algorithm to allow upper bound computation involves computing the best case, rather than the worst case, for \(y\) in for decision robustness (line 12 in Algorithm 3) and ensuring that the entire interval \([y^{L},y^{U}]\notin S\) (line 18) for probabilistic robustness. In Appendix D subsection B we also discuss how adversarial attacks can be leveraged to improve the upper bounds.

### _Computational Complexity_

Calculations for probabilistic robustness and decision robustness follow the same computational flow and include: bounding of the neural network output, sampling from the posterior distribution, and computation of integrals over boxes on the input and weight space.

Regarding Algorithm 1 (or equivalently Algorithm 2 for decision robustness), it is clear that the computational complexity scales linearly with the number of samples, \(N\), taken from the posterior distribution. Observe that, in order to obtain a tight bound on the integral computation, \(N\) needs to be large enough such that \(N\) samples of the posterior with width \(\gamma\) span an area of high probability mass for \(p(w|\mathcal{D})\). Unfortunately, this means that, for a given approximation error magnitude, \(N\) needs to scale quadratically on the number of hidden neurons. Given the sampling of the hyper-rectangles, computation of the integral over the weight boxes is done through Equations (13) and (14). The integration over the weight boxes is done in constant time for HMC (though a good quality HMC posterior approximation scales with the number of parameters) and \(\mathcal{O}(n_{w})\) for VI. The final step needed for the methodology is that of bound propagation, which clearly differs when using IBP or LBP. In particular, the cost of performing IBP is \(\mathcal{O}(K\hat{n}\hat{m})\), where \(K\) is the number of hidden layers and \(\hat{n}\times\hat{m}\) is the size of the largest weight matrix \(W^{(k)}\), for \(k=0,\ldots,K\). LBP is, on the other hand, \(\mathcal{O}(K^{2}\hat{n}\hat{m})\). Overall, the time complexity for certifying a VI BNN with IBP is therefore \(\mathcal{O}(Nn_{w}K\hat{n}\hat{m})\), and similar formulas can be obtained for alternative combinations of inference and propagation techniques that are employed. We remark that, while sound, the bounds we compute are not guaranteed to converge to the true values of \(P_{\text{safe}}(T,S)\) and
\(D_{\text{safe}}(T,S)\) in the limit of the number of sample \(N\) because of the introduction of over-approximation errors due to bound propagation.

## VIII Experiments

In this section, we empirically investigate the suitability of our certification framework for the analysis of probabilistic and decision robustness in BNNs. We focus our analysis on four different case studies. First, we provide a comprehensive evaluation of an airborne collision avoidance system [24]. To do so, we partition the entire input space into 1.8 million different specifications and bound \(P_{\text{safe}}(T,S)\) and \(D_{\text{safe}}(T,S)\) by computing the bounds for each specification. We then turn our attention to an industrial regression benchmark [25] and demonstrate how our analysis can provide tight characterization of the worst-case error of predictions in adversarial settings in relation to the magnitude of the maximum attack allowed. Next, we analyse the scalability of our method in the well-known MNIST dataset for handwritten digits recognition [55], along with its behaviour on out-of-distribution input samples. Finally, we we study a two-class subset of the German Traffic Sign Recognition Benchmark (GTSRB) dataset [26], whose input space is 1500 dimensions larger than what has previously been studied for BNN certification against adversarial examples, showcasing that we are still able to compute non-trivial guarantees in this setting. For each dataset, we first describe the problem setting and BNN used to solve it, along with its hyperparameters. We then discuss the properties of interest for each dataset. Finally, we provide discussion and illustration of our bounds performance. All the experiments have been run on 4 NVIDIA 2080Ti GPUs in conjunction with 4 24-core Intel Core Xeon 6230.

### _Airborne Collision Avoidance_

Our first case study is the Horizontal airborne Collision Avoidance System (HCAS) [24], a dataset composed of millions of labelled examples of intruder scenarios.

#### Vi-A1 Problem Setting

The task of the BNN is to predict a turn advisory for an aircraft given another oncoming aircraft, including clear of conflict (COC), weak left (WL), weak right (WR), strong left (SL), and strong right (SR). These are depicted in the top left of Figure 3. We follow the learning procedure described in [24], where encounter scenarios are partitioned into 40 distinct datasets. We then learn a BNN to predict the correct advisories for each dataset, resulting in 40 different BNNs which need to be analysed.

To analyze the system of 40 BNNs, we first discretize the entire state-space into 1.8 million mutually exclusive input specifications. The input specifications are sized according to the spacing of the ground truth labels supplied by [24]. Namely, we consider an \(\ell_{\infty}\) norm ball with different widths for each input dimension. Those widths are \([0.016,0.025,0.025,0.05]\). The output specification is taken to be the set of all softmax vectors such that the argmax of the softmax corresponds to the true label. We separate these output specifications into 5 different properties, which we termed \(\phi_{j}\) for \(j=0,\ldots,4\) corresponding to each of the possible advisories. For all properties in this section we use LBP with 5 samples with a weight margin of 2.5 standard deviations.

We train BNNs with Variational Online Gauss Newton (VOGN), where the posterior approximation is a diagonal covariance Gaussian, and wih Hamiltonian Monte Carlo (HMC). The BNN architecture has a single hidden layer with 125 hidden units, the same size as the original system proposed in [24]. We use a diagonal covariance Gaussian prior with variance 0.5 for VOGN and a prior variance of 2.5 for HMC.

#### Vi-A2 Analysis with \(P_{\text{safe}}(T,S)\) Certification

For each of the 1.8 million disjoint input specification, we compute both upper and lower bounds on \(P_{\text{safe}}(T,S)\). Given that probabilistic robustness is a real-valued probability and not a binary predicate, practitioners must select thresholds that reflect a strong belief that a value is safe or unsafe. We call these thresholds \(\tau_{\text{safe}}\)

Fig. 4: **Left:** Box plots showing the distribution of upper and lower bounds for VI (top) and HMC (bottom). **Right:** Histograms showing gap between upper and lower bounds for VI (top) and HMC (bottom).

Fig. 3: **Top Left:** Encounter geometry, ground truth and property labels: Clear of Conflict (COC), Strong Left/Right (SL/R), Weak Left/Right (WL/R), for a collision scenario. Diagrams modified from [24]. **Bottom Left:** Encounter geometry labelled with features used for collision avoidance prediction. **Right:** Bounds on decision robustness obtained for HMC and VI trained BNNs for each property.
and \(\tau_{\text{unsafe}}\). Once one has computed bounds on \(P_{\text{safe}}(T,S)\), the proportions of safe and unsafe states (as reported in Table I) can be computed by checking thresholds. We check our bounds against strict safety and unsafety thresholds \(\tau_{\text{safe}}=0.98,\tau_{\text{unsafe}}=0.05\).

For the selected threshold values, Table I reports the certified performance of the BNN system. Such a report can be used by regulators and practitioners to determine the if the system is safe for deployment. In this case, we find that across all properties 74% of the states are certified to be safe while 20% are certified to be unsafe. The remaining 6% fall somewhere in between the two safety thresholds. These statistics indicate that roughly 18% of the decisions issued by the system were correct but not robust, thus the systems accuracy of 92% does not paint the complete story of its performance. Moreover, we break down each of the properties of the system, represented by each row of Table I, to understand where the most common failure modes occur. We find that the most unsafe indicators are the strong left, \(\phi_{3}\), and for strong right, \(\phi_{4}\), the system has the lowest certified safety at 53.4% and 55.0% respectively. They also have the highest certified unsafety at 34.0% and 37.0%, respectively. We conjecture that these these specifications are less safe due to the fact that their is less labeled data representing them in the dataset. Less data has been shown to be correlated with less robustness for BNNs [10]. If the results in Table I are deemed to be insufficient for deployment, then practitioners can collect more data for unsafe properties e.g., \(\phi_{3}\) and \(\phi_{4}\), or could resort to certified safe training for BNNs as suggested in [27].

#### V-A3 Analysis with \(D_{\text{safe}}(T,S)\) Certification

In order to analyze the decision robustness of the BNNs, we again discretize the input space. For these results, we use a coarser discretization, with the input specification being an \(\ell_{\infty}\) ball radius of \(0.125\) over each input dimension, and for the sake of computational efficiency we allow some gaps between the input specifications. For each specification, we compute upper and lower bounds on \(D_{\text{safe}}(T,S)\). We plot the result of our bounds on decision robustness in Figures 3 and 4. In the right hand portion of Figure 3 we visualize the average lower bound on decision robustness for two BNNs, one trained with HMC (yellow) and the other trained with VI (green). We find that we are able to certify a higher lower bound, indicating heightened robustness, for the HMC-trained BNN. This corroborates previous robustness studies that highlight that HMC is more adversarially robust [10, 11]. In Figure 4 we analyze the tightness of our bounds in this scenario by comparing the lower and upper bounds provided by our method. For VI, the gap, plotted in purple in the upper right, is tightly centered around a mean of \(0.08\), with maximum gap observed in these experiments being approximately \(0.16\) and a minimum \(0.035\). For HMC, on the other hand, the mean gap is \(0.11\), which is higher than VI, but this mean is affected by a small proportion of inputs that have a very high gap between upper and lower bounds, with the highest gap being \(0.72\). We further highlight the higher variance bound distribution for the HMC-trained BNN (plotted as blue and red box plots). We hypothesize that this arises due to the higher uncertainty predictive of the HMC in areas of little data [7]

#### V-A4 Computational Requirements

For VI certification, we can compute upper and lower bounds in an average of 0.544 seconds. Thus, when run in serial mode, the 3.6 million probabilistic bound computations needed for Table I takes an estimated 11.347 computational days. However, our parallelized certification procedure produces Table I in under 3 days of computational time (61 hours). These computations include the 1.8 million lower bound runs and 1.8 million upper bound runs. For HMC, on the other hand, certification can be done in a fraction of the time, with bounds being computed in 0.073 seconds. This is due to the fact that weight intervals for HMC necessarily satisfy the pairwise disjoint precondition of Proposition 1, thus no Bonferroni correction is needed.

### _Industrial Benchmarks_

We now focus our analyses on two safety-critical industrial regression problems taken from the UCI database [25], and widely employed for benchmarking of Bayesian inference methods [53, 56, 57].

#### V-B1 Problem Setting

The _Concrete_ dataset involves predicting the compressive strength of concrete, based on 8 key factors including its ingredients and age. The _Powerplant_ dataset uses six years worth of observations from combined cycle power plants and poses the problem of predicting energy output from a plant given a range of environmental and plant parameters. For each dataset we learn a BNN by using the architecture (i.e., a single hidden layer with 100 hidden units) and inference settings proposed in [53]. The BNNs are inferred using VOGN with a diagonal covariance prior over the weights with variance 0.5 for the _Concrete_ dataset and 0.25 for the _Powerplant_ dataset. We use a Gaussian likelihood corresponding to a mean squarred error loss function. In this setting we use IBP with 10 samples and a weight margin of 2 standard deviations.
#### Vi-B2 Analysis with \(P_{\text{safe}}(T,s)\) Certification

In industrial applications it is useful to understand the maximum amount of adversarial noise that a learned system can tolerate, as failures can be costly and unsafe [58]. To this end, we introduce the maximum and minimum robust radius. Given a threshold \(\tau_{\text{safe}}\) (as before) the maximum robust radius (MaxRR) is the largest \(\ell_{\infty}\) radius for which we can certify the BNN satisfies \(P_{\text{safe}}(T,S)>\tau_{\text{safe}}\). Similarly, the minimum unrobust radius (MinUR) is the smallest radius such that we can certify \(P_{\text{safe}}(T,S)<\tau_{\text{ unsafe}}\). The MaxRR gives us a safe lower bound on the amount of adversarial noise a BNN is robust against, whereas the MinUR gives us a corresponding upper bound.

In our experiments on these datasets we considered \(\tau_{\text{safe}}=\tau_{\text{ unsafe}}=0.7\), meaning that we request that over \(70\%\) of the BNN probability mass is certifiably safe; however, we stress that similar results can be obtained for different values of \(\tau_{\text{safe}}\) similarly to what is discussed in our previous analysis of the HCAS dataset. In order to compute the MaxRR we start with \(\epsilon=0\), check that \(P_{\text{safe}}(T,S)>\tau_{\text{safe}}\) using our lower bound, and if the inequality is satisfied we increase epsilon by 0.01 and continue this process until the inequality to longer holds. Similarly for the MinUR, we start with \(\epsilon=0.5\) and iteratively decrease the value of \(\epsilon\) until the upper bound no longer certifies that \(P_{\text{safe}}(T,S)<\tau_{\text{ unsafe}}\); if the property does not hold at \(0.5\) one can increase the value of \(\epsilon\) until the bound holds.

The result of computing the MaxRR and MinUR over the test datasets for the _Concrete_ and _Powerplant_ datasets are plotted in Figure 5. We highlight that in the overwhelming majority of the cases our methods is able to return non-vacuous bounds on MinUR (i.e., strictly less than \(1\)) and on MaxRR (i.e., strictly greater than \(0\)). As expected we observe the MaxRR is strictly smaller than MinUR. Encouragingly, as MinUR grows, MaxRR tends to increase indicating that our bounds track the true value of \(P_{\text{safe}}(T,S)\). We see that the _Concrete_ dataset is typically guaranteed to be robust for radius \(\epsilon\approx 0.03\) and is typically guaranteed to be unsafe for \(\epsilon\approx 0.06\). For the _Powerplant_ posterior we compute a MaxRR of roughly \(0.18\) for most inputs and a MinUR lower than \(0.32\). Notice how the results for the _Concrete_ datasets systematically display more robustness than those for _Powerplant_ and the gap between MaxRR and MinUR is significantly smaller in the former datasets than in the latter.

#### Vi-B3 Computational Requirements

On average, it takes 1.484 seconds to compute a certified upper or lower bound on \(P_{\text{safe}}(T,S)\) for the _Powerplant_ dataset and 1.718 seconds for the _Concrete_ dataset. We use a linear search in order to compute the MaxRR and MinUR which require, on average, 5 certifications for both MaxRR and MinUR computations. We compute these values over the entire test datasets for both _Powerplant_ and _Concrete_, which requires tens of thousands of certifications and each input can be done in parallel.

### _Mnist_

We investigate the suitability of our methods in providing certifications for BNNs on larger input domains, specifically BNNs learnt for MNIST, a standard benchmark for the verification of deterministic neural networks whose inputs are 784-dimensional. In this setting we use IBP with 5 weight samples with a weight margin of 2.5 standard deviations.

#### Vi-B1 Problem Setting

MNIST poses the problem of handwritten digit recognition. Given handwritten digits encoded as a 28 by 28 black and white image, the task is to predict which digit - 0 through 9 - is depicted in the image (two images randomly sampled from the dataset are reproduced in the far left of Figure 6). We learn BNNs using the standard 50,000/10,000 train/test split that is provided in the original work [55]. For our experimental analysis, we use one-layer neural network with 128 hidden neurons, each of which uses rectified linear unit activation functions. The BNN has 10 output neurons that use a softmax activation function. We train the network using VOGN with a diagonal covariance Gaussian prior that has variance 2.0. We use a sparse categorical cross-entropy loss modified with the method presented in [27] to promote robustness in the BNN posteriors.

#### Vi-B2 Analysis using \(D_{\text{safe}}(T,S)\) Certification

We analyze the trained BNN using decision robustness on 1000 images taken from the MNIST test dataset. We compute bounds on \(D_{\text{safe}}(T,S)\) for increasing widths of an \(\ell_{\infty}\) input region \(\epsilon\). We plot the mean and standard deviation obtained for the upper (\(D_{\text{safe}}^{U}\), in red) and lower bound (\(D_{\text{safe}}^{L}\), in blue) on decision robustness for the ground truth label of each image in the left hand portion of Figure 6. As greater \(\epsilon\) implies a larger input

Fig. 5: Computation of the minimum (MinUR) and maximum (MaxRR) safe radius for _Concrete_ and _Powerplant_ datasets. **Left:** Boxplots for the empirical distribution of MinUR and MaxRR computed over all test inputs. **Centre:** Per-test-instance certified radii for the _Concrete_ dataset. **Right:** Per-test-instance certified radii for the _Powerplant_ dataset.
specification \(T\), increasing values of \(\epsilon\) leads to a widening of the gap between the lower and upper bounds, and hence an increased vulnerability of the network. Notice that even for \(\epsilon=0.25\), i.e., half of the whole input space, our method still obtains on average non-vacuous bounds (i.e., strictly within \((0,1)\)). In order to get a rough estimation of the adversarial robustness of the network, we observe that, for lower bound values above \(0.5\), the BNN is guaranteed to correctly classify all the inputs in the region \(T\) (however, as MNIST has 10 classes, even values of the lower bound lower than \(0.5\) could still result in correct classification). Using the \(0.5\) threshold, we notice that our method guarantees that the BNN is still robust on average for \(\epsilon=0.075\). Notice that this is on par with results obtained for verification of deterministic neural networks, where \(\epsilon=0.05\) leads to adversarial attack robustness of around \(70\%\)[17].

#### Vii-C3 Certification of Uncertainty Behaviour

In this section we study how to certify the uncertainty behavior of a BNN in the presence of adversarial noise. We assume we have an out-of-distribution input, i.e., an input whose ground-truth does not belong to any of the classes in the range of the learned model. As with previous specifications, we build the set \(T\) around such an input with an \(\ell_{\infty}\) ball of radius \(\epsilon\). Unlike for the previous specifications, we build \(S\) as the set of all softmax vectors such that no entry in the vector is larger than a specified value \(\tau_{\text{uncertain}}\). The function of \(\tau_{\text{uncertain}}\) is to determine the confidence at which a classification is ruled to be uncertain. For example, in Figure 6 we have set \(\tau_{\text{uncertain}}=0.4\), thus any classification that is made with confidence \(<0.4\) will be ruled uncertain. By certifying that all values of \(T\) are mapped into \(S\), we guarantee that the BNN is uncertain on all points around the out-of-distribution input.

In the right half of Figure 6, we plot two example images from the FashionMNIST dataset, which are considered out-of-distribution for the BNN trained on MNIST. In our experiments we use 1000 test set images from the FashionMNIST dataset. On the right of the out-of-distribution samples in Figure 6, we plot the bounds on decision robustness with various values of \(\epsilon\) for the \(\ell_{\infty}\) ball. We start by noticing that the BNN never outputs a confidence of more than \(\sim 0.25\) on the clean Fashion-MNIST dataset, which indicates that the network has good calibrated uncertainty on these samples. We notice that up to \(\epsilon=0.06\) we certify that no adversary can perturb the image to force a confident classification; however, at \(\epsilon=0.10\) no guarantees can be made.

#### Vii-C4 Architecture Width and Depth

We now analyse the behaviour of our method when computing bounds on the certified radius on MNIST while varying the width and depth of the BNN architecture. The results of this analysis are given in Figure 7. Notice that we are able to obtain non-vacuous bounds in all the cases analysed. However, as could be expected, we see that the gap between MinUR and MaxRR widens as we increase the depth and/or the width of the neural network. This inevitably arises from the fact that the tightness of bound propagation techniques decreases as we need to perform more boundings and/or propagations, and because increasing the number of weights in the network renders the bounds obtained by Proposition 2 more coarse, particularly as we increase the number of layers of the BNN, as explained in Section VII-C. In particular, we observe that MinRR increases drastically as we increase the number of layers in the BNN architecture, while, empirically, the bounding for MaxRR is more stable w.r.t. the architecture parameters.

#### Vii-C5 Computational Requirements

On average, it takes 24.765 seconds to verify an MNIST image on a single CPU core. Each of the images in our experiments is run in parallel across 96 cores which allows us to compute all of the results for Figure 6 in less than an hour.

### _German Traffic Sign Recognition_

In this section, we investigate the ability of our method to scale to a full-color image dataset, which represent safety

Fig. 6: **Left:** Mean and standard deviation of upper and lower bounds obtained on \(D_{\text{safe}}\) on 1000 images taken from the MNIST test set. **Right:** Mean and standard deviation on upper and lower bound on decision robustness on out-of-distribution samples taken from the FashionMNIST dataset.

Fig. 7: **Left:** Boxplots for the empirical distribution of the maximum safe radius and minimum unsafe radius for a BNN with 128 hidden units and a single hidden layer. **Right:** For a range of architectures, we plot the mean certified safe and unsafe radius.
critical tasks with high-dimensional inputs (\(2,352\) dimensions).

#### V-C1 Problem Setting

We study BNNs on a two-class subset of the German Traffic Sign Recognition Benchmark (GTSRB), consisting of the images that represent the 'construction ahead' and '50 Km/H speed limit' [26]. Though this dataset is only comprised of two classes, full-colour images stretch the capabilities of BNN training methods, especially robust Bayesian inference. The dataset is comprised of 5000 training images and 1000 test set images. We employ VOGN to train a Bayesian convolutional architecture, with 2 convolutional layers and one fully-connected layer first proposed in [18]. We employ the method of [27] in order to encourage robustness in the posterior. We find that this dataset poses a challenge to robust inference methods, with the BNN achieving 72% accuracy over the test set after 200 epochs. We found that, without robust training, we are able to achieve 98% accuracy over the test set, but were unable to certify robustness or uncertainty for any tested image (see discussion of limitations below). We verify these networks with 3 weight samples with a weight margin of 3.0 standard deviations.

#### V-C2 Analysis with \(D_{\text{safe}}(T,S)\)

For our analysis for GTSRB, we take \(T\) to be a \(\ell_{\infty}\) ball with radius \(2/255\). As in our previous analysis, for test set images we take \(S\) to be the set of all vectors such that the true class is the argmax. We study 250 images and find that 53.8% of the images are certified to be correct. We plot a visual sample of these images in the top row of Figure 8. We also study the out-of-distribution performance of various kinds of images with \(\tau_{\text{uncertain}}=0.55\). Of 400 images of random noise, visualized in the bottom row of Figure 8, we certified that the BNN was uncertain on 398 images, indicating that on that set the BNN has correctly calibrated uncertainty as it does not issue conflict predictions on random noise. We then turned our attention to two more realistic sets of out-of-distribution images: nonsense traffic signs and international traffic signs. We were limited to a small set of free-use images for these tests but found that for eight out of ten nonsense traffic signs we were able to certify the BNN's uncertainty, and for nine out of ten international traffic signs we were able to certify the uncertainty. On average these certifications took 34.2 seconds.

#### V-C3 Limitations

While this analysis represents an encouraging proof of concept for certification of BNNs, we find that datasets whose inputs are of this scale and complexity are not yet fully accessible to robust inference for BNNs, as 74% test set accuracy is not strong enough performance to warrant deployment. However, with approaches such as [59, 60] investigating more powerful methods for scaling Bayesian inference for neural networks, we are optimistic that future works will be able to apply our method to more advanced Bayesian approximate posteriors.

## IX Conclusion

In this work, we introduced a computational framework for evaluating robustness properties of BNNs operating under adversarial settings. In particular, we have discussed how probabilistic robustness and decision robustness - both employed in the adversarial robustness literature for Bayesian machine learning [42, 43] - can be upper- and lower-bounded via a combination of posterior sampling, integral computation over boxes and bound propagation techniques. We have detailed how to compute these properties for the case of HMC and VI posterior approximation, and how to instantiate the bounds for interval and linear propagation techniques, although the framework presented is general and can be adapted to different inference techniques and to most of the verification techniques employed for deterministic neural networks.

In an experimental analysis comprising 5 datasets (airborne collision avoidance, concrete, powerplant, MNIST, and GTSRB), we have showcased the suitability of our approach for computing effective robustness bounds in practice, and for various additional measures that can be computed using our technique including certified robust radius and analysis of uncertainty.

With verification of deterministic neural networks already being NP-hard, inevitably certification of Bayesian neural networks poses several practical challenges. The main limitation of the approach presented here arises directly from the Bayesian nature of the model analysed, i.e., the need to bound and partition at the weight space level (which is not needed for deterministic neural networks, with the weight fixed to a specific value). Unfortunately, this means that the computational complexity, and also the tightness of the bounds provided, scale quadratically with the number of neurons across successive layer connections. We have discussed methods for mitigating the resulting gap between the bounds, including adaptive partitioning based on weight variance and implementing a branch-and-bound refinement approach for the bound, which would, however, results in a sharp increase in computational time. Nevertheless, the methods presented here provide the first formal technique for the verification of robustness in Bayesian neural network systematically and across various robustness notions, and as such can provide a sound basis for for future practical applications in safety-critical scenarios.

Fig. 8: Certification of a Bayesian CNN on a two-class subset of the German Traffic Recognition (GTSRB) dataset. In the top row, we plot illustrative examples showing that we can verify correctness of test inputs. In the bottom three rows we visualize the the uncertainty guarantees on various out-of-distribution inputs including nonsense traffic signs (second row), international traffic signs (third row), and random noise (bottom row).