A Toolbox for Fast Interval Arithmetic in numpy with an Application to
  Formal Verification of Neural Network Controlled Systems
In this paper, we present a toolbox for interval analysis in numpy, with an
application to formal verification of neural network controlled systems. Using
the notion of natural inclusion functions, we systematically construct interval
bounds for a general class of mappings. The toolbox offers efficient
computation of natural inclusion functions using compiled C code, as well as a
familiar interface in numpy with its canonical features, such as n-dimensional
arrays, matrix/vector operations, and vectorization. We then use this toolbox
in formal verification of dynamical systems with neural network controllers,
through the composition of their inclusion functions.
Akash Harapanahalli, Saber Jafarpour, Samuel Coogan
A Toolbox for Fast Interval Arithmetic in numpy with an Application to Formal Verification of Neural Network Controlled Systems

###### Abstract

In this paper, we present a toolbox for interval analysis in numpy, with an application to formal verification of neural network controlled systems. Using the notion of natural inclusion functions, we systematically construct interval bounds for a general class of mappings. The toolbox offers efficient computation of natural inclusion functions using compiled C code, as well as a familiar interface in numpy with its canonical features, such as \(n\)-dimensional arrays, matrix/vector operations, and vectorization. We then use this toolbox in formal verification of dynamical systems with neural network controllers, through the composition of their inclusion functions.

Machine Learning, Neural Network, Neural Network, Neural Network, Neural Network, Neural Network

## 1 Introduction

Interval analysis is a classical field that provides a computationally efficient approach for propagating errors by computing function bounds (Jaulin et al., 2001). It has been successfully used for floating point error bounding in numerical and scientific analysis (Hickey et al., 2001). Interval bounds are often used for dynamical systems: (i) in reachability analysis, using methods such as Differential Inequalities (Scott and Barton, 2013) and Mixed Monotonicity (Meyer et al., 2019; Abate et al., 2021); (ii) for invariant set computation (Abate and Coogan, 2020). Recently, interval analysis has been increasingly used for the verification of learning algorithms: (i) standalone neural network verification approaches such as Interval Bound Propagation (Gowal et al., 2019); (ii) in-the-loop neural network control system verification approaches including simulation-guided approaches (Xiang et al., 2021), POLAR (Huang et al., 2022) and ReachMM (Jafarpour et al., 2023).

Since all of these techniques use a similar suite of tools, there is value in creating an efficient, user-friendly toolbox for general interval analysis. Python has become the standard for the learning community, and as such there are several existing tools for interval arithmetic. However, they come with key drawbacks: pyinterval does not natively support interval vectors and matrices, and portion supports lists of intervals, but not matrix and vector operations.

ContributionsIn this paper, we introduce a novel interval analysis framework called npinterval1, implemented in numpy(Harris et al., 2020), the computational backbone of most scientific Python packages. This framework is built upon the notion of inclusion functions, which provide interval bounds on the output of a given function. We first define tight inclusion functions for several elementary functions, then use Theorem 2.3 to build natural inclusion functions for a more general class of composed functions. The proposed package extends the prominent benefits of numpy directly to interval analysis, including its efficiency with compiled C implementations, versatility with \(n\)-dimensional arrays, matrix/vector operations, vectorization, and its familiar user interface. We then demonstrate its utility through an application in formal verification of neural network controlled systems, by composing CROWN (Zhang et al., 2018), a state-of-the-art neural network verification method with a natural inclusion functions of the system in Theorem 3.4. The proofs of all the Theorems are presented in Appendix B.

Footnote 1: The most recent code for npinterval can be viewed at [https://github.com/gtfactslab/npinterval](https://github.com/gtfactslab/npinterval); to reproduce the figures in this paper, see [https://github.com/gtfactslab/Harapanahalli_WFVML2023](https://github.com/gtfactslab/Harapanahalli_WFVML2023).

NotationWe denote the standard partial order on \(\mathbb{R}^{n}\) by \(\leq\), i.e., for \(x,y\in\mathbb{R}^{n}\), \(x\leq y\) if and only if \(x_{i}\leq y_{i}\) for all \(i\in\{1,\ldots,n\}\). A (bounded) _interval_ of \(\mathbb{R}^{n}\) is a set of form \(\{z:\underline{x}\leq z\leq\overline{x}\}=:[\underline{x},\overline{x}]\) for some endpoints \(\underline{x},\overline{x}\in\mathbb{R}^{n}\), \(\underline{x}\leq\overline{x}\). Let \(\mathbb{IR}^{n}\) denote the set of all intervals on \(\mathbb{R}^{n}\). We also use the notation \([x]\in\mathbb{IR}^{n}\) to denote an interval when its endpoints are not relevant or implicitly understood to be \(\underline{x}\) and \(\overline{x}\). For every two vectors \(v,w\in\mathbb{R}^{n}\) and every \(i\in\{1,\ldots,n\}\), we define the vector \(v_{\{i:w\}}\in\mathbb{R}^{n}\) by \(\left(v_{\{i:w\}}\right)_{j}=\left\{\begin{matrix}v_{j}&j\neq i\\ w_{j}&j=i\end{matrix}\right..\) For a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) and a set \(\mathcal{X}\subseteq\mathbb{R}^{n}\), define the set-valued extension
\(f(\mathcal{X})=\{f(x):x\in\mathcal{X}\}\). For two vectors \(x,y\in\mathbb{R}^{n}\), let \((x,y)\in\mathbb{R}^{2n}\) denote their concatenation.

## 2 Interval Analysis

### Interval Arithmetic

Interval analysis extends operations and functions to intervals (Jaulin et al., 2001). For example, if we know that some number \(a\in[\underline{a},\overline{a}]\), and \(b\in[\underline{b},\overline{b}]\), it is easy to see that the sum \((a+b)\in[\underline{a}+\underline{b},\overline{a}+\overline{b}]\).

**Definition 2.1** (Inclusion Function (Jaulin et al., 2001)).: Given a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\), the interval function \([f]:\mathbb{IR}^{n}\rightarrow\mathbb{IR}^{m}\) is called an

1. _inclusion function_ for \(f\) if, for every \([x]\in\mathbb{IR}^{n}\), \(f([x])\subseteq[f]([x])\);
2. \([y]\)_-localized inclusion function_ for \(f\) if for every \([x]\subseteq[y]\), we have \(f([x])\subseteq[f]([x])\).

Moreover, an inclusion function \([f]\) for \(f\) is

1. _monotone_ if \([x]\subseteq[y]\) implies that \([f]([x])\subseteq[f]([y])\).
2. _tight_ if, for every \([x]\), \([f]([x])\) is the smallest interval containing \(f([x])\).

In the next Theorem, we provide a closed-form expression for the tight inclusion function.

**Theorem 2.2** (Uniqueness and Monotonicity of the Tight Inclusion Function).: _Given a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\), the tight inclusion function can be characterized as_

\[[f]([x])=\bigg{[}\inf_{x\in[x]}f(x),\sup_{x\in[x]}f(x)\bigg{]},\]

_where the \(\inf\) and \(\sup\) are taken element-wise, which is unique and monotone._

For some common functions, the tight inclusion function is easily defined. For example, if a function is monotonic, the tight inclusion function is simply the interval created by the function evaluated at its endpoints.

More generally, the tight inclusion function can be alternatively computed using the fact that on a closed and bounded interval, a continuous function will achieve its maximal and minimal values at either the endpoints or a critical value in the interval. For any bounded interval, the tight inclusion function can be evaluated by taking the maximum and minimum of the function on each critical point within and the endpoints of the input interval. Tight inclusion functions for some elementary functions such as \(\sin\) and \(\cos\) can be defined in this manner (see Table 1). However, when considering general functions, finding the tight inclusion function is often not computationally viable. The following theorem shows a more computational approach, by chaining known inclusion functions.

**Theorem 2.3** (Natural Inclusion Functions).: _Given a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) defined by a composition of functions with known monotone inclusion functions, i.e., \(f=e_{\ell}\circ e_{\ell-1}\circ\cdots\circ e_{1}\), an inclusion function for \(f\) is formed by replacing each composite function with its inclusion function, i.e. \([f]=[e_{\ell}]\circ[e_{\ell-1}]\circ\cdots\circ[e_{1}]\) and is called the natural inclusion function. Additionally, if each of the \([e_{j}]\) are monotone inclusion functions, the natural inclusion function is also monotone._

Note that two different decompositions of the function \(f\) can lead to two different natural inclusion functions for \(f\). Thus, the natural inclusion function is not guaranteed to be the tight inclusion function. For example, consider the function

\[(x+1)^{2}=x^{2}+2x+1,\]

on the interval \([-1,1]\). With the natural inclusion function for the first expression (LHS), the output interval is \(([-1,1]+1)^{2}=[0,4]\). With the natural inclusion function for the second expression (RHS), the output interval is \([-1,1]^{2}+2*[-1,1]+1=[0,1]+[-2,2]+1=[-1,4]\). Figure 1 demonstrates this phenomenon in further detail.

### Automated Interval Analysis using npinterval

The main contribution of this paper is to introduce the open source npinterval package, an extension of numpy to allow native support for interval arithmetic. npinterval defines a new interval data-type, internally represented as a tuple of two doubles, \([a]=(a.l,a.u)\). The interval type is fully implemented in C, and therefore the standard operations from Table 1 are all compiled into efficient machine-code executed as needed at runtime. Additionally, since interval is implemented as a dtype in numpy, all of numpy's prominent features, including \(n\)-dimensional arrays, fast matrix multiplication, and vectorization, are available to use. In particular, each function from Table 1 is registered as a numpy universal function, allowing for quick, element-wise operation over an \(n\)-dimensional array. While there are existing interval arithmetic toolboxes, none plug directly into numpy, opting instead to rewrite every operation in Python. While these packages do support the same standard operations from Table 1, they lose the flexibility and utility of numpy, as well as the efficiency of compiled C code.

## 3 Interval Reachability of Neural Network Controlled Systems

One application of interval analysis is in reachability analysis of neural network controlled systems. This section
revisits and extends the framework considered in (Jafarpour et al., 2023).

### Problem Statement

Consider a dynamical system of the following form

\[\dot{x}=f(x,u,w), \tag{1}\]

where \(x\in\mathbb{R}^{n}\) is the system state, \(u\in\mathbb{R}^{p}\) is the control input, \(w\in\mathcal{W}\subseteq\mathbb{R}^{q}\) is a unknown disturbance input in a compact set \(\mathcal{W}\), and \(f:\mathbb{R}^{n}\times\mathbb{R}^{p}\times\mathbb{R}^{q}\to\mathbb{R}^{n}\) is a parameterized vector field. We assume that a feedback control policy for the system (1) is given by a \(k\)-layer fully connected feed-forward neural network \(N:\mathbb{R}^{n}\to\mathbb{R}^{p}\) as follows:

\[\begin{split}\xi^{(i)}&=\sigma^{(i-1)}\left(W^{(i-1 )}\xi^{(i-1)}+b^{(i-1)}\right),\,i=1,\ldots k\\ &\xi^{(0)}=x,\quad N(x)=W^{(k)}\xi^{(k)}+b^{(k)}\end{split} \tag{2}\]

where \(m_{i}\) is the number of neurons in the \(i\)-th layer, \(W^{(i)}\in\mathbb{R}^{m_{i}\times m_{i-1}}\) is the weight matrix on the \(i\)-th layer, \(b^{(i)}\in\mathbb{R}^{m_{i}}\) is the bias vector on the \(i\)-th layer, \(\xi^{(i)}\in\mathbb{R}^{m_{i}}\) is the \(i\)-th layer hidden variable and \(\sigma_{i}\) is the activation function for the \(i\)-th layer. Thus, we consider the closed-loop system

\[\dot{x}=f(x,N(x),w)=f^{c}(x,w). \tag{3}\]

Given an initial time \(t_{0}\), an initial state \(x_{0}\), and a piecewise continuous mapping \(\mathbf{w}:\mathbb{R}\to\mathcal{W}\), denote the trajectory of the system for any \(t\geq t_{0}\) as \(\phi_{f^{c}}(t,t_{0},x_{0},\mathbf{w})\). Given an initial set \(\mathcal{X}_{0}\), we denote the reachable set of \(f^{c}\) at some \(t\geq t_{0}\):

\[\mathcal{R}_{f}(t,t_{0},\mathcal{X},\mathcal{W})=\left\{\begin{split}&\phi_{f^{c}}(t,t_{0},x_{0}, \mathbf{w}),\,\forall x_{0}\in\mathcal{X}_{0},\\ & w:\mathbb{R}\to\mathcal{W}\text{ piecewise cont.}\end{split}\right\} \tag{4}\]

One can use the reachable set of the system to verify safety specifications, e.g. by ensuring an empty intersection with unsafe states for all time. However, in general, computing the reachable set exactly is not computationally tractable--instead, approaches typically compute an over-approximation \(\overline{\mathcal{R}}_{f^{c}}(t,t_{0},\mathcal{X},\mathcal{W})\supseteq \mathcal{R}_{f^{c}}(t,t_{0},\mathcal{X},\mathcal{W})\). The main challenge addressed in this section is to develop an approach for providing tight over-approximations of reachable sets while remaining computationally tractable for runtime computation.

### Open-loop System Interval Reachability

Previously, (Jafarpour et al., 2023) consider a known decomposition function of the open-loop system. The interval analysis framework from Section 2 allows us to extend this theory to remove the need to define a decomposition function _a priori_.

**Assumption 3.1**.: For the dynamical system (1), there exists a known monotone inclusion function \([f]\) for \(f\).

Using Theorem 2.3, this assumption reduces to knowing a particular form \(f=e_{\ell}\circ\cdots\circ e_{1}\) with known monotone inclusion functions for each \(e_{j}\), thus removing the need to manually define a decomposition function for \(f\).

The open-loop embedding function is defined by:

\[\begin{split}&\underline{\mathsf{E}}_{i}([x],[u],[w])=\underline{f }_{i}([\underline{x},\overline{x}_{(i:\underline{x})}],[u],[w]),\\ &\overline{\mathsf{F}}_{i}([x],[u],[w])=\overline{f}_{i}([ \underline{x}_{(i:\underline{x})},\overline{x}],[u],[w]),\end{split} \tag{5}\]

for every \(i\in\{1,\ldots,n\}\), where \([f]=[f,\overline{f}]\) is the monotone inclusion function of \(f\), and \(\underline{\mathsf{E}},\overline{\mathsf{F}}:\mathbb{IR}^{n}\times\mathbb{IR} ^{p}\times\mathbb{IR}^{q}\to\mathbb{R}^{n}\). Using this embedding function, the following embedding dynamics can be defined

\[\begin{bmatrix}\dot{x}\\ \dot{\widehat{x}}\end{bmatrix}=\begin{bmatrix}\underline{\mathsf{E}}([x, \widehat{x}],[u],[w])\\ \overline{\mathsf{F}}([x,\widehat{x}],[u],[w])\end{bmatrix}, \tag{6}\]

where \((x,\widehat{x})\in\mathbb{R}^{2n}\), \(x\leq\widehat{x}\).

**Theorem 3.2** (Open-Loop System Interval Reachability).: _Let \(t\mapsto(\underline{x}(t),\overline{x}(t))\) denote the trajectory of the system (6) with initial condition \((\underline{x}_{0},\overline{x}_{0})\), under interval control mapping \([\mathbf{u}]:\mathbb{R}\to\mathbb{IR}^{p}\) and interval disturbance mapping \([\mathbf{w}]:\mathbb{R}\to\mathbb{IR}^{q}\). Let \(t\mapsto x(t)\) denote the trajectory of the system (1) with initial condition \(x_{0}\), under control \(\mathbf{u}:\mathbb{R}\to\mathbb{R}^{p}\) and disturbance \(\mathbf{w}:\mathbb{R}\to\mathcal{W}\). If \(x_{0}\in[\underline{x}_{0},\overline{x}_{0}]\), \(\mathbf{u}(t)\in[\mathbf{u}](t)\), and \(\mathbf{w}(t)\in[\mathbf{w}](t)\) for every \(t\geq t_{0}\), then_

\[x(t)\in[\underline{x}(t),\overline{x}(t)],\quad\text{for every $t\geq t_{0}$.}\]

Figure 1: **Left:**npinterval is used to generate interval approximations for a function \(f\) using two different natural inclusion functions. Blue: \(f(x_{1},x_{2})=[(x_{1}+x_{2})^{2},4\sin((x_{1}-x_{2})/4)]^{T}\) Green: \(f(x_{1},x_{2})=[x_{2}^{2}+2x_{1}x_{2}+x_{2}^{2},4\sin(x_{1}/4)\cos(x_{2}/4)-4 \cos(x_{1}/4)\sin(x_{2}/4)]^{T}\). The approximations are generated using the initial set \([-1,1]\times[-1,1]\), and 2000 uniformly sampled outputs are shown in red. **Right:** The same function is analyzed, with the same two natural inclusion functions, but the initial set is partitioned into 1024 uniform sections, and the union of the interval approximations are shown.
### Interconnected Closed-Loop System Interval Reachability

While Theorem 3.2 provides a method of over-approximating the system (3) using global bounds on the control input, it disregards all interactions between the neural network controller and the dynamical system.

**Assumption 3.3**.: For the neural network (2), there exists a known monotone inclusion function \([N]\) for \(N\).

For example, one can use CROWN (Zhang et al., 2018) to obtain these bounds. Using CROWN, given an interval \([y]\), we can obtain affine upper and lower bounds of the following form

\[\underline{C}_{[y]}x+\underline{d}_{[y]}\leq N(x)\leq\overline{C}_{[y]}x+ \overline{d}_{[y]}, \tag{7}\]

valid for any \(x\in[y]\), which can be used to create the following monotone \([y]\)-localized inclusion function \([N]_{[y]}=[\underline{N}_{[y]},\overline{N}_{[y]}]\), with

\[\underline{N}_{[y]}([x]) =\underline{C}_{[y]}^{+}\underline{x}+\underline{C}_{[y]}^{-} \overline{x}+\underline{d}_{[y]}, \tag{8}\] \[\overline{N}_{[y]}([x]) =\overline{C}_{[y]}^{+}\overline{x}+\overline{C}_{[y]}^{-} \underline{x}+\overline{d}_{[y]},\]

valid for any \([x]\subseteq[y]\). One can then construct the following "hybrid" closed-loop embedding function by interconnecting the open-loop embedding function with the monotone inclusion function for the neural network as follows

\[\underline{\mathsf{E}}_{i}^{c}([x],[w]) =\underline{f}_{i}([\underline{x},\overline{x}_{\{i:\underline{ x}\}}],[N]_{[x]}([\underline{x},\overline{x}_{\{i:\underline{x}\}}]),[w]), \tag{9}\] \[\overline{\mathsf{F}}_{i}^{c}([x],[w]) =\overline{f}_{i}([\underline{x}_{\{i:\overline{x}\}},\overline {x}],[N]_{[x]}([\underline{x}_{\{i:\overline{x}\}},\overline{x}]),[w]),\]

for every \(i\in\{1,\ldots,n\}\), and \(\underline{\mathsf{E}}^{c},\overline{\mathsf{F}}^{c}:\mathbbm{R}^{n}\times \mathbbm{R}^{q}\rightarrow\mathbb{R}^{n}\). Using this embedding function, the following embedding dynamics can be defined

\[\begin{bmatrix}\dot{x}\\ \dot{\dot{x}}\end{bmatrix}=\begin{bmatrix}\underline{\mathsf{F}}^{c}([x, \widehat{x}],[w])\\ \overline{\mathsf{F}}^{c}([x,\widehat{x}],[w])\end{bmatrix}, \tag{10}\]

where \((x,\widehat{x})\in\mathbb{R}^{2n}\), \(x\leq\widehat{x}\).

**Theorem 3.4** (Closed-Loop System Interval Reachability).: _Let \(t\mapsto(\underline{x}(t),\overline{x}(t))\) denote the trajectory of the system (10) with initial condition \((\underline{x}_{0},\overline{x}_{0})\), under interval disturbance mapping \([\mathbf{w}]:\mathbb{R}\rightarrow\mathbbm{R}^{q}\). Let \(t\mapsto x(t)\) denote the trajectory of the closed-loop system (3) with initial condition \(x_{0}\), under disturbance \(\mathbf{w}:\mathbb{R}\rightarrow\mathcal{W}\). If \(x_{0}\in[\underline{x}_{0},\overline{x}_{0}]\) and \(\mathbf{w}(t)\in[\mathbf{w}](t)\) for every \(t\geq t_{0}\), then_

\[x(t)\in[\underline{x}(t),\overline{x}(t)],\quad\text{ for every }t\geq t_{0}.\]

### Experiments

Vehicle ModelConsider the nonlinear dynamics of a vehicle adopted from (Polack et al., 2017):

\[\dot{p_{x}} =v\cos(\phi+\beta(u_{2})) \dot{\phi} =\frac{v}{\ell_{r}}\sin(\beta(u_{2})) \tag{11}\] \[\dot{p_{y}} =v\sin(\phi+\beta(u_{2})) \dot{v} =u_{1}. \tag{12}\]

where \([p_{x},p_{y}]^{\top}\in\mathbb{R}^{2}\) is the displacement of the center of mass, \(\phi\in[-\pi,\pi)\) is the heading angle in the plane, and \(v\in\mathbb{R}_{\geq 0}\) is the speed of the center of mass. Control input \(u_{1}\) is the applied force, input \(u_{2}\) is the angle of the front wheels, and \(\beta(u_{2})=\arctan\left(\frac{\ell_{f}}{\ell_{r}+\ell_{r}}\tan(u_{2})\right)\) is the slip slide angle. Let \(x=[p_{x},p_{y},\phi,v]^{\top}\) and \(u=[u_{1},u_{2}]^{\top}\). We use the neural network controller (\(4\times 100\times 100\times 2\) ReLU) defined in (Jafarpour et al., 2023), which is applied at evenly spaced intervals of \(0.25\) seconds apart. This neural network is trained to mimic an MPC that stabilizes the vehicle to the origin while avoiding a circular obstacle centered at \((4,4)\) with a radius of \(2\).

The natural inclusion function is constructed using npinterval with Table 1, and the monotone inclusion function (8) is computed using autoLiRPA (Xu et al., 2020). The closed-loop embedding function (10) is then used to over-approximate the reachable set of the system using Theorem 3.4. The dynamics are simulated using Euler integration with a step size of \(0.05\). The results are shown in Figure 2.

## 4 Conclusion

In this paper, we introduced a framework for interval analysis implemented directly in numpy, called npinterval. The framework provides an automatic way to generate provide interval bounds on the output of a general class of functions. We use this framework to formally verify the output of nonlinear neural network controlled systems. In the future, npinterval will be updated to account for floating point error, as well as to support a wider array of standard inclusion functions.

Figure 2: The over-approximated reachable set of the nonlinear vehicle model in the \(p_{x}\)-\(p_{y}\) coordinates are shown in blue for the initial set \([7.95,8.05]^{2}\times[-\frac{2\pi}{3}-0.005,-\frac{2\pi}{3}+0.005]\times[1.995,2.005]\) over the time interval \([0,1.25]\). 100 true trajectories of the system are shown in red, and the average runtime and standard deviation over 100 runs is shown.
## Acknowledgements

This work was supported in part by the National Science Foundation under grant #2219755 and by the Ford Motor Company.
