BNN-DP: Robustness Certification of Bayesian Neural Networks via Dynamic
  Programming
In this paper, we introduce BNN-DP, an efficient algorithmic framework for
analysis of adversarial robustness of Bayesian Neural Networks (BNNs). Given a
compact set of input points $T\subset \mathbb{R}^n$, BNN-DP computes lower and
upper bounds on the BNN's predictions for all the points in $T$. The framework
is based on an interpretation of BNNs as stochastic dynamical systems, which
enables the use of Dynamic Programming (DP) algorithms to bound the prediction
range along the layers of the network. Specifically, the method uses bound
propagation techniques and convex relaxations to derive a backward recursion
procedure to over-approximate the prediction range of the BNN with piecewise
affine functions. The algorithm is general and can handle both regression and
classification tasks. On a set of experiments on various regression and
classification tasks and BNN architectures, we show that BNN-DP outperforms
state-of-the-art methods by up to four orders of magnitude in both tightness of
the bounds and computational efficiency.
Steven Adams, Andrea Patane, Morteza Lahijanian, Luca Laurenti
# BNN-DP: Robustness Certification of Bayesian Neural Networks via Dynamic Programming

###### Abstract

In this paper, we introduce BNN-DP, an efficient algorithmic framework for analysis of adversarial robustness of Bayesian Neural Networks (BNNs). Given a compact set of input points \(T\subset\mathbb{R}^{n}\), BNN-DP computes lower and upper bounds on the BNN's predictions for all the points in \(T\). The framework is based on an interpretation of BNNs as stochastic dynamical systems, which enables the use of Dynamic Programming (DP) algorithms to bound the prediction range along the layers of the network. Specifically, the method uses bound propagation techniques and convex relaxations to derive a backward recursion procedure to over-approximate the prediction range of the BNN with piecewise affine functions. The algorithm is general and can handle both regression and classification tasks. On a set of experiments on various regression and classification tasks and BNN architectures, we show that BNN-DP outperforms state-of-the-art methods by up to four orders of magnitude in both tightness of the bounds and computational efficiency.

Machine Learning, BNN-DP, BNN-DP, BNN-DP: Robustness Cer
precision and computational time. For instance, on the Fashion MNIST dataset, our approach achieves an average 93% improvement in certified lower bound compared to Berrada et al. (2021), while being around 3 orders of magnitude faster. In summary, this paper makes the following main contributions:

* we introduce a framework based on stochastic dynamic programming and convex relaxation for the analysis of adversarial robustness of BNNs,
* we implement an efficient algorithmic procedure of our framework for BNNs trained with Gaussian variational inference (VI) in both regression and classification settings,1 and Footnote 1: Our code is available at [https://github.com/sjladams/BNN_DP](https://github.com/sjladams/BNN_DP).
* we benchmark the robustness of a variety of BNN models on five datasets, empirically demonstrating how our method outperforms state-of-the art approaches by orders of magnitude in both tightness and efficiency.

Related WorksMany algorithms have been developed for certification of deterministic (i.e., non Bayesian) neural networks (NNs) (Katz et al., 2017; Weng et al., 2018; Wong and Kolter, 2018; Bunel et al., 2020). However, these methods cannot be employed to BNNs because they all assume the weights of the network have a fixed value, whereas in the Bayesian setting they are distributed according to the BNN posterior. Methods for certification of BNNs have recently presented in (Wicker et al., 2020; Berrada et al., 2021; Lechner et al., 2021). Wicker et al. (2020) consider a different notion of robustness than the one in this paper, not directly related to adversarial attacks on the BNN decision. Furthermore, that work considers a partitioning procedure in weight space that makes it applicable only to small networks and/or with small variance. The method proposed in (Berrada et al., 2021) is based on dual optimization. Hence, it is restricted to distributions with bounded support and needs to solve non-convex problems at large computational costs for classification tasks. Separately, Lechner et al. (2021) aims to build an intrinsic safe BNN by truncating the posterior in the weight space. Cardelli et al. (2019); Wicker et al. (2021) introduced statistical approaches to quantify the robustness of a BNN, which however, does not return formal guarantees, which are necessary in safety-critical settings. Empirical methods that use the uncertainty of BNNs to flag adversarial examples are introduced in (Rawat et al., 2017; Smith and Gal, 2018). These, however, consider only point-wise uncertainty estimates, specific to a particular test point and do not account for worst-case adversarial perturbations.

Various recent works have proposed formal methods to compute adversarial robustness for Gaussian Processes (GPs) (Cardelli et al., 2019; Smith and Gal, 2018; Patane et al., 2022; Smith et al., 2022). In BNNs, however, due to the non-linearity of activation functions, the distribution over the space of functions induced by a BNN is generally non-Gaussian, even if a Gaussian distribution in weight space is assumed. Hence, the techniques that are developed for GPs cannot be directly applied to BNNs.

## 2 Robust Certification of BNNs Problem

### Bayesian Neural Networks (BNNs)

For an input vector \(x\in\mathbb{R}^{n_{0}}\), we consider fully connected neural networks \(f^{w}:\mathbb{R}^{n_{0}}\rightarrow\mathbb{R}^{n_{K+1}}\) of the following form for \(k=0,\ldots,K\):2

Footnote 2: Note that the formulation of neural networks considered in Eqn (1) also includes convolutional neural networks (CNNs). In fact, the convolutional operation can be interpreted as a linear transformation into a larger space; see, e.g., Chapter 3.4.1 in (Gal and Ghahramani, 2016). This allows us to represent convolutional layers equivalently as fully connected layers, and do verification for CNNs as we show in Section 7.

\[\begin{split} z_{0}&=x,\qquad\qquad\zeta_{k+1}=W_{k} (z_{k}^{T},1)^{T},\\ z_{k}&=\phi_{k}(\zeta_{k}),\quad f^{w}\left(x\right) =\zeta_{K+1},\end{split} \tag{1}\]

where \(K\) is the number of hidden layers, \(n_{k}\) is the number of neurons of layer \(k\), \(\phi_{k}:\mathbb{R}^{n_{k}}\rightarrow\mathbb{R}^{n_{k}}\) is a vector of continuous activation functions (one for each neuron) in layer \(k\), and \(W_{k}\in\mathbb{R}^{n_{k}\times n_{k+1}}\) is the matrix of weights and biases that correspond to the \(k\)th layer of the network. We denote the vector of parameters by \(w=(W_{1}^{T},\ldots,W_{K}^{T})^{T}\) and the mapping from \(\zeta_{k_{1}}\) to \(\zeta_{k_{2}}\) by \(f^{w}_{k_{1}:k_{2}}:\mathbb{R}^{n_{k_{1}}}\rightarrow\mathbb{R}^{n_{k_{2}}}\) for \(k_{1},k_{2}\in\{0,...,K\}\). \(\zeta_{K+1}\) is the final output of the network (or the logit in the case of classification problems).

In the Bayesian setting, one starts by assuming a prior distribution \(p(w)\) over the parameters \(w\) and a likelihood function \(p(y|x,w)\). We adopt bold notation to denote random variables and write \(f^{w}\) to denote a BNN defined according to Eqns. (1). The likelihood is generally assumed to be Gaussian in case of regression and categorical for classification, where the probability for each class is given as the softmax of the neural network final logits (MacKay, 1992). Then, given a training dataset \(\mathcal{D}=\{(x_{i},y_{i})\}_{i=1}^{N_{\mathcal{D}}}\), learning amounts to computing the posterior distribution \(p(w|\mathcal{D})\) via the Bayes rule (MacKay, 1992). The posterior predictive distribution over an input \(x^{*}\) is finally obtained by marginalising the posterior over the likelihood, i.e., \(p(y^{*}|x^{*},\mathcal{D})=\int p(y^{*}|x^{*},w)p(w|\mathcal{D})dw\). The final output (decision) of the BNN, \(\hat{y}(x^{*})\), is then computed by minimising a loss function \(\mathcal{L}\) averaged over the predictive distribution, i.e.,

\[\hat{y}(x^{*})=\arg\min_{y}\int\mathcal{L}(y,y^{*})p(y^{*}|x^{*},\mathcal{D})dy ^{*}.\]
In this paper, we focus on both regression and classification problems. In regression, an \(l_{2}\) loss is generally used, which leads to an optimal decision \(\hat{y}(x^{*})\) given by the mean of the predictive posterior distribution (Neal, 2012), i.e., \(\hat{y}(x^{*})=\mathbb{E}_{y\sim p(y|x^{*},\mathcal{D})}\left[\mathbf{y}\right].\)3 For classification, \(\ell_{0-1}\) loss is typically employed, which results in

Footnote 3: In the remainder, we may omit the probability measure of an expectation or probability when it is clear from the context.

\[\hat{y}(x^{*})=\operatorname*{arg\,max}_{i\in\{1,...,n_{K+1}\}}\mathbb{E}_{ \boldsymbol{w}\sim p(w|\mathcal{D})}\left[\text{softmax}^{(i)}(f^{\boldsymbol {w}}\left(x^{*}\right))\right],\]

where \(\text{softmax}^{(i)}\) is the \(i\)th component of the \(n_{K+1}\)-dimensional softmax function.4 Unfortunately, because of the non-linearity introduced by the neural network architecture, the computation of the posterior distribution and consequently of \(\hat{y}(x^{*})\) cannot be done analytically. Therefore, approximate inference methods are required. In what follows, we focus on mean-field Gaussian Variational Inference (VI) approximations (Blundell et al., 2015). Specifically, we fit an approximating multivariate Gaussian \(q(w)=\mathcal{N}\left(w\mid\mu_{w};\,\Sigma_{w}\right)\approx p(w\mid\mathcal{ D})\) with mean \(\mu_{w}\) and block diagonal covariance matrix \(\Sigma_{w}\) such that for \(k\in\{0,\ldots,K\}\) and \(i\in\{1,\ldots,n_{k}\}\), the approximating distribution of the parameters corresponding to the \(i\)th node of the \(k\)th layer is

Footnote 4: Analogous formulas can be obtained for the weighted classification loss by factoring in misclassification weights in the argmax.

\[q(W_{k}^{(i,:)})=\mathcal{N}\left(W_{k}^{(i,:)}\mid\mu_{w,k,i};\,\Sigma_{w,k,i }\right) \tag{2}\]

with mean \(\mu_{w,k,i}\) and covariance matrix \(\Sigma_{w,k,i}\).5

Footnote 5: BNN-DP can be extended to Gaussian approximation distributions with inter-node or inter-layer correlations. In that case, to solve the backward iteration scheme of Theorem 1, the value functions need to be marginalized over partitions in weight space.

**Remark 1**.: _While our primary focus is on VI, the techniques presented in this paper can be applied to other approximate inference methods, such as HMC (Neal, 2012) and Dropout (Gal & Ghahramani, 2016). In these cases, the prediction of the BNN is obtained by averaging over a finite ensemble of NNs. For this setting, the dynamic programming problem in Theorem 1 reduces to computing piecewise linear relaxations for each layer of a weighted sum, i.e., an average, of deterministic neural networks, and propagating the resulting relaxations backward._

### Problem Statement

Given a BNN \(f^{\boldsymbol{w}}\) trained on a dataset \(\mathcal{D}\), as common in the literature (Madry et al., 2017), for a generic test point \(x^{*}\), we represent the possible adversarial perturbations by defining a compact neighbourhood \(T\) around \(x^{*}\) and measure the changes in the BNN output caused by limiting the perturbations to lie within \(T\).

**Definition 1** (Adversarial Robustness).: _Consider a BNN \(f^{\boldsymbol{w}}\), a compact set \(T\subset\mathbb{R}^{n_{0}}\), and input point \(x^{*}\in T\). For a given threshold \(\gamma>0\), \(f^{\boldsymbol{w}}\) is adversarially robust in \(x^{*}\) iff_

\[\forall x\in T,\quad\|\hat{y}(x)-\hat{y}(x^{*})\|_{p}\leq\gamma,\]

_where \(\|\cdot\|_{p}\) is an \(\ell_{p}\) norm._

Definition 1 is analogous to the standard notion of adversarial robustness employed for deterministic neural networks (Katz et al., 2017) and Bayesian models (Patane et al., 2022). As discussed in Section 2.1, the particular form of a BNN's output depends on the specific application considered. Below, we focus on regression and classification problems.

**Problem 1**.: _Let \(T\subset\mathbb{R}^{n_{0}}\) be a compact subset. Define functions \(I(y)=y\) and \(\text{softmax}(y)=\left[\text{softmax}^{(1)}(y),...,\text{softmax}^{(n_{K+1}) }(y)\right]\). Then, for a BNN \(f^{\boldsymbol{w}}\), \(h\in\{1,\text{softmax}\}\), and \(i\in\{1,...,n_{K+1}\}\), compute:_

\[\begin{split}\pi_{\min}^{(i)}(T)&=\min_{x\in T} \mathbb{E}_{\boldsymbol{w}\sim q(\cdot)}\left[h^{(i)}(f^{\boldsymbol{w}}\left( x\right))\right],\\ \pi_{\max}^{(i)}(T)&=\max_{x\in T}\mathbb{E}_{ \boldsymbol{w}\sim q(\cdot)}\left[h^{(i)}(f^{\boldsymbol{w}}\left(x\right)) \right].\end{split} \tag{3}\]

In the regression case (\(h=I\)), Problem 1 seeks to compute the ranges of the expectation of the BNN for all \(x\in T\). Similarly, in the classification case (\(h=\text{softmax}\)), Eqns. (3) define the ranges of the expectation of the softmax of each class for \(x\in T\). It is straightforward to see that these quantities are sufficient to check whether \(f^{\boldsymbol{w}}\) is adversarially robust for \(x\in T\); that is, if \(\sup_{x\in T}||\hat{y}(x)-\hat{y}(x^{*})||_{p}\leq\gamma\).

**Remark 2**.: _Our method can be extended to other losses, i.e., other forms of \(h\) in Eqns. (3), as long as affine relaxations of \(h\) can be computed._

Approach OutlineDue to the non-convex nature of \(f^{\boldsymbol{w}}\) and possibly \(h\), the computation of \(\mathbb{E}_{\boldsymbol{w}\sim q(\cdot)}\left[h(f^{\boldsymbol{w}}\left(x \right))\right]\) is analytically infeasible. To solve this problem, in Section 4, we view BNNs as stochastic dynamical systems evolving over the layers of the neural network. Through this, we show that adversarial robustness can be characterized as the solution of a dynamic programming (DP) problem. This allows us to break its computation into \(K\) simpler optimization problems, one for each layer. Each problem essentially queries a back-propagation of the uncertainty of the BNN through \(h\) and from one layer of the neural network to the next. Due to the non-convex nature of the layers of the BNN, these problems still cannot be solved exactly. We overcome this problem by using convex relaxations. Specifically, in Section 5, we show that efficient PWA relaxations can be obtained by recursively bounding the DP problem. In Section 6, we combine the theoretical results into a general algorithm called BNN-DP that solves Problem 1 efficiently.
## 3 Preliminaries on Relaxations of Functions

To propagate the uncertainty of the BNN from one layer to the other, we rely on upper and lower approximations of the corresponding Neural Network (NN), also known as _relaxations_. For vectors \(\hat{x},\hat{x}\in\mathbb{R}^{n}\), we denote by \([\hat{x},\hat{x}]\) the \(n\)-dimensional hyper-rectangle defined by \(\hat{x}\) and \(\hat{x}\), i.e., \([\hat{x},\hat{x}]=[\hat{x}^{(1)},\hat{x}^{(1)}]\times[\hat{x}^{(2)},\hat{x}^{(2 )}]\times\ldots\times[\hat{x}^{(n)},\hat{x}^{(n)}]\). We consider two types of relaxations, interval and affine.

**Definition 2** (Interval Relaxation).: _An interval relaxation of a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) over a set \(T\subseteq\mathbb{R}^{n}\) are two vectors \(\hat{b},\hat{b}\in\mathbb{R}^{m}\) such that \(f(x)\in[\hat{b},\hat{b}]\) for all \(x\in T\)._

**Definition 3** (Affine Relaxation).: _An affine relaxation of a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) over a set \(T\subseteq\mathbb{R}^{n}\) are two affine functions \(\hat{A}x+\hat{b}\) and \(\hat{A}x+\hat{b}\) with \(\hat{A},\hat{A}\in\mathbb{R}^{m\times n}\) and \(\hat{b},\hat{b}\in\mathbb{R}^{m}\) such that \(f(x)\in[\hat{A}x+\hat{b},\hat{A}x+\hat{b}]\) for all \(x\in T\)._

Interval and symbolic arithmetic can be used to propagate relaxations through the layers of a NN. Let, \([\alpha]_{+}\coloneqq\max\{\alpha,0\}\) and \([\alpha]_{-}\coloneqq\min\{\alpha,0\}\) represent the saturation operators on \(\alpha\). For a vector or matrix, \([\cdot]_{+}\) and \([\cdot]_{-}\) represent element-wise max and min, respectively. We adopt the notation of Liu et al. (2021) and write interval arithmetic w.r.t. a linear mapping \(M\) compactly as \(\otimes\) where \(M\otimes[\hat{b},\hat{b}]\coloneqq[[M]_{+}\hat{b}+[M]_{-}\hat{b},\,[M]_{+} \hat{b}+[M]_{-}\hat{b}]\), and use the similar notation for symbolic arithmetic.

## 4 BNN Certification via Dynamic Program

As observed in Marchi et al. (2021), NNs and consequently BNNs can be viewed as dynamical systems evolving over the layers of the network. In particular, for \(k\in\{0,...,K\}\), Eqn. (1) can be rewritten as:

\[z_{k+1}=\phi_{k+1}(\mathbf{W}_{k}(z_{k}^{T},1)^{T}) \tag{4}\]

with initial condition \(z_{0}=x\). Since, in a BNN, weights and biases are random variables sampled from the approximate posterior \(q(\cdot)\), Eqn. (4) describes a non-linear stochastic process evolving over the layers of the NN. This observation leads to the following theorem, which shows that \(\mathbb{E}_{\mathbf{w}\sim q(\cdot)}\left[h(\mathbf{W}_{k}(z^{T},1)^{T})\right]\) can be characterized as the solution to a backward recursion DP problem.

**Theorem 1**.: _Let \(f^{\mathbf{w}}\left(x\right)\) be a fully connected BNN with \(K\) hidden layers and \(h:\mathbb{R}^{n_{K+1}}\rightarrow\mathbb{R}^{l}\) be an integrable function. For \(k=0,...,K\), define functions \(V_{k}:\mathbb{R}^{n_{k}}\rightarrow\mathbb{R}^{l}\) backwards-recursively as:_

\[V_{K}(z)=\mathbb{E}_{\mathbf{W}_{K}\sim q(\cdot)}\left[h(\mathbf{W}_{K}( z^{T},1)^{T})\right], \tag{5a}\] \[V_{k-1}(z)=\mathbb{E}_{\mathbf{W}_{k-1}\sim q(\cdot)}\left[V_{k}( \phi_{k}(\mathbf{W}_{k-1}(z^{T},1)^{T}))\right]. \tag{5b}\]

_Then, it holds that \(\mathbb{E}_{\mathbf{w}\sim q(\cdot)}\left[h(f^{\mathbf{w}}\left(x\right))\right]=V_{0}(x)\)._

The proof of Theorem 1 is reported in Appendix A.1 and obtained by induction over the layers of the NN by relying on the law of total expectation and independence of the parameters distribution at different layers.6 Figure 1 illustrates the backward-iteration scheme of Theorem 1 for a two-hidden-layer BNN. Starting from the last layer, value functions \(V_{k}\) are constructed according to Eqns. (5a) and (5b) describing how the output of layer \(k\) is transformed in the previous layers. Theorem 1 is a central piece of our framework as it allows one to break the computation of \(\mathbb{E}_{\mathbf{w}\sim q(\cdot)}\left[h(f^{\mathbf{w}}\left(x\right))\right]\) into \(K+1\) (simpler) sub-problems, one for each layer of the BNN. In fact, note that \(V_{k}\) is a deterministic function. Hence, all the uncertainty in \(V_{k-1}\) depends only on the weights of layer \(k-1\). This is a key point that we use to derive efficient methods to solve Problem 1. Nevertheless, we stress that since \(V_{k}(z)\) is obtained by propagating \(z\) over \(K-k\) layers of the BNN, this is still generally a non-convex function, whose exact optimisation is infeasible in practice. Consequently, we employ the following corollary, which guarantees that, to solve Problem 1, it suffices to recursively bound \(V_{k}\) following Eqns. (5a) and (5b).

Footnote 6: While the vast majority of VI algorithms make this assumption, Theorem 1 can be generalized to the case where there is inter-layer correlation by marginalizing Eqn. (5a) and (5b) over partitions in correlated weight-spaces.

**Corollary 1**.: _For \(k\in\{1,\ldots,K\}\), let functions \(\hat{V}_{k},\hat{V}_{k}:\mathbb{R}^{n_{k}}\rightarrow\mathbb{R}^{n_{l}}\) be relaxations of \(V_{k}(z_{k})\), i.e, \(\forall z_{k}\in\mathbb{R}^{n_{k}},\hat{V}_{k}(z_{k})\leq V_{k}(z_{k})\leq \hat{V}_{k}(z_{k})\). Then_

\[\mathbb{E}_{\mathbf{W}_{k-1}\sim q(\cdot)}\left[\hat{V}_{k}(\phi_{k} (\mathbf{W}_{k-1}(z^{T},1)^{T}))\right]\leq V_{k-1}(z_{k})\leq\] \[\mathbb{E}_{\mathbf{W}_{k-1}\sim q(\cdot)}\left[\hat{V}_{k}(\phi_{k }(\mathbf{W}_{k-1}(z^{T},1)^{T}))\right].\]

Figure 1: Illustration of the DP algorithm in Theorem 1 for a BNN with two hidden layers. Value functions \(V_{k}\) are mappings from the latent and input spaces of the BNN to the mean of the output distribution. For each mapping, the distribution and mapping of a single point is displayed in orange. Starting from the last hidden layer, we recursively compute PWA approximations of the mappings. The true mean of the BNN for all \(z_{2}\in\mathcal{Z}_{2}\) is in the green oval, which we over-approximate by the blue hexagon.
_Further, for \(i\in\{1,\ldots,l\}\), it holds that \(\pi^{(i)}_{\min}(T)\geq\min_{x\in T}\hat{V}^{(i)}_{0}(x)\) and \(\pi^{(i)}_{\max}(T)\leq\max_{x\in T}\hat{V}^{(i)}_{0}(x)\)._

Corollary 1.1 allows us to recursively find relaxations of \(V_{k}\) via Theorem 1. In what follows, we focus on finding PWA relaxations \(\hat{V}_{k}\) and \(\hat{V}_{k}\). To achieve that, there are two basic steps: (i) initialization of \(\hat{V}_{k},\hat{V}_{k}\) via Eqn. (5a), and (ii) backward propagation of \(\hat{V}_{k},\hat{V}_{k}\) through a hidden layer of the BNN via Eqn. (5b). In Section 5, we first show an efficient method to perform step (ii) and then focus on (i).

## 5 PWA Relaxation for Dynamic Programming

Our goal in this section is to find PWA relaxations of \(V_{k}\). To do that, we first show how to propagate affine relaxations of the value function backwards through a single hidden layer of the BNN via Eqn. (5b) and then generalize this result to PWA relaxations. Note that, because the support of a BNN is generally unbounded, affine relaxations of Eqn. (5b) and (5a) lead to overly conservative results (an affine function should over-approximate a non-linear function over an unbounded set). Thus, PWA relaxations are necessary to obtain tight approximations. Finally, in Subsection 5.3 we show how to compute relaxations for Eqn. (5a).

### Affine Value Functions

For the sake of presentation, we focus on upper bound \(\hat{V}_{k-1}\); the lower bound case follows similarly. Let \(\hat{V}_{k}:\mathbb{R}^{n_{k}}\rightarrow\mathbb{R}^{l}\) be an affine upper bound on \(V_{k}\). Then, by Corollary 1.1 and the linearity of expectation, it holds that

\[\hat{V}_{k-1}(z)=\hat{V}_{k}(\mathbb{E}_{\mathbf{W}_{k-1}\sim q(\cdot)}\left[ \phi_{k}(\mathbf{W}_{k-1}(z^{T},1)^{T})\right]). \tag{6}\]

Recall that here \(q\) is a Gaussian distribution (see Section 2.1). Hence, due to the closure of Gaussian random variables w.r.t. linear transformations, we can rewrite Eqn. (6) as:

\[\hat{V}_{k-1}(z)=\hat{V}_{k}(\mathbb{E}_{\mathbf{\zeta}\sim\mathcal{N}(m_{k}(z); \text{ diag}(s_{k}(z)))}\left[\phi_{k}(\mathbf{\zeta})\right]), \tag{7}\]

where \(m_{k}:\mathbb{R}^{n_{k-1}}\rightarrow\mathbb{R}^{n_{k}}\) and \(s_{k}:\mathbb{R}^{n_{k-1}}\rightarrow\mathbb{R}^{n_{k}}_{\geq 0}\) are defined component-wise as

\[\begin{split}& m^{(i)}_{k}(z)=\mu_{w,k-1,i}(z^{T},1)^{T},\\ & s^{(i)}_{k}(z)=(z^{T},1)\Sigma_{w,k-1,i}(z^{T},1)^{T},\end{split} \tag{8}\]

for all \(i\in\{1,\ldots,n_{k}\}\) with \(\mu_{w,k-1,i},\Sigma_{w,k-1,i}\) being mean and covariance of the \(i\)th node of the \(k\)th layer. \(\text{diag}(s)\) is a diagonal matrix with the elements of \(s\) on the main diagonal. Note that Eqn (7) reduces the propagation of the value function to the propagation of a Gaussian random variable (\(\mathbf{\zeta}\)) through an activation function (\(\phi_{k}\)). In Proposition 2, we show how this propogation can be achieved analytically for ReLU activation functions. Generalization to other activation functions is discussed in Remark 3.

**Proposition 2**.: _For \(k\in\{1,\ldots,K\}\), let \(\hat{V}_{k}\) be an affine function and \(Z\subset\mathbb{R}^{n_{k-1}}\) be a compact set. Define function \(r_{k}:\mathbb{R}^{n_{k-1}}\rightarrow\mathbb{R}^{n_{k}}_{\geq 0}\) as \(r_{k}(z)=\sqrt{s_{k}(z)}\), and let \(\check{r}_{k},\check{r}_{k}:\mathbb{R}^{n_{k-1}}\rightarrow\mathbb{R}^{n_{k}}_ {\geq 0}\) be an affine-relaxation of \(r_{k}\) w.r.t. \(Z\). Further, define \(g:\mathbb{R}^{2}\rightarrow\mathbb{R}\) as_

\[g(\mu,\sigma)=\frac{\mu}{2}\left[1-\text{erf}\!\left(\frac{-\mu}{\sigma\sqrt{2 }}\right)\right]+\frac{\sigma}{\sqrt{2\pi}}\,e^{-(\mu/\sigma\sqrt{2})^{2}},\]

_and, let \(\check{g}_{i},\hat{g}_{i}:\mathbb{R}^{2}\rightarrow\mathbb{R}\) be an affine-relaxation of \(g\) w.r.t. \(\{(m^{(i)}_{k}(z),r^{(i)}_{k}(z))\mid\forall z\in Z\}\), Then, for \(\check{A},\hat{A}\in\mathbb{R}^{n_{k-1}\times n_{k}}\) and \(\check{b},\hat{b}\in\mathbb{R}^{n_{k}}\) defined as, \(\forall i\in\{1,\ldots,n_{k}\}\),_

\[\check{A}^{(i,\cdot)}=[\nabla_{z}g(m^{(i)}_{k}(z),r^{(i)}_{k}(z)) ]_{z=z^{*}},\] \[\check{b}^{(i)}=g(m^{(i)}_{k}(z^{*}),r^{(i)}_{k}(z^{*}))-\check{A} ^{(i,\cdot)}z^{*},\] \[[\cdot,\check{A}^{(i,z)}z+\check{b}^{(i)}]=(m^{(i)}_{k},\check{r} ^{(i)}_{k})^{T}\otimes[\check{g},\hat{g}],\]

_with \(z^{*}\in Z\) and \(\nabla_{z}\) being the gradient w.r.t. \(z\), it holds that \(\forall z\in Z\), \(\mathbb{E}_{\mathbf{\zeta}\sim\mathcal{N}(m_{k}(z);\text{ }s_{k}(z))}\left[\hat{V}_{k}(\text{ReLU}(\mathbf{\zeta}))\right]\in\hat{V}_{k} \otimes[\check{A}z+\check{b},\check{A}z+\hat{b}]\)._

The proof of Proposition 2 is based on the convexity of the expected value of a rectified Gaussian w.r.t. its mean and variance. The proof and detailed procedures for obtaining affine relaxations of \(g\) and \(r\) are reported in Appendix B.1. Next, we show how the result of Proposition 2 can be extended to PWA relaxations of the value functions.

**Remark 3**.: _The results of Proposition 2 (as well as Propositions 4 and 5 below) extend to any continuous activation function \(\phi_{k}\). That is, as shown in (Beunissi et al., 2022), every continuous activation function can be under and over-approximated by PWA functions \(\hat{\phi}_{k},\hat{\phi}_{k}:\mathbb{R}^{n_{k}}\rightarrow\mathbb{R}^{n_{k}}\) such that \(\hat{\phi}_{k}\leq\phi_{k}\leq\hat{\phi}_{k}\). Consequently, \(\mathbb{E}\!\left[\hat{\phi}_{k}(\mathbf{\zeta})\right]\leq\mathbb{E}\left[\phi_{k}( \mathbf{\zeta})\right]\leq\mathbb{E}\!\left[\hat{\phi}_{k}(\mathbf{\zeta})\right]\), which allows the extension of Proposition 2 from ReLU to general continuous \(\phi_{k}\)._

### Piecewise Affine Value Functions

For \(N\in\mathbb{N}\), let \(\mathcal{Z}_{k}=\{Z_{k,1},\ldots,Z_{k,N}\}\subseteq\mathbb{R}^{n_{k}}\) be a partition of the support of \(f^{\mathbf{\alpha}_{k}}_{0:k}\), and let \(\hat{V}_{k,j},\hat{V}_{k,j}:\mathbb{R}^{n_{k}}\rightarrow\mathbb{R}^{l}\) be an affine relaxation of \(V_{k}(z_{k})\) w.r.t. \(Z_{k,j}\) for all \(j\in\{1,\ldots,N\}\), i.e., \(\forall z_{k}\in Z_{k,j}\)\(V_{k}(z_{k})\leq\hat{V}_{k,j}\) with \(\hat{V}_{k,j}\coloneqq\hat{A}_{k,j}z_{k}+\hat{b}_{k,j}\). Then, by Eqn. (5b) and the law of total expectation, we obtain an upper bound on \(V_{k-1}\):

\[V_{k-1}(z)\leq\sum_{j=1}^{N}\hat{b}_{k,j}\underbrace{\mathbb{P}_{\mathbf{\zeta}\sim \mathcal{N}(m_{k}(z);\text{ diag}(s_{k}(z)))}\left[\mathbf{\zeta}\in Z_{k,j}\right]}_{9 \text{a}}+ \tag{9}\] \[\hat{A}_{k,j}\underbrace{\mathbb{E}_{\mathbf{\zeta}\sim\mathcal{N}(m_{k} (z);\text{ diag}(s_{k}(z)))}\left[\phi_{k}(\mathbf{\zeta})\mid\mathbf{\zeta}\in Z_{k,j}\right]}_{9 \text{b}},\]

where \(\mathbb{E}_{\mathbf{\zeta}\sim p}\left[\mathbf{\zeta}\mid\mathbf{\zeta}\in Z\right]:= \mathbb{E}_{\mathbf{\zeta}\sim p}\left[\mathbf{\zeta}\mid\mathbf{\zeta}\in Z\right]\mathbb{P}_{ \mathbf{\zeta}\sim p}\left[\mathbf{\zeta}\in Z\right]\). The lower bound on \(V_{k-1}\) follows similarly. Term 9a is
simply the probability that a Gaussian random variable (\(\zeta\)) is in a given set (partition \(Z_{k,j}\)). If the partition is hyperrectangular, in Lemma 3 we express Term 9a in closed-form.

**Lemma 3**.: _For \(k\in\{1,\ldots,K\}\), \(\tilde{\zeta},\hat{\zeta}\in\mathbb{R}^{n_{k}}\), it holds that 7_

Footnote 7: A similar result holds for unbounded regions defined by vector \(\tilde{z}\), that is, \(\mathbf{\zeta}\in[\tilde{z},\infty)\) or \(\mathbf{\zeta}\in(\infty,\tilde{z}]\), as shown in Appendix B.2.

\[\mathbb{P}_{\mathbf{\zeta}\sim\mathcal{N}(m_{k}(z);\;\text{diag}(s_{ k}(z)))}\left[\mathbf{\zeta}\in[\tilde{\zeta},\hat{\zeta}]\right]= \tag{10}\] \[\frac{1}{2^{n_{k}}}\prod_{i=1}^{n_{k}}\text{erf}\left(\frac{ \hat{\zeta}^{(i)}-m_{k}^{(i)}(z)}{\sqrt{2s_{k}^{(i)}(z)}}\right)-\text{erf} \left(\frac{\hat{\zeta}^{(i)}-m_{k}^{(i)}(z)}{\sqrt{2s_{k}^{(i)}(z)}}\right)\]

Term 9b is the conditional expectation of the random variable propagated through an activation function. The following shows that we can decompose this term in expectations, which we can bound using the result of Proposition 2, and probabilities for which Lemma 3 can be applied.

**Proposition 4**.: _For \(k\in\{1,\ldots,K\}\), vectors \(\tilde{\zeta},\hat{\zeta}\in\mathbb{R}^{n_{k}}\), and \(\mathbf{\zeta}\sim\mathcal{N}\left(m_{k}(z);\;\text{diag}\;(s_{k}(z))\right)\), it holds that8_

Footnote 8: A similar relation can be obtained for \(\phi_{k}\) being the identity function, as shown in Appendix B.3.

\[\mathbb{\mathbb{E}}\left[\text{ReLU}\left(\mathbf{\zeta}\right)\mid\mathbf{\zeta}\in[ \tilde{\zeta},\hat{\zeta}]\right]=\] \[\mathbb{\mathbb{E}}\left[\text{ReLU}\left(\mathbf{\zeta}+[\tilde{ \zeta}]_{+}\right)\right]-\mathbb{\mathbb{E}}\left[\text{ReLU}\left(\mathbf{\zeta }+[\tilde{\zeta}]_{+}\right)\right]+\] \[[\tilde{\zeta}]_{+}\mathbb{P}\left[\mathbf{\zeta}\in[[\tilde{\zeta}]_ {+},\infty)\right]-[\hat{\zeta}]_{+}\mathbb{P}\left[\mathbf{\zeta}\in[[\hat{\zeta}] _{+},\infty)\right].\]

Next, we show how these results can be extended to unbounded sets in partition \(\mathcal{Z}_{k}\), i.e., unbounded support \(f_{0:k}^{\mathbf{w}}\).

Unbounded SupportIf \(f_{0:k}^{\mathbf{w}}\) has an unbounded support, then there must necessarily be at least a region that is unbounded in the partition \(\mathcal{Z}_{k}\). While for this region we can still apply Lemma 3 to compute Term 9a, we cannot use Proposition 4 to compute a bound for Term 9b. Instead, we rely on Proposition 5 (below), which derives relaxations based on the fact that Gaussian distributions decay exponentially fast (thus, faster than a linear function).

**Proposition 5**.: _For \(k\in\{1,\ldots,K\}\), \(i\in\{1,\ldots,n_{k}\}\), and vector \(\tilde{\zeta}\in\mathbb{R}^{n_{k}}\), it holds that9_

Footnote 9: A similar relation can be obtained for \(\phi_{k}\) being the identity function, as shown in Appendix B.4.

\[\frac{1}{2}[m_{k}^{(i)}(z)]_{-}\leq\] \[\leq\frac{1}{2}[m_{k}^{(i)}(z)]_{+}+\sqrt{\frac{s_{k}^{(i)}(z)}{ 2\pi}}.\]

**Algorithm 1** Adversarial Robustness for Classification

\[\mathbb{\mathbb{E}}_{\mathbf{\zeta}\sim\mathcal{N}(m_{k}(z);\;\text{diag}(s_{ k}(z)))}\left[\text{ReLU}\left(\mathbf{\zeta}\right)\mid\mathbf{\zeta}\in[\tilde{ \zeta},\infty)\right]\] \[\leq\frac{1}{2}[m_{k}^{(i)}(z)]_{+}+\sqrt{\frac{s_{k}^{(i)}(z)}{ 2\pi}}.\]

\(\mathbb{\mathbb{E}}\left[\text{ReLU}\left(\mathbf{\zeta}\right)\mid\mathbf{\zeta} \in[\tilde{\zeta},\infty)\right]\)

### Relaxation of the Last Layer of BNN

We show how to compute interval relaxations of Eqn (5a). For the regression case (\(h=I\)), the process is simple since Eqn. (5a) becomes an affine function. That is, \(V_{K}(z)=m_{K}(z)\), where \(m_{K}(z)\) as defined in Eqn. (8), and hence no relaxation is required. For classification, however, further relaxations are needed because \(h=\text{softmax}\), i.e., the output distribution of the BNN (the logit) is propagated through the softmax. The following proposition shows that an interval relaxation can be obtained by relaxing the distribution of \(\mathbf{W}_{K}(z^{T},1)^{T}\) by Dirac delta functions on the extremes of \(h\) for each set in the partition of the BNN's output.

**Proposition 6**.: _For \(N\in\mathbb{N}\), let \(\{Z_{1},\ldots,Z_{N}\}\subseteq\mathbb{R}^{n_{K+1}}\) be a partition of \(\text{supp}(f^{\mathbf{w}}\left(x\right))\). Then, for \(i\in\{1,\ldots,n_{K+1}\}\) and \(\mathbf{w}\sim q(\cdot)\), it holds that_

\[\sum_{j=1}^{N}[\underset{\zeta\in Z_{j}}{\min}h^{(i)}(\zeta)] \mathbb{P}\left[f^{\mathbf{w}}\left(x\right)\in Z_{j}\right]\leq\mathbb{E}\left[ h^{(i)}(f^{\mathbf{w}}\left(x\right)\right]\leq\] \[\sum_{j=1}^{N}[\underset{\zeta\in Z_{j}}{\max}h^{(i)}(\zeta)] \mathbb{P}\left[f^{\mathbf{w}}\left(x\right)\in Z_{j}\right].\]

A particularly simple case is when there are only two sets in the partition of the BNN's output layer. Then, the following corollary of Proposition 6 guarantees that, similarly to deterministic NNs (Zhang et al., 2018), we can determine adversarial robustness by simply looking at the logit.

**Corollary 6**.: _Let \(\{[\tilde{\zeta},\hat{\zeta}],Z\}\subseteq\mathbb{R}^{n_{K+1}}\) be a partition of \(\text{supp}(f^{\mathbf{w}})\). Then, for \(i,j\in\{1,\ldots,n_{K+1}\}\) and \(\mathbf{w}\sim q(\cdot)\), it holds that_

\[\hat{e}^{\hat{\zeta}^{(j)}}-e^{\hat{\zeta}^{(i)}}+(\frac{1}{ \mathbb{P}\left[f^{\mathbf{w}}(x)\in[\tilde{\zeta},\hat{\zeta}]\right]}-1)\sum_{l=1 }^{n_{K+1}}e^{\hat{\zeta}^{(l)}}\leq 0\] \[\implies\mathbb{E}\left[\text{softmax}^{(j)}(f^{\mathbf{w}}(x))- \text{softmax}^{(i)}(f^{\mathbf{w}}(x))\right]\leq 0.\]
## 6 BNN-DP Algorithm

We summarize our overall procedure to solve Problem 1 in an algorithm called BNN-DP. Algorithm 1 presents BNN-DP for the classification setting; the procedure for the regression setting follows similarly and is provided in Appendix D. Algorithm 1 consists of a forward pass to partition the latent space of the BNN (Lines 2-4), and a backward pass to recursively approximate the value functions via Eqns. (5a) and (5b) (Lines 7-10). The last layer of the BNN (Eqn. 5a) is handled by the IBPSoftmax function in Line 7 using the results of Proposition 6. The BP function in Line 9 performs the back-propagation over the hidden layers of the BNN (Eqn. 5b) using the results of Lemma 3 and Proposition 4 and 5. The detailed procedures of IBPSoftmax and BP can be found in Appendix D. In what follows, we describe how we partition the support of the latent space of the BNN, and discuss the computational complexity of BNN-DP.

PartitioningRecall that our results rely on hyperrectangular partitions. Hence, for each layer \(k\), we employ the following proposition to find a hyper-rectangular subset of the support of each layer that captures at least \(1-\epsilon\) of the probability mass of \(\text{supp}(f_{0:k}^{u})\).

**Proposition 7**.: _For \(k\in\{1,\ldots,K\}\), let \(\epsilon\in[0,1]\) be a constant, and \(Z\subset\mathbb{R}^{n_{k-1}}\) be a compact set. Then, for vectors \(\hat{\zeta}_{k},\hat{\zeta}_{k}\in\mathbb{R}^{n_{k}}\) defined such that \(\forall i\in\{1,\ldots,n_{k}\}\),_

\[\hat{\zeta}_{k}^{(i)} =\max_{z\in Z}\left[\text{erf}^{-1}\left(-\eta\right)\sqrt{2s_{k} ^{(i)}(z)}+m_{k}^{(i)}(z)\right], \tag{11}\] \[\hat{\zeta}_{k}^{(i)} =\min_{z\in Z}\left[\text{erf}^{-1}\left(\eta\right)\sqrt{2s_{k} ^{(i)}(z)}+m_{k}^{(i)}(z)\right], \tag{12}\]

_where \(\eta=(1-\epsilon)^{\frac{1}{n_{k}}}\), it holds that, \(\forall z\in Z\),_

\[\mathbb{P}_{\boldsymbol{\zeta}\sim\mathcal{N}(m_{k}(z);\text{ diag}(s_{k}(z)))}\left[\boldsymbol{\zeta}\in[\tilde{\zeta}_{k},\hat{\zeta}_{k}] \right]\geq 1-\epsilon.\]

Here, Eqns (11) and (12) are convex minimization problems, which can be efficiently solved via, e.g., the gradient descent algorithm. We denote the resulting region obtained via Proposition 7 as \(Z_{k,main}\subset\text{supp}(f_{0:k}^{u})\). Then, \(Z_{k,main}\) can be further refined by interval splitting.

Computational ComplexitySimilarly as for linear bounding procedures for deterministic neural networks, see e.g. [Zhang et al.2018], the cost of computing piecewise-affine relaxations of a BNN with \(K\) layers and \(n\) neurons per layer is polynomial in both \(K\) and \(n\). Refinement, which is not part of the main algorithm, has exponential cost in \(n\). In practice, however, in NNs and consequently in BNNs, only a few neurons are generally active, and those are the ones that most influence the posterior [Frankle and Carbin2018]. Therefore, the refining procedure can focus only on these neurons. Because of this, in almost all the experiments in Section 7, only \(2\) regions in the partition per hidden layer were required to certify robustness, even in cases where the BNN had large posterior variance and thousands of neurons.

## 7 Experimental Results

We empirically evaluated BNN-DP on various regression and classification benchmarks. We ran our experiments on an AMD EPYC 7252 8-core CPU and train the BNNs using Noisy Adam [Zhang et al., 2018] and variational online Gauss-Newton [Khan et al., 2018]. We first validate the bounds obtained by BNN-DP for BNNs trained on samples from an 1D sine with additive noise (referred to as the 1D Noisy Sine). We then analyse a set of BNNs with various architectures trained on the 2D dimensional equivalent of 1D Noisy Sine and the Kinsnm dataset.10 The latter dataset contains state-space readings for the dynamics of an 8 link robot arm, and is commonly used as a regression task to benchmark BNNs [Hernandez-Lobato and Adams, 2015, Gal and Ghahramani, 2016]. Last, we turn our attention to classification and evaluate BNNs trained on the MNIST, Fashion MNIST and CIFAR-10 datasets. 11

Footnote 10: Available at [http://www.cs.toronto.edu/~delve](http://www.cs.toronto.edu/~delve).

Footnote 11: Our code is available at [https://github.com/sjladams/BNN_DP](https://github.com/sjladams/BNN_DP).

As a baseline for our experiments, we consider the state-of-the-art approach of berrada2021, to which we refer as "FL". In fact, FL is the only existing method that can provide robustness certification for BNNs in similar settings as our BNN-DP. Nevertheless, we must remark that even

Figure 2: Certified affine bounds on the mean of BNNs trained on the 1D Noisy Sine dataset w.r.t. the grey marked interval of input.
FL is not fully formal; it works by truncating the Gaussian posterior distribution associated to each weight at a given multiple of its standard deviation (std), disregarding a large part of the posterior distribution. Hence, the returned bound is not sound over the full posterior but only a subset of it. More importantly, the disregarded portion of the posterior grows exponentially with the number of weights of the networks. Already for a two hidden layer BNN with 48 neurons per layer, FL verifies only \(0.1\%\) of the BNN posterior when truncated at 3 std. Thus, the bounds computed by FL are optimistic and not mathematically guaranteed to hold. In contrast, not only BNN-DP returns formal bounds accounting for the whole posterior, but also the benchmark results show that BNN-DP bounds are much tighter than FL ones.

### Bound Validation

We validate and qualitatively compare the bounds obtained by BNN-DP and FL on BNNs with 1 and 2 hidden layers trained on 1D Noisy Sine. The results of these analyses are reported in Figure 2. Visually, we see that BNN-DP is able to compute tight affine relaxations (blue lines) on the mean of the BNNs over the grey shaded intervals. In contrast, already in this simple scenario, and even when truncating the posterior distribution at just 1 std, FL returns as guaranteed output intervals \([-1.68,0.59]\) and \([0.23,0.71]\) for the 1 and 2 hidden layer BNN, respectively. Hence, even though FL disregards most of the BNNs posterior, BNN-DP still produces tighter bounds. When using 3 std, the FL interval bounds become even wider, that is \([-7.07,9.31]\) and \([-0.65,1.54]\), for the 1 and 2 hidden layer BNN, respectively. Intuitively, the major improvement of the bounds can be explained by the fact that, while BNN-DP directly averages the uncertainty of each layer by solving the DP in Theorem 1, FL solves an overall optimisation problem that at each layer considers the worst combination of parameters in the support of the (truncated) distribution, leading to conservative bounds. In fact, the bound computed by FL is looser in the one-hidden layer case than in the two-hidden layers one by one order of magnitude, precisely because of the higher variance of the former BNN compared to the second. In what follows, we see that analogous observations apply to more complex deep learning benchmarks.

### Regression Benchmarks

We consider a set of BNNs with various architectures trained on the 2D Noisy Sine and Kin8nm regression datasets. To asses the certification performance of BNN-DP, we compute the difference between the upper and lower bounds on the expectation of the BNNs, referred to as the \(\gamma\)-robustness, for its input in a \(\ell_{\infty}\)-norm ball of radius \(\epsilon\) centered at a sampled data point. Clearly, a smaller value of \(\gamma\) implies a tighter bound computation. Results, averaged over 100 randomly sampled test points, are reported in Table (a)a. For all experiments, BNN-DP greatly improves the value of \(\gamma\)-robustness provided by the FL-baseline by 1 to 4 orders of magnitude with similar computation times. We also note that the larger the BNN is, the larger the improvement in the (tightness of the) bounds are, which empirically demonstrates the superior scalability of BNN-DP. Figure 3 explicitly shows the impact of the model size and variance on the certified \(\gamma\)-robustness. For BNNs with 1 hidden layer, BNN-DP guarantees small \(\gamma\)-robustness (and hence tighter bounds) irrespective of the number of neurons as well as the amount of uncertainty. In contrast, as already observed for the 1D Noisy Sine case, FL is particularly impacted by the variance of the posterior distribution. For BNNs with two hidden layers, BNN-DP requires partitioning the latent space, which leads to a positive correlation with the value of \(\gamma\)-robustness and the number of hidden neurons. A similar, but more extreme, trend is also observed for FL.

\begin{table}

\end{table}
Table 1: Comparison between BNN-DP and FL on various fully connected BNN architectures, with \(K\) being the number of hidden layers, and \(n_{hid}\) the number of neurons per layer. The results are the average over \(100\) test point, and the computation times are averaged over all architectures. The best values for each comparison are reported in bold.
### Classification Benchmarks

We now evaluate BNN-DP on the MNIST, Fashion MNIST and CIFAR-10 classification benchmarks. In order to quantitatively measure the robustness of an input point \(x^{*}\), we consider the maximum radius \(\epsilon\) for which the decisions on \(\ell_{\infty}\)-norm perturbations of \(x^{*}\) with radii \(\epsilon\) is invariant. That is, any perturbation of \(x^{*}\) smaller than \(\epsilon\) does not change the classification output; hence, the larger \(\epsilon\) in \(x^{*}\), the more robust the BNN in the specific point. Results are reported in Table 0(b) and 2. For the fully connected BNN architectures, BNN-DP not only is able to certify a substantially larger \(\epsilon\) compared to the baseline, but also it does so by orders of magnitude smaller computation time. This is because our approach uses interval relaxations (Proposition 6) to bound the softmax, whereas FL explicitly considers a non-convex optimization problem, which is computationally demanding. For the Bayesian CNN architectures, FL is able to certify a slightly larger \(\epsilon\), at the costs of magnitudes of orders increase of computation time. This can be explained by the decreasing support of the BNN posterior certified by FL for increasing network size, whereas, the \(\epsilon\) certified by BNN-DP holds for the whole posterior.

## 8 Conclusion

We introduced BNN-DP, an algorithmic framework to certify adversarial robustness of BNNs. BNN-DP is based on a reformulation of adversarial robustness for BNNs as a solution of a dynamic program, for which efficient relaxations can be derived. Our experiments on multiple datasets for both regression and classification tasks show that our approach greatly outperforms state-of-the-art competitive methods, thus paving the way for applications of BNNs in safety-critical applications.

## Acknowledgements

This work was supported in part by the NSF grant 2039062.
