# Vertex-based reachability analysis for verifying ReLU deep neural networks

João<EMAIL_ADDRESS>Eduardo
<EMAIL_ADDRESS>Eric
<EMAIL_ADDRESS>
(January, 2023)

###### Abstract

Neural networks achieved high performance over different tasks, i.e. image
identification, voice recognition and other applications. Despite their
success, these models are still vulnerable regarding small perturbations,
which can be used to craft the so-called adversarial examples. Different
approaches have been proposed to circumvent their vulnerability, including
formal verification systems, which employ a variety of techniques, including
reachability, optimization and search procedures, to verify that the model
satisfies some property. In this paper we propose three novel reachability
algorithms for verifying deep neural networks with ReLU activations. The first
and third algorithms compute an over-approximation for the reachable set,
whereas the second one computes the exact reachable set. Differently from
previously proposed approaches, our algorithms take as input a V-polytope. Our
experiments on the ACAS Xu problem show that the Exact Polytope Network
Mapping (EPNM) reachability algorithm proposed in this work surpass the state-
of-the-art results from the literature, specially in relation to other
reachability methods.

##  1 Introduction

Regardless of the success of deep neural networks in computer vision and
natural language processing, these models are susceptible to small
perturbations applied to their inputs, i.e. it is possible to misguide the
model output by applying a designed perturbation to a given input. For
instance, the ACAS Xu model [1] (explained later in detail), that responded
differently from expected while under specific circumstances. The inputs
purposefully designed to force a misbehavior of the neural network are denoted
as adversarial examples [2].

To overcome such vulnerability many different approaches have been previously
employed. One proposed the application of algorithms that were able to
generate adversarial examples, the so called adversarial attacks, and
subsequently applied these inputs in the training process of the network [2,
3, 4, 5]. There were also those approaches that aim to identify the
adversarial examples before feeding them as input to the neural network [6, 7,
8, 9]. Even though these procedures helped to reduce the vulnerability of the
neural networks, these models remained vulnerable to adversarial attacks.

Formal methods were also applied to certify or guarantee that the model
behaves as expected under some circumstances or within a specified domain
region, nevertheless [10] showed that the verification problem is NP-hard,
leaving the process of certifying large models still an open problem.

The existing formal procedures can be classified into three different
categories: 1) reachability methods; 2) optimization methods; and 3) search
methods. The first one relies on calculating the output mapping of an input
set [11, 12, 13], the second one comprises the application of mathematical
optimization (Mixed Integer-Linear Programming or Convex Optimization) to
identify counter-examples [14, 15, 16, 17], and the third makes use of both
reachability and optimization approaches in conjunction with search methods
for identifying counter-examples [10, 18].

In this paper, we propose three novel reachability algorithms: APNM and PAPNM
algorithms that compute an over-approximation for the output, while EPNM which
computes the exact mapping. We present demonstrations on the behavior and
correctness of these algorithms and case studies of their applications for
comparison with existing algorithms from the literature. We also show that the
algorithms proposed in this work are highly parallelizable.

The rest of this paper is organized into five sections: Section 2 gives an
overview of the existing algorithms and related works; Section 3 describes the
proposed algorithms; Section 4 addresses the demonstrations regarding the
completeness and soundness of the proposed algorithms; Section 5 presents a
study case for application and comparison; and finally Section 6 discusses the
outcome of the procedures presented in this work.

##  2 Related Works

Formal verification of neural networks is receiving a huge amount of attention
mainly because of its importance in security sensitive tasks such as the
application of neural networks in autonomous vehicles, systems controllers,
aeronautics, and several other applications that can possibly involve
financial, human or environmental injury [18, 1].

As previously presented, there are three main research areas developed to
provide guarantees for neural networks: 1) reachability methods; 2)
optimization methods; and 3) search methods. In this work, we propose three
novel approaches regarding reachability analysis, which consist of two major
steps: computing the output set by mapping a subset of the neural network
domain and comparing the output set with the specification.

One of the algorithms developed in the literature to verify neural networks by
means of reachability analysis is denoted as MaxSens [12], which computes an
over-approximation to the output reachable set and compares it to the desired
specification. As this algorithm computes an over-approximation it does not
satisfy the completeness property of a verification algorithm, although it is
sound (if the property is verified then the algorithm guarantees that it is
satisfied). To compute such an over-approximation, MaxSens divides the input
set into several smaller hyperrectangles and computes the maximum sensitivity
for each of them layer-by-layer. The output of this algorithm consists of the
union of several hyperrectangles, which approximates the exact reachable set.
Note that their approach can be applied to neural networks with any activation
function.

Another approach from the literature, denoted as ExactReach [11], computes the
exact reachable set given an input set. This procedure takes as inputs a
convex H-polytope (a polytope represented by its inequalities) and computes
the exact mapping layer-by-layer. The authors proposed this approach
specifically for the ReLU activation function, which is reasonable as this
function has achieved promising results for convolutional and fully connected
neural networks, which are widely applied in different applications. Due to
the nature of the ReLU activation function, the authors separated the non-
linear mapping process in three cases: 1) all the elements of the input are
positive; 2) all the elements of the input are negative; and 3) the input has
positive, negative or null elements. These cases cover all the mapping
possibilities regarding the ReLU activation.

Similarly to the aforementioned approaches, we propose in this paper three
novel algorithms for reachability analysis. However, instead of using the
H-polytope representation, each of our proposed procedures take as input set a
V-polytope (a polytope represented by its vertices). We demonstrate the
correctness of both algorithms and compare them with the literature (not only
with those that make use of reachability analysis) by using a case study (ACAS
Xu [1]).

##  3 Algorithms

In this section we will state the verification problem and the three
algorithms proposed in this work. The demonstrations regarding the correctness
of each procedure will be presented in the following section.

###  3.1 Problem statement

Let $F:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$ represent a mapping given by a
neural network composed of $L$ layers. Then, for a given
$\mathbf{x}\in\mathbb{R}^{n}$, $F(\mathbf{x})=(F_{L}\circ
F_{L-1}\circ\cdots\circ F_{2}\circ F_{1})(\mathbf{x})$, where $F_{l}$ is the
mapping given by the $l$-th layer. Further, $F_{l}$ consists of the
composition of an affine mapping and a non-linear mapping (in this case a ReLU
function):
$F_{l}(\mathbf{x}_{l})=ReLU(\mathbf{W}_{l}\mathbf{x}_{l}+\boldsymbol{\theta}_{l})$,
where $\mathbf{W}_{l}$ and $\boldsymbol{\theta}_{l}$ denote the weight matrix
and bias vector of the $l$-th layer, respectively, and $\mathbf{x}_{l}$ is the
input to the same layer.

Suppose that $\mathcal{X}\subseteq\mathbb{R}^{n}$ is the input set that we
want to verify and
$\mathcal{R}=\\{F(\mathbf{x})\mid\mathbf{x}\in\mathcal{X}\\}$ is the exact
output set associated with $\mathcal{X}$, which resulted from the application
of the neural network $F$ to each input in $\mathcal{X}$. Moreover, the set
$\mathcal{Y}$ comprises the expected output set for the inputs from
$\mathcal{X}$. Then, the verification problem consists of assuring that:

$\mathcal{R}\cap\lnot\mathcal{Y}=\emptyset.$ (1)

In other words, we want to guarantee that there will not exist an input of $F$
in $\mathcal{X}$ that will cause the network to generate an undesired output
(an output that is not in $\mathcal{Y}$).

To perform such a verification, one needs to start by calculating the output
set associated with $\mathcal{X}$. As we presented in previous sections, there
are approaches that compute the exact reachable set and other approaches that
over-approximate it. For those that compute an approximation for the output
set, denoted by $\widehat{\mathcal{R}}$, we will have that
$\mathcal{R}\subseteq\widehat{\mathcal{R}}$. Then, we can see that if
$\widehat{\mathcal{R}}\cap\lnot\mathcal{Y}=\emptyset$, then the property given
by the Equation 1 will still be assured. In the following sections we will
present the proposed algorithms.

###  3.2 Approximate Polytope Network Mapping (APNM)

The first algorithm proposed in this work is called Approximate Polytope
Network Mapping (APNM) and consists of a procedure to compute an over-
approximation for the reachable set. Let $\mathcal{P}$ be a convex closed
polytope, defined as a convex combination of its vertices, namely

$\mathcal{P}=\left\\{\sum_{i=1}^{o}\lambda_{i}\mathbf{v}_{i}\mathrel{\bigg{|}}\sum_{i=1}^{o}\lambda_{i}=1,\lambda_{i}\geq
0\text{ and }\mathbf{v}_{i}\in\mathcal{V},\forall
i\in\\{1,\ldots,o\\}\right\\},$ (2)

where $\mathcal{V}$ is the set of vertices of $\mathcal{P}$
($\mathcal{V}=\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{o}\\}$). To achieve its
goal, the algorithm has five parts:

  1. 1.

Affine map of the vertices with weights and biases;

  2. 2.

Adjacent vertex identification;

  3. 3.

Polytope intersection with an orthant’s hyperplanes;

  4. 4.

ReLU mapping; and

  5. 5.

Removing non-vertices.

In what follows, the way the algorithm works will be explained with a simple
$2$-dimensional problem for visualization purposes. The input set
$\mathcal{P}$ is presented in Figure 1: the set of vertices is
$\mathcal{V}=\\{\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3},\mathbf{v}_{4}\\}$,
where $\mathbf{v}_{1}=(1.0,1.0)$, $\mathbf{v}_{2}=(-1.0,1.0)$,
$\mathbf{v}_{3}=(-1.0,-1.0)$ and $\mathbf{v}_{4}=(1.0,-1.0)$; and the set of
edges is
$\mathcal{E}=\\{(\mathbf{v}_{1},\mathbf{v}_{2}),(\mathbf{v}_{1},\mathbf{v}_{4}),(\mathbf{v}_{2},\mathbf{v}_{3}),(\mathbf{v}_{3},\mathbf{v}_{4})\\}$.

Figure 1: The input set for visualizing each step of the algorithm.

#### Affine map of the vertices with weights and biases

The affine map comprises the product of each vertex of $\mathcal{P}$ by the
weights associated with layer $l$ plus the biases of the same layer. Algorithm
1 contains the steps to compute the affine map of all vertices from
$\mathcal{V}$.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$,
$\mathbf{V}=\left[\mathbf{v}[1],\ldots,\mathbf{v}[o]\right]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$,
$\mathbf{W}\in\mathbb{R}^{m\times n}$ and
$\boldsymbol{\theta}\in\mathbb{R}^{m}$

Function  _$\mathbf{AM}$(_$\mathbf{V},\mathbf{W},\boldsymbol{\theta}$ _)_:

let $\mathbf{Z}[1..o]$ be a new array, where each $\mathbf{Z}[i]\text{ is an
empty array},\,\forall i\in\\{1,\ldots,o\\}$

for  _$i\in\\{1,\ldots,o\\}$_ do

$\mathbf{Z}[i]=\mathbf{W}\mathbf{v}[i]+\boldsymbol{\theta}$ ;

end for

return $\mathbf{Z}$

Algorithm 1 Affine Map (AM)

By applying the affine map to the input set, presented in Figure 1, we have
the output of this step as shown in Figure 1. The weight matrix $\mathbf{W}$
and biases $\boldsymbol{\theta}$ employed in this toy example are as follows:

$\mathbf{W}=\begin{bmatrix}0.492693&-1.29232\\\
0.925861&0.675146\end{bmatrix},\quad\boldsymbol{\theta}=\begin{bmatrix}-0.18857972\\\
-0.14839205\end{bmatrix}$ (3)

Note that, as expected, the output of the current operation is a simple affine
transformation of the vertices of $\mathcal{P}$.

Figure 2: Representation of the output of the affine map operation.

#### Adjacent vertex identification

Following the affine map, the edge identification process takes place. Let
$\mathbf{v}_{i},\mathbf{v}_{j}\in\mathcal{V}$ be two distinct vertices, such
that $i\neq j$. If there is at least one combination of the vertices of
$\mathcal{V}$, except from $\mathbf{v}_{i}$ or $\mathbf{v}_{j}$, that allow us
to compute the middle point of $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$, given by
$(\mathbf{v}_{i}+\mathbf{v}_{j})/2$, then $\mathbf{v}_{i}$ and
$\mathbf{v}_{j}$ are not adjacent. Equation 4 presents this feasibility
problem as a MILP (Mixed Integer-Linear Programming), where
$\boldsymbol{\lambda}\in\mathbb{R}^{o}$ is the vector used to represent the
convex combination of the vertices of $\mathcal{P}$ and $\eta\in\\{0,1\\}$ is
the binary variable that enables evaluating $\mathbf{v}_{i}$ or
$\mathbf{v}_{j}$ separately.

$\displaystyle\max$ $\displaystyle 0\,\text{,}$ (4a) s.t.
$\displaystyle\frac{(\mathbf{v}_{i}+\mathbf{v}_{j})}{2}\,=\,\sum_{k=1}^{o}\mathbf{v}_{k}\lambda_{k}\,\text{,}$
(4g) $\displaystyle\sum_{k=1}^{o}\lambda_{k}\,=\,1\,\text{,}$
$\displaystyle\lambda_{k}\,\geq\,0\,\text{,}\forall
k\in\\{1,\ldots,o\\}\,\text{,}$
$\displaystyle\lambda_{i}\,\leq\,\eta\,\text{,}$
$\displaystyle\lambda_{j}\,\leq\,1-\eta\,\text{,}$
$\displaystyle\eta\,\in\,\\{0,1\\}$

Algorithm 2 presents the identification of adjacent vertices. It creates an
undirected graph using an adjacency list as the data structure to represent
the $1$-skeleton of the polytope, which is the edge structure of some polytope
[19, 20]. We denote the adjacency list of the undirected graph by
$\mathbf{E}$, were each element is a set containing the indices of the
adjacent vertices.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$,
$\mathbf{V}=[\mathbf{v}[1],\ldots,\mathbf{v}[o]]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$

Function  _$\mathbf{EI}$(_$\mathbf{V}$ _)_:

let $\mathbf{E}[1\ldots o]$ be a new array, where each
$\mathbf{E}[i]\,,\,\forall i\in\\{1,\ldots,o\\}$, corresponds to an empty set

for  _$i\in\\{1,\ldots,o\\}$_ do

for  _$j\in\\{i+1,\ldots,o\\}$_ do

if
_$\lnot\exists\boldsymbol{\lambda}\in\mathbb{R}^{o}:\frac{1}{2}(\mathbf{v}[i]+\mathbf{v}[j])=\sum_{k=1}^{o}\mathbf{v}[k]\lambda[k]\land\sum_{k=1}^{o}\lambda[k]=1\land\lambda[k]\geq
0\,,\,\forall k\in\\{1,\ldots,o\\}\land(\lambda[i]=0\lor\lambda[j]=0)$_ then

$\mathbf{E}[i]\leftarrow\mathbf{E}[i]\cup\\{j\\}$

end for

end for

return $\mathbf{E}$

Algorithm 2 Edge-skeleton Identification (EI)

Following the toy problem, where the affine map is shown in Figure 2, the
undirected graph that represents the edge structure of the affine map of
$\mathcal{P}$ is presented in Figure 3. This graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ consists of the vertex set
$\mathcal{V}=\\{\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3},\mathbf{v}_{4}\\}$
and edge set
$\mathcal{E}=\\{(\mathbf{v}_{1},\mathbf{v}_{2}),(\mathbf{v}_{1},\mathbf{v}_{4}),(\mathbf{v}_{2},\mathbf{v}_{3}),(\mathbf{v}_{3},\mathbf{v}_{4})\\}$.

Notice that each edge in the undirected graph indicates that the associated
vertices are connected by an edge of the corresponding polytope.

Figure 3: Undirected graph that represents the edge structure of the output of
the affine map of $\mathcal{P}$.

#### Polytope intersection with orthant’s hyperplanes

The next step for the algorithm concerns the identification of those points at
the intersection of the hyperplanes that define each orthant with the edges of
$\mathcal{P}$. As the edges of the polytope were previously computed, only the
edges whose extreme points belong to different orthants need to be identified,
i.e., extreme points that have at least one component with a different sign.
So, we compute the difference of the sign for each element of both extreme
points for a given edge, given by Equation 5:

$\mathbf{a}={\tt sign}(\mathbf{v}_{i})-{\tt sign}(\mathbf{v}_{j})$ (5)

where $i$ denotes the vertex under consideration, $j\in\mathbf{E}_{i}$
represents the index of the vertices that are adjacent to $\mathbf{v}_{i}$,
and ${\tt sign}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the signal mapping
that associates $1$ for the positive elements, $-1$ for the negative elements,
and $0$ for the null ones. We also define a function
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$, given by Equation 6. For the cases
where a component of $\mathbf{a}$ is greater than or equal to $2$, we have
that the sign indeed changed. For those cases where the difference is equal to
$1$, there was a null component in one of the vertices, which does not
indicate that they are in different orthants.

$\sigma(\alpha)=\begin{cases}1,&\text{if }|\alpha|\geq 2\\\
0,&\text{otherwise}\end{cases}$ (6)

Then, for those cases where there is a sign change, or, in other words, for
those cases that satisfy the condition:

$\sum_{k=1}^{n}\sigma(a_{k})\neq 0$ (7)

we compute those points at which the intersection between the edge and the
hyperplane occurred. Each sign change generates an intersection point.

Given two vertices $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$, suppose that there
is a $k\in\\{1,\ldots,n\\}$ for which $\sigma(a_{k})\neq 0$, indicating that
$\mathbf{v}_{i}$ and $\mathbf{v}_{j}$ belong to different orthants. So, there
exists a value of $\lambda$, such that $0\leq\lambda\leq 1$, defining a convex
combination of the vertices that lies on the orthant hyperplane, namely a
$\lambda$ that satisfies:

$\lambda
v_{j,k}+(1-\lambda)v_{i,k}=0\iff\lambda=\frac{-v_{i,k}}{v_{j,k}-v_{i,k}}$ (8)

Figure 4 contains a visual representation of this process for the case of
vertices $\mathbf{v}_{1}=(-0.988207,1.45261)$ and
$\mathbf{v}_{4}=(1.59643,0.102323)$. Notice that $\mathbf{a}={\tt
sign}(\mathbf{v}_{1})-{\tt sign}(\mathbf{v}_{4})=(-1,1)-(1,1)=(-2,0)$, thus
$\boldsymbol{\sigma}(\mathbf{a})=(1,0)$ implying that $\mathbf{v}_{1}$ and
$\mathbf{v}_{4}$ lie in different orthants, and the edge connecting them
intercepts the hyperplane defining the respective orthant at a point
$\mathbf{p}$. By using the above equation, we compute the convex combination
parameter $\lambda$ to be $0.382338$.

Figure 4: Representation of the intersection identification process. The blue
dots represent two adjacent vertices and the line segment is the edge
connecting them. The red dot is the intersection point of the edge with
respect to the orthant hyperplane supporting plane.

Then, we compute each component $s\in\\{1,\ldots,n\\}$ of the point
$\mathbf{p}\in\mathbb{R}^{n}$, as presented by Equation 9,

$p_{s}=(v_{j,s}-v_{i,s})\lambda+v_{i,s},$ (9)

which represents the intersection between an edge of $\mathcal{P}$ and one of
the hyperplanes of an orthant from $\mathbb{R}^{n}$. For the example in Figure
4, the intercept of the edge $(\mathbf{v}_{1},\mathbf{v}_{4})$ with the
supporting hyperplane $\\{\mathbf{x}=(x_{1},x_{2}):x_{1}=0\\}$ of the orthant
$\\{\mathbf{x}:\mathbf{x}\geq 0\\}$ is $\mathbf{p}=(0,0.93635)$.

The procedure, formalized by Algorithm 3, computes the intersection points
between each edge of $\mathcal{P}$ and each orthant’s supporting hyperplane
that are intercepted by the edge.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$ where
$\mathbf{V}=[\mathbf{v}[1],\ldots,\mathbf{v}[o]]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$,

$\mathbf{E}=[\mathbf{E}[1],\ldots,\mathbf{E}[o]]$

Function  _$\mathbf{II}$(_$\mathbf{V},\mathbf{E}$_)_:

for  _$i\in\\{1,\ldots,o\\}$_ do

/* For each edge of vertice i */

for  _$j\in\mathbf{E}[i]$ _ do

$\mathbf{a}\leftarrow{\tt sign}(\mathbf{v}[i])-{\tt sign}(\mathbf{v}[j])$

for  _$k\in\\{1,\ldots,n\\}$_ do

if  _$\sigma(a[k])\neq 0$_ then

$\lambda\leftarrow\frac{-v[i,k]}{v[j,k]-v[i,k]}$

let $\mathbf{p}$ be a point in $\mathbb{R}^{n}$

for  _$s\in\\{1,\ldots,n\\}$_ do

$p[s]\leftarrow\frac{v[j,s]-v[i,s]}{v[i,s]}\lambda$

end for

$\mathbf{V}\leftarrow[\mathbf{V}\,|\,\mathbf{p}]$ /* Append $\mathbf{p}$ to
matrix $\mathbf{V}$ */

end for

end for

end for

return $\mathbf{V}$

Algorithm 3 Orthant intersection identification (II)

The result of the intersection identification for the output of the affine map
of $\mathcal{P}$ is presented in Figure 5. There we present all the vertices
of the polytope in addition to the computed intersection points.

Figure 5: Intersection identification results representation.

#### ReLU Mapping

After calculating the intersection points between edges and orthant’s
hyperplanes, the ReLU mapping of the resulting points is computed. We apply
the mapping $ReLU:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, defined by
Equation 10, to all the points given by the previous procedure. Algorithm 4
summarizes the process where $\mathbf{V}$ is the output vertex set from
Algorithm 3.

$ReLU(\mathbf{x})=\max(0,\mathbf{x})$ (10)

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$,
$\mathbf{V}=[\mathbf{v}[1],\ldots,\mathbf{v}[o]]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$

Function  _$\mathbf{ReLU}$(_$\mathbf{V}$ _)_:

for  _$i\in\\{1,\ldots,o\\}$_ do

for  _$j\in\\{1,\ldots,n\\}$_ do

$v[i,j]\leftarrow\max(0,v[i,j])$

end for

end for

return $\mathbf{V}$

Algorithm 4 ReLU map

We can see the result of the application of the $ReLU$ mapping as presented in
Figure 6. As expected, all the vertices were projected to the positive
orthant, whereby vertices $\mathbf{v}_{2}$, $\mathbf{v}_{5}$, and
$\mathbf{v}_{6}$ are projected onto the origin. Notice that the output set is
not convex.

Figure 6: Representation of the output of the $ReLU$ mapping of the vertices
and intersection points.

#### Removing Non-vertices

To simplify the output of a single layer, algorithm APNM computes the convex
hull on the output of the ReLU mapping. To compute the convex hull, the
algorithm removes the points that are not vertices. To create a generalized
process for eliminating non-vertices, and due to the high dimensionality of
the internal layers of the neural network, we propose the application of a
feasibility analysis to identify the points that can be expressed as a convex
combination from the remaining points.

The feasibility problem is formalized in Equation 11, as follows:

$\displaystyle\min$ $\displaystyle 0\,\text{,}$ (11a) s.t.
$\displaystyle\sum_{i=1}^{o}\lambda_{i}\mathbf{v}_{i}\,=\,\mathbf{v}_{k}\,\text{,}$
(11e) $\displaystyle\sum_{i=1}^{o}\lambda_{i}\,=\,1\,,$
$\displaystyle\lambda_{k}\,=\,0\,,$
$\displaystyle\lambda_{i}\,\geq\,0\,,\,\forall i\in\\{1,\ldots,o\\},$

where, for each vertex candidate, denoted by $\mathbf{v}_{k}$, we search for a
$\boldsymbol{\lambda}$ vector such that $\mathbf{v}_{k}$ can be described as a
convex combination of the remaining vertices. If that is the case, then
$\mathbf{v}_{k}$ is not a vertex and it can be removed. Otherwise, the point
is a vertex and must remain in the set of vertices. The non-vertex elimination
is formalized in Algorithm 5.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$,
$\mathbf{V}=[\mathbf{v}[1],\ldots,\mathbf{v}[o]]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$

Function  _$\mathbf{RP}$(_$\mathbf{V}$ _)_:

for  _$k\in\\{1,\ldots,o\\}$_ do

if
_$\lnot\exists\boldsymbol{\lambda}\in\mathbb{R}^{o}:\sum_{i=1}^{o}\lambda_{i}\mathbf{v}[i]=\mathbf{v}[k]\land\sum_{i=1}^{o}\lambda_{i}=1\land\lambda_{k}=0\land\lambda_{i}\geq
0\,,\,\forall i\in\\{1,\ldots,o\\}$_ then

$\mathbf{V}\leftarrow[\mathbf{v}[1],\ldots,\mathbf{v}[k-1],\mathbf{v}[k+1],\ldots,\mathbf{v}[o]]$

end for

return $\mathbf{V}$

Algorithm 5 Removing Internal Points (RP)

The output of the process of removing non-vertices can be viewed in Figure 7.
As we can see, instead of returning the exact non-convex set, this algorithm
computes an over-approximation for the output set which is simpler than the
exact output as it is a unique set.

Figure 7: Representation of the convex hull of the output of the $ReLU$
mapping.

#### Approximate Polytope Network Mapping

Finally, we show the complete process, presented by the Algorithm 6, to
compute an approximation for the output reachable set for a given input
polyhedron $\mathcal{P}$.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$, $\mathbf{W}_{l}$ and
$\boldsymbol{\theta}_{l}$ for $l\in\\{1,\ldots,L\\}$

Function  _$\mathbf{SAPNM}$(_$\mathbf{V},\mathbf{W},\boldsymbol{\theta}$ _)_:

$\mathbf{Z}[0]\leftarrow\mathbf{V}$

for  _$l\in\\{1,\ldots,L\\}$_ do

$\hat{\mathbf{Z}}[l]\leftarrow\textbf{AM}(\mathbf{Z}[l-1],\mathbf{W}_{l},\boldsymbol{\theta}_{l})$
/* Compute affine mapping for layer $l$ */

if  _$l <L$_ then

$\mathbf{E}[l]\leftarrow\textbf{EI}(\hat{\mathbf{Z}}[l])$ /* Compute edge-
skeleton of polyhedron given by vertices $\hat{\mathbf{Z}}[l]$ */

$\hat{\mathbf{Z}}[l]\leftarrow\textbf{II}(\hat{\mathbf{Z}}[l],\mathbf{E}[l])$
/* Obtain intercept of edges with horthant hyperplanes */

$\hat{\mathbf{Z}}[l]\leftarrow\textbf{ReLU}(\hat{\mathbf{Z}}[l])$ /* Apply
ReLU mapping to the vertex set $\hat{\mathbf{Z}}[l]$ */

$\mathbf{Z}[l]\leftarrow\textbf{RP}(\hat{\mathbf{Z}}[l])$ /* Remove non-vertex
points from vertex set $\mathbf{Z}[l]$ */

end for

return $\mathbf{Z}[L]$

Algorithm 6 Simple Approximated Politope Network Mapping

###  3.3 Exact polytope network mapping (EPNM)

We present in this section the second proposed approach, which similarly to
the first one computes the reachable set by using the vertices of a given
input polytope, but computes the exact output instead. We consider the same
input closed convex polytope $\mathcal{P}$ and the same set of vertices
$\mathcal{V}$ as presented in the previous section. This approach comprises
six parts:

  1. 1.

Affine map of the vertices with weights and biases;

  2. 2.

Adjacent vertex identification;

  3. 3.

Origin verification;

  4. 4.

Polytope intersection with orthant’s hyperplanes;

  5. 5.

Separate points according to the orthant to which they belong;

  6. 6.

ReLU mapping; and

  7. 7.

Removing non-vertices;

In this section we only present parts 3 and 5, which were not previously
introduced.

#### Origin verification

After we compute the edges of $\mathcal{P}$, we need to verify if the origin
is inside it. In case it is true, then we need to insert it into the set of
vertices that will be partitioned with respect to the orthant that they belong
in the next part of the algorithm. This verification is compulsory due to the
separation of vertices according to the orthant to which they belong (for
instance, if the origin is in $\mathcal{P}$ and is not in the set of the
vertices of the partitions of the polytope, then the union of the partitions
will not be equal $\mathcal{P}$), as presented in Figure 8. Notice that, each
partition represent the portion of $\mathcal{P}$ that is inside an orthant,
i.e., the portion of $\mathcal{P}$ that is in the positive orthant, in Figure
8, is the union of the red area and the shaded area. We can see that if the
origin is not part of the set of vertices of each partition of $\mathcal{P}$,
after the orthant separation process, the shaded area will be missing for the
portion of $\mathcal{P}$ inside the positive orthant. Thus, we need to include
the origin point in the respective set so that the portion of $\mathcal{P}$ in
the positive orthant also includes the shaded area.

Figure 8: Representation of the issue that will occur if the origin is not
inserted to the vertex set for the orthant separation process.

To perform such a verification, we state the search problem as a feasibility
analysis, in which we try to verify if there exists a vector
$\boldsymbol{\lambda}\in\mathbb{R}^{o}$ such that the origin is a convex
combination of the vertices in $\mathcal{V}$. This problem is formalized by
Equation (12).

$\displaystyle\min$ $\displaystyle 0\,\text{,}$ (12a) s.t.
$\displaystyle\sum_{i=1}^{o}\lambda_{i}\mathbf{v}_{i}\,=\,0\,\text{,}$ (12d)
$\displaystyle\sum_{i=1}^{o}\lambda_{i}\,=\,1\,,$
$\displaystyle\lambda_{i}\,\geq\,0\,,\,\forall i\in\\{1,\ldots,o\\}\text{.}$

Algorithm 7 formalizes the search process associated with the verification of
the origin.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$,
$\mathbf{V}=[\mathbf{v}[1],\ldots,\mathbf{v}[o]]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$

Function  _$\mathbf{OS}$(_$\mathbf{V}$ _)_:

if
_$\exists\boldsymbol{\lambda}\in\mathbb{R}^{o}:\sum_{i=1}^{o}\lambda_{i}\mathbf{v}_{i}=0\land\sum_{i=1}^{o}\lambda_{i}=1\land\lambda_{i}\geq
0\,,\,\forall i\in\\{1,\ldots,o\\}$_ then

$\mathbf{V}\leftarrow[\mathbf{v}[1],\ldots,\mathbf{v}[o]\,|\,\mathbf{0}]$ /*
Adding the origin to the vertex set */

return $V$

Algorithm 7 Origin Search (OS)

#### Separating points to their respective orthants

We present in this section the process of splitting the vertex set
$\mathcal{V}$ into several sets, each associated with a different orthant. In
other words, we split $\mathcal{V}$ in such a way that each resulting set
contains points located in a single orthant (observe that the convex hull of
each set of vertices represents the partition of $\mathcal{P}$ that is inside
an orthant). We begin by computing:

$\mathbf{b}=\frac{1}{2}({\tt
sign}(\mathbf{v}_{i})+\mathbf{1}_{|\mathbf{v}_{i}|})$ (13)

where, $\mathbf{v}_{i}$ is the vertex under analysis, ${\tt
sign}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the sign mapping previously
presented and $\mathbf{1}_{|\mathbf{v}_{i}|}$ is a vector with all the
components equal to one and has size $|\mathbf{v}_{i}|$.

The elements of $\mathbf{b}$ with value equal to $1$ are associated with a
positive element of $\mathbf{v}_{i}$, the elements equal to $0$ are associated
with negative components of the vertex, and those components of $\mathbf{b}$
equal to $1/2$ are associated with the null components of the vertex. Notice
that a null component of a vertex means that this point belongs to at least
two different orthants (the origin being a special case inside all the
orthants).

For those cases where the component of $\mathbf{b}$ is equal to $1/2$, we must
guarantee that the associated vertex $\mathbf{v}_{i}$ is properly inserted
into each of those sets that are associated with the orthants it belongs to
(i.e., if $\mathbf{v}_{i}$ has one null component, it must be inserted in two
sets). Then, the algorithm starts a verification process to identify those
component of $\mathbf{b}$ that are equal $1/2$. In case it is equal true, two
copies of $\mathbf{v}_{i}$ must be created. This verification process is
repeated until all components of $\mathbf{b}$ have been checked. This process
is detailed in Algorithm 8.

Input: $\mathbf{b}\in[0,1]^{n}$

Function  _$\mathbf{ZV}$(_$\mathbf{b}$ _)_:

$\mathcal{A}\leftarrow\\{\mathbf{b}\\}$

$\mathcal{B}\leftarrow\\{\mathbf{b}\\}$

while  _$|\mathcal{B}|\neq 0$_ do

$\mathcal{B}\leftarrow\emptyset$

for  _$\mathbf{a}\in\mathcal{A}$_ do

for  _$i\in\\{1,\ldots,|\mathbf{a}|\\}$_ do

/* Check for those components that are equal $1/2$ */

if  _$a[i]=1/2$_ then

$\mathbf{b}\leftarrow\mathbf{a}$

$b[i]\leftarrow 1$

$\mathcal{B}\leftarrow\mathcal{B}\cup\mathbf{b}$ /* Insert the first copy into
$\mathcal{B}$ */

$\mathbf{b}\leftarrow\mathbf{a}$

$b[i]\leftarrow 0$

$\mathcal{B}\leftarrow\mathcal{B}\cup\mathbf{b}$ /* Insert the second copy
into $\mathcal{B}$ */

break

end for

end for

if  _$|\mathcal{B}|\neq 0$_ then

$\mathcal{A}\leftarrow\mathcal{B}$

end while

return $\mathcal{A}$ /* return a set of vectors */

Algorithm 8 Zeros Verification (ZV)

After checking for null values in $\mathbf{v}_{i}$, the sets in which it must
be placed must be defined. It should be noted that at the end of this process,
$\mathbf{b}$ is a binary vector that represents the index in which
$\mathbf{v}_{i}$ must be placed (there might be more than one $\mathbf{b}$ for
a single vertex). Equation (14) presents the process for converting
$\mathbf{b}$ into a decimal number $q$. The calculation of the index is
presented in Algorithm 9.

$q=\sum_{i=1}^{n}b_{i}2^{i-1}$ (14)

Input: $\mathbf{b}\in\\{0,1\\}^{m}$

Function  _$\mathbf{AP}$(_$\mathbf{b}$ _)_:

let $q$ be an integer

$q\leftarrow 0$

for  _$i\in\\{1,\ldots,n\\}$_ do

$q\leftarrow q+b[i]2^{i-1}$

end for

return $q$

Algorithm 9 Get Array Position (AP)

Figure 9 illustrates a toy example of how this process works for a single
vertex $\mathbf{v}_{i}$. In this example, $\mathbf{v}_{i}=(-2,0)$, with the
second component being null. As a result, $\mathbf{b}=(0,\frac{1}{2})$. Next,
two copies of $\mathbf{b}$ are created. The components of the first (second)
copy that are equal to $\frac{1}{2}$ are replaced by $0$ ($1$). In this case,
as there was only one null component, the generated vectors are $(0,0)$ and
$(0,1)$. The vector generated from $\mathbf{b}$ contains the binary
representation of the indices of the sets that $\mathbf{v}_{i}$ must be placed
in. In this example, $\mathbf{v}_{i}$ must be placed in the orthants with
indices $0$ and $2$. It is important to note that the indices used by the
algorithm to identify the destination sets are different from the traditional
orthant enumeration (2nd and 3rd quadrants in this case).

Figure 9: Vertex separation process for a single vertex. Given a vertex
$\mathbf{v}_{i}=(-2,0)$, the associated $\mathbf{b}$ is given by
$(0,\frac{1}{2})$. As we have that there is a component of $\mathbf{b}$ equal
to $\frac{1}{2}$, we split it into two vectors in such a way that the
component with value equal to $\frac{1}{2}$ is substituted by $1$ and $0$ in
the newly created vectors. These new vertices contain the binary
representation of the orthant’s index that $\mathbf{v}_{i}$ belong to. In the
presented example, $\mathbf{v}_{i}$ belong to 2 different orthants (the blue
and the green ones).

The orthant separation complete process is computed accordingly by Algorithm
10.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$,
$\mathbf{V}=[\mathbf{v}[1],\ldots,\mathbf{v}[o]]$,
$\mathbf{v}[i]\in\mathbb{R}^{n}\,,\,\forall i\in\\{1,\ldots,o\\}$

Function  _$\mathbf{SP}$(_$\mathbf{V}$ _)_:

let $\mathbf{Z}[1\ldots 2^{n}]$ be an array

for  _$i\in\\{1,\ldots,n\\}$_ do

let $\mathbf{b}[1,\ldots,n]$ be an array

$\mathbf{b}\leftarrow\frac{1}{2}({\tt
sign}(\mathbf{v}[i])+\mathbf{1}_{|\mathbf{v}[i]|})$

% identify if there are null coordinates

$\mathcal{B}\leftarrow\textbf{ZV}(\mathbf{b})$ /* Identify the orthants to
which $v[i]$ belongs */

for  _$\mathbf{b^{\prime}}\in\mathcal{B}$_ do

$q\leftarrow\textbf{AP}(\mathbf{b^{\prime}})$ /* Calculate index position of
orthant $\mathbf{b}^{\prime}$ to which $\mathbf{v}[i]$ belongs */

if  _$\mathbf{Z}[q]=\emptyset$_ then

$\mathbf{Z}[q]\leftarrow\mathbf{v}[i]$

else

$\mathbf{Z}[q]\leftarrow\left[\mathbf{Z}[q]\,|\,\mathbf{v}[i]\right]$ /*
Appending the new vertex */

end for

end for

return $\mathbf{Z}$

Algorithm 10 Separate points per orthant (SP)

Continuing the example from the previous subsection, the orthant separation
performed by Algorithm 10 will divide the polytope illustrated in Figure 5
into four separate polytopes, one for each orthant, as shown in Figure 10. It
is worth noting that vertices $\mathbf{v}_{5}$, $\mathbf{v}_{6}$,
$\mathbf{v}_{7}$, $\mathbf{v}_{8}$, and $\mathbf{v}_{9}$ have been placed in
more than one set, as they are located in more than one orthant.

Figure 10: Representation of the separation of the polytope $\mathcal{P}$
after the application of the affine map.

#### Exact Polytope Network Mapping

The layer mapping and the complete process for the Exact Polytope Network
Mapping is formalized in Algorithm 11. The algorithm takes as input the
vertices of the input polytope aimed to be verified. Then, it computes the
affine map of these vertices (calculated by the algorithm AM) and identifies
the edges between adjacent vertices (with the algorithm EI). For the vertices
connected by an edge, the algorithm computes the intersection of the
corresponding edge with each orthant’s supporting hyperplane, when those
vertices are not in the same orthant (by means of the algorithm II), and
verifies whether or not the origin belongs to the polytope under analysis
(checked by the algorithm OS). Next, it separates all of these points in
different sets, where each set contains those vertices that belong to a single
orthant (computed by the algorithm SP). Finally, the algorithm performs the
$ReLU$ mapping and removes non-vertices for each partition generated. As later
presented, the $ReLU$ mapping preserves the convexity inside a given orthant,
though there are points that are not vertices after the application of the
non-linear mapping (check vertex $\mathbf{v}_{3}$ in Figure 6).

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$, /* vertices of the polytope */

$\mathbf{W}_{l}$ and $\boldsymbol{\theta}_{l}$ for $l\in\\{1,\ldots,L\\}$ /*
network layers’ weights and biases */

Function  _$\mathbf{EPNM}$(_$\mathbf{V},\mathbf{W},\boldsymbol{\theta}$ _)_:

let $\mathcal{P}$ be an empty set

$\mathcal{P}\leftarrow\\{\mathbf{V}\\}$ /* The list of polytopes’s vertices;
starting with the input polytope */

for  _$l\in\\{1,\ldots,L\\}$_ do

$\hat{\mathcal{P}}\leftarrow\emptyset$ /* $\hat{\mathcal{P}}$ is the auxiliary
set that represents $\mathcal{P}$ in the next layer */

for  _$\mathbf{Z}\in\mathcal{P}$_ do

/* iterate for each set $\mathcal{Z}$ of vertices in $\mathcal{P}$ */

if  _$l >1$_ then

$\mathbf{Z}\leftarrow\textbf{ReLU}(\mathbf{Z})$ /* Apply ReLU mapping to the
vertex set $\mathbf{Z}$ */

$\mathbf{Z}\leftarrow\textbf{RP}(\mathbf{Z})$ /* Remove non-vertex points from
vertex set $\mathbf{Z}$ */

$\mathbf{Z}\leftarrow\textbf{AM}(\mathbf{Z},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$
/* Perform affine mapping,
$\mathbf{Z}\leftarrow\mathbf{W}_{l}\mathbf{Z}+\boldsymbol{\theta}_{l}$ */

if  _$l <L$_ then

$\mathbf{E}\leftarrow\textbf{EI}(\mathbf{Z})$ /* Compute edge-skeleton
$\mathbf{E}$ of the polyhedron given by vertices $\mathbf{Z}$ */

$\mathbf{Z}\leftarrow\textbf{II}(\mathbf{Z},\mathbf{E})$ /* Add edge
intercepts with orthant’s supporting hyperplanes to vertex set */

$\mathbf{Z}\leftarrow\textbf{OS}(\mathbf{Z})$ /* Add origin to the vertex set
if needed */

$\mathbf{Z}\leftarrow\textbf{SP}(\mathbf{Z})$ /* Transform the vertex set into
a list of vertices for each orthant, considering the orthants to which each
vertex belongs */

for  _$k\in\\{1,\ldots,|\mathbf{Z}|\\}$_ do

if  _$\mathbf{Z}[k]\neq\emptyset$_ then

$\hat{\mathcal{P}}\leftarrow\hat{\mathcal{P}}\cup\mathbf{Z}[k]$

end for

end for

$\mathcal{P}\leftarrow\hat{\mathcal{P}}$

end for

return $\mathcal{P}$

Algorithm 11 Exact Polytope Network Mapping

Finally, after the application of the $ReLU$ mapping and the removal of non-
vertices (RP) for each partition in the set, as illustrated in Figure 10, the
result of the EPNM algorithm for a single layer mapping is shown in Figure 11.
Observe that this output, which is a union of sets, represents the mapping of
a unique layer. For each of the sets that result from the previous layer, the
algorithm computes the exact mapping for the next layer. This process will
repeat until the algorithm maps all the sets to the last layer of the neural
network. For instance, considering that the visual problem has one extra
layer, each of the four output sets (one set for each orthant, as the
intersection of $\mathcal{P}$ with each orthant is not empty) will be mapped
to the next layer following the same described process.

Figure 11: Representation of the output reachable set for a single layer of
the EPNM algorithm.

###  3.4 Partially approximated polytope network mapping (PAPNM)

We have proposed a third approach, which involves a slight modification of the
EPNM algorithm. Instead of an exact mapping between layers, this approach
allows for the merging of some of the resulting sets. By computing the convex
hull of the union for some of the output sets, we aim to achieve a more
accurate approximation of the output set (compared to APNM) while also
reducing the execution time of the algorithm (compared to EPNM). The result is
presented in Algorithm 12.

Input: $\mathbf{V}\in\mathbb{R}^{o\times n}$, /* vertices of the polytope */

$\mathbf{W}_{l}$ and $\boldsymbol{\theta}_{l}$ for $l\in\\{1,\ldots,L\\}$ /*
network layers’ weights and biases */

/* $d$ is the size of each group of sets that will be merged */

Function  _$\mathbf{PAPNM}$(_$\mathbf{V},\mathbf{W},\boldsymbol{\theta},d$
_)_:

let $\mathcal{P}$ be an empty set

$\mathcal{P}\leftarrow\\{\mathbf{V}\\}$

for  _$l\in\\{1,\ldots,L\\}$_ do

$\hat{\mathcal{P}}\leftarrow\emptyset$

for  _$\mathbf{Z}\in\mathcal{P}$_ do

if  _$l >1$_ then

$\mathbf{Z}\leftarrow\textbf{ReLU}(\mathbf{Z})$

$\mathbf{Z}\leftarrow\textbf{RP}(\mathbf{Z})$

$\mathbf{Z}\leftarrow\textbf{AM}(\mathbf{Z},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$

if  _$l <L$_ then

$\mathbf{E}\leftarrow\textbf{EI}(\mathbf{Z})$

$\mathbf{Z}\leftarrow\textbf{II}(\mathbf{Z},\mathbf{E})$

$\mathbf{Z}\leftarrow\textbf{OS}(\mathbf{Z})$

$\mathbf{Z}\leftarrow\textbf{SP}(\mathbf{Z})$

$\mathbf{Z}\leftarrow\textbf{MS}(\mathbf{Z},d)$ /* new function for merging
sets */

for  _$k\in\\{1,\ldots,|\mathbf{Z}|\\}$_ do

$\hat{\mathcal{P}}\leftarrow\hat{\mathcal{P}}\cup\mathbf{Z}_{k}$

end for

end for

$\mathcal{P}\leftarrow\hat{\mathcal{P}}$

end for

return $\mathcal{P}$

Algorithm 12 Partially approximate polytope network mapping

The main difference between Algorithm 11 and Algorithm 12 is the application
of the $MS$ procedure after the separation process. The $MS$ algorithm
comprises the merging procedure, where a set of sets of vertices is given as
input. As each set of vertices represents a single polytope, this procedure
merges some of these sets of vertices to produce fewer sets compared to the
exact mapping. It is worth noting that the output of $MS$ is an
overapproximation of the input set it receives, and that the extreme case in
which a single set of vertices is computed is exactly the case of APNM.

Here, we propose a basic approach for the merging process in Algorithm 13. The
$MS$ algorithm takes the set $\mathcal{P}$ which contains the sets of vertices
and an integer $d$ as inputs. The sets of vertices are then grouped into sets
of size $d$ and each group is merged (i.e, for the case where $\mathcal{P}$
has 13 sets and $d=3$, there will be 4 groups of 3 sets and 1 group of 1 set).
It is important to note that this is just a simple example of the merging
procedure and that different heuristics can be implemented to improve this
process.

Input:
$\mathcal{P}=\\{\mathbf{V}_{1},\ldots,\mathbf{V}_{q}\\},\mathbf{V}_{i}\in\mathbb{R}^{o_{i}\times
n},i\in\\{1,\ldots,q\\}$ /* set of polytopes’ vertices */

$d\in\mathbb{N}$ /* $d$ is the size of each group of sets that will be merged
*/

Function  _$\mathbf{MS}$(_$\mathcal{P},d$ _)_:

let $\mathcal{A}$ be an empty set

for  _$i\in\\{1,\ldots,\left\lceil\frac{|\mathcal{P}|}{d}\right\rceil\\}$_ do

let $\mathbf{B}[1\ldots n]$ be an array

for  _$j\in\\{(i-1)\times d+1,\ldots,\min(i\times d,|\mathcal{P}|)\\}$_ do

if  _$\mathbf{B}=\emptyset$_ then

$\mathbf{B}\leftarrow\mathbf{V}[j]$

else

$\mathbf{B}\leftarrow(\mathbf{B}|\mathbf{V}[j])$

end for

$\mathcal{A}\leftarrow\mathcal{A}\cup\mathbf{B}$

end for

return $\mathcal{A}$

Algorithm 13 Merge Sets (MS)

The use of this merging procedure allows for different levels of
approximation, resulting in a range of possible approximations from the exact
mapping (EPNM) to the coarser case (APNM). It is important to ensure that the
merging procedure returns an overapproximation of the exact mapping of each
layer, which is a necessary condition to ensure the soundness of PAPNM.

##  4 Demonstrations

We present in this section demonstrations that provide theoretical guarantees
for the correctness of each proposed algorithm.

###  4.1 Identification of adjacent vertices

Let
$\mathcal{P}=\left\\{\sum_{i=1}^{o}\lambda_{i}\mathbf{v}_{i}\mid\sum_{i=1}^{o}\lambda_{i}=1,\lambda_{i}\geq
0\text{ e }\mathbf{v}_{i}\in\mathcal{V},\forall i\in\\{1,\ldots,o\\}\right\\}$
be a convex V-polytope defined as the convex combination of its vertices,
where $\mathcal{V}=\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{o}\\}$ is the set of
vertices of a polyhedron $\mathcal{P}$.

######  Definition 1.

Given a vector $\mathbf{c}\in\mathbb{R}^{n}$ and
$\delta=\max\\{\mathbf{c}^{T}\mathbf{x}\mid\mathbf{x}\in\mathcal{P}\\}$, we
have that $\mathcal{H}=\\{\mathbf{x}\mid\mathbf{c}^{T}\mathbf{x}=\delta\\}$ is
a supporting hyperplane of $\mathcal{P}$.

######  Definition 2.

$\mathcal{F}$ is a face of $\mathcal{P}$ if $\mathcal{F}=\mathcal{P}$ or
$\mathcal{F}=\mathcal{P}\cap\mathcal{H}$ for some supporting hyperplane
$\mathcal{H}$. In other words, $\mathcal{F}$ is a face of $\mathcal{P}$ if,
and only if, $\mathcal{F}$ is the set of optimal solutions for
$\max\\{\mathbf{c}^{T}\mathbf{x}\mid\mathbf{x}\in\mathcal{P}\\}$ for a given
$\mathbf{c}\in\mathbb{R}^{n}$ [21].

######  Definition 3.

$\mathbf{v}_{i}$ e $\mathbf{v}_{j}$ are adjacent vertices if there is a vector
$\mathbf{c}\in\mathbb{R}^{n}$ such that
$\mathbf{c}^{T}\mathbf{v}_{i}=\mathbf{c}^{T}\mathbf{v}_{j}=\max\\{\mathbf{c}^{T}\mathbf{x}\mid\mathbf{x}\in\mathcal{P}\\}>\mathbf{c}^{T}\mathbf{v}_{k},\forall
k\in\\{1,\ldots,o\\},k\neq i\text{ and }k\neq j$ [21].

Notice that the face $\mathcal{F}$ is an edge in the particular case stated by
Definition 3. We propose that:

######  Proposition 1.

Two extreme point $\mathbf{v}_{i},\mathbf{v}_{j}\in\mathcal{V}$ are adjacent
if, and only if, it is not possible to compute the median point,
$\bar{\mathbf{v}}=(\mathbf{v}_{i}+\mathbf{v}_{j})/2$, as a convex combination
of the vertices in $\mathcal{V}\setminus\\{\mathbf{v}_{i}\\}$ and
$\mathcal{V}\setminus\\{\mathbf{v}_{j}\\}$.

###### Proof.

Given two vertices $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$, we have two
possibilities regarding their adjacency: 1) they are adjacent, or else 2) they
are not adjacent.

For the first case, as $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$ are adjacent, we
have by definition that there is a hyperplane
$\mathcal{H}_{a}=\\{\mathbf{x}\mid\mathbf{c}^{T}\mathbf{x}=d\\}$ that contains
both $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$, such that
$\mathcal{H}_{a}\cap\mathcal{P}=\max\\{\mathbf{c}^{T}\mathbf{x}\mid\mathbf{x}\in\mathcal{P}\\}>\mathbf{c}^{T}\mathbf{v}_{k},\forall
k\in\\{1,\ldots,o\\},k\neq i\text{ and }k\neq j$ for a given
$\mathbf{c}\in\mathbb{R}^{n}$. Therefore, all the remaining extreme points
from $\mathcal{P}$, except from $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$, are in
the open half-space
$\mathcal{H}_{b}=\\{\mathbf{x}\mid\mathbf{c}^{T}\mathbf{x}<d\\}$.
Consequently, the median point of $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$,
denoted by $\bar{\mathbf{v}}={(\mathbf{v}_{i}+\mathbf{v}_{j})}/{2}$, can not
be expressed as a convex combination of $\mathbf{v}_{k}$, for
$k\in\\{1,\ldots,o\\}\setminus\\{i\\}$ or
$k\in\\{1,\ldots,o\\}\setminus\\{j\\}$.

Considering that $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$ are not adjacent, let
$\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ be two polytopes given by the convex
combination of the extreme points $\mathcal{V}\setminus\\{\mathbf{v}_{i}\\}$
and $\mathcal{V}\setminus\\{\mathbf{v}_{j}\\}$, respectively. Hence, there are
two distinct possibilities: $\bar{\mathbf{v}}\in\mathcal{P}_{j}$ or
$\bar{\mathbf{v}}\not\in\mathcal{P}_{j}$. For the first case, as
$\bar{\mathbf{v}}\in\mathcal{P}_{j}$, then we can compute this point as a
convex combination of the vertices in
$\mathcal{V}\setminus\\{\mathbf{v}_{j}\\}$. For the second case, as
$\bar{\mathbf{v}}\not\in\mathcal{P}_{j}$ and $\bar{\mathbf{v}}\in\mathcal{P}$,
then $\bar{\mathbf{v}}\in\mathcal{P}\setminus\mathcal{P}_{j}$. However, since
the vertices that are adjacent to $\mathbf{v}_{j}$ are in
$\mathcal{V}\setminus\\{\mathbf{v}_{i},\mathbf{v}_{j}\\}$, then
$\mathcal{P}\setminus\mathcal{P}_{j}\subseteq\mathcal{P}_{i}$. From this fact
it follows that $\bar{\mathbf{v}}\in\mathcal{P}_{i}$ and, consequently,
$\bar{\mathbf{v}}$ can be expressed as a convex combination of the vertices
$\mathcal{V}\setminus\\{\mathbf{v}_{i}\\}$. ∎

###  4.2 ReLU convexity inside an orthant

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ be a function that denotes the
$ReLU$ mapping, defined by the Equation (15):

$f(\mathbf{x})_{i}=\max(0,x_{i})$ (15)

######  Proposition 2.

Given $\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}$, if
$\sup\\{\gamma\mid\gamma={\tt sign}(\mathbf{x})_{i}-{\tt
sign}(\mathbf{y})_{i},\forall i\in\\{1,\ldots,n\\}\\}\leq 1$, or, in other
words, if $\mathbf{x}$ and $\mathbf{y}$ belong to the same orthant, we have
that $f(\mathbf{x}+\mathbf{y})=f(\mathbf{x})+f(\mathbf{y})$ and that
$f(\alpha\mathbf{x})=\alpha f(\mathbf{x})$.

###### Proof.

Denoting each orthant of $\mathbb{R}^{n}$ as $\mathcal{O}_{j},\forall
j\in\\{1,\ldots,2^{n}\\}$, we can rewrite the $ReLU$ mapping, given by $f$, as
$f_{j}:\mathcal{O}_{j}\rightarrow\mathbb{R}^{n}$, such that:

$f_{j}(\mathbf{x})=\boldsymbol{\Lambda}_{j}\mathbf{x}$ (16)

where $\boldsymbol{\Lambda}_{j}\in\mathbb{R}^{n\times n}$ is the matrix in
which the elements of the principal diagonal associated with a negative
component of $\mathbf{x}\in\mathcal{O}_{j}$ are equal $0$, while the remaining
elements are equal $1$ ($\lambda_{k,l}$ for $k=l$). The elements off the
principal diagonal are all equal to zero ($\lambda_{k,l}=0$ for all $k\neq
l$).

Therefore, $f_{j}$ preserves convexity since it consists of a linear mapping
that satisfies both properties stated in the proposition, the addition
property:

$\displaystyle f_{j}(\mathbf{x}+\mathbf{y})$
$\displaystyle=\boldsymbol{\Lambda}_{j}(\mathbf{x}+\mathbf{y})$
$\displaystyle=\boldsymbol{\Lambda}_{j}\mathbf{x}+\boldsymbol{\Lambda}_{j}\mathbf{y}$
$\displaystyle=f_{j}(\mathbf{x})+f_{j}(\mathbf{y})$

and the product by a scalar:

$\displaystyle f_{j}(\alpha\mathbf{x})$
$\displaystyle=\boldsymbol{\Lambda}_{j}(\alpha\mathbf{x})$
$\displaystyle=\alpha\boldsymbol{\Lambda}_{j}\mathbf{x}$ $\displaystyle=\alpha
f_{j}(\mathbf{x})$

where $\alpha\in\mathbb{R}$.

Hence, as $f_{j}$ satisfies the linearity conditions within each orthant
$\mathcal{O}_{j}$ and the linear mapping preserves convexity, it follows that:

$\displaystyle f_{j}(\theta\mathbf{x}+(1-\theta)\mathbf{y})$
$\displaystyle=\boldsymbol{\Lambda}_{j}(\theta\mathbf{x}+(1-\theta)\mathbf{y})$
$\displaystyle=\theta\boldsymbol{\Lambda}_{j}\mathbf{x}+(1-\theta)\boldsymbol{\Lambda}_{j}\mathbf{y}$
$\displaystyle=\theta f_{j}(\mathbf{x})+(1-\theta)f_{j}(\mathbf{y})$

for all $\mathbf{x},\mathbf{y}\in\mathcal{O}_{j}$ e $\theta\in[0,1]$.
Therefore, the $ReLU$ mapping preserves the convexity inside a given orthant.
∎

###  4.3 V-polytope and half-space intersection

Let
$\mathcal{P}=\left\\{\sum_{i=1}^{o}\lambda_{i}\mathbf{v}_{i}\mid\sum_{i=1}^{o}\lambda_{i}=1,\lambda_{i}\geq
0\text{ and }\mathbf{v}_{i}\in\mathcal{V},\forall
i\in\\{1,\ldots,o\\}\right\\}$ be a closed convex polyhedron defined in terms
of the convex combination of its vertices (or extreme points), where
$\mathcal{V}=\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{n}\\}$ is the set of the
vertices of $\mathcal{P}$. Such a polyhedron can also be defined as a set of
inequalities
$\mathcal{P}=\\{\mathbf{x}\mid\mathbf{C}\mathbf{x}\leq\mathbf{d}\\}$, such
that $\mathbf{C}\in\mathbb{R}^{m\times n}$ and $\mathbf{d}\in\mathbb{R}^{m}$.

(a) $\mathcal{P}\cap\mathcal{H}=\mathcal{P}$

(b) $\mathcal{P}\cap\mathcal{H}\subset\mathcal{P}$

(c) $\mathcal{P}\cap\mathcal{H}=\emptyset$

Figure 12: Representation of all the three different possible cases with
respect to the intersection of $\mathcal{P}$ with a half-space $\mathcal{H}$.

There are three different possible cases regarding the intersection of
$\mathcal{P}$ with a half-space
$\mathcal{H}=\\{\mathbf{x}\mid\mathbf{a}^{T}\mathbf{x}\leq b\\}$, where
$\mathbf{a}\in\mathbb{R}^{n}$ and $b\in\mathbb{R}$ (Figure 12 presents a
visual representation for each case):

  1. 1.

$\mathcal{P}\cap\mathcal{H}=\mathcal{P}$;

  2. 2.

$\mathcal{P}\cap\mathcal{H}\subset\mathcal{P}$;

  3. 3.

$\mathcal{P}\cap\mathcal{H}=\emptyset$.

The first and the third cases are trivial. For the first case, we have that
all the extreme points of $\mathcal{P}$ are in $\mathcal{H}$. For the third
case none of the extreme points of $\mathcal{P}$ belongs to $\mathcal{H}$.

For the second case, a new face is generated for $\mathcal{P}$, given by
$\mathcal{P}_{s}=\mathcal{P}\cap\mathcal{H}_{s}$, where
$\mathcal{H}_{s}=\\{\mathbf{x}\mid\mathbf{a}^{T}\mathbf{x}=b\\}$ such that
$\mathbf{a}\in\mathbb{R}^{n}$ and $b\in\mathbb{R}$. Note that
$\mathcal{H}_{s}$ is the supporting hyperplane of $\mathcal{H}$. Therefore,
those vertices of $\mathcal{P}$ that are not in $\mathcal{H}$, are not extreme
points of $\mathcal{P}\cap\mathcal{H}$. Thus, it is necessary to find those
extreme points of the new polytope
$\mathcal{P}_{h}=\mathcal{P}\cap\mathcal{H}$. Figure 13 presents a
visualization of the elements previously defined.

Figure 13: Representation of the elements of interest from the intersection
between $\mathcal{P}$ and $\mathcal{H}$ for the second case. The supporting
hyperplane $\mathcal{H}_{s}$ of $\mathcal{H}$ is presented in blue. In green
we have the intersection between $\mathcal{H}_{s}$ and $\mathcal{P}$, and in
red the result of the intersection between $\mathcal{P}$ with $\mathcal{H}$.
Notice that, for this example,
$\mathcal{V}_{h}=\\{\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3},\mathbf{v}_{4}\\}$.

We denote by $\mathcal{V}_{h}$ the subset of vertices of $\mathcal{P}$ that
belong to $\mathcal{H}$, given by
$\mathcal{V}_{h}=\mathcal{V}\cap\mathcal{H}$. Let
$\mathcal{E}=\\{\mathcal{E}_{1},\ldots,\mathcal{E}_{p}\\}$ be the set of edges
of $\mathcal{P}$. $\mathcal{E}_{i}$ denotes the set of points that belong to
the $i$-th edge of $\mathcal{P}$. Finally, we denote by $\mathcal{V}_{p}$ the
set of vertices from $\mathcal{P}_{h}$. Observe that
$\mathcal{V}_{h}\subseteq\mathcal{V}_{p}$.

Let $\mathbf{c}$ be a non-zero vector and
$\delta=\max\\{\mathbf{c}^{T}\mathbf{x}\mid\mathbf{C}\mathbf{x}\leq\mathbf{d}\\}$.
The affine hyperplane
$\mathcal{H}_{a}=\\{\mathbf{x}\mid\mathbf{c}^{T}\mathbf{x}=\delta\\}$ is a
supporting hyperplane of $\mathcal{P}$.

######  Definition 4.

A subset $\mathcal{F}$ of $\mathcal{P}$ is called a face of $\mathcal{P}$ if
$\mathcal{F}=\mathcal{P}$ or $\mathcal{F}=\mathcal{P}\cap\mathcal{H}_{a}$
[22].

######  Definition 5.

$\mathcal{F}$ is a face of $\mathcal{P}$ if, and only if, $\mathcal{F}$ is not
empty and
$\mathcal{F}=\\{\mathbf{x}\in\mathcal{P}\mid\mathbf{C}^{\prime}\mathbf{x}=\mathbf{d}^{\prime}\\}$,
for a subsystem $\mathbf{C}^{\prime}\mathbf{x}\leq\mathbf{d}^{\prime}$ of
$\mathbf{C}\mathbf{x}\leq\mathbf{d}$ [22].

######  Proposition 3.

Considering that $\mathcal{P}\cap\mathcal{H}\subset\mathcal{P}$ and given that
$\mathcal{V}_{p}$ is the set of extreme points of $\mathcal{P}_{h}$, if
$\mathbf{v}\in\mathcal{V}_{p}$ and $\mathbf{v}\not\in\mathcal{V}_{h}$ then
$\mathbf{v}\in\mathcal{P}_{s}$.

###### Proof.

By definition, an extreme point of a polytope is a $0$-dimensional face, thus:

$\mathcal{F}=\\{\mathbf{x}\mid\mathbf{C}^{\prime}\mathbf{x}=\mathbf{d}^{\prime}\\}$
(17)

where $\mathbf{C}^{\prime}\mathbf{x}\leq\mathbf{d}^{\prime}$ is a subsystem of
$\mathbf{C}\mathbf{x}\leq\mathbf{d}$. Remember that
$\mathcal{P}=\\{\mathbf{x}\mid\mathbf{C}\mathbf{x}\leq\mathbf{d}\\}$ and that
$\mathcal{H}=\\{\mathbf{x}\mid\mathbf{a}^{T}\mathbf{x}\leq b\\}$. Furthermore,
since a vertex $\mathbf{v}$ is the particular case of a face given by the
intersection of $n$ hyperplanes, for a $n$-dimensional space, then:

$\mathbf{v}=\\{\mathbf{x}\mid\mathbf{C}^{\prime\prime}\mathbf{x}=\mathbf{d}^{\prime\prime}\\}$
(18)

for some $\mathbf{C}^{\prime\prime}\in\mathbb{R}^{n\times n}$ and
$\mathbf{d}^{\prime\prime}\in\mathbb{R}^{n}$, such that ${\tt
det}(\mathbf{C}^{\prime\prime})\neq 0$ and $\mathbf{v}\neq\emptyset$.

However, for $\mathbf{x}\in\mathcal{H}\setminus\mathcal{H}_{s}$, no new
solution for the subsystem
$\mathbf{C}^{\prime\prime}\mathbf{x}=\mathbf{d}^{\prime\prime}$ subsystem of:

$\begin{bmatrix}\mathbf{C}\\\
\mathbf{a}^{T}\end{bmatrix}\mathbf{x}\leq\begin{bmatrix}\mathbf{d}\\\
b\end{bmatrix}$ (19)

such that ${\tt det}(\mathbf{C}^{\prime\prime})\neq 0$ and
$\mathbf{v}\neq\emptyset$, is possible. Put in other words, as no constraint
was placed in $\mathcal{H}\setminus\mathcal{H}_{s}$, no new vertex was
generated in $\mathcal{P}\cap(\mathcal{H}\setminus\mathcal{H}_{s})$.
Consequently, those vertices generated by the intersection with $\mathcal{H}$,
in case they exist, must belong to $\mathcal{H}_{s}$. Hence, if
$\mathbf{v}\in\mathcal{V}_{p}$ and $\mathbf{v}\not\in\mathcal{V}_{h}$, then
$\mathbf{v}\in\mathcal{P}_{s}$. ∎

######  Proposition 4.

Considering that $\mathcal{P}\cap\mathcal{H}\subset\mathcal{P}$, if
$\mathcal{V}_{p}$ is the set of extreme points of $\mathcal{P}_{h}$, then
$\mathcal{V}_{p}\setminus\mathcal{V}_{h}=\mathcal{E}\cap\mathcal{H}_{s}$.

###### Proof.

As $\mathcal{P}\cap\mathcal{H}\subset\mathcal{P}$, then $\mathcal{P}_{s}$ is a
face of $\mathcal{P}_{h}$. Furthermore, if $\mathcal{F}$ is a face of
$\mathcal{P}_{s}$, then $\mathcal{F}$ is a face of $\mathcal{P}_{h}$.
Therefore, the 0-dimensional faces of $\mathcal{P}_{s}$ (vertices) are also
faces of $\mathcal{P}_{h}$.

Recall that a vertex is defined by the intersection of $n$ hyperplanes, for a
$n$-dimensional space, that
$\mathcal{P}=\\{\mathbf{x}\mid\mathbf{C}\mathbf{x}\leq\mathbf{d}\\}$ and that
$\mathcal{H}=\\{\mathbf{x}\mid\mathbf{a}^{T}\mathbf{x}\leq b\\}$. Therewith,
there are two possible cases for the vertices of $\mathcal{P}_{s}$:

  * •

$\mathbf{v}$ is given by $\mathbf{C}^{\prime}\mathbf{x}=\mathbf{d}^{\prime}$,
where $\mathbf{C}^{\prime}x\leq\mathbf{d}^{\prime}$ is a subsystem of
$\mathbf{C}\mathbf{x}\leq\mathbf{d}$;

  * •

$\mathbf{v}$ is given by
$\mathbf{C}^{\prime\prime}\mathbf{x}=\mathbf{d}^{\prime\prime}$, where:

$\mathbf{C}^{\prime\prime}=\begin{bmatrix}\mathbf{C}^{\prime}\\\
\mathbf{a}^{T}\end{bmatrix},\
\mathbf{d}^{\prime\prime}=\begin{bmatrix}\mathbf{d}^{\prime}\\\
b\end{bmatrix}$ (20)

For the first case, we trivially observe that $\mathbf{v}\in\mathcal{V}_{h}$.
For the second case, as $\mathbf{C}^{\prime\prime}\in\mathbb{R}^{n\times n}$
then $\mathbf{C}^{\prime}\mathbf{x}=\mathbf{d}^{\prime}$ represents the
intersection of $n-1$ hyperplanes, resulting in an edge of $\mathcal{P}$,
given by
$\mathcal{E}_{i}=\\{\mathbf{x}\mid\mathbf{C}^{\prime}\mathbf{x}=\mathbf{d}^{\prime}\\}$.
Hence, for the case where $\mathcal{P}\cap\mathcal{H}\subset\mathcal{P}$, if
$\mathcal{V}_{p}$ is the set of vertices of $\mathcal{P}_{h}$, then
$\mathcal{V}_{p}\setminus\mathcal{V}_{h}=\mathcal{E}\cap\mathcal{H}_{s}$. ∎

###  4.4 Origin regarding polytope intersection

As presented in the previous section, considering that $\mathcal{P}$ is a
closed polytope defined by the convex combination of its vertices, the
intersection between $\mathcal{P}$ and a half-space $\mathcal{H}$ is given by
the convex combination of the vertices from $\mathcal{P}$ in the intersection
with $\mathcal{H}$, along with the vertices obtained by the intersection of
the supporting hyperplane of $\mathcal{H}$ with the edges of $\mathcal{P}$.

Consequently, the intersection between $\mathcal{P}$ and $\mathcal{O}_{j}$,
where $\mathcal{O}_{j}$ represents one of the $2^{n}$ orthants in a
$n$-dimensional space, can be rewrite as
$\mathcal{P}\cap\mathcal{H}_{1}\cap\cdots\cap\mathcal{H}_{n}$. Note that
$\mathcal{O}_{j}=\mathcal{H}_{1}\cap\cdots\cap\mathcal{H}_{n}$ for suitable
half-spaces $\mathcal{H}_{1},\dots,\mathcal{H}_{n}$. By induction, it can be
shown that the vertices of $\mathcal{P}\cap\mathcal{O}_{j}$ consist of the
union of vertices of $\mathcal{P}$ that also belong to $\mathcal{O}_{j}$ with
the vertices given by the intersection of the edges of $\mathcal{P}$ with the
supporting hyperplanes of $\mathcal{O}_{j}$ (denoted as
$\mathcal{H}_{i},\forall i\in\\{1,\ldots,n\\}$). We define $\mathcal{O}_{j}$
as:

$\mathcal{O}_{j}=\\{\mathbf{x}\mid\boldsymbol{\Phi}_{j}\mathbf{x}\leq\mathbf{0}\\}$
(21)

where $\mathbf{\Phi}_{j}$ is the suitable matrix given by:

$\mathbf{\Phi}_{j}=\begin{bmatrix}\phi_{j,1,1}&0&\ldots&0&\\\
0&\phi_{j,2,2}&\ldots&0&\\\ \vdots&\vdots&\ddots&\vdots&\\\
0&0&\ldots&\phi_{j,n,n}&\end{bmatrix}$ (22)

such that $\phi_{j,i,i}\in\\{-1,1\\},\forall i\in\\{1,\ldots,n\\}$.

We can see that the equation system $\mathbf{\Phi}_{j}\mathbf{x}=\mathbf{0}$
has the origin as its sole solution for any orthant $j$. Therefore, the unique
extreme point of $\mathcal{O}_{j}$ is the origin for all
$j\in\\{1,\ldots,2^{n}\\}$.

######  Proposition 5.

If the origin belongs to $\mathcal{P}$, then it is an extreme point of
$\mathcal{P}\cap\mathcal{O}_{j}$, for all $j\in\\{1,\ldots,2^{n}\\}$.

###### Proof.

Let $\mathcal{P}$ be a convex closed polytope and $\mathcal{O}_{j}$ a convex
cone. It is known that the intersection between $\mathcal{P}$ and
$\mathcal{O}_{j}$ is a closed convex polytope. We have that the vertices from
$\mathcal{P}$ that belong to $\mathcal{O}_{j}$ are also vertices of
$\mathcal{P}\cap\mathcal{O}_{j}$, denoted by $\mathcal{V}_{\mathcal{O}_{j}}$.
Also, we have that
$\mathcal{P}\cap\mathcal{O}_{j}=\mathcal{P}\cap\mathcal{H}_{1}\cap\cdots\cap\mathcal{H}_{n}$.

Based on Proposition 4, we have by induction that the intersection of the
edges of $\mathcal{P}$ with the supporting hyperplanes of $\mathcal{O}_{j}$
are also vertices of $\mathcal{P}\cap\mathcal{O}_{j}$, denoted by
$\mathcal{V}_{h}$.

For the origin, which is the single vertex of $\mathcal{O}_{j}$, there are two
possible cases:

  1. 1.

the origin is not in $\mathcal{P}$: this is the trivial case in which the
origin can not be a vertex of $\mathcal{P}\cap\mathcal{O}_{j}$, as it is not
in $\mathcal{P}$;

  2. 2.

the origin is in $\mathcal{P}$: in this case, by contradiction we suppose that
the origin is not a vertex of $\mathcal{P}\cap\mathcal{O}_{j}$, denoted by
$\mathbf{0}$. Then, there must exist two points
$\mathbf{v},\mathbf{u}\in\mathcal{P}\cap\mathcal{O}_{j}$, such that:

$\mathbf{0}=\lambda\mathbf{v}+(1-\lambda)\mathbf{u}$ (23)

given that $0<\lambda<1$, since $\mathbf{0}\neq\mathbf{v}$ and
$\mathbf{0}\neq\mathbf{u}$. As the origin is given by
$\mathbf{0}=(0,0,\ldots,0,0)$, there must exist a solution for the equation
$0=\lambda v_{i}+(1-\lambda)u_{i}$ for all $i\in\\{1,\ldots,n\\}$, such that:

$v_{i}=\frac{(\lambda-1)u_{i}}{\lambda}$ (24)

Since $\frac{(\lambda-1)}{\lambda}<0$, the sign of $v_{i}$ and $u_{i}$ must by
different if $v_{i}\neq 0$ and $u_{i}\neq 0$. However, inside a given orthant
there must not exist a point with components that have opposite signs. Hence,
by contradiction, if the origin is in $\mathcal{P}$, it must be a vertex of
$\mathcal{P}\cap\mathcal{O}_{j}$.

∎

###  4.5 Correctness of APNM algorithm

Let $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ denote a neural network
mapping, where $L$ is the number of layers, the weights of the layer $l$ are
given by
$\mathbf{W}_{l}\in\mathbb{R}^{|\mathbf{x}_{l}|\times|\mathbf{x}_{l-1}|}$ and
the biases by $\boldsymbol{\theta}_{l}\in\mathbb{R}^{|\mathbf{x}_{l}|}$.
Thereby, $F(\mathbf{x})=(F_{L}\circ F_{L-1}\circ\cdots\circ
F_{1})(\mathbf{x})$ where
$F_{l}(\mathbf{x})=ReLU(\mathbf{W}_{l}\mathbf{x}+\boldsymbol{\theta}_{l})$ is
the mapping for each layer $l\in\\{1,\ldots,L-1\\}$. The only difference for
the last layer is the activation (usually sigmoid for binary problems, or
softmax for multi-class problems).

For a given layer $l$, we have that the set $\mathcal{X}_{l}$ denotes the
inputs that we aim to map regarding $F_{l}$, such that
$\mathcal{X}_{l}=\\{\sum_{i=1}^{o}\lambda_{i}\mathbf{v_{i}}\mid\sum_{i=1}^{o}\lambda_{i}\land\lambda_{i}\geq
0\land\mathbf{v}_{i}\in\mathcal{V}_{l},\forall i\in\\{1,\ldots,o\\}\\}$ and
$\mathcal{V}_{l}$ is the set of vertices of $\mathcal{X}_{l}$. Associated with
$\mathcal{X}_{l}$, there is
$\mathcal{Y}_{l}=\\{F_{l}(\mathbf{x}_{l})\mid\mathbf{x}_{l}\in\mathcal{X}_{l}\\}$,
which corresponds to the output set for layer $l$.

Finally, let $\widehat{F}_{l}$ denote the output mapping for a given layer by
the Algorithm 6, given by:

$\widehat{F}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})=RP\left(ReLU\bigl{(}II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})))\bigr{)}\right)$
(25)

where $AM$ corresponds the affine map given by Algorithm 1, $EI$ is the edge
identification given by the Algorithm 2, $II$ is the intersection
identification defined by Algorithm 3, $ReLU$ is given by Algorithm 4, and
$RP$ stands for removing non-vertices defined by Algorithm 5. We denote the
convex hull mapping for a given set of vertices by $CH$.

######  Proposition 6.

Given a closed convex polytope $\mathcal{X}_{l}$ as input set and
$\mathcal{V}_{l}$ as the set of its vertices, then it implies that
$\mathcal{Y}_{l}\subseteq
CH(\widehat{F}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$.
In other words, every output of a layer $l$ associated with an input in
$\mathcal{X}_{l}$ is in the convex hull of
$\widehat{F}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$

###### Proof.

The mapping of a layer $l$ from $F$ is composed of an affine map
($\mathbf{W}_{l}\mathbf{x}_{l}+\boldsymbol{\theta}_{l}$) and a non-linear map
($ReLU$). Hence, since $\mathcal{X}_{l}$ is a closed convex polytope and
$\mathcal{V}_{l}$ its vertices, the affine map of $\mathcal{X}_{l}$ is given
by the convex hull of the affine map of its vertices, as computed by Algorithm
1.

As previously presented, the $ReLU$ map is non-linear and therefore does not
necessarily preserve the convexity of a given input set. However, as
established by Proposition 2, the $ReLU$ mapping preserves the convexity of a
convex set inside a given orthant.

Therefore, we divide the resulting set of the affine map in such a way that
each partition is inside a single orthant, so that we can apply the $ReLU$ map
to the set
$\widehat{\mathcal{X}}_{l}=CH(AM(\mathcal{V},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$.
As assured by Proposition 4, the intersection of a given orthant
$\mathcal{O}_{j}$ with the polytope $\widehat{\mathcal{X}}_{l}$ consists of
the convex hull of the union of the vertices of $\widehat{\mathcal{X}}_{l}$
that are in $\mathcal{O}_{j}$, with the vertices from the intersection of the
edges of $\widehat{\mathcal{X}}_{l}$ with the supporting hyperplanes that
define $\mathcal{O}_{j}$. Note that $\mathcal{O}_{j}$ represents a given
orthant, for all $j\in\\{1,\ldots,2^{|\mathbf{x}_{l}|}\\}$.

Thus, we first need to compute the edges of $\widehat{\mathcal{X}}_{l}$,
calculated by $EI(AM(\mathcal{V},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$, as
stated by Proposition 1, followed by the determination of the intersection of
these edges with the supporting hyperplanes of $\mathcal{O}_{j}$, given by
$II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})))$
and computed by Algorithm 3.

Therewith, we have that $\widehat{\mathcal{X}}_{l}$ was divided in such a way
that each partition is inside a single orthant. Then, the $ReLU$ mapping is
applied to the vertices of each partition of $\widehat{\mathcal{X}}_{l}$,
resulting in
$ReLU(II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))))$.
Proposition 2 assures that the $ReLU$ mapping preserves the convexity inside a
single orthant, which allows its previous application. Note that all the
vertices will be in the non-negative orthant after applying the $ReLU$
mapping.

Finally, we compute the convex hull of all the output sets of the $ReLU$
mapping (at most $2^{|\mathbf{x}_{l}|}$) by removing those points that are not
vertices of
$CH(ReLU(II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})))))$.
This final step is computed by Algorithm 5 (RP) as stated in Equation (25).

Consequently,
$CH(\widehat{F}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$
is the convex hull of the exact map $\mathcal{Y}_{l}$ of $\mathcal{X}_{l}$. ∎

We denote by $\widehat{F}$ the composition of $\widehat{F}_{l}$ for all layers
of the neural network, except the last one, where the non-linear map is not
applied, given by
$\widehat{F}(\mathcal{V},\mathbf{W},\boldsymbol{\theta})=(\widehat{F}_{L}\circ\widehat{F}_{L-1}\circ\cdots\circ\widehat{F}_{1})(\mathcal{V},\mathbf{W},\boldsymbol{\theta})$.
Furthermore, we define
$\mathcal{Y}=\\{F(\mathbf{x})\mid\mathbf{x}\in\mathcal{X}\\}$ as the exact
output set of the network, regarding the closed convex input polytope
$\mathcal{X}$ and its corresponding set of vertices $\mathcal{V}$.

######  Proposition 7.

Given a closed convex polytope $\mathcal{X}$ as the input set and
$\mathcal{V}$, the set of its vertices, then we have that
$\mathcal{Y}\subseteq
CH(\widehat{F}(\mathcal{V},\mathbf{W},\boldsymbol{\theta}))$. Put in different
terms, each output of the network associated with an input in $\mathcal{X}$ is
in $\widehat{F}(\mathcal{V},\mathbf{W},\boldsymbol{\theta})$.

###### Proof.

As established by Proposition 6, $\mathcal{Y}_{l}\subseteq
CH(\widehat{F}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$
for a given layer $l$ of the neural network $F$. Then, for $l=1$, it follows
that:

$\mathcal{Y}_{1}\subseteq
CH(\widehat{F}_{1}(\mathcal{V}_{1},\mathbf{W}_{1},\boldsymbol{\theta}_{1}))$
(26)

where $\mathcal{V}_{1}$ is the set of vertices from $\mathcal{X}_{1}$ and
$\mathcal{X}_{1}=\mathcal{X}$. As the output of the first layer is the input
of the second one, we have that
$\mathcal{X}_{2}=CH(\widehat{F}_{1}(\mathcal{V}_{1},\mathbf{W}_{1},\boldsymbol{\theta}_{1}))$
and that
$\mathcal{V}_{2}=\widehat{F}_{1}(\mathcal{V}_{1},\mathbf{W}_{1},\boldsymbol{\theta}_{1})$.

For $l=2$:

$\mathcal{Y}_{2}\subseteq
CH(\widehat{F}_{2}(\mathcal{V}_{2},\mathbf{W}_{2},\boldsymbol{\theta}_{2}))$
(27)

Then, by replacing $\mathcal{V}_{2}$ with the output of layer 1, results in:

$\mathcal{Y}_{2}\subseteq
CH(\widehat{F}_{2}(\widehat{F}_{1}(\mathcal{V}_{1},\mathbf{W}_{1},\boldsymbol{\theta}_{1}),\mathbf{W}_{2},\boldsymbol{\theta}_{2}))$
(28)

Now, for layers $k$ and $k+1$, it follows by induction that:

$\mathcal{Y}_{k+1}\subseteq
CH(\widehat{F}_{k+1}(\widehat{F}_{k}(\mathcal{V}_{k},\mathbf{W}_{k},\boldsymbol{\theta}_{k}),\mathbf{W}_{k+1},\boldsymbol{\theta}_{k+1}))$
(29)

Consequently, the mapping $\hat{F}$ in fact computes an over-approximation for
the exact output set $\mathcal{Y}$. ∎

###  4.6 Correctness of EPNM algorithm

The set $\mathcal{Y}_{l}$ denotes the exact output associated with the input
set $\mathcal{X}_{l}$, regarding the layer $l$ of the neural network $F$. The
mapping of a given layer $l$ implemented by Algorithm 11, denoted here by
$\widehat{E}_{l}$, is stated as:

$\widehat{E}_{l,k}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})=ReLU\left(SP\left(OS\bigl{(}II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})))\bigr{)}\right)_{k}\right)$
(30)

where $AM$ is the affine map computed by Algorithm 1, $EI$ is the edge
identification implemented by Algorithm 2, $II$ is the intersection
identification given by Algorithm 3, $ReLU$ is computed by Algorithm 4, $OS$
is implemented by Algorithm 7 to verify if the origin is in the polytope, and
finally $SP$, which separates the vertices in orthants, is computed by
Algorithm 10. Notice that $k\in\\{1,\ldots,K\\}$ is the index that represents
each
$\widehat{E}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$. We
denote the convex hull mapping for a given input set of vertices by $CH$.

######  Proposition 8.

Given a convex closed polytope $\mathcal{X}_{l}$ as input set and
$\mathcal{V}_{l}$, the set of its vertices, we have that
$\mathcal{Y}_{l}=\bigcup_{k=1}^{K}CH(\widehat{E}_{l,k}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$,
where $K$ is the number of sets that comprise
$\widehat{E}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$. Put
another way, the set of outputs of the layer $l$, resulting from all inputs in
$\mathcal{X}_{l}$, consists of the union of the convex hull of each set
$\widehat{E}_{l,k}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$.

###### Proof.

As presented previously, the mapping of a given layer $l$ of the neural
network $F$ is a composition of two different functions: one affine mapping
($\mathbf{W}_{l}\mathbf{x}_{l}+\boldsymbol{\theta}_{l}$) with one non-linear
mapping ($ReLU$). As $\mathcal{X}_{l}$ is a closed convex polytope, the affine
mapping is obtained by the convex hull of the affine map of each of its
vertices, implemented by Algorithm 1.

$ReLU$ does not necessarily preserve the convexity of a given input set.
However, as shown by Proposition 2, the $ReLU$ mapping preserves the convexity
of an input set if the input is inside a single orthant.

Note that the affine mapping of the input set is computed by
$AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$, where
$\widehat{\mathcal{X}}_{l}=CH(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$
denotes the application of the affine mapping in $\mathcal{X}_{l}$. Therefore,
it is necessary to split $\widehat{\mathcal{X}}_{l}$ in such a way that each
partition lies inside a single orthant, so that it becomes possible to apply
the $ReLU$ mapping to the vertices of $\widehat{\mathcal{X}}_{l}$. As assured
by Proposition 4, the intersection of an orthant $\mathcal{O}_{j}$ with a
polytope $\widehat{\mathcal{X}}_{l}$ consists of the union of the vertices of
$\widehat{\mathcal{X}}_{l}$ that are in $\mathcal{O}_{j}$, with the vertices
in the intersection of the edges of $\widehat{\mathcal{X}}_{l}$ with the
supporting hyperplanes of $\mathcal{O}_{j}$ and the origin, if the latter lies
inside $\widehat{\mathcal{X}}_{l}$, according to Proposition 5.

Thus, we firstly compute the edges of $\widehat{\mathcal{X}}_{l}$, a step
denoted by $EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}))$
and computed by Algorithm 2, as stated by Proposition 1. Then, the
intersection of these edges with the supporting hyperplanes of orthant
$\mathcal{O}_{j}$ is obtained with Algorithm 3, and finally Algorithm 7
verifies whether the origin belongs to $\widehat{\mathcal{X}}_{l}$. The result
of such an operation is given by:

$OS\bigl{(}II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})))\bigr{)}$
(31)

In the next step, we separate the vertices of $\widehat{\mathcal{X}}_{l}$ in
different sets, such that those vertices in the same set represent the portion
of $\widehat{\mathcal{X}}_{l}$ that is inside a single orthant. This process
takes place to enable the application of the $ReLU$ mapping, as the separation
allows the non-linear mapping to be applied in each partition of
$\widehat{\mathcal{X}}_{l}$ while ensuring convexity, as established by
Proposition 2. Algorithm 10 performs the partitioning operation, given by:

$SP\left(OS\bigl{(}II(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l}),EI(AM(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})))\bigr{)}\right)$
(32)

The result of $SP$ is the set of sets of vertices, where the convex hull of
each set represents a partition of $\widehat{\mathcal{X}}_{l}$ restricted to a
single orthant. Thus,
$\widehat{E}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$
denotes the result of the application of the $ReLU$ mapping to each of these
sets. Finally, we have that
$\widehat{E}_{l}(\mathcal{V}_{l},\mathbf{W}_{l},\boldsymbol{\theta}_{l})$
represents the exact mapping of $\mathcal{X}_{l}$, as we applied both the
affine and the non-linear mapping without any over-approximation.
Consequently, Algorithm 11 in fact returns the exact output set for a given
layer $l$ of the neural network $F$. ∎

Finally, we denote by $\widehat{E}$ the composition of $\widehat{E}_{l}$, for
each $l\in\\{1,\ldots,L\\}$. To avoid repetition, the Proposition 9 is not
presented. However, it follows the same inductive process as presented for the
APNM (Proposition 7).

######  Proposition 9.

Given a closed convex polytope $\mathcal{X}$ as input set and $\mathcal{V}$,
the set of its vertices, we have that
$\mathcal{Y}=\bigcup_{k=1}^{K}CH(\widehat{E}_{k}(\mathcal{V},\mathbf{W},\boldsymbol{\theta}))$,
where $K$ is the number of sets in
$\widehat{E}(\mathcal{V},\mathbf{W},\boldsymbol{\theta})$. In other terms, the
set of outputs of the neural network $F$, associated with the input set
$\mathcal{X}$, is in the union of the convex hull of each set
$\widehat{E}_{k}(\mathcal{V},\mathbf{W},\boldsymbol{\theta})$.

##  5 Application

In this section we present a comparative analysis between the proposed vertex-
based reachability approach and representative algorithms from the literature.
This comparison is carried out by verifying one of the properties from ACAS XU
[1].

###  5.1 ACAS XU

The Aircraft Collision Avoidance System (ACAS XU) [1] comprises a set of fully
connected neural networks that aim to eliminate the possibility of collisions
between two aircraft. The system comprises 45 trained models, where each model
receives as inputs five properties of the ownship and the intruder (the
aircraft that is invading the space of the ownship): the distance from the
ownship to the intruder ($\rho$), the angle to the intruder regarding the
ownship heading direction ($\theta$), the heading angle of intruder relative
to ownship heading direction ($\psi$), the speed of the ownship
($v_{ownship}$), and the speed of the intruder ($v_{intruder}$). We can see a
representation of the inputs in Figure 14.

Figure 14: Representation of the ACAS XU inputs. The black dot represents the
position of the ownship and the red dot the position of the intruder. $\rho$
represents the Euclidean distance between both aircrafts, $\theta$ the angle
between the ownship heading direction and the vector that connects both
aircraft, and $\psi$ the angle between the ownship and the intruder heading
direction.

There are two extra parameters that are not used as inputs to the neural
network. The first one is the time until loss of vertical separation ($\tau$),
whereas the second one is the previous prediction advice ($a_{previous}$).
These $\tau$ and $a_{previous}$ parameters are discretized such that for each
possible combination of values for $\tau$ and $a_{previous}$, a different
model is trained. This process resulted in a total of 45 trained models, as
mentioned before. Each trained model associates an input pattern to five
possible categories: clear of conflict ($y_{0}$), weak left ($y_{1}$), weak
right ($y_{1}$), strong left ($y_{3}$), and strong right ($y_{4}$). Figure 15
contains a simple representation of a single the neural network classifier
associated with a value for $\tau$ and $a_{previous}$. Each neural network
comprises 6 hidden layers, containing 50 neurons each. For further details on
the training and prediction process of these models, please refer to [1].

Figure 15: Representation of the inputs and the outputs of a single neural
network from ACAS XU. On the left of the neural network we present each of the
five expected inputs: $\rho$, $\theta$, $\psi$, $\mathbf{v}_{own}$ and
$\mathbf{v}_{int}$. The neural network is defined according to a given $\tau$
and a $a_{prev}$. To the right of the neural network we have the expected
probability for each output class. $CC$ stands for clear of conflict, $WL$ for
weak left, $WR$ for weak right, $SL$ for strong left, and $SR$ for strong
right.

Based on the trained models of the ACAS XU, ten desired properties have been
proposed so that this system works correctly according to its crafted design.
In this paper, to compare with existing verification approaches, we verified
Property 1, formally stated as:

######  Property 1.

The conditions are established as follows:

  * •

input constraint: $\rho\geq 55947,691$, $v_{ownship}\geq 1145$ and
$v_{intruder}\leq 60$;

  * •

output constraint: $y_{0}\leq 1500$;

where $y_{0}$ is the output associated with the clear of conflict output
class.

Recall from the problem statement that the reachability analysis aims to
verify a reachable set $\mathcal{R}$, obtained from $\mathcal{X}$, such that
$\mathcal{R}\cap\lnot\mathcal{Y}=\emptyset$ holds. $\mathcal{Y}$ denotes the
expected output for the inputs in $\mathcal{X}$. For the Property 1, we have
that:

$\displaystyle\mathcal{X}=\\{(\rho,\theta,\phi,v_{ownship},v_{intruder})\mid$
$\displaystyle 55947.691\leq\rho\leq 60760,$ (33)
$\displaystyle-\pi\leq\theta\leq\pi,$ $\displaystyle-\pi\leq\phi\leq\pi,$
$\displaystyle 1145\leq v_{ownship}\leq 1200,$ $\displaystyle 0\leq
v_{intruder}\leq 60\\}$

and that:

$\displaystyle\mathcal{Y}=\\{(y_{0},y_{1},y_{2},y_{3},y_{4})\mid y_{0}\leq
1500\\}$ (34)

Property 1 is the only one presented in this paper because only this property
is evaluated by means of comparison with other formal verification algorithms.
However, all of the remaining properties are formally described in [10].

###  5.2 Experimental description

We present in this section the procedures of our experimental setup. The
experimental results are divided into two main parts: the first part aiming to
validate and compare the results with algorithms from the literature, and the
second part to evaluate the features of our approaches. For the first part, we
conducted the experiments by:

  1. 1.

Generate the vertices of the input polytope, based on the input constraint
imposed by Property 1;

  2. 2.

Set a timeout of 24 hours for each algorithm verification;

  3. 3.

Compute the output reachable set and the verification status for each
reachability algorithm (Exact polytope network mapping (EPNM), MaxSens [12],
Ai2 [13] and ExactReach [11]);

  4. 4.

Compute counterexamples and the verification status for the search algorithm
(Reluplex [10], Duality [16], MIPVerify [15] and NSVerify [14]); and

  5. 5.

Compare results.

For the second experiment, which aims to analyze the parallelism behavior of
the algorithm, the experimental setup consists of:

  1. 1.

Generate the vertices of the input polytope, based on the input constraint
imposed by Property 1;

  2. 2.

Set the number of parallel processes, denoted by $p$, such that
$p\in\\{1,4,8,12,16,20,24,28,32\\}$; and

  3. 3.

Verify Property 1 for each $p$.

###  5.3 Hardware and software specification

For comparative matters, we provide the specification for both hardware
resources and software language. The comparative experiment previously
presented was performed in an Intel Xeon CPU E5-2630 V4 of 2.20GHz, with 40
available CPU. Those algorithms from the literature and the algorithms
proposed in this work were developed in Julia language. The implementation of
the algorithms from the literature were available on [23]. The implementation
of our algorithms is available on [24].

###  5.4 Comparative results

The validation and comparative results of the verification of ACAS Xu models
for Property 1 are presented in two different perspectives: firstly among
those verification procedures that follow reachability approaches, then
comparing with approaches that make use of different techniques (search or
optimization). Finally, we present some useful features of the proposed
approach.

####  5.4.1 EPNM versus reachability approaches

By comparing the EPNM approach with those verification algorithms that follow
a reachability approach, as presented in Table 1, it can be seen that the
proposed exact approach verified most of the neural networks within the
stipulated timeout time (43 out of 45 neural networks). As we can see, none of
the other existing exact approaches were able to verify a single model within
a day of execution (24 hours). Notice also that the approximate approaches
(MaxSens and Ai2) finished their execution, though, due to their over-
approximation, these procedures did not estimate the correct status well
(which is acceptable, as these approaches are sound but not complete).

Table 1: Comparative results between the proposed approach (EPNM) regarding existing reachability approaches from the literature. These results refer to the verification of Property 1 of ACAS XU models. | Proposed | Approaches from the Literature
---|---|---
Status | EPNM | ExactReach | MaxSens | Ai2
holds | 43 | - | - | -
violated | - | - | 45 | 45
timeout | 2 | 45 | - | -

####  5.4.2 EPNM versus search and optimization approaches

In comparison to optimization and search approaches, the proposed EPNM
approach also reached interesting results. Compared to Reluplex results from
the literature, EPNM could verify more neural networks within the specified
timeout. However, by executing Reluplex in the same hardware conditions of
EPNM, the verification was not completed within 24 hours. The same occurred
for the NSVerify procedure.

Table 2: Comparative results between the proposed approach regarding existing search and optimization approaches from the literature. These results refer to the verification of Property 1 of ACAS XU models. | Proposed | Approaches from the Literature
---|---|---
Status | EPNM | Reluplex111 | Reluplex (24 hours)222 | Duality | NSVerify
holds | 43 | 41 | - | - | -
violated | - | - | - | - | -
timeout | 2 | 4 | 45 | - | 45
unknown | - | - | - | 45 | -
11footnotetext: Results extracted from [10].22footnotetext: Results from the
authors for a $24$-hour execution.

####  5.4.3 Parallel Computation

Due to the characteristics of the proposed approaches, their implementation
allows the parallelization in a procedural level. We chose to implement the
parallelization in two of the procedures (EI and II), because the remaining
algorithms did not respond well due to the tradeoff between the overhead and
the speed-up of the parallelization.

The first one is the edge identification (EI) algorithm. For this procedure,
we created a pool for the execution with the size equal to the number of
available threads. The pool guarantees that, after each thread ends its
execution, a new thread is started and takes the empty space. For each vertex,
a new thread was initiated, which verified if the adjacency property holds for
the current and each of the remaining vertices. As sharing memory was not
necessary, because each thread has its own adjacency list, no synchronization
approach was implemented. By the end of the execution, the algorithm
concatenated the adjacency list calculated for each vertex.

The intersection identification (II) was the second parallelized procedure.
Following the same idea, a new thread was created for each vertex. After the
initialization is completed, the procedure identifies those intersection
points between the current and each of its adjacent vertices. In this case,
similarly to the previous one, no synchronization was necessary, as each
thread carried its own list of intersections. At the end of the pool
execution, the procured concatenated those intersection points associated with
each vertex.

The results of the second part of the experiments are depicted in Figure
16(a). As can be seen in this figure, as the number of available threads for
the execution increases, the running time decreases significantly. This
characteristic of the algorithm can be explored for huge problems.

Figure 16(b) presents the speedup behavior for the algorithm EPNM. As
expected, the algorithm indeed reduce the runtime as there is an increment on
the available threads, though the difference between the real and the ideal
curve indicates that the parallel processes are not ideally balanced (there
are threads waiting for some execution to end).

(a) The graph illustrates the execution speed-up obtained from parallelizing
the proposed algorithm. As the number of available threads increases, there is
a significant reduction in the running time (duration).

(b) Comparison between the actual speed up of the EPNM and the ideal speed up.

Figure 16: Illustrative visualization of the EPNM parallelization behavior. We
present both, the runtime and the speed up for the algorithm EPNM execution on
a single ACAS Xu model for the property 1.

(a) Total of number vertices processed at each layer of a single neural
network. The y-axis is in logarithmic scale. The APNM algorithm presents a
significantly higher increment in the number of vertices after layer $3$, in
comparison to PAPNM and EPNM. Both PAPNM and EPNM executions have the same
value at layer $2$, as expected, though EPNM has a lower number of vertices at
layer 3.

(b) Total number of sets generated at each layer processed by each of the
algorithms. As expected, APNM keeps $1$ set at the end of every layer
execution. The EPNM has the greatest increment on the number of sets, which is
an expected behavior.

Figure 17: Illustration of the complexity behavior of APNM, EPNM and PAPNM
with $d=2$ and $d=4$ for the first $3$ layers of a single neural network from
ACAS XU model.

####  5.4.4 Complexity behavior of the proposed approaches

Our experiments showed that, differently from the expected, the algorithm EPNM
has a shorter running time in comparison to APNM, for the ACAS XU model. This
behavior can be explained considering the total number of vertices that are
processed in both algorithms. Figure 17 presents the behavior of both, the
total number of vertices and sets processes at each layer of a single neural
network from ACAS XU.

Table 3: Results on the running time of EPNM, APNM and PAPNM with $d=2$ and
$d=4$ for the first three layers of a single neural network.

(a) Behavior of the total number of vertices processed for each layer of a neural network by APNM, EPNM and PAPNM. Layer | APNM | PAPNM | PAPNM | EPNM
---|---|---|---|---
|  | $(d=2)$ | $(d=4)$ |
1 | 32 | 32 | 32 | 32
2 | 323 | 580 | 580 | 580
3 | 313286 | 85404 | 21977 | 2502

(b) Behavior of the total number of sets processed for each layer of a neural network by APNM, EPNM and PAPNM. Layer | APNM | PAPNM | PAPNM | EPNM
---|---|---|---|---
|  | $(d=2)$ | $(d=4)$ |
1 | 1 | 1 | 1 | 1
2 | 1 | 11 | 21 | 42
3 | 1 | 78 | 117 | 149

Figure 17(a) shows that, from the same input set, the algorithm APNM generates
the greater set of vertices for representing its approximation of the output
set compared to both, EPNM and PAPNM (i.e., $313286$ vertices at the third
layer). Figure 17(b) depicts the increment on the number of sets for EPNM and
PAPNM. Table 3(b) reports the data used to create Figure 17.

Both results lead to a lower average number of vertices across each of the
sets for EPNM. On the other hand, the opposite occurs to APNM which has a
single set with a strongly increasing number of vertices. Table 4 shows the
average number of vertices per set for each algorithm at each layer of the
neural network. This means that the merging process (performed partially by
PAPNM and completely by APNM) induces a simplification on the total number of
sets, though as a side effect it significantly increases the total number of
vertices to be processed.

Table 4: Average of the total number of vertices within each set for each algorithm. The results are presented layer by layer. Layer | APNM | PAPNM | PAPNM | EPNM
---|---|---|---|---
|  | $(d=2)$ | $(d=4)$ |
1 | 32 | 32 | 32 | 32
2 | 323 | 52.7 | 27.6 | 13.8
3 | 313286 | 1094.9 | 187.8 | 16.8

The number of vertices within a set is directly related to the running time
for each algorithm because the computational complexity of each of the
procedures that comprise APNM, PAPNM and EPNM is directly related to the total
of vertices processed.

##  6 Conclusion

In this work, we proposed two vertex-based reachability algorithms for formal
verification of deep neural networks. These algorithms compute a reachable
output set, which may consist of a set of polyhedral sets, for a given input
polyhedral set, satisfying different properties: the first one (APNM) computes
an approximation for the output reachable set, while the second one (EPNM)
computes the exact output reachable set.

Supported by formal demonstrations of correctness, the proposed algorithms
were shown to correctly verify properties of neural networks with the $ReLU$
activation function. More specifically, it was shown that APNM yields an
overestimation of the output reachable set, while EPNM computes the exact
reachable set.

Our proposal was applied to a benchmark problem for neural network
verification and compared to some of the algorithms previously proposed in the
literature. The results showed that among the verification algorithms that
make use of reachability analysis, the presented EPNM approach concluded most
of the verifications (43 out of 45 neural networks), differently from the
ExactReach, which is another exact approach from the literature that could not
verify any neural network within the specified timeout. Compared to those
algorithms that are not complete, despite the fact that these approaches were
able to finish their executions, the expected output was not reached in any
case (see Table 1).

In comparison to those methods that make use of optimization and search
strategies, the result reported in the literature for Reluplex surpassed ours
in terms of running time. However, under the same hardware conditions, our
approach was able to overcome their results.

Finally, we argue that the algorithms proposed in this paper are strong
candidates for the verification of neural networks with $ReLU$ activation.
Additionally, the running time can be reduced drastically by using multiple
processing cores. Future work includes the investigation of a different
construction for the search approach in EPNM and of a third approach that can
be designed to improve performance by using heuristics for the reduction of
sets and vertices during the search process.

## Acknowledgment

This work was funded in part by Fundação de Amparo à Pesquisa e Inovação do
Estado de Santa Catarina (FAPESC) under grant 2021TR2265.

## References

  * [1] K. D. Julian, M. J. Kochenderfer, and M. P. Owen, “Deep neural network compression for aircraft collision avoidance systems,” Journal of Guidance, Control, and Dynamics, vol. 42, no. 3, pp. 598–608, 2019.
  * [2] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. J. Goodfellow, and R. Fergus, “Intriguing properties of neural networks,” CoRR, vol. abs/1312.6199, 2014.
  * [3] I. J. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” CoRR, vol. abs/1412.6572, 2015.
  * [4] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, “Towards deep learning models resistant to adversarial attacks,” in International Conference on Learning Representations, ICLR 2018, 2018.
  * [5] F. Tramèr, A. Kurakin, N. Papernot, I. Goodfellow, D. Boneh, and P. McDaniel, “Ensemble adversarial training: Attacks and defenses,” in International Conference on Learning Representations, ICLR 2017, 2017.
  * [6] Y. Song, T. Kim, S. Nowozin, S. Ermon, and N. Kushman, “Pixeldefend: Leveraging generative models to understand and defend against adversarial examples,” CoRR, vol. abs/1710.10766, 2017.
  * [7] K. Grosse, P. Manoharan, N. Papernot, M. Backes, and P. McDaniel, “On the (statistical) detection of adversarial examples,” CoRR, vol. abs/1702.06280, 2017.
  * [8] J. H. Metzen, T. Genewein, V. Fischer, and B. Bischoff, “On detecting adversarial perturbations,” in Proceedings of 5th International Conference on Learning Representations (ICLR), 2017.
  * [9] J. Lu, T. Issaranon, and D. Forsyth, “Safetynet: Detecting and rejecting adversarial examples robustly,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 446–454, 2017.
  * [10] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer, “Reluplex: An efficient SMT solver for verifying deep neural networks,” in International Conference on Computer Aided Verification, pp. 97–117, Springer, 2017.
  * [11] W. Xiang, H.-D. Tran, and T. T. Johnson, “Reachable set computation and safety verification for neural networks with ReLU activations,” arXiv preprint arXiv:1712.08163, 2017.
  * [12] W. Xiang, H.-D. Tran, and T. T. Johnson, “Output reachable set estimation and verification for multilayer neural networks,” IEEE Transactions on Neural Networks and Learning Systems, vol. 29, no. 11, pp. 5777–5783, 2018.
  * [13] T. Gehr, M. Mirman, D. Drachsler-Cohen, P. Tsankov, S. Chaudhuri, and M. Vechev, “Ai2: Safety and robustness certification of neural networks with abstract interpretation,” in IEEE Symposium on Security and Privacy (SP), pp. 3–18, IEEE, 2018.
  * [14] A. Lomuscio and L. Maganti, “An approach to reachability analysis for feed-forward relu neural networks,” arXiv preprint arXiv:1706.07351, 2017\.
  * [15] V. Tjeng, K. Xiao, and R. Tedrake, “Evaluating robustness of neural networks with mixed integer programming,” arXiv preprint arXiv:1711.07356, 2017\.
  * [16] K. Dvijotham, R. Stanforth, S. Gowal, T. A. Mann, and P. Kohli, “A dual approach to scalable verification of deep networks,” in the Conference on Uncertainty in Artificial Intelligence (UAI), vol. 1, p. 3, 2018.
  * [17] E. Wong and Z. Kolter, “Provable defenses against adversarial examples via the convex outer adversarial polytope,” in International Conference on Machine Learning, pp. 5286–5295, PMLR, 2018.
  * [18] X. Huang, M. Kwiatkowska, S. Wang, and M. Wu, “Safety verification of deep neural networks,” in International Conference on Computer Aided Verification, pp. 3–29, Springer, 2017.
  * [19] P. McMullen and E. Schulte, Abstract regular polytopes, vol. 92.  Cambridge University Press, 2002.
  * [20] I. Z. Emiris, V. Fisikopoulos, and B. Gärtner, “Efficient edge-skeleton computation for polytopes defined by oracles,” Journal of Symbolic Computation, vol. 73, pp. 139–152, 2016.
  * [21] A. Schrijver et al., Combinatorial Optimization: Polyhedra and Efficiency, vol. 24.  Springer, 2003.
  * [22] A. Schrijver, Theory of Linear and Integer Programming.  John Wiley & Sons, 1998.
  * [23] C. Liu, T. Arnon, C. Lazarus, C. Strong, C. Barrett, M. J. Kochenderfer, et al., “NeuralVerification.jl,” 2019.
  * [24] J. Zago, E. Camponogara, and E. Antonelo, “vertexBasedRechabilityAnalysis.jl,” 2023.