Make Interval Bound Propagation great again

## Abstract

Abstract In various scenarios motivated by real life, such as medical data analysis, autonomous driving, and adversarial training, we are interested in robust deep networks. A network is robust when a relatively small perturbation of the input cannot lead to drastic changes in output (like change of class, etc.). This falls under the broader scope field of Neural Network Certification (NNC).
Two crucial problems in NNC are of profound interest to the scientific community: how to calculate the robustness of a given pre-trained network and how to construct robust networks. The common approach to constructing robust networks is Interval Bound Propagation (IBP).
This paper demonstrates that IBP is sub-optimal in the first case due to its susceptibility to the wrapping effect. Even for linear activation, IBP gives strongly sub-optimal bounds. Consequently, one should use strategies immune to the wrapping effect to obtain bounds close to optimal ones. We adapt two classical approaches dedicated to strict computations – Dubleton Arithmetic and Affine Arithmetic – to mitigate the wrapping effect in neural networks. These techniques yield precise results for networks with linear activation functions, thus resisting the wrapping effect. As a result, we achieve bounds significantly closer to the optimal level than IBPs.

## 1 Introduction

Deep neural networks find application in medical data analysis, autonomous driving, and adversarial training (Zhang et al., ) where safety-critical and robustness guarantees against adversarial examples (Biggio et al., ; Szegedy et al., ) are extremely important.
The rapid development of artificial intelligence models does not correspond to their robustness (Luo et al., ).
Therefore, certifiable robustness (Zhang et al., ; Ferrari et al., ) becomes an important task in deep learning. The aim of Neural Network Certification lies in rigorous validation of a classifier’s robustness within a specified input region.

Most commonly applied certification method is interval bound propagation (IBP) (Gowal et al., ; Mirman et al., ). It is based on application of interval arithmetic, which allows to propagate the input intervals through a neural network. If in the case of classification tasks such propagation gives an unambiguous output then all elements of the interval inputs are guaranteed to have identical prediction. Therefore, we can control the behavior of predictions of a neural network in an explicit neighborhood of the input data.
Among the approaches studied most extensively in robustness of neural networks is Certified Training (Singh et al., ; Mao et al., ).
These certified training methods try to estimate and optimize the worst-case loss approximations of a network across an input domain defined by adversary specifications. They achieve this by computing an over approximation of the network reachable set through symbolic bound propagation techniques (Singh et al., ; Gowal et al., ). Interestingly, training techniques based on the least accurate bounds derived from interval-bound propagation (IBP) have delivered the best empirical performance (Shi et al., ; Mao et al., ).

The certification process uses the network’s upper bound of the propagated input interval. Although classical IBP gives reasonable estimations in robust training, it ultimately fails in the certification of classically trained neural networks. In practice, for a given pre-trained networks, intervals that store intermediate values in a neural network evaluation increase exponentially with respect to number of layers, see Theorem 2.1. This phenomena is known as *the wrapping effect* (Neumaier, ), which in the context of neural networks applications was previously an unexplored area.
Such a growth of obtained bounds makes them often dramatically sub-optimal in practical applications.

The aim of this paper is to analyse and adapt two existing methods for reduction of the wrapping effect to the context of neural networks applications(^1^11[https://github.com/gmum/Make-Interval-Bound-Propagation-great-again](https://github.com/gmum/Make-Interval-Bound-Propagation-great-again)): Doubleton Arithmetics (DA) (Mrozek & Zgliczyński, ) and Affine Arthmetic (AA) (de Figueiredo & Stolfi, ).

Doubletons is very special family of subsets of $\mathbb{R}^{n}$, which has been extensively used to reduce and control the wrapping effect in validated solvers to initial value problems of ODEs (Kapela et al., ). Direct application of interval arithmetics to propagate sets along trajectories, that is enclosing them into the Cartesian product of closed intervals, leads to accumulated overestimation known as the wrapping effect. This overestimation becomes larger and larger when we use smaller time steps $h$ of the underlying ODE solver.
Lohner (Lohner, ) observed, that one can propagate coordinate system (approximate space derivative of the flow) between subsequent time steps along trajectories of the flow. This method proved to be very efficient and one of the reasons is that the for small time steps the mapping defined as a time shift along trajectories is close to identity.

Figure: Figure 1: The figure presents how the interval is propagated throughout linear layers. By red color we marked wrapping obtain by IBP and by green by Affine Arithmetic. As we can see, Affine Arithmetic produces significantly lower wrapping effects. In the case of linear transformations, Affine Arithmetic gives an exact approximation. We can work with more complex objects than hyper-cubes from IBP and obtain bounds close to optimal ones. In Fig. 3 we present the procedure used in Affine Arithmetic to obtain $\mathrm{ReLU}(I^{1})$.
Refer to caption: extracted/5901384/img/in_2.png

In this paper we adopt Doubleton Arithmetics to the context of neural networks. The main difference in comparison to ODEs is that the dimensions of subsequent layers in a network are usually different, while in ODEs we have a fixed dimension of the phase space. Moreover, in neural networks we often deal with non-smooth functions, such as $\mathrm{ReLU}$. In Section 3 we will formally define doubleton representation and give algorithm for propagation of non-smooth functions in this arithmetics.

The second method, called Affine Arthmetic (AA) (de Figueiredo & Stolfi, ), is a special case of Taylor Models by (Berz & Makino, ; Makino & Berz, ). Here subsets of $\mathbb{R}^{n}$ are represented as a range of an affine map (often sparse) over a cube $[-1,1]^{m}$, where $m$ in general is not related to $n$. Similarly to DA, evaluation of affine layers in AA causes no wrapping effect and it is sharp. In Section 3 and Appendix we show how to implement $\mathrm{ReLU}$ and $\mathrm{softmax}$ functions over a set represented in this way. The main numerical drawback of AA is that its computational cost is non-constant and depends on actual input arguments. Our experiments show that AA outperforms IBP in obtained bounds, see Fig. 1. Although AA and DA provide bounds of comparable sizes, AA is much faster (orders of magnitude) than DA on large networks. Thus, we recommend AA for evaluation of interval inputs through neural networks.

To make our approach completely certifiable, we need to have the full control over the numerical and rounding errors appearing in floating point arithmetics (Kahan, ). To obtain this we have decided to switch from Python-based networks to interval arithmetics (IEEE Std 1788.1-2017, ) in C++ with the use of the CAPD library (Kapela et al., ), which gives us certifiable control over the rounding errors. Consequently, to the best of the authors’ knowledge, the presented approach is the first model that deals with rounding errors and obtains guaranteed boundaries.

Our contributions can be summarized as follows:

- •
We theoretically analyze the wrapping effect in the Neural Network certification task and show that classical IBP is sub-optimal even for linear transformations.
- •
We adapt two approaches, Doubleton Arithmetic and Affine Arithmetic, with full control over numerical and rounding errors to the neural network certification task.
- •
Using empirical evaluation, we show that Affine Arithmetic gives the best bounds of the neural network output and significantly outperforms classical IBP.

## 2 IBP and wrapping effect

In this section, we examine how interval bounds propagate through a linear layer. We show that the appearance of wrapping effect, even in the case of isometric transformations, leads to an exponential growth of interval bounds.
Wrapping effect is typically studied in the context of strict estimations for solutions of dynamical systems, where the propagated set is at each iteration “wrapped” in the minimal interval bound (Neumaier, ).
Therefore, applying the standard interval bound propagation layer after layer leads to an exponential increase of bounds.

Given a bounded set $X\subset\mathbb{R}^{n}$, by $\mbox{$\mathrm{IB}$}(X)$ (interval bounds) we denote the smallest interval bounding box for $X$. The aim of IBP (Interval Bound Propagation) lies in obtaining the $\mathrm{IB}$ for the processing of $X$ through a network, i.e. a series of possibly nonlinear maps.
In the case of linear map $A=[a_{ij}]$, the optimal bounds are given by

$$ $\mbox{$\mathrm{IB}$}(A(x+[-r,r]))=Ax+[-|A|r,|A|r],$ (1) $$

where $x+[-r,r]=\prod_{i}[x_{i}-r_{i},x_{i}+r_{i}]$ and $|A|=[|a_{ij}|]$.
To propagate an interval through the $\mathrm{ReLU}$ activation, we propagate the lower and upper bound separately:
$\mathrm{ReLU}([x,y])=[\mathrm{ReLU}(x),\mathrm{ReLU}(y)].$
We can propagate intervals through the standard network $\Phi$, which is represented as a sequence of mappings corresponding to the successive layers
$y=\Phi(x)=\phi_{k}\circ\ldots\circ\phi_{1}(x)$.
The aim of IBP is to obtain the estimate of
$\mbox{$\mathrm{IB}$}(\Phi(x+[-r,r])),$
where commonly we restrict to the case when $r=\varepsilon\mathbbm{1}$:

$$ $\mbox{$\mathrm{IB}$}(\Phi(x+\varepsilon[-\mathbbm{1},\mathbbm{1}])),\text{ where }\mathbbm{1}=(1,\ldots,1)\in\mathbb{R}^{n}.$ $$

The standard classical approach used for IBP in the networks uses the naive iterative approach, where we process through each layer the interval bounds obtained from the previous one:

$$ $[I=x+[-r,r]]\to[I^{1}=\mbox{$\mathrm{IB}$}(\phi_{1}(I))]\to[I^{2}=\mbox{$ \mathrm{IB}$}(\phi_{2}(I^{1}))]\to\ldots\to[y=\mbox{$\mathrm{IB}$}(\phi_{k}(I^ {k}))].$ $$

Since we compute interval bound in each stage, the estimations are far from optimal; see Fig. 1. In practice, wrapping effects appear in neural networks. We will show that intervals grow exponentially, even for linear networks. We consider linear orthogonal ones, as they can be seen as the natural initialization of the deep network ().

We will need the following lemma which proof is given in the Appendix.

###### Lemma 2.1 .

Let $V=(V_{1},\ldots,V_{n})$ be a random vector uniformly chosen from the unit sphere in $\mathbb{R}^{n}$. Let $R$ be a random variable given by $R=|V_{1}|+\ldots+|V_{n}|.$
Then

$$ $\mathbb{E}(R)=\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}+O(1/\sqrt{n}),\,\mathbb{V}(R )=1+\frac{1}{\pi}+O(1/n).$ $$

Now we will show how a uniform interval bound is processed through an orthogonal map (isometry).

###### Proposition 2.1 .

Let $U$ be a randomly chosen orthogonal map in $\mathbb{R}^{n}$. Then

$$ $\mbox{$\mathrm{IB}$}(U([-\mathbbm{1},\mathbbm{1}]))\approx\frac{\sqrt{2}}{ \sqrt{\pi}}\sqrt{n}\cdot([-\mathbbm{1},\mathbbm{1}]+O(1/\sqrt{n})).$ $$

###### Proof.

For each fixed $i=1\ldots,n$ the $i$-th row of $U$ is a random vector uniformly chosen from the unit sphere. Thus $U(x)=[U_{1}(x),\ldots,U_{n}(x)]^{T}$.
By (1), $\mbox{$\mathrm{IB}$}(U_{i}([-\mathbbm{1},\mathbbm{1}]))=[-R_{i},R_{i}]$,
where $R_{i}$ is a random variable given by
$R_{i}=|U_{i1}|+\ldots+|U_{in}|$.
Now by Lemma 2.1,
$\mathbb{E}(R_{i})=\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}+O(1/n)$, $\mathbb{V}(R_{i})=1+\frac{1}{\pi}+O(1/n)$.
By the Chebyshev inequality,

$$ $\operatorname{P}(|R_{i}-\mathbb{E}[R_{i}]|\geq a)\leq{\frac{\mathbb{V}([R_{i}] )}{a^{2}}},$ $$

and consequently asymptotically for large $n$

$$ $\operatorname{P}(|R_{i}-\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}|\geq a)\leq{\frac{ 1+\frac{1}{\pi}+O(1/n)}{a^{2}}}.$ $$

Consequently, with an arbitrary large probability

$$ $U_{i}\approx\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}\cdot[-1,1]+O(1)=\frac{\sqrt{2} }{\sqrt{\pi}}\sqrt{n}\cdot\left([-1,1]+O(1/\sqrt{n})\right).$ $$

∎

The following theorem shows that the standard IBP leads to an exponential increase of the bound with respect to the number of layers, even when the true optimal bound does not increase.
We obtain the formal proof for the linear layers with orthogonal activations.

###### Theorem 2.1 .

Let $U_{1},\ldots,U_{k}$ be a sequence of orthogonal maps in $\mathbb{R}^{n}$, and let $U=U_{k}\circ\ldots\circ U_{1}$.

Let $B_{0}=[-\mathbbm{1},\mathbbm{1}]$. Then

$$ $\mbox{$\mathrm{IB}$}(U(B_{0}))\approx\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}([- \mathbbm{1},\mathbbm{1}]+O(1/\sqrt{n}).$ $$

Let $B_{i}$ be defined iteratively by
$B_{i}=\mbox{$\mathrm{IB}$}(U_{i}(B_{i-1})).$
Then

$$ $B_{k}\approx(\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n})^{k}([-\mathbbm{1},\mathbbm{1 }]+O(1/\sqrt{n}))$ $$

###### Proof.

The proof follows from the recursive use of the previous proposition.
∎

Observe, that the above theorem says, that applying standard interval bounds propagation layer after layer leads to exponential increase in the bound, as compared to the true optimal bound.
This paper modifies two Doubleton and Affine Arithmetic models, which provide optimal bounds for linear transformations.

###### Lemma 2.1 .

###### Proposition 2.1 .

###### Proof.

###### Theorem 2.1 .

###### Proof.

## 3 Doubleton and Affine Arithmetics

As shown in the previous section, classical interval bound propagation leads to an exponential increase in the bounds, even for the case of most superficial linear networks, which implies that it is suboptimal for pre-trained networks. Consequently, we postulate that we should develop methods that obtain strict estimation in the case of linear networks. In this paper, we propose adapting two Doubleton and Affine Arithmetics models for deep neural networks.

### Doubleton Arithmetics

Doubleton is a class of subsets of $X\subset\mathbb{R}^{n}$ that are represented in the following form

$$ $X=\{x+Cr+Qq:r\in\mathbf{r},q\in\mathbf{q}\},$ $$

for some $x\in\mathbb{R}^{n}$, $C\in\mathbb{R}^{n\times m}$, $Q\in\mathbb{R}^{n\times k}$ and $\mathbf{r}\subset\mathbb{R}^{m}$, $\mathbf{q}\subset\mathbb{R}^{k}$ are interval vectors (product of intervals) containing zero. For a possibly nonlinear function $f:\mathbb{R}^{m}\to\mathbb{R}^{n}$ and a compact set $W\subset\mathbb{R}^{m}$ we use doubletons to enclose range $f(W)$. The component $x+C\mathbf{r}$ is supposed to store linear approximation to $f$, while $Q\mathbf{q}$ stores accumulated errors (usually bounds on nonlinear terms) in certain (often orthogonal) coordinate system.

Such a family is one of the most frequently used in validated integration of ODEs (Lohner, ; Mrozek & Zgliczyński, ) and provides a good balance between accuracy (size of overestimation) add time complexity of operations on such objects. Here we would like to adopt it to the special case of neural networks. We have to extend doubleton arithmetics to functions with different dimensions of domain and codomain and also for non-smooth $\mathrm{ReLU}$ frequently used as an activate function in neural networks.

In the context of neural networks Doubleton Arithmetics is promising since we obtain sharp bound when mapping a doubleton by an affine transformation.

###### Theorem 3.1 .

Evaluation of an affine function $A(t)=x_{0}+Lt$ over a doubleton $X=x+C\mathbf{r}+Q\mathbf{q}$ is exact, that is

$$ $A(X)=\left\{\tilde{x}+\tilde{C}r+\tilde{Q}q:q\in\mathbf{q},r\in\mathbf{r} \right\}\text{, where }\tilde{x}=x_{0}+Lx,\quad\tilde{C}=LC,\quad\tilde{Q}=LQ.$ $$

Consequently, we can process sets described by doubletons through linear layers without any wrapping effect. Enclosing a classical neural network activation function in Doubleton Arithmetics is more challenging. Below we proceed with the formulation how it can be done for a general nonlinear map.

Assume that $f:\mathbb{R}^{n}\to\mathbb{R}^{d}$ is a nonlinear function (even not continuous) and assume that for $z\in X=x+C\mathbf{r}+Q\mathbf{q}$ there holds

$$ $f(z)=x_{0}+L(z-x)+e(z)$ $$

and let us assume that we have computed a bound $e(z)\in\mathbf{e}$ for $z\in X$. Then we have

$$ $\displaystyle f(z)$ $\displaystyle=$ $\displaystyle x_{0}+(LC)r+(LQ)q+e(z)=\left(x_{0}+\mathrm{mid}(\mathbf{e}) \right)+(LC)r+(LQ)q+(e(z)-\mathrm{mid}(\mathbf{e}))$ $\displaystyle\in$ $\displaystyle\tilde{x}+\tilde{C}\mathbf{r}+\tilde{Q}\tilde{\mathbf{q}},$ $$

where
$\tilde{x}=x_{0}+\mathrm{mid}(\mathbf{e})\in\mathbb{R}^{d},\quad\tilde{C}=LC\in
\mathbb{R}^{d\times m}$ and the term $\tilde{Q}\tilde{\mathbf{q}}$ is computed as follows. To simplify notation put $\Delta=\mathbf{e}-\mathrm{mid}(\mathbf{e})$.
Let $\tilde{Q}\in\mathbb{R}^{d\times n}$ and $A\in\mathbb{R}^{n\times d}$ be arbitrary matrices so that $\tilde{Q}A=\mathrm{Id}_{d}$. Then for $q\in\mathbf{q}$ and $\delta\in\Delta$ we have

$$ $(LQ)q+\delta=\tilde{Q}A(LQq+\delta)=\tilde{Q}\left((ALQ)q+A\delta\right)$ $$

Now we define $\tilde{\mathbf{q}}:=(ALQ)\mathbf{q}+A\Delta\subset\mathbb{R}^{d}$. It should be emphasized, that it is very important to first evaluate the product of matrices $(ALQ)$ and then multiply the result by the interval vector $\mathbf{q}$. Here is the place when we can reduce wrapping effect provided we make a good choice of $\tilde{Q}$ and $A$. There are various strategies for that.

Strategy 1. If $n=d$ and $LQ\in\mathbb{R}^{n\times n}$ is nonsingular then we may set $\tilde{Q}=LQ$ and $A=\tilde{Q}^{-1}$. Then $\tilde{\mathbf{q}}:=\mathbf{q}+\tilde{Q}^{-1}\Delta\subset\mathbb{R}^{n}$.

Strategy 2 - $QR$-decomposition. If $d\geq n$ then we can first compute $QR$-decomposition of $LQ=\tilde{Q}R$, where $\tilde{Q}\in\mathbb{R}^{d\times d}$ is orthonormal. Then $A=\tilde{Q}^{T}$ and $R=ALQ\in\mathbb{R}^{d\times k}$. We can set $\tilde{\mathbf{q}}:=R\mathbf{q}+\tilde{Q}^{T}\Delta\subset\mathbb{R}^{d}$. If $d<n$ then we may compute $QR$-decomposition of the leading $d\times d$ block of $LQ$ and proceed as before.

Strategy 3 - $QR$-decomposition with pivots. In Strategy 2 we add a preconditioning step – that is permutation of columns of $LQ$ before $QR$-factorization. The permutation should take into account widths of components $\mathbf{q_{i}}$ and $\Delta_{i}$.

We can also try hybrid strategies in the case $n=d$. For instance, we can start from Strategy 1 and if $LQ$ is singular or close to singular we switch to Strategy 3.

Figure: Figure 2: Graphs of $\mathrm{ReLU}$ over a hyperplane crossing zero. (Left) first affine approximation of $\mathrm{ReLU}$ with $\tau=1$, that is $\widetilde{b}_{0}+c\sum{a_{i}t_{i}}$ and (right) its final affine approximation $b_{0}+c\sum{a_{i}t_{i}}$.
Refer to caption: extracted/5901384/img/relu1.png

We can implement $\mathrm{ReLU}$ and $\mathrm{softmax}$ in Doubleton Arithmetics thanks to the above methodology. Consequently, we can track how the input interval is propagated through the neural network. Since the transformation by linear mapping does not cause the wrapping effect, we get much better estimates than the classical IBP.

###### Theorem 3.1 .

### Affine Arithmetic

The main drawback of the Doubleton Arithmetics is that it is expensive, because it involves multiplication of full dimensional non-sparse matrices. Affine arithmetics (de Figueiredo & Stolfi, ) is a concept of reducing overestimation in evaluating an expression in interval arithmetics coming from dependency, that is multiple occurrence of a variable in an expression.

Affine arithmetics keeps track of linear dependencies between variables through evaluation of an expression.
Affine Arithmetic gives sharp bounds for linear transformations.

In affine aritmetics we represent subsets of $\mathbb{R}^{n}$ as a range of an affine functions $A([-1,1]^{m})$ for some affine map $A:\mathbb{R}^{m}\to\mathbb{R}^{n}$. Clearly the composition of $A$ with another affine map $B:\mathbb{R}^{n}\to\mathbb{R}^{k}$ is again an affine map $B\circ A:\mathbb{R}^{m}\to\mathbb{R}^{k}$ and thus the image of $[-1,1]^{m}$ via $B\circ A$ is represented as an affine expression with no overestimation.

Let us present on an easy example the main property of affine arithmetics, which shows its superiority over interval arithemtics. Assume we have two expressions $A(x,y)=1+x+2y$ and $B(x,y,z)=1-x-2y+z$, where $x,y,z\in I:=[-1,1]$ and we would like to compute a bound on $A(I,I)+B(I,I,I)$. Evaluation in interval arithmetics gives

$$ $\displaystyle A(I,I)+B(I,I,I)\subset\left(1+[-1,1]+2[-1,1]\right)+\left(1-[-1, 1]-2[-1,1]+[-1,1]\right)=[-5,9].$ $$

We see that multiple occurrence of a variable in an expression leads to large overestimation. In affine arithmetics we keep linear track of variables and only in the end we evaluate expression in interval arithmetics. This gives the following (sharp) bound

$$ $\displaystyle A(x,y)+B(x,y,z)$ $\displaystyle=$ $\displaystyle 2+z,$ $\displaystyle A(I,I)+B(I,I,I)$ $\displaystyle=$ $\displaystyle 2+[-1,1]=[1,3].$ $$

Because different affine functions may have different number of arguments it is convenient to treat them (formally) as functions $A:\ell^{0}\to\mathbb{R}$, where $\ell^{0}$ is a set of sequences with all but finite number of non-zero elements. Then we have a straightforward interpretation of addition of such functions and multiplication of affine function by a scalar.

To adapt Affine Arithmetic to neural networks we need to implement $\mathrm{ReLU}$ and $\mathrm{softmax}$ functions. In the case of nonlinear transformation in Affine Arithmetic we approximate our nonlinear mapping by an affine function with known precision, see Fig. 2. Then we propagate our input interval throught this affine transformation and add a new interval equal to upper bound of the difference between linear approximation and original function, see Fig. 3.

Figure: Figure 3: Affine Arithmetic works with more complicated shapes than hypercubes from IBP. In the example, we take Interval $I=[-1,1]^{2}$ and see how AA produces an approximation of output from the linear layer with ReLu activation. We use affine transformation from Fig. 1. To approximate AA output from $\mathrm{ReLU}(A(I))$, we first approximate nonlinear function $\mathrm{ReLU}(A(\cdot))$ by linear $B(\cdot)$ . Then, we propagate input interval $I$ through $B(\cdot)$. Then we add interval correction, which is equal to the maximal error between $\mathrm{ReLU}(A(I))$ and $B(\cdot)$ denoted by $b_{n+1}$. Finally, we obtain bound in Affine Arithmetic in the case of mapping interval through a linear layer with linear activation.
Refer to caption: extracted/5901384/img/I_ex.png

In order to implement $\mathrm{ReLU}$ in Affine Arithmetic let us consider an affine function

$$ $A(t_{1},\cdots,t_{n})=a_{0}+\sum_{i=1}^{n}a_{i}t_{i}$ $$

defined on the cube $t=(t_{1},\ldots,t_{n})\in[-1,1]^{n}=:I^{n}$ and assume $0\in A(I^{n})$. Clearly the composition $\mathrm{ReLU}\circ A$ is nonlinear and the set $\mathrm{ReLU}(A(I^{n}))$ cannot be represented exactly as a range of an affine function.

Our strategy is to find an affine map $B:\mathbb{R}^{n}\to\mathbb{R}$, which approximates well the composition $\mathrm{ReLU}\circ A$ on the hypercube $I^{n}$. Bound on the difference

$$ $\max_{t\in I^{n}}|\mathrm{ReLU}(A(t))-B(t)|$ $$

will be treated as a new variable and finally the range $\mathrm{ReLU}(A(I^{n}))$ will be covered by a range of an affine function but with $n+1$ variables. This scenario is visualised in Fig. 3.

We impose that $B$ is of the form

$$ $B(t)=b_{0}+\sum_{i=1}^{n}(ca_{i})t_{i}$ $$

for some $c\in\mathbb{R}$, that is $b_{i}=ca_{i}$ for $i>0$. Put $S:=\sum_{i=1}^{n}|a_{i}|$ and let $M:=\sup_{t\in I^{n}}A(t)=a_{0}+S$. Let $\tau\in[0,1]$ be a parameter to be specified later. We impose that $B$ vanishes for its all arguments being $-1$, while it reaches maximum value in $I^{n}$ equal to $\tau M$ for all arguments equal to $1$ – see Fig.2 left panel. This gives the following system of equations with two unknowns

$$ $\widetilde{b}_{0}-cS=0,\qquad\widetilde{b}_{0}+cS=\tau\cdot M.$ $$

The solution is $\widetilde{b}_{0}=\frac{1}{2}\tau M$ and $c=\frac{1}{2}\tau M/S$.
The graph of the first affine approximation $\widetilde{B}(t)=\widetilde{b}_{0}+c\sum_{i=1}^{n}a_{i}t_{i}$ of $\mathrm{ReLU}(A(t))$ is shown in Fig. 2 left panel.

Now, we have to bound the difference between $\widetilde{B}$ and $\mathrm{ReLU}\circ A$ on $[-1,1]^{n}$. By the choice of $c$ and $\widetilde{b}_{0}$ the maximal value of $\widetilde{B}(t)$ in the cube $[-1,1]^{n}$ is $\tau M$. Hence

$$ $D_{+}:=\max_{t\in I^{n}}\left(\mathrm{ReLU}(A(t))-\widetilde{B}(t)\right)=\max _{t\in I^{n}}\left(A(t)-\widetilde{B}(t)\right)=M-\tau M=M(1-\tau).$ (2) $$

The minimal value of this difference is achieved, when $(a_{0}+\sum a_{i}t_{i})=0$, that is $\sum a_{i}t_{i}=-a_{0}$ – see Fig. 2 (left panel). This minimal value is then

$$ $D_{-}=\min_{t\in I^{n}}\left(\mathrm{ReLU}(A(t))-\widetilde{B}(t)\right)=\min_ {t\in I^{n}}\left(0-\widetilde{B}(t)\right)=-\widetilde{b}_{0}+ca_{0}=ca_{0}- \frac{1}{2}\tau M.$ (3) $$

Gathering (2)-(3) we obtain

$$ $\left(\mathrm{ReLU}(A(t))-\widetilde{B}(t)\right)\in[D_{-},D_{+}],\qquad\text{ for }t\in[-1,1]^{n}.$ $$

The above considerations lead to an algorithm for computation of $\mathrm{ReLU}$ in the affine arithmetics. Given coefficients $(a_{0},\ldots,a_{n})$ of an affine function $A(t_{1},\ldots,t_{n})$ we compute coefficients $(b_{0},\ldots,b_{n+1})$ of an affine function $B(t_{1},\ldots,t_{n+1})$ so that the range of $B(t)$, $t\in[-1,1]^{n+1}$ covers the range of $\mathrm{ReLU}(A(t))$, $t\in[-1,1]^{n}$ in the following way

$$ $\displaystyle S$ $\displaystyle=$ $\displaystyle\sum_{i=1}^{n}|a_{i}|,\quad M=a_{0}+S,\quad c=\frac{1}{2}\tau M/S ,\quad D_{+}=M(1-\tau),\quad D_{-}=ca_{0}-\frac{1}{2}\tau M,$ $\displaystyle b_{0}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\left(\tau M+D_{+}+D_{-}\right),\quad b_{n+1}=\frac{1} {2}(D_{+}-D_{-}),\quad b_{i}=ca_{i},\quad i=1,\ldots,n.$ $$

There remains to explain how we choose the parameter $\tau\in[0,1]$. Set $U=\max_{t\in[-1,1]^{n}}A(t)=a_{0}+S$ and $L=\min_{t\in[-1,1]^{n}}A(t)=a_{0}-S$. Experimentally we have found that the choice $\tau\approx\frac{U}{U-L}=\frac{a_{0}+S}{2S}$ gives reasonable small overestimation of $\mathrm{ReLU}\circ A$ and it is very fast to compute.

Figure: Figure 4: The average maximal diameter of the NN output measured for points near the classification boundary. The X axis represents the perturbation size applied to the data points, while the Y axis shows the average maximal diameter of the NN output in the logarithmic scale. As we can see, the AA and DA methods give better approximation of interval bounds than the IBP method. Note that the DA cannot be calculated for large CNN architectures according to CPU constraints. We can see that IBP training in relation to standard training allows to reduce wrapping effect.
Refer to caption: extracted/5901384/img/mnist_the_most_uncertain_points_plot_with_doubletons_eps_0_01.png

The computation of softmax in affine arithmetics is presented in Appendix.

## 4 Experiments

In this section we present the results obtained by our two proposed methods: Affine and Doubleton Arithmetics. For a fixed point $x$ from the dataset, we define a box $B=x+\varepsilon[-\mathbbm{1},\mathbbm{1}]$ and then we compare bounds on the output of a neural network $\Phi(B)$ obtained by means of IBP, DA and AA methods. Additionally we compute Lower Bound (LB) on $\Phi(B)$ as the smallest box (interval hull) containing the set $\{\Phi(\xi_{k})\}_{k=1}^{1000}$, where $\{\xi_{k}\}_{k=1}^{1000}\subset B$ are randomly chosen points.
All the experiments are
implemented in C++ with the full control over numerical and rounding errors obtained due to the use of CAPD library (Kapela et al., ).

The results presented for the MNIST, CIFAR-10, and SVHN datasets are shown only for the small CNN architecture unless stated otherwise. The results for the Digits dataset are presented using an MLP architecture. For more details about the training hyperparameters, architectures, datasets, and hardware used – see Section C in Appendix.

Figure: Figure 5: The average maximal diameter of the NN output measured for points near the classification boundary for the medium and large CNN architectures. The X axis represents the perturbation size applied to the data points, while the Y axis shows the average maximal diameter of the NN output in the logarithmic scale.
Refer to caption: extracted/5901384/img/mnist_cnn_medium_the_most_uncertain_points_plot_eps_0_01.png

### Interval bounds for points sampled near the decision boundary

We compare our methods on points sampled near the decision boundary. We start by selecting one point from each class. From these points, we sample one point and connect it to the remaining points with line segments. For each of these line segments, we sample a point that lies near the decision boundary. We then calculate the interval neural network output for each of these selected points and average maximal diameters of these intervals to assess the model’s uncertainty. It is important to emphasize that these selected points may not be actual points from the real dataset, as they are, in fact, convex combinations of points from the real dataset. The experiments are conducted on the MNIST and CIFAR-10. We provide results for neural networks with weights obtained through a classical training procedure (without IBP training) and weights obtained through IBP training as well, for comparison.

As shown in Fig. 4, for a neural network trained using the IBP method ($\epsilon_{\text{train}}=0.01$) and standard training, the AA and DA methods perform significantly better compared to the IBP method. It is important to highlight that both the AA and DA methods produce nearly identical results. We would like to emphasize, that the AA method gives useful answer even for large perturbation size $\epsilon$, while the IBP even for not very large perturbations gives useless outputs of length $1$, which means that the probability is somewhere between 0 0 and $1$. Such phenomenon is well visible in each experiment we conducted – see Figs. 4, 5 and 6.

### Influence of network size on interval bounds

It is a fair question how interval bounds change depending on the size of a neural network architecture. We address this question by using medium, and large CNN architectures for the MNIST dataset. The architectures were trained using the IBP method with $\epsilon_{\text{train}}=0.01$, as well as without the IBP method for comparison.

For the medium and large CNN architectures trained on the MNIST dataset (Fig. 5), the AA method produces results close to those of the LB method, significantly outperforming the IBP method. These differences are particularly noticeable when IBP training is not applied. In this case, the AA method once again outperforms the IBP method, while the differences between the LB and AA methods are only slightly worse compared to the scenario when IBP training is used. It is also worth emphasizing that when medium and large CNN architectures are used without IBP training, the neural network becomes extremely uncertain about the investigated data points.

### Influence of various perturbation sizes used in IBP-based training on the resulting interval bounds

Figure: Figure 6: The average maximal diameter of the NN output measured for points near the classification boundary for network trained with different interval lengths $\epsilon_{\text{train}}$for the Digits dataset. The X axis represents the perturbation size applied to the data points, while the Y axis shows the average maximal diameter of the NN output in the logarithmic scale.
Refer to caption: extracted/5901384/img/digits_the_most_uncertain_points_plot_eps_0_0001.png

Generally, the presented plots in Fig. 6 show that the larger the perturbation size applied during IBP training, the smaller the difference between the results obtained using the IBP and AA/DA methods. However, it is important to emphasize that increasing the perturbation size during IBP training makes training a neural network more difficult, leading to challenges in achieving satisfactory accuracy. Therefore, our proposed methods offer a much easier way to reduce the wrapping effect, and they can be applied regardless of whether IBP training is used, making them both more practical and efficient for real-world applications.

## 5 Conclusion

This paper analyzes wrapping effect in a neural network. We show that for linear models, interval bounds can grow exponentially. Such effects have a strong influence on the IBP certification of neural networks. To solve such a problem, we propose adapting two models from strict numerical calculations: Doubleton and Affine Arithmetics. Both models give sharp bounds for linear transformations. The experimental section shows that Affine Arithmetic returns bounds close to optimal within reasonable computational time.

### Limitations

Doubleton Arithmetics provides near-optimal bounds, but the computational complexity is $O(n^{3})$, where $n$ is the largest dimension of the hidden layers. Even for the small CNN architecture on the CIFAR-10 and SVHN datasets, the computation time was unacceptably high.

## Appendix A Integral computations

We will show the following lemma, which
gives a detailed estimations for Lemma 2.1.

###### Lemma A.1 .

Let $V=(V_{1},\ldots,V_{n})$ be a random vector uniformly chosen from the unit sphere in $\mathbb{R}^{n}$. Let $R$ be a random variable given by

$$ $R=|V_{1}|+\ldots+|V_{n}|.$ $$

Then

$$ $ER=\frac{2n}{n-1}\frac{\Gamma(\frac{n}{2})}{\sqrt{\pi}\Gamma(\frac{n-1}{2})}, \,ER^{2}=1+\frac{2}{\pi}(n-\frac{1}{n-2}),\,VR=ER^{2}-(ER)^{2}.$ $$

Moreover, we have the asymptotics

$$ $ER=\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}+\frac{1}{2\sqrt{2\pi n}}+O(n^{-3/2}),\, VR=1+\frac{1}{\pi}+O(1/n).$ $$

###### Proof.

To calculate $ER$, we compute

$$ $ER=\frac{1}{S_{n-1}}\int_{x:\|x\|=1}|x|_{1}dS(x)=\frac{1}{S_{n-1}}\int_{x:\|x \|=1}x_{1}+\ldots+x_{n}dS(x)=$ $$

$$ $\frac{2^{n}}{S_{n-1}}\int_{x:\|x\|=1,x_{1},\ldots,x_{n}\geq 0}x_{1}+\ldots+x_{ n}dS(x)=\frac{n2^{n}}{S_{n-1}}\int_{x:\|x\|=1,x_{1},\ldots,x_{n}\geq 0}x_{1}dS (x)=$ $$

$$ $\frac{2n}{S_{n-1}}\int_{x:\|x\|=1,x_{1}\geq 0}x_{1}dS(x)=\frac{2n}{S_{n-1}} \int_{0}^{1}x_{1}S_{n-2}(\sqrt{1-x_{1}^{2}})^{n-2}\cdot\frac{1}{\sqrt{1-x_{1}^ {2}}}dx_{1}=$ $$

$$ $\frac{2nS_{n-2}}{S_{n-1}}\int_{0}^{1}x_{1}(1-x_{1}^{2})^{\frac{n-3}{2}}dx_{1}= \frac{nS_{n-2}}{S_{n-1}}\int_{0}^{1}t^{\frac{n-3}{2}}dt=$ $$

$$ $=\frac{2nS_{n-2}}{(n-1)S_{n-1}}=\frac{2n}{n-1}\frac{\Gamma(\frac{n}{2})}{\sqrt {\pi}\Gamma(\frac{n-1}{2})}$ $$

Finally we obtain the assymptotic expansion

$$ $\approx\frac{2}{\sqrt{\pi}}[\frac{\sqrt{n}}{\sqrt{2}}+\frac{1}{4\sqrt{2n}}+O(1 /n^{3/2})].=\frac{\sqrt{2}}{\sqrt{\pi}}\sqrt{n}+\frac{1}{2\sqrt{2\pi n}}+O(n^{ -3/2}).$ $$

We proceed to computation of $R^{2}$. We have

$$ $ER^{2}=\frac{1}{S_{n-1}}\int_{x:\|x\|=1}(|x|_{1}+\ldots+|x_{n}|)^{2}dS(x)$ $$

$$ $=\frac{1}{S_{n-1}}\int_{x:\|x\|=1}|x|_{1}^{2}+\ldots+|x_{n}|^{2}dS(x)+\frac{n( n-1)}{S_{n-1}}\int_{x:\|x\|=1}|x|_{1}|x_{2}|dS(x)=$ $$

$$ $=1+\frac{4n(n-1)}{S_{n-1}}\int_{x:\|x\|=1,x_{1},x_{2}\geq 0}x_{1}x_{2}dS(x)$ $$

Now

$$ $\int_{x:\|x\|=1,x_{1},x_{2}\geq 0}x_{1}x_{2}dS(x)=\int_{x_{1},x_{2}\geq 0,x_{1 }^{2}+x_{2}^{2}\leq 1}x_{1}x_{2}S_{n-3}\sqrt{1-(x_{1}^{2}+x_{2}^{2})}^{n-3} \frac{1}{\sqrt{1-(x_{1}^{2}+x_{2}^{2})}}dx_{1}dx_{2}$ $$

Now

$$ $\int_{x_{1},x_{2}\geq 0,x_{1}^{2}+x_{2}^{2}\leq 1}x_{1}x_{2}(1-(x_{1}^{2}+x_{2 }^{2}))^{n/2-2}dx_{1}dx_{2}=\int_{0}^{\pi/2}\int_{0}^{1}r\sin\phi r\cos\phi(1- r^{2})^{n/2-2}rdrd\phi$ $$

$$ $=\int_{0}^{\pi/2}\frac{\sin 2\phi}{2}d\phi\cdot\frac{1}{2}\int_{0}^{1}(1-t)t^{ n/2-2}dt=\frac{1}{4}\cdot(\frac{2}{n-2}-\frac{2}{n})=\frac{1}{n(n-2)}.$ $$

Finally

$$ $ER^{2}=1+4\frac{n-1}{n-2}\frac{S_{n-3}}{S_{n-1}}=1+4\frac{n-1}{n-2}\frac{n/2}{ \pi}=1+\frac{2}{\pi}(n-\frac{1}{n-2}).$ $$

Clearly $VR=ER^{2}-(ER)^{2}$, which trivially yields the asymptotic expansion

$$ $VR=ER^{2}-(ER)^{2}=1+\frac{1}{\pi}+O(1/n).$ $$

∎

###### Lemma A.1 .

###### Proof.

## Appendix B Doubleton and affine arithmetics

### Softmax in affine arithmetics

Assume we have an affine function $At=x+Lt$ defined on the cube $t=(t_{1},\ldots,t_{n})\in[-1,1]^{n}=:I^{n}$, with $x\in\mathbb{R}^{m}$, $L\in\mathbb{R}^{m\times n}$. Our goal is to find an affine map

$$ $Bt=\tilde{x}+\tilde{L}t$ $$

and a vector $e\in\mathbb{R}^{m}$, so that for $i=1,\ldots,m$ there holds

$$ $\max_{t\in I^{n}}|(\mathrm{softmax}(A(t))-B(t))_{i}|\leq e_{i}.$ $$

The vector $\tilde{x}$ and the matrix $\tilde{L}$ will be computed from first order Taylor expansion of $\mathrm{softmax}$. A bound on error term $e$ will be computed from second derivatives.

Recall, that for $z\in\mathbb{R}^{m}$

$$ $\mathrm{softmax}(z)=(s_{1},\ldots,s_{m}):=\left(\frac{\exp(z_{i})}{\sum_{j=1}^ {m}\exp(z_{j})},\ldots,\frac{\exp(z_{m})}{\sum_{j=1}^{m}\exp(z_{j})}\right).$ $$

In order to avoid numerical instabilities in evaluation of the above expression we take $R=\|z\|_{\infty}$ and compute $\mathrm{softmax}(z)=\mathrm{softmax}(z_{1}-R,\ldots,z_{m}-R)$.

It is well known that the Jacobian of $\mathrm{softmax}$ is given by

$$ $J(z):=D\mathrm{softmax}(z)=\begin{bmatrix}s_{1}(1-s_{1})&-s_{1}s_{2}&\ldots&-s _{1}s_{m}\\ -s_{1}s_{2}&s_{2}(1-s_{2})&\ldots&-s_{2}s_{m}\\ \vdots&\vdots&\ddots&\vdots\\ -s_{m}s_{1}&\cdots&-s_{m-1}s_{m}&s_{m}(1-s_{m})\end{bmatrix}.$ $$

Thus, we can compute a linear approximation of $\mathrm{softmax}$ by

$$ $B(t)=\tilde{x}+\tilde{L}t=\mathrm{softmax}(x)+\left(J(x)L\right)t.$ $$

In the above $\mathrm{softmax}(x)$ and $J(x)$ are evaluated at a single point and therefore neither dependency error nor wrapping effect is present.

The error term $e_{i}$ can be bounded using second order Taylor expansion. We would like to find a bound

$$ $\left|(I^{n})^{T}D^{2}g_{i}(I^{n})I^{n}\right|\leq e_{i},\quad i=1,\ldots,m,$ $$

where $g(t)=\mathrm{softmax}(x+Lt)$. Differentiation of $Dg(t)=J(x+Lt)L$ gives

$$ $\displaystyle D^{2}g_{i}(t)$ $\displaystyle=$ $\displaystyle L^{T}DJ_{i}(x+Lt)L.$ $$

There remains to derive formula for $DJ_{i}(z)$. Differentiation gives

$$ $\displaystyle\frac{\partial J_{ij}(z)}{\partial z_{c}}$ $\displaystyle=$ $\displaystyle\frac{\partial}{\partial z_{c}}\left(\delta_{ij}s_{i}-s_{i}s_{j}\right)$ $\displaystyle=$ $\displaystyle\left(\delta_{ijc}s_{i}-\delta_{ij}s_{i}s_{c}\right)-(\delta_{ic} s_{i}-s_{i}s_{c})s_{j}-s_{i}(\delta_{jc}s_{j}-s_{j}s_{c})$ $\displaystyle=$ $\displaystyle\delta_{ijc}s_{i}-\delta_{ij}s_{i}s_{c}-\delta_{ic}s_{i}s_{j}- \delta_{jc}s_{i}s_{j}+2s_{i}s_{j}s_{c}.$ $$

Evaluation of the above formula in interval arithmetics leads to a rough bound on the error term $e_{i}$. We will show, however, that increasing time complexity we can significantly reduce dependency problem in this expression.

### Evaluation of g i subscript 𝑔 𝑖 g_{i} italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and products g i ⁢ g c subscript 𝑔 𝑖 subscript 𝑔 𝑐 g_{i}g_{c} italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_g start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and g i ⁢ g c ⁢ g r subscript 𝑔 𝑖 subscript 𝑔 𝑐 subscript 𝑔 𝑟 g_{i}g_{c}g_{r} italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_g start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT italic_g start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT .

We have

$$ $g_{i}(t)=\frac{\exp\left(x_{i}+L_{i1}t_{1}+\ldots+L_{in}t_{n}\right)}{\sum_{j= 1}^{m}\exp\left(x_{j}+L_{j1}t_{1}+\ldots L_{jn}t_{n}\right)}$ (4) $$

Dependency in (4) can be reduced using equivalent formula

$$ $g_{i}(t)=\frac{\exp\left(y_{i}\right)}{\sum_{j=1}^{m}\exp\left(y_{j}+(L_{j1}-L _{i1})t_{1}+\ldots+(L_{jn}-L_{in})t_{n}\right)},$ (5) $$

where $R=\max_{i=1,\ldots,m}x_{i}$ and $y_{i}=x_{i}-R$.

Let us recall an important in this context property of interval arithmetics. It is well known that multiplication is not distributive, that is for intervals $a,b,c$ there holds $a(b+c)\subset ab+ac$. However, if all intervals are nonnegative then we have equality. Such situation appears in evaluation of the product

$$ $\displaystyle g_{i}(t)g_{c}(t)$ $\displaystyle=$ $\displaystyle\left(\frac{\exp\left(y_{i}\right)}{\sum_{j=1}^{m}\exp\left(y_{j} +\sum_{p=1}^{n}(L_{jp}-L_{ip})t_{p}\right)}\right)\left(\frac{\exp\left(y_{c} \right)}{\sum_{k=1}^{m}\exp\left(y_{k}+\sum_{p=1}^{n}(L_{kp}-L_{cp})t_{p} \right)}\right)$ $\displaystyle=$ $\displaystyle\frac{\exp\left(y_{i}+y_{c}\right)}{\sum_{j,k=1}^{m}\exp\left(y_{ j}+y_{k}+\sum_{p=1}^{n}(L_{jp}+L_{kp}-L_{ip}-L_{cp})t_{p}\right)}.$ $$

The above two expressions, when evaluated in interval arithmetics, may lead to different bounds (of course they can be intersected). Time complexity of the second evaluation is $O(M^{2}N)$ while direct evaluation (first expression) is of order $O(MN)$. However, $\mathrm{softmax}$ is applied to the output of last layer in a neural network, which is usually of low dimension and therefore this should not be a serious additional cost.

Similarly, we can evaluate the product of three functions as

$$ $g_{i}(t)g_{c}(t)g_{r}(t)=\frac{\exp\left(y_{i}+y_{c}+y_{r}\right)}{\sum_{j,k,s =1}^{m}\exp\left(y_{j}+y_{k}+y_{s}+\sum_{p=1}^{n}(L_{jp}+L_{kp}+L_{sp}-L_{ip}- L_{cp}-L_{rp})t_{p}\right)}$ $$

and intersect the result with direct multiplication of three intervals $g_{i}(t)$, $g_{c}(t)$ and $g_{r}(t)$.

## Appendix C Experimental setting

### Datasets

We use the following publicly available datasets: 1) MNIST dataset, consisting of 60,000 training and 10,000 testing $28\times 28$ pixel gray-scale images of 10 classes of digits; 2) CIFAR-10 dataset, consisting of 50,000 training and 10,000 testing $32\times 32$ colour images in 10 classes; 3) SVHN dataset, consisting of 600,000 $32\times 32$ pixel colour images of 10 classes of digits; 4) Digits dataset, consisting of 1797 $8\times 8$ pixel gray-scale images of 10 classes of digits.

### Architectures

We use three CNN architectures (small, medium and large) as defined in Table 1 in Gowal et al. (). Additionally, we consider an MLP architecture consisting of four hidden layers with 100 neurons per layer. A classification head is added on top of these layers.

### Training parameters

During training, we use the Adam optimizer Kingma & Ba () with the default configuration of $\beta_{1}=0.9$ and $\beta_{2}=0.999$, but with different learning rates ($lr$) across all datasets. We consistently use the ReLU activation function. Whenever a scheduler is mentioned, we apply the MultiStepLR scheduler with a default multiplicative learning rate decay factor set to 0.1. The scheduler steps are applied twice: once after $\frac{1}{3}$ of the total number of iterations and once after $\frac{2}{3}$ of the total number of iterations. Additionally, there is a parameter $\kappa$ scheduled over the entire training process as $\kappa_{i}=\max\{1-0.00005\cdot i,\kappa_{max}\}$, where $i$ denotes the current training iteration and $\kappa_{max}$ is set to 0.5. A perturbation value $\epsilon$ grows linearly from 0 at the beginning of training to the $\epsilon_{max}$ hyperparameter value at the midpoint of the total number of iterations. The considered $\epsilon_{max}$ values are from the set $\{0.0001,0.001,0.01,0.05,0.1\}$ and remain the same regardless of the architecture used. We use 10% of the training samples as the validation set.

- •
For the MNIST, SVHN, and CIFAR-10 datasets, we train small, medium, and large CNNs using the best set of hyperparameters identified in Gowal et al. (). We apply the same normalization and augmentation scheme. The only differences are in the epsilons $\epsilon$ used during training and the number of epochs. We decreased the total number of epochs for the CIFAR-10 and SVHN datasets to 100 for the large CNN.
- •
For Digits, we train the MLP for 50 epochs with batch sizes of 32. No normalization or augmentation is applied. The rest of the hyperparameters remain the same as for the MNIST, SVHN, and CIFAR-10 datasets.

### Hardware and software resources used

The implementation is done in Python 3.10.13, utilizing libraries such as PyTorch 2.3.1 with CUDA support, NumPy 1.26.4, Pandas 2.1.1, and others. Most computations are performed on an NVIDIA GeForce RTX 4090 GPU, with some training sessions also conducted on NVIDIA GeForce RTX 3080 and NVIDIA DGX GPUs. The experiments involving Affine and Doubleton Arithmetics were implemented using the CAPD library (Kapela et al., ).

## Appendix D Experimental results

### Interval bounds for partially masked data

In this subsection, we aim to present the interval bounds obtained through a neural network for data from the Digits and SVHN datasets. We sample 10 points, each belonging to a single class, and apply a mask where 50% of the values are masked (replaced by zero) and the remaining values stay unchanged. We then measure the average diameter of the neural network output.

Figure: Figure 7: The average maximal diameter of the NN output measured for points near the classification boundary in the case where parts of the images was masked. The X axis represents the perturbation size applied to the data points, while the Y axis shows the average maximal diameter of the NN output in the logarithmic scale.
Refer to caption: extracted/5901384/img/digits_masked_points_plot_eps_0_1.png

Even for partially masked data, the output interval bounds obtained using the AA and DA methods are very close to those of the LB methods (Fig.  7), significantly outperforming the IBP method.

These results indicate that the AA and DA methods, compared to the IBP method, effectively minimizes the wrapping effect in neural networks. Consequently, these methods can be regarded as a viable approach for quantifying the uncertainty of a neural network’s output when some pixels of an image are masked with a value of 0.

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: Could the authors give a clarification of proposition 2.1 in the context of IBP?
soundness: 2
strengths: 1. Precisely certifying the robustness of neural networks remains a well-motivated problem.

2. This paper is well-written.
summary: This study investigates methods to enhance the precision of Interval Bound Propagation (IBP) in the context of neural network certification, specifically under $\ell_\infty$-perturbations. The research demonstrates IBP's vulnerability to the wrapping effect and proposes the implementation of Dubleton Arithmetic and Affine Arithmetic as mitigation strategies.
weaknesses: 1. Overall, the setting is ad-hoc. The motivation of this work is the wrapping effect, specifically, the precision loss for a sequence randomly orthogonal maps is exponential. This is not a realistic setting. Moreover, the interpretation of Proposition 2.1 in the context of IBP is questionable.

2. Novelty is limited. The mitigation strategies, dubleton arithmetic and affine arithmetic, were already known. This work adapts these two known techniques to the problem.

3. The baseline is weak. The empirical evaluation only employs IBP as a baseline while overlooking state-of-the-art CROWN variants for $\ell_\infty$ certification.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: Do you have any suggestions for practical use?
soundness: 3
strengths: 1. Certification of pre-trained NNs is an important task for many different applications. The paper suggests novel theoretical approaches to achieve tighter bounds.

2. The intuition is well-explained, and the experiments are well-organized and thoroughly analyzed.
summary: The paper suggests a novel method, Doubleton and Affine Arithmetic, to control rounding errors in NN certification tasks by removing the wrapping effect in IBP. The theoretical findings are supported by experimental evidence.
weaknesses: 1. The method is only applicable for the linear transformation. However, there are various practical models with non-linear activation functions, which make its application is limited.

2. The time-complexity is too high for real applications such as ResNet and transformer-based models.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: Please see the weaknesses above.
soundness: 3
strengths: 1. Considering the approximation error in IBP from a different perspective is very interesting to me.
2. The experiment results are strong compared to the original IBP method.
summary: This paper first studies the wrapping effect of IBP, which causes large approximation errors in forward propagation. The authors then propose to apply two existing methods, Doubletons and Affine Arithmetic to mitigate the wrapping effect, inspired by a similar problem in ODE. Experiments show that the proposed method can achieve a much tighter bounds compared to original IBP.
weaknesses: 1. The conclusions about IBP wrapping effect has been discussed in many existing works such as [1], with a very similar theoretical result.
2. The proposed method looks very similar to the forward CROWN [2], which uses affine functions to bound the neural networks.

Therefore, the novelty of this paper is questionable given these relevant existing works.



[1] Fast certified robust training with short warmup. 
[2] Automatic Perturbation Analysis for Scalable Certified Robustness and Beyond.

---

confidence: 5
contribution: 1
initial_rating: 3
presentation: 2
questions: 1. In the definition of DA, why is the term $C\textbf{r}$ needed? Why can’t this term just be combined with $Q\textbf{q}$?
2. I would like to see a comparison between the proposed methods and different SOTA certification methods regarding certification power and time (both convex-relaxation-based and multi-neuron or branch-and-bound-based ones). 
3. A better lower bound might be obtained by using adversarial samples (i.e. generated by PGD adversarial attacks) instead of just randomly sampled points. Could the authors provide some comparisons in this regard?
4. Is there any significant improvement gained from using the CAPD C++ library for having full control over numerical and rounding errors? Without a proper discussion and/or experimental results, the computational overhead seems hard to justify.
soundness: 3
strengths: The paper is well-structured and easy to read. The proofs are (mostly) easy to follow. 

Using bounding methods from the field of ODE solvers in Neural Network Certification is an interesting direction.

The experiments are well organized and the results are presented and interpreted in insightful ways.
summary: The paper explores the performance of bound propagation methods in the context of neural network certification. The paper demonstrates theoretically and empirically that bounds obtained by IBP (Interval Bound Propagation) produce exponentially growing over-approximations due to the wrapping effect. The authors propose two alternatives to IBP, namely Dubleton Arithmetic (DA) and Affine Arithmetic (AA) that produce tighter and closer to optimal over-approximations when compared to IBP.
weaknesses: **Lack of discussion on relevant related work**

* The first claimed contribution of the paper (theoretical analysis of the wrapping effect in the case of IBP bounds) has been widely discussed in multiple recent works: Shi et al. (2021)[1] show that IBP bounds grow by a factor of $O(\sqrt n)$ after each affine layer, leading to exponential growth for deep networks. Various works [2,3,4,5,6] have proposed methods to compute tighter bounds for certification because IBP was not good enough. Other works [7,8,9] also acknowledge and experimentally verify the concept of exploding IBP bounds in the setting of Certified Training. Consequently, none of the theoretical results presented in section 2 of the paper can be considered novel.

* The Affine Arithmetic bound propagation method has been extensively used for Neural Network Certification in the literature under different names: Zonotope Abstraction [2], Symbolic Linear Relaxation [3], Fast-Lin [4], none of which have been discussed in the paper. Moreover, [2] is wrongly cited in Introduction (line 041) as a Certified Training method, when it is in fact a Certification method.

* Other works [5,6] focus on propagating linear constraints through the network. Most SOTA certification methods [10,11] combine this with further optimization such as multi-neuron analysis and branch-and-bound for analyzing unstable neurons.

* For the claim that their proposed methods are novel and better than existing methods, the authors should include discussions and justifications about the novelty and significance of their work when compared to previous work.

**Lack of experimental comparisons to related work**

* The authors only compare the results obtained by AA and DA with results obtained by IBP and with the Lower Bound obtained by sampling random points, but provide no comparisons with other existing certification methods. I recommend including some figures or tables with such comparisons.

* The paper presents plots of the maximal interval length for different networks and methods, but there are no results concerning the actual certification power of the methods (e.g. number of certified samples). I recommend including some figures or tables with such evaluations.

**Minor points**

* The notation mid() is not defined. Although a reader might infer from context that it refers to the middle point of the interval, it would be better to offer a formal definition or at least a reference to its definition from [12].

* Line 234: It is not clear if the $n$ considered for the dimensions of $\tilde Q$ and $A$ is the same as the dimension of the input space. In strategy 2 the matrices $\tilde Q$ and $A$ have dimensions $d \times d$ instead of $n \times d$ even if $d \neq n$.

* Figure 6 contains IBP training and Standard Training as titles for the two columns, but from the caption I understand there is only IBP training involved under different perturbation sizes.

* It would be nice to have a table with network architectures and hyperparameters used in the appendix instead of just referencing prior work. 

**References**

[1] Shi et al., Fast certified robust training with short warmup, NeurIPS 2021

[2] Singh et al., Fast and effective robustness certification, NeurIPS 2018

[3] Wang et al., Efficient Formal Safety Analysis of Neural Networks, NeurIPS 2018

[4] Weng et al., Towards fast computation of certified robustness for relu networks, ICML 2018

[5] Zhang et al., Efficient neural network robustness certification with general activation functions, NeurIPS 2018

[6] Singh et al., An abstract domain for certifying neural networks, POPL 2019

[7] Muller et al., Certified training: Small boxes are all you need, ICLR 2023

[8] De Palma et al., Expressive losses for verified robustness via convex combinations, ICLR 2024

[9] Mao et al., Understanding certified training with interval bound propagation, ICLR 2024

[10] Wang et al., Beta-CROWN: Efficient Bound Propagation with Per-neuron Split Constraints for Neural Network Robustness Verification, NeurIPS 2021

[11] Ferrari et al., Complete Verification via Multi-Neuron Relaxation Guided Branch-and-Bound, ICLR 2022

[12] Mrozek and Zgliczyński, Set arithmetic and the enclosing problem in dynamics, Annales Polonici Mathematici 2000