\documentclass[journal]{IEEEtran}

\usepackage[utf8]{inputenc} %
\usepackage[T1]{fontenc}    %
\usepackage{hyperref}       %
\usepackage{url}            %
\usepackage{booktabs}       %
\usepackage{amsfonts}       %
\usepackage{nicefrac}       %
\usepackage{microtype}      %
\usepackage{amsmath, amsfonts}       %
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{bbm}
\usepackage{lipsum}
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage{lipsum} 
\usepackage[noadjust]{cite}

\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\N}{\mathbb{N}}

\newcommand{\orderof}[1]{\mathcal{O}\left(#1\right)}

\renewcommand{\Re}[1]{\operatorname{Re}\left\{#1\right\}}
\renewcommand{\Im}[1]{\operatorname{Im}\left\{#1\right\}}
\newcommand{\conj}[1]{\mkern 1.5mu\overline{\mkern-1.5mu#1\mkern-1.5mu}\mkern 1.5mu}

\renewcommand{\P}[1]{\operatorname{P}\left(#1\right)}
\newcommand{\E}{\operatorname{E}}
\newcommand{\var}{\operatorname{var}}
\newcommand{\cov}{\operatorname{cov}}
\newcommand{\normal}{\mathcal{N}}

\renewcommand{\d}[1]{d#1}

\newcommand{\e}{e}
\renewcommand{\j}{j}

\newcommand{\vct}[1]{\boldsymbol{#1}}
\newcommand{\mtx}[1]{\boldsymbol{#1}}
\newcommand*{\vertbar}{\rule[-1ex]{0.5pt}{2.5ex}}
\newcommand*{\horzbar}{\rule[.5ex]{2.5ex}{0.5pt}}
\newcommand{\bvct}[1]{\mathbf{#1}}
\newcommand{\bmtx}[1]{\mathbf{#1}}
\newcommand{\<}{\langle}
\renewcommand{\>}{\rangle}
\renewcommand{\H}{\mathrm{H}}
\newcommand{\T}{\mathrm{T}}
\newcommand{\pinv}{\dagger}
\newcommand{\Null}{\operatorname{Null}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Span}{\operatorname{Span}}
\newcommand{\trace}{\operatorname{trace}}
\newcommand{\rank}{\operatorname{rank}}

\newcommand{\set}[1]{\mathcal{#1}}
\newcommand{\closure}{\operatorname{cl}}	%
\newcommand{\interior}{\operatorname{int}}
\newcommand{\boundary}{\operatorname{bd}}
\newcommand{\diameter}{\operatorname{diam}}

\newcommand{\domain}{\operatorname{dom}}
\newcommand{\epigraph}{\operatorname{epi}}
\newcommand{\hypograph}{\operatorname{hypo}}
\newcommand{\linop}[1]{\mathscr{#1}}	%

\DeclareMathOperator*{\minimize}{\text{minimize}}
\DeclareMathOperator*{\maximize}{\text{maximize}}
\newcommand{\argmin}[1]{\underset{#1}{\operatorname{arg}\,\operatorname{min}}\;} %
\newcommand{\argmax}[1]{\underset{#1}{\operatorname{arg}\,\operatorname{max}}\;} %
\newcommand{\va}{\vct{a}}
\newcommand{\vb}{\vct{b}}
\newcommand{\vc}{\vct{c}}
\newcommand{\vd}{\vct{d}}
\newcommand{\ve}{\vct{e}}
\newcommand{\vf}{\vct{f}}
\newcommand{\vg}{\vct{g}}
\newcommand{\vh}{\vct{h}}
\newcommand{\vi}{\vct{i}}
\newcommand{\vj}{\vct{j}}
\newcommand{\vk}{\vct{k}}
\newcommand{\vl}{\vct{l}}
\newcommand{\vm}{\vct{m}}
\newcommand{\vn}{\vct{n}}
\newcommand{\vo}{\vct{o}}
\newcommand{\vp}{\vct{p}}
\newcommand{\vq}{\vct{q}}
\newcommand{\vr}{\vct{r}}
\newcommand{\vs}{\vct{s}}
\newcommand{\vt}{\vct{t}}
\newcommand{\vu}{\vct{u}}
\newcommand{\vv}{\vct{v}}
\newcommand{\vw}{\vct{w}}
\newcommand{\vx}{\vct{x}}
\newcommand{\vy}{\vct{y}}
\newcommand{\vz}{\vct{z}}
\newcommand{\valpha}{\vct{\alpha}}
\newcommand{\vbeta}{\vct{\beta}}
\newcommand{\vdelta}{\vct{\delta}}
\newcommand{\vepsilon}{\vct{\epsilon}}
\newcommand{\vgamma}{\vct{\gamma}}
\newcommand{\vlambda}{\vct{\lambda}}
\newcommand{\vmu}{\vct{\mu}}
\newcommand{\vnu}{\vct{\nu}}
\newcommand{\vphi}{\vct{\phi}}
\newcommand{\vpsi}{\vct{\psi}}
\newcommand{\vsigma}{\vct{\sigma}}
\newcommand{\vtau}{\vct{\tau}}
\newcommand{\vtheta}{\vct{\theta}}
\newcommand{\vzero}{\vct{0}}
\newcommand{\vone}{\vct{1}}
\newcommand{\mA}{\mtx{A}}
\newcommand{\mB}{\mtx{B}}
\newcommand{\mC}{\mtx{C}}
\newcommand{\mD}{\mtx{D}}
\newcommand{\mE}{\mtx{E}}
\newcommand{\mF}{\mtx{F}}
\newcommand{\mG}{\mtx{G}}
\newcommand{\mH}{\mtx{H}}
\newcommand{\mJ}{\mtx{J}}
\newcommand{\mK}{\mtx{K}}
\newcommand{\mL}{\mtx{L}}
\newcommand{\mM}{\mtx{M}}
\newcommand{\mN}{\mtx{N}}
\newcommand{\mO}{\mtx{O}}
\newcommand{\mP}{\mtx{P}}
\newcommand{\mQ}{\mtx{Q}}
\newcommand{\mR}{\mtx{R}}
\newcommand{\mS}{\mtx{S}}
\newcommand{\mT}{\mtx{T}}
\newcommand{\mU}{\mtx{U}}
\newcommand{\mV}{\mtx{V}}
\newcommand{\mW}{\mtx{W}}
\newcommand{\mX}{\mtx{X}}
\newcommand{\mY}{\mtx{Y}}
\newcommand{\mZ}{\mtx{Z}}
\newcommand{\mDelta}{\mtx{\Delta}}
\newcommand{\mLambda}{\mtx{\Lambda}}
\newcommand{\mPhi}{\mtx{\Phi}}
\newcommand{\mPsi}{\mtx{\Psi}}
\newcommand{\mSigma}{\mtx{\Sigma}}
\newcommand{\mUpsilon}{\mtx{\Upsilon}}
\newcommand{\mId}{{\bf I}}
\newcommand{\mEx}{{\bf J}}
\newcommand{\mzero}{{\bf 0}}
\newcommand{\mone}{{\bf 1}}

\newcommand{\mAbar}{\underline{\mtx{A}}}
\newcommand{\mRbar}{\underline{\mtx{R}}}
\newcommand{\vebar}{\underline{\vct{e}}}
\newcommand{\vxbar}{\underline{\vct{x}}}
\newcommand{\vybar}{\underline{\vct{y}}}
\newcommand{\loF}{\linop{F}}
\newcommand{\setA}{\set{A}}
\newcommand{\setB}{\set{B}}
\newcommand{\setC}{\set{C}}
\newcommand{\setD}{\set{D}}
\newcommand{\setE}{\set{E}}
\newcommand{\setF}{\set{F}}
\newcommand{\setG}{\set{G}}
\newcommand{\setH}{\set{H}}
\newcommand{\setI}{\set{I}}
\newcommand{\setJ}{\set{J}}
\newcommand{\setK}{\set{K}}
\newcommand{\setL}{\set{L}}
\newcommand{\setM}{\set{M}}
\newcommand{\setN}{\set{N}}
\newcommand{\setO}{\set{O}}
\newcommand{\setP}{\set{P}}
\newcommand{\setQ}{\set{Q}}
\newcommand{\setR}{\set{R}}
\newcommand{\setS}{\set{S}}
\newcommand{\setT}{\set{T}}
\newcommand{\setU}{\set{U}}
\newcommand{\setV}{\set{V}}
\newcommand{\setW}{\set{W}}
\newcommand{\setX}{\set{X}}
\newcommand{\setY}{\set{Y}}
\newcommand{\setZ}{\set{Z}}

\newtheorem{assumption}{Assumption}[section]
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}{Corollary}[theorem]
\newtheorem{lemma}[theorem]{Lemma}
\newenvironment{proof}{\paragraph{Proof:}}{\hfill$\square$}


\hyphenation{op-tical net-works semi-conduc-tor}

\begin{document}
\title{Traversing the Local Polytopes of ReLU Neural Networks: A Unified Approach for \\ Network Verification}

\author{Shaojie~Xu,
            Joel~Vaughan,
            Jie~Chen,
            Aijun~Zhang,
	       Agus~Sudjianto%
\thanks{The authors are with Wells Fargo \& Company. The views expressed in the paper are those of the authors and do not represent the views of Wells Fargo.}%
}
\maketitle

\begin{abstract}
Although neural networks (NNs) with ReLU activation functions have found success in a wide range of applications, their adoption in risk-sensitive settings has been limited by the concerns on robustness and interpretability. Previous works to examine robustness and to improve interpretability partially exploited the piecewise linear function form of ReLU NNs. In this paper, we explore the unique topological structure that ReLU NNs create in the input space, identifying the adjacency among the partitioned local polytopes and developing a traversing algorithm based on this adjacency. Our polytope traversing algorithm can be adapted to verify a wide range of network properties related to robustness and interpretability, providing an unified approach to examine the network behavior. As the traversing algorithm explicitly visits all local polytopes, it returns a clear and full picture of the network behavior within the traversed region. The time and space complexity of the traversing algorithm is determined by the number of a ReLU NN's partitioning hyperplanes passing through the traversing region.
\end{abstract}

\begin{IEEEkeywords}
ReLU NNs, Piecewise-Linear NNs, Adversarial Attack, Robustness, Interpretability, Network Verification
\end{IEEEkeywords}

\IEEEpeerreviewmaketitle

\section{Introduction \& Related Work} \label{sec:intro}
Neural networks with rectified linear unit activation functions (ReLU NNs) are arguably the most popular type of neural networks in deep learning. This type of network enjoys many appealing properties including better performance than NNs with sigmoid activation \cite{glorot2011deep}, universal approximation ability \cite{arora2018understanding, lu2017expressive, montufar2014number, schmidt2020nonparametric}, and fast training speed via scalable algorithms such as stochastic gradient descent (SGD) and its variants \cite{zou2020gradient}. 

Despite their strong predictive power, ReLU NNs have seen limited adoption in risk-sensitive settings \cite{bunel2018unified}. These settings require the model to make robust predictions against potential adversarial noise in the input \cite{athalye2018synthesizing, carlini2017towards, goodfellow2014explaining, szegedy2014intriguing}. The alignment between model behavior and human intuition is also desirable \cite{liu2019algorithms}: prior knowledge such as monotonicity may be incorporated into model design and training \cite{daniels2010monotone, gupta2019incorporate, liu2020certified, sharma2020testing}; users and auditors of the model may require a certain degree of explanations of the model predictions \cite{gopinath2019property, chu2018exact}.

The requirements in risk-sensitive settings has motivated a great amount of research on verifying certain properties of ReLU NNs. These works often exploit the piecewise linear function form of ReLU NNs. In \cite{bastani2016measuring} the robustness of a network is verified in very small input region via linear programming (LP). To consider the non-linearity of ReLU activation functions, \cite{ehlers2017formal, katz2017reluplex, pulina2010abstraction, pulina2012challenging} formulated the robustness verification problem as a satisfiability modulo theories (SMT) problem. A more popular way to model ReLU nonlinearality is to introduce a binary variable representing the on-off patterns of ReLU neurons. Property verification can then be solved using mixed-integer programming (MIP) \cite{anderson2020strong, fischetti2017deep, liu2020certified, tjeng2018evaluating, weng2018towards}.

The piecewise linear functional form of ReLU NNs also creates distinct topological structures in the input space. Previous studies have shown that a ReLU NN partitions the input space into convex polytopes and has one linear model associated with each polytope \cite{montufar2014number, serra2018bounding, croce2019provable, robinson2019dissecting, sudjianto2020unwrapping, yang2020reachability}. Each polytope can be coded by a binary activation code, which reflects the on-off patterns of the ReLU neurons. The number of local polytopes is often used as a measure of the model's expressivity \cite{hanin2019deep, lu2017expressive}. Built upon this framework, multiple studies \cite{sudjianto2020unwrapping, yang2020enhancing, zhao2021self} tried to explain the behavior of ReLU NNs and to improve their interpretability. They viewed ReLU NN as a collection of linear models. However, the relationship among the local polytopes and their linear models was not fully investigated.

When the network's behavior within some specific region in the input space is of interest, one can collect all the local polytopes overlapped with the region to conduct analysis. The methods to collect these polytopes can be categorized into top-down and bottom-up approaches. The top-down approaches in \cite{xiang2017reachable, yang2020reachability} pass the entire region of interest through a ReLU NN and calculate how the hyperplanes corresponding to the neurons partition the region into local polytopes. The major drawback of the top-down approach is that the analysis must start after the computationally expensive forward passing is fully finished. 

One the contrary, the bottom-up approaches start from a point of interest inside the region, moving from one local polytope to another while running the analysis, and can be stopped at any time. \cite{croce2018randomized, croce2020scaling} achieved the movement among polytopes by generating a sequence of samples in the input space using randomized local search. Although being computationally simple, this sample-based method does not guarantee covering all polytopes inside the region of interest. The most recent work and also the closest to ours is \cite{vincent2021reachable}, where polytope boundaries and adjacency are identified using LP, and the traversing is done directly on the polytopes.

In this paper, we explore the topological relationship among the local polytopes created by ReLU NNs. We propose algorithms to identify the adjacency among these polytopes, based on which we develop traversing algorithms to visit all polytopes within a bounded region in the input space. Compared with \cite{vincent2021reachable}, our polytope traversing algorithm exploits ReLU NNs' hierarchical partitioning of the input space to reduce computational overhead and accelerates the discovering of adjacent polytopes. The thoroughness of our traversing algorithm is proved. Our paper has the following major contributions:
\begin{enumerate}
\item The polytope traversing algorithm provides a unified framework to examine the network behavior. Since each polytope contains a linear model whose properties are easy to verify, the full verification on a bounded domain is achieved after all the covered polytopes are visited and verified. We provide theoretical guarantees on the thoroughness of the traversing algorithm.
\item Property verification based on the polytope traversing algorithm can be easily customized. Identifying the adjacency among the polytopes is formulated as LP. Within each local polytope, the user has the freedom to choose the solver most suitable for the verification sub-problem. We demonstrate that many common applications can be formulated as convex problems within each polytope.
\item Because the polytope traversing algorithm explicitly visits all the local polytopes, it returns a full picture of the network behavior within the traversed region and improves interpretability.
\end{enumerate}
Although we focus on ReLU NN with fully connected layers through out this paper, our polytope traversing algorithm can be naturally extended to other piecewise linear networks such as those containing convolutional and maxpooling layers.

The rest of this paper is organized as follows: Section \ref{sec:llpolytopes} reviews how polytopes are created by ReLU NNs. Section \ref{sec:boundary} introduces two related concepts: the boundaries of a polytope and the adjacency among the polytopes. Our polytope traversing algorithm is described in Section \ref{sec:polytope_traversing}. Section \ref{sec:apps} demonstrates several applications of adapting the traversing algorithm for network property verification. Two specific cases studies are shown in Section \ref{sec:casestudies}. The paper is concluded in Section \ref{sec:conclusion}.

\section{The Local Polytopes in ReLU NNs} \label{sec:llpolytopes}
\subsection{The case of one hidden layer} \label{sec:llpolytopesI}
A ReLU NN partitions the input space $\R^P$ into several polytopes and forms a linear model within each polytope. To see this, we first consider a simple NN with one hidden layer of $M$ neurons. It takes an input $\vx \in \R^P$ and outputs $\vo \in \R^Q$ by calculating:
\small{
\begin{equation}
\begin{split}
\vo = \mW^o\vh + \vb^o &= \mW^o\left(\sigma(\mW\vx + \vb)\right) + \vb^o \\
\text{where}\ 
\sigma(\vx)_m &= 
\begin{cases} 
0,\ & \vx_m < 0 \\
\vx_m,\ & \vx_m \geq 0
\end{cases}
\ .
\end{split} \label{eq:relu_nn_I}
\end{equation}
}%
For problems with a binary or categorical target variable (i.e. binary or multi-class classification), a sigmoid or softmax layer is added after $\vo$ respectively to convert the convert the NN outputs to proper probabilistic predictions.

The ReLU activation function $\sigma({\cdot})$ inserts non-linearity into the model by checking a set of linear inequalities: $\vw_m^T\vx + b_m \geq 0, \ m = 1 , 2, \ldots, M$, where $\vw_m^T$ is the $m$th row of matrix $\mW$ and $b_m$ is the $m$th element of $\vb$. Each neuron in the hidden layer creates a \textbf{partitioning hyperplane} in the input space with the linear equation $\vw_m^T\vx + b_m = 0$. The areas on two sides of the hyperplane are two \textbf{halfspaces}. The entire input space is, therefore, partitioned by these $M$ hyperplanes. We define a \textbf{local polytope} as a set containing all points that fall on the same side of each and every hyperplane. The polytope encoding function (\ref{eq:polytope_encode}) uses an element-wise indicator function $\mathbbm{1}(\cdot)$ to create a unique binary code $\vc$ for each polytope. Since the $m$th neuron is called ``ON'' for some $\vx$ if $\vw_m^T\vx + b_m \geq 0$, the code $\vc$ also represents the on-off pattern of the neurons. Using the results of this encoding function, we can express each polytope as an intersection of $M$ halfspaces as in (\ref{eq:polytope}), where the binary code $\vc$ controls the directions of the inequalities.
{\small
\begin{align}
C(\vx) = &\mathbbm{1}(\mW\vx + \vb \geq 0) \ . \label{eq:polytope_encode} \\
\setR_{\vc} = \{ \vx\ |\ (-1)^{c_m} (\vw_m^T\vx &+ b_m \leq 0),\ \forall m=1,\ldots,M \} \ .  \label{eq:polytope}
\end{align}
}%

Figure \ref{fig:grid_nets}.(b) shows an example of ReLU NN trained on a two-dimensional synthetic dataset (plotted in Figure \ref{fig:grid_nets}.(a)). The bounded input space is $[-1, 1]^2$ and the target variable is binary. The network has one hidden layer of 20 neurons. The partitioning hyperplanes associated with these neurons are plotted as the blue dashed lines. They form in total 91 local polytopes within the bounded input space.

For a given $\vx$, if $\vw_m^T\vx + b_m \geq 0$, the ReLU neuron turns on and passes through the value. Otherwise, the neuron is off and suppresses the value to zero. Therefore, if we know the $m$th neuron is off, we can mask the corresponding $\vw_m$ and $b_m$ by zeros and create $\tilde{\mW}_{\vc}$ and $\tilde{\vb}_{\vc}$ that satisfy (\ref{eq:zero_masking_locally_linear}). The non-linear operation, therefore, can be replaced by the a locally linear operation after zero-masking. Because each local polytope $\setR_{\vc}$ has a unique neuron activation pattern encoded by $\vc$, the zero-masking process in (\ref{eq:zero_masking}) is also unique for each polytope. Here, $\mathbf{1}$ is a vector of 1s of length $p$ and $\otimes$ denotes element-wise product.
{\small
\begin{align}
\tilde{\mW}_{\vc} = \mW \otimes (\vc\mathbf{1}^T) \ ,\ \tilde{\vb}_{\vc} = \vb \otimes \vc \ , \label{eq:zero_masking} \\
\sigma(\mW\vx + \vb) = \tilde{\mW}_{\vc} \vx + \tilde{\vb}_{\vc},\quad \forall \vx \in \setR_{\vc} \ . \label{eq:zero_masking_locally_linear}
\end{align}
}%

Within each polytope, as the non-linearity is taken out by the zero-masking process, the input $\vx$ and output $\vo$ have a linear relationship:
{\small
\begin{equation}
\begin{split}
\vo = \mW^o(\sigma(\mW\vx + \vb)) + \vb^o &= \hat{\mW}_{\vc}^o\vx + \hat{\vb}_{\vc}^o \ ,\ \forall \vx \in \setR_{\vc} \ , \\
\text{where}\ \hat{\mW}_{\vc}^o =\mW^o\tilde{\mW}_{\vc} \ &,\ \hat{\vb}_{\vc}^o = \mW^o\tilde{\vb}_{\vc} + \vb^o
\end{split}
\end{equation}
}%
The linear model associated with polytope $\setR_{\vc}$ has the weight matrix $\hat{\mW}_{\vc}^o$ and the bias vector $\hat{\vb}_{\vc}^o$. The ReLU NN is now represented by a collection of linear models, each defined on a local polytope $\setR_{\vc}$. 

In Figure \ref{fig:grid_nets}.(b), we represent the linear model in each local polytopes as a red solid line indicating $\left(\hat{\vw}^o_{\vc}\right)^T\vx + \hat{b}^o_{\vc} = 0$. In this binary response case, the two sides of this line have the opposite class prediction. We only plot the line if it passes through its corresponding polytope. For other polytopes, the entire polytopes fall on one side of their corresponding class-separating lines and the predicted class is the same within the whole polytope. The red lines all together form the decision boundary of the ReLU NN and are continuous when passing through one polytope to another. This is a direct result of ReLU NN being a continuous model.

\begin{figure*}[t]
\center
\includegraphics[width=1.75\columnwidth]{fig_grid_nets}
\caption{\small Examples of trained ReLU NNs and their local polytopes. (a) The grid-like training data with binary target variable. (b) A trained ReLU NN with one hidden layer of 20 neurons. The heatmap shows the predicted probability of a sample belong to class 1. The blue dashed lines are the partitioning hyperplanes associated with the ReLU neurons, which form 91 local polytopes in total. The red solid lines represent the linear model within each polytope where class separation occurs. (c) A trained ReLU NN with two hidden layers of 10 and 5 neurons respectively. The blue dashed lines are the partitioning hyperplanes associated with the first 10 ReLU neurons, forming 20 level-1 polytopes. The orange dashes lines are the partitioning hyperplanes associated with the second 5 ReLU neurons within each level-1 polytope. There are in total 41 (level-2) local polytopes. The red solid lines represent the linear model within each level-2 polytope where class separation occurs.}
\label{fig:grid_nets}
\end{figure*} 

\subsection{The case of multiple layers} \label{sec:hierarchical_polytopes}
We can generalize the results to ReLU NNs with multiple hidden layers. A ReLU NN with $L$ hidden layers hierarchically partitions the input space and is locally linear in each and every \textbf{level-$L$ polytope}. Each level-$L$ polytope $\setR^L$ has a unique binary code $\vc^1\vc^2\ldots\vc^L$ representing the activation pattern of the neurons in all $L$ hidden layers. The corresponding partitioning hyperplanes of each level, $\hat{\mW}^{l} \vx + \hat{\vb}^{l} = 0$, $l=1,2,\ldots,L$, can be calculated recursively level by level, using the zero masking procedure:
{\small
\begin{align}
&\hat{\mW}^1 = \mW^1 \ , \ \hat{\vb}^1 = \vb^1 \label{eq:cal_ieq_begin} \\ 
&\tilde{\mW}^{l} = \hat{\mW}^{l} \otimes (\vc^{l}\mathbf{1}^T) \ ,\ \tilde{\vb}^{l} = \hat{\vb}^{l} \otimes \vc^{l} \label{eq:zero_masking_level_l} \\
&\hat{\mW}^{l+1} = \mW^{l+1}\tilde{\mW}^{l}\ , \  \hat{\vb}^{l+1} = \mW^{l+1}\tilde{\vb}^{l} + \vb^{l+1} \label{eq:coeffs_level_l} \ .
\end{align}
}%
We emphasis that $\tilde{\mW}^l$, $\tilde{\vb}^l$, $\hat{\mW}^{l+1}$, and $\hat{\vb}^{l+1}$ depend on all polytope code up to level $l$: $\vc^1\vc^2\ldots\vc^l$. The subscription $\vc$ is dropped to simplify the notations.

At each level $l$, the encoding function $C^l(\cdot)$ and the polytope $\setR^l$ expressed as an intersection of $\sum_{t=1}^l M_t$ halfspaces can be written recursively as:
{\small
\begin{align}
&C^1(\vx) = \mathbbm{1}(\mW^1\vx + \vb^1 \geq 0)  \\
\begin{split}
&\setR^1 = \{ \vx\ |\ (-1)^{c_{m}} \left((\vw^1)_{m}^T\vx + (b^1)_{m} \leq 0 \right),\\ 
&\quad\quad\quad\quad\forall m=1,2,\ldots,M_1 \}
\end{split}\\
&C^{l+1}(\vx) = \mathbbm{1}(\hat{\mW}^{l+1}\vx + \hat{\vb}^{l+1} \geq 0) \ ,\ \forall \vx \in \setR^{l} \label{eq:polytope_encoding_l} \\
\begin{split}
&\setR^{l+1} = \{ \vx\ |\ (-1)^{c_{m}} \left( (\hat{\vw}^{l+1})_{m}^T\vx + (\hat{b}^{l+1})_{m} \leq 0 \right),\\
&\quad\quad\quad\quad\forall m=1,2,\ldots,M_{l+1} \}\ \cap\ \setR^{l} \ . 
\end{split} \label{eq:polytope_level_l}
\end{align}
}%
Finally, the linear model in a level-$L$ polytope is:
{\small
\begin{equation}
\begin{split}
\vo = \hat{\mW}^o\vx + \hat{\vb}^o \ &,\ \forall \vx \in \setR^L \ , \\ 
\text{where}\ \hat{\mW}^o =\mW^o\tilde{\mW}^L \ &,\ \hat{\vb}^o = \mW^o\tilde{\vb}^L + \vb^o \ . \label{eq:local_model}
\end{split}
\end{equation}
}%

Figure \ref{fig:grid_nets}.(c) shows an example of ReLU NN with two hidden layers of size 10 and 5 respectively. The partitioning hyperplanes associated with the first 10 neuron are plotted as the blue dashed lines. They form 20 level-1 polytopes within the bounded input space. Within each of the level-1 polytope, the hyperplanes associated with the second 5 neurons further partition the polytope. In many cases, some of the 5 hyperplanes are outside the level-1 polytope and, therefore, not creating a new sub-partition. The hyperplanes do create new partitions are plotted as the orange dashed lines. The orange lines are only straight within a level-1 polytope but are continuous when passing through one polytope to another, which is also a result of ReLU NN being a continuous model. In total, this ReLU NN creates 41 (level-2) local polytopes. As in Figure \ref{fig:grid_nets}.(b), the linear model within each level-2 polytope is represented as a red solid line if class separation occurs within the polytope.

\section{Polytope Boundaries and Adjacency} \label{sec:boundary}
Beyond viewing ReLU NNs as a collection of linear models defined on local polytopes, we explore the topological relationship among these polytopes. A key concept is the \textbf{boundaries} of each polytope. As shown in (\ref{eq:polytope_level_l}), each level-$l$ polytope $\setR_{\vc}$ with corresponding binary code $\vc=\vc^1\vc^2\ldots\vc^l$ is an intersection of $\sum_{t=1}^l M_t$ halfspaces induced by a set of inequality constraints. Two situations can rise among these inequalities. First, an arbitrary $\vc$ may lead to conflicting inequalities and makes $\setR_{\vc}$ an empty set. This situation can be common when the number of neurons is much larger than the dimension of the input space. Second, there can be \textbf{redundant inequalities} which means removing them does not affect set $\setR_{\vc}$. We now show that the non-redundant inequalities are closely related to the boundaries of a polytope.

\begin{definition}
Let $\setR$ contains all $\vx\in\R^P$ that satisfy $M$ linear inequalities:  $\setR = \{ \vx | g_1(\vx) \leq 0, g_2(\vx) \leq 0,\ldots, g_M(\vx) \leq 0 \}$. Assume that $\setR \neq \emptyset$. Let $\tilde{\setR}$ contains all $\vx$'s that satisfy $M-1$ linear inequalities:  $\tilde{\setR} = \{ \vx | g_1(\vx) \leq 0, \ldots, g_{m-1}(\vx) \leq 0 ,g_{m+1}(\vx) \leq 0, \ldots, g_M(\vx) \leq 0 \}$. Then the inequality $g_m(\vx) \leq 0$ is a \textbf{redundant inequality} with respect to (w.r.t.) $\setR$ if $\setR = \tilde{\setR}$. \label{def:redundant_ieq}
\end{definition}

With the redundant inequality defined above, the following lemma provides an algorithm to identify them. The proof of this lemma is in the Appendix.
\begin{lemma}
Given a set $\setR = \{ \vx | g_1(\vx) \leq 0,\ldots, g_M(\vx) \leq 0 \} \neq \emptyset$, then $g_m(\vx)$ is a redundant inequality if the new set formed by flipping this inequality is empty: $\hat{\setR} = \{ \vx | g_1(\vx) \leq 0, \ldots, g_{m}(\vx) \geq 0, \ldots, g_M(\vx) \leq 0 \} = \emptyset$. \label{them:redundant_ieq}
\end{lemma}

We can now define the boundaries of a polytope formed by a set of linear inequalities using a similar procedure in Lemma\ref{them:redundant_ieq}. The concept of polytope boundaries also leads to the definition of adjacency. Intuitively, we can move from one polytope to its adjacent polytope by crossing a boundary. 
\begin{definition}
Given a non-empty set formed by $M$ linear inequalities: $\setR = \{ \vx | g_1(\vx)\leq0,\ldots, g_M(\vx)\leq0 \} \neq \emptyset$, then the hyperplane $g_m(\vx) = 0$ is a \textbf{boundary} of $\setR$ if the new set formed by flipping the corresponding inequality is non-empty: $\hat{\setR} = \{ \vx | g_1(\vx) \leq 0, \ldots, g_{m}(\vx) \geq 0, \ldots, g_M(\vx) \leq 0 \} \neq \emptyset$. Polytope $\hat{\setR}$ is called \textbf{one-adjacent} to $\setR$. \label{def:boundary_adj}
\end{definition}

Since for each polytope the directions of its linear inequalities are reflected by the binary code, two one-adjacent polytopes must have their code differ by one bit. Figure \ref{fig:polytope_traversing}.(a) demonstrates the adjacency among the local polytopes. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(b). Using the procedure in Definition \ref{def:boundary_adj}, 4 out of the 20 partitioning hyperplanes are identified as the boundaries of polytope No.0 and marked in red. The 4 one-adjacent neighbors to polytope No.0 are No.1, 2, 3, and 4; each can be reached by crossing one boundary. 

As we have shown in the Section \ref{sec:hierarchical_polytopes}, ReLU NNs create polytopes level by level. We follow the same hierarchy to define the polytope adjacency. Assume two non-empty level-$l$ polytopes, $\setR$ and $\hat{\setR}$, are inside the same level-$(l-1)$ polytope, which means their corresponding code $\vc=\vc^1\vc^2\ldots\vc^l$ and $\hat{\vc}=\vc^1\vc^2\ldots\hat{\vc}^l$ only differs at level-$l$. We say that polytope $\hat{\setR}$ is a \textbf{level-$l$ one-adjacent neighbor} of $\setR$ if $\hat{\vc}^l$ and $\vc^l$ only differs in one bit.

The condition that $\vc=\vc^1\vc^2\ldots\vc^l$ and $\hat{\vc}=\vc^1\vc^2\ldots\hat{\vc}^l$ only differs at level-$l$ is important. In this way, the two linear inequalities associated with each pair of bits in $\vc$ and $\hat{\vc}$ have the same coefficients, and the difference in $\vc^l$ and $\hat{\vc}^l$ only changes the direction of the linear inequality. On the other hand, if the two codes differ at a level $l' < l$, then according to the recursive calculation in (\ref{eq:zero_masking_level_l}) and (\ref{eq:coeffs_level_l}), the codes starting from level $l'+1$ will correspond to linear inequalities of different coefficients, leaving our Definition \ref{def:boundary_adj} of adjacency not applicable.

Figure \ref{fig:polytope_traversing}.(b) demonstrates the hierarchical adjacency among the local polytopes. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(c). Level-1 polytopes $(1,\cdot)$ and $(2,\cdot)$ are both (level-1) one-adjacent to $(0,\cdot)$. Within the level-1 polytope $(0,\cdot)$, Level-2 polytopes $(0,0)$ and $(0,1)$ are (level-2) one-adjacent to each other. Similarly, we can identify the level-2 adjacency of the other two pairs $(1,0)-(1,1)$ and $(2,0)-(2,1)$. Note that in the plot, even thought one can move from polytope $(2,1)$ to $(0,1)$ by crossing one partitioning hyperplane, we do not define these two polytopes as adjacent, as they lie into two different level-1 polytopes.

\section{Polytope Traversing} \label{sec:polytope_traversing}

\begin{figure*}[t]
\center
\includegraphics[width=1.68\columnwidth]{fig_polytope_traversing}
\caption{\small Demonstration of the BFS-base polytope traversing algorithm. (a) Traversing the 8 local polytopes within the bounded regions. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(b). The lines marked in red are the boundaries of polytope No.0. (b) Traversing the 6 local polytopes within the bounded region. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(c). The polytopes are indexed as ``(level-1, level-2)''. (c) The evolution of the BFS queue for traversing the local polytopes in (a). The gray arrows show the traversing order. The colored arrows at the bottom indicate the one-adjacent neighbors added to the queue. (d) The evolution of the hierarchical BFS queue for traversing the local polytopes in (b). The level-1 BFS queue is shown vertically while the level-2 BFS queue is shown horizontally.}
\label{fig:polytope_traversing}
\end{figure*} 

\subsection{The case of one hidden layer} \label{sec:polytope_traversing_I}
The adjacency defined in the previous section provides us an order to traverse the local polytopes: starting from an initial polytope $\setR$, visiting its all one-adjacent neighbors, then visiting all the neighbors' neighbors and so on.

This algorithm can be viewed as breath-first search (BFS) on a \textbf{polytope graph}. To create this graph, we turn each polytope created by the ReLU NN into a node. An edge is added between each pair of polytopes that are one-adjacent to each other. The BFS algorithm uses a queue to keep track the traversing progress. At the beginning of traversing, the initial polytope is added to an empty queue and is marked as visited afterwards. In each iteration, we pop out the first polytope from the queue and identify all of its one-adjacent neighbors. Among these identified polytopes, we add those that have not been visited to the back of the queue and mark them as visited. The iteration stops when the queue is empty.

The key component of the polytope traversing algorithm is to identify a polytope's one-adjacent neighbors. For a polytope $\setR_{\vc}$ coded by $\vc$ of $M$ bits, there are at most $M$ one-adjacent neighbors with codes corresponding to flipping one of the bits in $\vc$. Each valid one-adjacent neighbor must be non-empty and can be reached by crossing a boundary. Therefore, we can check each linear inequality in (\ref{eq:polytope}) and determine whether it is a boundary or redundant. Some techniques of identifying redundant inequalities are summarized in \cite{telgen1983identifying}. By flipping the bits corresponding to the identified boundaries, we obtain the codes of the one-adjacent polytopes.

Equivalently, we can identify the one-adjacent neighbors by going through all $M$ candidate codes and selecting those corresponding to non-empty sets. Checking the feasibility of a set constrained by a set of linear inequalities is often referred to as the ``Phase-I Problem'' of LP and can be solved efficiently by modern LP solvers. During BFS iterations, we can hash the checked codes to avoid checking them repetitively. The BFS-based polytope traversing algorithm is summarized in Algorithm \ref{algo:traverseI}. We now state the correctness of this algorithm with its proof in Appendix.
\begin{theorem}
Given a ReLU NN with one hidden layer of $M$ neurons as specified in (\ref{eq:relu_nn_I}), Algorithm \ref{algo:traverseI} covers all non-empty local polytopes created by the neural network. That is, for all $\vx \in \R^P$, there exists one $\setR_{\vc}$ as defined in (\ref{eq:polytope}) such that $\vx \in \setR_{\vc}$ and $\vc \in \setS_R$, where $\setS_R$ is the result returned by Algorithm \ref{algo:traverseI}.
\label{them:traverseI}
\end{theorem}

Algorithm \ref{algo:traverseI} visits all the local polytopes created by a ReLU NN within $\R^P$. The time complexity is exponential to the number of neurons, as all $2^M$ possible activation patterns are checked once in the worst-case scenario. The space complexity is also exponential to the number of neurons as we hash all the checked activation patterns. Furthermore, for each activation pattern, we solve a phase-I problem of LP with $M$ inequalities in $\R^P$. Traversing all local polytopes in $\R^P$, therefore, becomes intractable for neural networks with a large number of neurons. 

Fortunately, traversing in $\R^P$ is usually undesirable. Firstly, a neural network may run into extrapolation issues for points outside the sample distribution. The polytopes far away from the areas covered by the samples are often considered unreliable. Secondly, many real-life applications, to be discussed in Section \ref{sec:apps}, only require traversing within small bounded regions to examine the local behavior of a model. In the next section, we introduce a technique to improve the efficiency when traversing within a bounded region. 

\begin{algorithm}[thb]
\small
\caption{BFS-Based Polytope Traversing} \label{algo:traverseI}
\begin{algorithmic}[1]
\Require A ReLU NN with one hidden layer of $M$ neurons as specified in (\ref{eq:relu_nn_I}).
\Require An initial point $\vx\in\R^P$.
\State Initialize an empty queue $\setQ$ for BFS.
\State Initialize an empty set $\setS_R$ to store the codes of all visited polytopes.
\State Initialize an empty set $\setS_{\vc}$ to store all checked codes. 
\State Calculate $\vx$'s initial polytope code $\vc$ using (\ref{eq:polytope_encode}).
\State Append $\vc$ to the end of the $\setQ$.
\State Add $\vc$ to both $\setS_R$ and $\setS_{\vc}$.
\While {$\setQ$ is not empty}
	\State Pop out the first element in the front of BFS queue: $\vc = \setQ.\text{pop}()$.
	\For {$m=1,2,\ldots,M$}
\State Create a candidate polytope code $\hat{\vc}$ by flipping one bit in $\vc$: $\hat{c}_m = 1-c_m$ and $\hat{c}_k = c_k \forall k \neq m$.
\If {$\hat{\vc} \notin \setS_{\vc}$}
\State Check if $\setR_{\hat{\vc}} = \{ \vx|(-1)^{\hat{c}_k}\left(\vw_k^T\vx + b_k\right) \leq 0,\ k=1,2\ldots,M \}$ is empty using LP.
\State Add $\hat{\vc}$ to $\setS_{\vc}$.
\If {$\setR_{\hat{\vc}} \neq \emptyset$}
\State Append $\hat{\vc}$ to the end of the $\setQ$.
\State Add $\hat{\vc}$ to $\setS_R$.
\EndIf
\EndIf
	\EndFor
\EndWhile
\State Return $\setS_R$.
\end{algorithmic}
\end{algorithm}

\subsection{Polytope traversing within a bounded region} \label{sec:bounded_polytope_traversing}
We first consider a region with each dimension bounded independently: $l_j \leq x_j \leq u_j$, $j=1,2,\ldots,P$. These $2\times P$ linear inequalities creates a hypercube denoted as $\setB$. During the BFS-based polytope traversing, we repetitively flip the direction of one of the $M$ inequalities to identify the one-adjacent neighbors. When the bounded region is small, it is likely that only a small number of the $M$ hyperplanes cut through the hypercube. For the other hyperplanes, the entire hypercube falls onto only one side. Flipping to the other sides of these hyperplanes would leave the bounded region. Therefore, at the very beginning of polytope traversing, we can run through the $M$ hyperplanes to identify those cutting through the hypercube. Then in each neighbor identifying step, we only flip these hyperplanes.

To identify the hyperplanes cutting through the hypercube, we denote the two sides of a hyperplane $\setH$ and $\bar{\setH}$:  $\setH=\{\vx | \vw_m^T\vx + b_m \leq 0 \}$ and $\bar{\setH}=\{\vx | \vw_m^T\vx + b_m \geq 0 \}$. If neither $\setH\cap\setB$ nor $\hat{\setH}\cap\setB$ is empty, we say the hyperplane $\vw_m^T\vx + b_m = 0$ cuts through $\setB$. $\setH\cap\setB$ and $\hat{\setH}\cap\setB$ are both constrained by $2\times P + 1$ inequalities, checking their feasibility can again be formulated as a phase-I problem of LP. We name this technique \textbf{hyperplane pre-screening} and summarize it in algorithm \ref{algo:prescreening}.

\begin{algorithm}[thb]
\small
\caption{Hyperplane Pre-Screening} \label{algo:prescreening}
\begin{algorithmic}[1]
\Require A set of hyperplanes $\vw_m^T\vx + b_m \leq 0$, $m=1,2,\ldots,M$.
\Require A bounded traversing region $\setB$, e.g. $\{\vx | l_j \leq x_j \leq u_j$, $j=1,2,\ldots,P\}$.
\State Initialize an empty set $\setT$ to store all hyperplanes cutting through $\setB$. 
\For {$m=1,2,\ldots,M$}
	\State Get two halfspaces $\setH=\{\vx | \vw_m^T\vx + b_m \leq 0 \}$ and $\bar{\setH}=\{\vx | \vw_m^T\vx + b_m \geq 0 \}$.
	\If {$\setH\cap\setB\neq\emptyset$ and $\hat{\setH}\cap\setB\neq\emptyset$}
\State Add $m$ to $\setT$.
	\EndIf
\EndFor
\State Return $\setT$.
\end{algorithmic}
\end{algorithm}

Hyperplane pre-screening effectively reduces the complexity from $\orderof{2^M}$ to $\orderof{2^{|\setT|}}$, where $|\setT|$ is the number of hyperplanes cutting through the hypercube. The number $2^{|\setT|}$ corresponds to the worst-case scenario. Since the BFS-based traversing only checks non-empty polytopes and their potential one-adjacent neighbors, the number of activation patterns actually checked can be less than $2^{|\setT|}$. In general, the fewer hyperplanes go through $\setB$ the faster polytope traversing finishes.

Figure \ref{fig:polytope_traversing}.(a) shows traversing the 8 local polytopes within the bounded region. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(b). The lines marked in red are the hyperplanes cutting through the bounded region and are identified by the pre-screening algorithm. The evolution of the BFS queue is shown in Figure \ref{fig:polytope_traversing}.(c). The gray arrows show the traversing order. The colored arrows at the bottom indicate the one-adjacent neighbors added to the queue. When polytope No.0 is popped from the queue, its one-adjacent neighbors, No.1, 2, 3, and 4, are added to the queue. Next, when polytope No.1 is popped, its one-adjacent neighbors, No.5 and 6, are added. Polytope No.0, although as a one-adjacent neighbor to No.1, is ignored since it has been visited. Similarly, when polytope No.2 is popped, only one of its one-adjacent neighbors, No. 7, is added, since all others have been visited (including those in the queue). The algorithm finished after popping Polytope No.8 as no new polytopes can be added and the queue is empty. All 8 local polytopes in the bounded region are traversed. 

Because $\setB$ is bounded by a set of linear inequalities, the correctness of BFS-based polytope traversing as stated in Theorem \ref{them:traverseI} can be easily extended to this bounded traversing case. Following similar steps of the proof of Theorem \ref{them:traverseI} in the Appendix, we can show that for any two non-empty polytopes overlapped with $\setB$, we can move from one to another by repetitively finding a one-adjacent neighbor within $\setB$. We emphasis that the correctness of BFS-based polytope traversing can be proved for any traversing region bounded by a set of linear inequalities. This realization is critical to generalize our results to the case of ReLU NNs with multiple hidden layers. Furthermore, as any closed convex set can be represented as the intersection of a set of (possibly infinite) halfspaces, the correctness of BFS-based polytope traversing is true for any closed convex $\setB$. 

\subsection{Hierarchical polytope traversing in the case of multiple hidden layers} \label{sec:hierarchical_polytope_traversing}
The BFS-based polytope traversing algorithm can be generalized to ReLU NNs with multiple hidden layers. In section \ref{sec:hierarchical_polytopes}, we described how a ReLU NN with $L$ hidden layers hierarchically partition the input space into polytopes of $L$ different level. Then in section\ref{sec:boundary}, we showed the adjacency of level-$l$ polytopes is conditioned on all of them belonging to the same level-$(l-1)$ polytope. Therefore, to traverse all level-$L$ polytopes, we need to traverse all level-$(L-1)$ polytopes and within each of them traversing the sub-polytopes by following the one-adjacent neighbors.

The procedure above leads us to a recursive traversing scheme. Assume a ReLU NN with L hidden layers and a closed convex traversing region $\setB$. Starting from a sample $\vx \in \setB$, we traverse all level-1 polytopes using the BFS-based algorithm. Inside each level-1 polytope, we traverse all the contained level-2 polytopes, so on and so forth until we reach the level-L polytopes. As shown in (\ref{eq:polytope_level_l}), each level-$l$ polytope is constrained by $\sum_{t=1}^l M_t$ linear inequalities, the way to identify level-$l$ one-adjacent neighbors is largely the same as what we have described in Section \ref{sec:polytope_traversing_I}. Two level-$l$ one-adjacent neighbors must have the same $\sum_{t=1}^{l-1} M_t$ linear inequalities corresponding to $\vc^1\vc^2\ldots\vc^{l-1}$, and have one of the last $M_l$ inequalities differ in direction, so there are $M_l$ cases to check.

We can use hyperplane pre-screening at each level of traversing. When traversing the level-$l$ polytopes within in a level-$(l-1)$ polytope $\setR^{l-1}$, we update the bounded traversing region by taking the intersection of $\setR^{l-1}$ and $\setB$. We then screen the $M_l$ partitioning hyperplanes and only select those passing through this update traversing region.

The BFS-based hierarchical polytope traversing algorithm is summarized in Algorithm \ref{algo:hierarchical_traverse}. The correctness of this algorithm can be proved based on the results in Section \ref{sec:bounded_polytope_traversing}, which guarantees the thoroughness of traversing the level-$l$ polytopes within in any level-$(l-1)$ polytope. Then the overall thoroughness is guaranteed because each level of traversing is thorough. We state the result in the following theorem.
\begin{theorem}
Given a ReLU NN with $L$ hidden layers and a closed convex traversing region $\setB$. Algorithm \ref{algo:hierarchical_traverse} covers all non-empty level-$L$ polytopes created by the neural network that overlap with $\setB$. That is, for all $\vx \in \setB$, there exists one $\setR_{\vc}$ as defined in (\ref{eq:polytope_level_l}) such that $\vx \in \setR_{\vc}$ and $\vc \in \setS_R$, where $\setS_R$ is the result returned by Algorithm \ref{algo:hierarchical_traverse}.
\label{them:hierarchical_traverse}
\end{theorem}

Figure \ref{fig:polytope_traversing}.(b) shows traversing the 6 local polytopes within the bounded region. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(c).  The evolution of the hierarchical BFS queue is shown in Figure \ref{fig:polytope_traversing}.(d). The level-1 BFS queue is shown vertically while the level-2 BFS queue is shown horizontally. Starting from level-1 polytope $(0,\cdot)$, the algorithm traverses the two level-2 polytopes inside it (line 10 in Algorithm \ref{algo:hierarchical_traverse}). It then identifies the two (level-1) one-adjacent neighbors of $(0,\cdot)$:  $(1,\cdot)$ and $(2,\cdot)$. Every time a level-1 polytope is identified, the algorithm goes into it to traverse all the level-2 polytopes inside (line 36). At the end of the recursive call, all 6 local polytopes in the bounded region are traversed. 

\begin{algorithm}[thb]
\small
\caption{BFS-Based Hierarchical Polytopes Traversing in a Bounded Region} \label{algo:hierarchical_traverse}
\begin{algorithmic}[1]
\Require A ReLU NN with $L$ hidden layers.
\Require A closed convex traversing region $\setB$.
\Require An initial point $\vx\in\setB$.
\State Initialize an empty set $\setS_R$ to store the codes of all visited polytopes.
\State
\Function{HIERARCHICAL\_TRAVERSE}{$\vx, l$} 
\State Initialize an empty queue $\setQ^l$ for BFS at level $l$.
\State Initialize an empty set $\setS_{\vc}^l$ to store all checked level-$l$ codes. 
\State Calculate $\vx$'s initial polytope code $\vc$ recursively using (\ref{eq:polytope_encoding_l}).
\If {$l == L$}
	\State Add $\vc$ to $\setS_R$
\Else
	\State HIERARCHICAL\_TRAVERSE($\vx$,$l$+1) 
\EndIf
\If {$l>1$}
\State Get the level-$(l-1)$ polytope code specified by the front segment of $\vc$: $\vc^{1:l-1}=\vc^1\vc^2\ldots\vc^{l-1}$.
	\State Use $\vc^{1:l-1}$ to get the level-$(l-1)$ polytope $\setR_{\vc}^{l-1}$ as in (\ref{eq:polytope_level_l}).
\Else
	\State $\setR_{\vc}^0 = \R^P$
\EndIf
\State Form the new traversing region  $\setB^{l-1} = \setB\cap\setR_{\vc}^{l-1}$.
\State Append the code segment $\vc^l$ to the end of the $\setQ^l$.
\State Add the code segment $\vc^l$ to $\setS_{\vc}$.
\State Get the $M_l$ hyperplanes associated with $\vc^l$.
\State Pre-Screen the hyperplanes associated with $\vc^l$ using Algorithm \ref{algo:prescreening} with bounded region $\setB^{l-1}$.
\State Collect the pre-screening results $\setT$.
\While {$\setQ^l$ is not empty}
	\State Pop the first element in the front of BFS queue: $\vc^l = \setQ^l.\text{pop}()$.
	\For {$m\in\setT$}
\State Create a candidate polytope code $\hat{\vc}^l$ by flipping one bit in $\vc^l$: $\hat{c}_m^l = 1-c_m^l$ and $\hat{c}_k^l = c_k^l \forall k \neq m$.
\If {$\hat{\vc}^l \notin \setS_{\vc}$}
\State Get set $\setR_{\hat{\vc}} =  \{ \vx|(-1)^{\hat{c}_k}\left(\langle\hat{\vw}_k^l,\vx\rangle + \hat{b}_k^l \right) \leq 0,\ k=1,2\ldots,M_l \}$ 
\State Check if $\setR_{\hat{\vc}} \cap \setB^{l-1}$ is empty using LP.
\State Add $\hat{\vc}^l$ to $\setS_{\vc}$.
\If {$\setR_{\hat{\vc}} \cap \setB^{l-1} \neq \emptyset$}
\State Append $\hat{\vc}^l$ to the end of the $\setQ^l$.
\If {$l == L$}
\State Add $\hat{\vc}=\vc^1\vc^2\ldots\hat{\vc}^l$ to $\setS_R$
\Else
\State Find a point $\hat{\vx} \in \setR_{\hat{\vc}} \cap \setB^{l-1}$
\State HIERARCHICAL\_TRAVERSE($\hat{\vx}$,$l$+1) 
\EndIf
\EndIf
\EndIf
	\EndFor
\EndWhile
\EndFunction
\State
\State HIERARCHICAL\_TRAVERSE($\vx$,1) 
\State Return $\setS_R$.
\end{algorithmic}
\end{algorithm}

\section{Network Property Verification Based on Polytope Traversing} \label{sec:apps}
The biggest advantage of the polytope traversing algorithm is its ability to be adapted to solve many different problems of practical interest. Problems such as local adversarial attacks, searching for counterfactual samples, and local monotonicity verification can be solved easily when the model is linear. As we have shown in Sections \ref{sec:hierarchical_polytopes}, the local model within each level-$L$ polytope created by a ReLU NN is indeed linear. The polytope traversing algorithm provides a way to analyze not only the behavior of a ReLU NN at one local polytope but also the behavior within the neighborhood, and therefore enhances our understanding of the overall model behavior. In this section, we describe the details of adapting the polytope traversing algorithm to verify several properties of ReLU NNs. 

\begin{figure*}[t]
\center
\includegraphics[width=1.75\columnwidth]{fig_apps}
\caption{\small Demonstration of different applications of the polytope traversing algorithm. We use the ReLU NN in Figure \ref{fig:grid_nets}.(b) as an example. (a) Conducting local adversarial attack by finding the maximum (green) and minimum (red) model predictions within a bounded region. (b) Creating counterfactual samples that are closest to the original sample. The distance are measured in $L_1$ (green) and $L_2$ (red) norms. (c) Monotonicity verification in a bounded region. The polytope in red violates condition of model prediction monotonically decreasing along the horizontal axis.}
\label{fig:apps}
\end{figure*} 

\subsection{Local Adversarial Attacks}
We define the local adversarial attack problem as finding the perturbation within a bounded region such that the model output can be changed most adversarially. Here, we assume the model output to be a scalar in $\R$ and consider three regression cases with different types of response variable: continuous, binary, and categorical. The perturbation region is a convex set around the original sample. For example, we can allow certain features to increase or decrease by certain amount; or we can use a norm ($L_1$, $L_2$, $L_\infty$) ball centered at the original sample. 

In the continuous response case, the one-dimensional output after the last linear layer of a ReLU NN is directly used as the prediction of the target variable. Denote the model function as $f(\cdot)$, the original sample as $\vx_0$, and the perturbation region as $\setB$. The local adversarial attack problem can be written as:
{\small
\begin{equation}
\begin{split}
\max_{\vx\in\setB} |f(\vx) - f(\vx_0)| = \max\Big( \max_{\vx\in\setB} f(\vx) - f(\vx_0), \\
f(\vx_0) - \min_{\vx\in\setB} f(\vx) \Big) \ , \label{eq:local_adversarial_attack}
\end{split}
\end{equation}
}%
which means we need to find the range of the model outputs on $\setB$. We can traverse all local polytopes covered by $\setB$, finding the model output range within each intersection $\setB\cap\setR$, then aggregating all the local results to get the final range. Finding the output range within each $\setB\cap\setR$ is a convex problem with linear objective function, so the optimality can be guaranteed within each polytope. Because our traversing algorithm covers all polytopes overlapped with $\setB$, the final solution also has guaranteed optimality.

In the case of binary response, the one-dimensional output after the last linear layer of a ReLU NN is passed through a logistic/sigmoid function to predict the probability of a sample belonging to class 1. To conduct adversarial attack, we minimize the predicted probability $f(\vx)$ if the true response $y$ is 1, and maximize the prediction if the true response is 0:
{\small
\begin{equation}
\begin{cases}
\max_{\vx\in\setB} f(\vx), \quad y = 0 \\
\min_{\vx\in\setB} f(\vx), \quad y = 1 \ .
\end{cases}
\end{equation}
}%
Because of the monotonicity of the logistic function, the minimizer and maximizer of the probabilistic output are the same minimizer and maximizer of the output after the last linear layer (i.e. the predicted log odds), making it equivalent to the case of continuous response.

In the case of categorical response with levels 1 to $Q$, the output after the last linear layer of a ReLU NN is in $\R^Q$ and is passed through a softmax layer to be converted to probabilistic predictions of a sample belonging to each class. The adversarial sample is generated to minimize the predicted probability of the sample being in its true class. Within each local polytope, the linear models are given by (\ref{eq:local_model}), and the predicted probability of class $q$ can be minimized by finding the maximizer of the following optimization problem:
{\small
\begin{equation}
\max_{\vx\in\setB\cap\setR} \sum_{i=1, i\neq q}^Q e^{(\hat{\vw}_i^o - \hat{\vw}_q^o )^T\vx+ (\hat{b}_i^o - \hat{b}_q^o )} \ , \label{eq:multiclass_adversarial_attack}
\end{equation} 
}%
where $\left(\hat{\vw}_i^o\right)^T$ is the $i$th row of the matrix $\hat{\mW}^o$ and $\hat{b}_i^o$ is the $i$th element in $\hat{\vb}^o$. Since the objective function in (\ref{eq:multiclass_adversarial_attack}) is convex, the optimality of local adversarial attack with polytope traversing is guaranteed.

Figure \ref{fig:apps}.(a) demonstrates local adversarial attack in the case of regression with binary response. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(b), which predicts the probability of a sample belong to class 1. The predictions across the whole domain are shown as the heat map. Within the region bounded by the black box, we find the minimum and maximum predictions and marked them by red and green respectively. Due to the nature of linear models, the minimizer and maximizer always fall on the intersections of partitioning hyperplanes and/or region boundaries.

\subsection{Counterfactual sample generation}
In classification problems, we are often interested in finding the smallest perturbation on a sample such that the model changes its class prediction. The magnitude of the perturbation is often measured by $L_1$, $L_2$, or $L_\infty$ norm. The optimization problem can be written as:
{\small
\begin{equation}
\min_{\vx} ||\vx-\vx_0||_p \quad \text{s.t.} f_{\setC}(\vx) \neq f_{\setC}(\vx_0) \ , \label{boundary_proj}
\end{equation}
}%
where $\vx_0$ is the original sample, $p$ indicates a specific type of norm, and $f_{\setC}(\cdot)$ is a ReLU NN outputting class predictions.

We can adapt the polytope traversing algorithm to solve this problem. In the case of binary response, each local polytope has an associated hyperplane separating the two classes: $(\hat{\vw}^o)^T\vx + \hat{b}^o=\gamma$, where $\hat{\vw}^o$ and $\hat{b}^o$ are given in (\ref{eq:local_model}), and $\gamma$ is the threshold converting predicted log odds to class. Finding the counterfactual sample within a local polytope $\setR$ can be written as a convex optimization problem:
{\small
\begin{equation}
\min_{\vx} ||\vx-\vx_0||_p \quad \text{s.t.} (-1)^{\hat{y}_0} \left((\hat{\vw}^o)^T\vx + \hat{b}^o\right) > \gamma,\ \vx\in\setR \ . \label{binary_boundary_proj}
\end{equation}
}%
where $\hat{y}_0$ is the original class (0 or 1) predicted by the model.

We start the traversing algorithm from the polytope where $\vx_0$ lies. In each polytope, we solve  (\ref{binary_boundary_proj}). It is possible that the entire polytope fall on one side of the class separating hyperplane and (\ref{binary_boundary_proj}) does not have any feasible solution. If a solution can be obtained, we compare it with the solutions in previously traversed polytopes and keep the one with the smallest perturbation. Furthermore, we use this perturbation magnitude to construct a new bounded traversing region around $\vx_0$. Because no points outside this region can have a smaller distance to the original points, once we finish traversing all the polytopes inside this region, the algorithm can conclude. In practice we often construct this dynamic traversing region as $\setB = \{ \vx\ |\ ||\vx-\vx_0||_{\infty} < d^* \}$, where $d^*$ is the smallest perturbation magnitude so far. When solving for  (\ref{binary_boundary_proj}) in the proceeding polytopes, we add $x\in\setB$ to the constraints. $\setB$ is updated whenever a smaller $d^*$ is found. Because the new traversing region is always a subset of the previous one, our BFS-based traversing algorithm covers all polytopes within the final traversing region under this dynamic setting. The final solution to (\ref{boundary_proj}) is guaranteed to be optimal, and the running time depends on how far the original point is away from a class boundary.

In the case of categorical response with levels 1 to $Q$, the output after the last linear layer of a ReLU NN has $Q$ dimensions and the dimension of the largest value is the predicted class. We ignore the softmax layer at the end because it does not change the rank of the dimensions. Assuming the original example is predicted to belong to class $\hat{q}_0$, we generate counterfactual samples in the rest of $Q-1$ classes. 

We consider one of these classes at a time and denote it as $q$. Within each ReLU NN's local polytope,  the linear models are given by (\ref{eq:local_model}). The area where a sample is predicted to be in class $q$ is enclosed by the intersection of $Q-1$ halfspaces:
{\small
\begin{equation}
\setC_q = \{ \vx|\left(\hat{\vw}_q^o - \hat{\vw}_i^o\right)^T\vx + (\hat{b}_q^o - \hat{b}_i^o ) > 0, \forall i=1,\ldots,Q, i\neq q \}.
\end{equation}
}%
Therefore, within each local polytope, we solve the convex optimization problem:
{\small
\begin{equation}
\min_{\vx} ||\vx-\vx_0||_p \quad \text{s.t.}\  \vx\in\setC_q \cap \setR \ . \label{multi_boundary_proj}
\end{equation}
}%
We compare all feasible solutions of (\ref{multi_boundary_proj}) under different $q$ and keep the one counterfactual sample that is closest to $\vx_0$. The traversing procedure and the dynamic traversing region update is the same as in the binary response case. Since (\ref{multi_boundary_proj}) is convex, the final solution to (\ref{boundary_proj}) is guaranteed to be optimal.

Figure \ref{fig:apps}.(b) demonstrates counterfactual sample generation in the case of binary classification. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(b) whose class decision boundaries are plotted in red. Given an original sample plotted as the black dot, we generate two counterfactual samples on the decision boundaries. The red dot has the smallest $L_2$ distance to the original point while the green dot has the smallest $L_1$ distance. 

\begin{figure*}[t]
\center
\includegraphics[width=2\columnwidth]{fig_acasxu}
\caption{\small Network verification results of all 45 ACAS Xu networks for property II (a), III (b), and IV (c). The blue lines and markers show the number of local polytopes traversed during verification. The red lines and markers show the time (in seconds) used. A dot marker indicates the corresponding network satisfies the property while a cross marker indicates the property is violated in at least one of the local polytopes.}
\label{fig:acasxu}
\end{figure*}

\subsection{Local monotonicity verification}
We can adapt the polytope traversing algorithm to verify if a trained ReLU NN is monotonic w.r.t. certain features. We consider the regression cases with continuous and binary response. In both cases, the output after the last linear layer is a scalar. Since the binary response case uses a logistic function at the end which is monotonically increasing itself, we can ignore this additional function. The verification methods for the two cases, therefore, are equivalent.

To check whether the model is monotonic w.r.t. a specific feature within a bounded convex domain, we traverse the local polytopes covered by the domain. Since the model is linear within each polytope, we can easily check the monotonicity direction (increasing or decreasing) by checking the sign of the corresponding coefficients. After traversing all local polytopes covered by the domain, we check their agreement on the monotonicity direction. Since a ReLU NN produces a continuous function, if the local models are all monotonically increasing or all monotonically decreasing, the network is monotonic on the checked domain. If there is a disagreement in the direction, the network is not monotonic. The verification algorithm based on polytope traversing not only provides us the final monotonicity result but also tells us in which part of the domain monotonicity is violated.

Figure \ref{fig:apps}.(c) demonstrates local adversarial attack in the case of regression with binary response. The ReLU NN is the same as in Figure \ref{fig:grid_nets}.(b), which predicts the probability of a sample belong to class 1. The predictions across the whole domain are shown as the heat map. We check if the model is monotonically increasing w.r.t. $x_1$ along the horizontal axis. The domain to check is bounded by the black box. Among the 5 polytopes overlapped with the domain, one of them violates the monotonically increasing condition and is marked in red.

\subsection{Comparison with algorithms based on mixed-integer programming}
The three applications above have been traditionally solved using MIP \cite{anderson2020strong, fischetti2017deep, liu2020certified, tjeng2018evaluating, weng2018towards}. Our algorithms based on polytope traversing have several advantages. First, our method exploits the topological structure created by ReLU NNs and fully explains the model behavior in small neighborhoods. For the $2^M$ cases created by a ReLU NN with $M$ neurons, MIP eliminates the searching branches using branch-and-bound. Our method, on the other hand, eliminates the searching branches by checking the feasibility of the local polytopes and their adjacency. Since a small traversing region often covers a limited number of polytopes, our algorithm has short running time when solving local problems. 

Second, since our algorithm explicitly identifies and visits all the polytopes, the final results contain not only the optimal solution but also the whole picture of the model behavior, providing explainability to the often-so-called black-box model.

Third, our method requires only linear and convex programming solvers and no MIP solvers. Identifying adjacent polytopes requires only linear programming. Convex programming may be used to solve the sub-problem within a local polytope. Our algorithm allows us to incorporate any convex programming solver that is most suitable for the sub-problem, providing much freedom to customize.

Last but probably the most important, our algorithm is highly versatile and flexible. Within each local polytope, the model is linear, which is often the simplest type of model to work with. Any analysis that one runs on a linear model can be transplanted here and wrapped inside the polytope traversing algorithm. Therefore, our algorithm provides a unified framework to verify different properties of piecewise linear networks.

\section{Case Studies} \label{sec:casestudies}

\begin{figure*}[t]
\center
\includegraphics[width=1.85\columnwidth]{fig_mnist}
\caption{\small Adversarial testing of a MNIST digit classification network w.r.t. 50 testing samples (5 samples per digit). The maximum change of an individual pixel value is (a) $+/-0.01$ or (b) $+/-0.05$. The blue lines and markers show the number of local polytopes traversed during verification. The red lines and markers show the time (in seconds) used. A dot marker indicates the network is robust w.r.t. the corresponding sample, while a cross marker indicates at least one adversarial sample can be found. Two adversarial samples are shown in (c).}
\label{fig:mnist}
\end{figure*}

\subsection{ACAS Xu}
We applied the polytope traversing algorithm to verify the safety of ACAS Xu networks \cite{julian2016policy}. The ACAS Xu networks contain an array of 45 ReLU NNs to issue advisories to avoid mid-air collisions for unmanned aircraft. This array of networks was developed to approximate a large lookup table traditionally used in an Airborne Collision Avoidance System so that the massive memory occupation by the table and the lookup time can be reduced. Each network takes five inputs: distance from ownship to intruder, angle from ownship to intruder, heading angle of the intruder w.r.t. ownship, speed of ownship, and speed of intruder. The five possible advisories output from each network are: Clear-of-Conflict (COC), week right, strong right, week left, and strong left. Each network contains six hidden layers with 50 neurons in each layer, resulting in a total of 300 neurons.

The appendix of \cite{katz2017reluplex} listed 10 desired properties that each network should satisfy. In our case study, we selected property II, III, and VI. Given a bounded set in the input space, these properties impose constraints on the rank of the networks' multi-class outputs. The verification of these properties can be formulated as a set of LP within each local polytope. We coded the polytope traversing algorithm in Python and used the LP solver in the SciPy package. Figure \ref{fig:acasxu} shows the verification results. The blue lines and markers show the number of local polytopes traversed during the verification. The red lines and markers show the total verification time in seconds. A dot marker indicates the corresponding network satisfies the property while a cross marker indicates the property is violated in at least one of the local polytopes. For property III and IV, the violating networks are identified after traversing only one of their local polytopes. For property II, most of the violating networks can be identified after traversing 10,000 local polytopes.

\subsection{MNIST}
We also applied the polytope traversing algorithm to verify the robustness of a MNIST digits classifier. The neural network \footnote{https://github.com/vtjeng/MIPVerify\_data/blob/master/weights/mnist/n1.mat} we tested takes a vectorized image as input which contains 784 pixels. It has two hidden layers of sizes 40 and 20 respectively. The output has a dimension of 10, corresponding to each possible digit. This network is trained using traditional techniques without special enforcement on robustness.

The robustness property requires the network's prediction to remain the same when a small perturbation is applied to the pixels of the original sample. In our test, we scaled all the pixel values to the range of 0 to 1. We tested two budget levels: 0.01 and 0.05, which means the maximum change each pixel can take is plus or minus the budget level while remaining inside the 0-1 range. We selected 50 samples from the testing dateset, five for each digit. We ran the polytope traversing algorithm until an adversarial sample was found or the network was verified w.r.t. the testing sample.

Figure \ref{fig:mnist} shows the robustness test results. As in previous experiments, we use blue and red lines to show the number of traversed local polytoeps and computational time (in seconds) respectively. A dot marker indicates the network is robust w.r.t. the corresponding sample, while a cross marker indicates at least one adversarial sample is found. Even under the small budget of 0.01, the network is not robust w.r.t. 11 out of the 50 testing samples. The number of local polytopes covered within the budget varies significantly with different original samples. When the budget is increased to 0.05, many adversarial samples can be found in the same polytope that the original samples fall into. While in some other cases more than 10,000 local polytopes are traversed before finding the first adversarial sample. Two adversarial samples are shown in Figure \ref{fig:mnist}.(c). The small perturbations hardly perceivable to us fool the neural network.

\section{Conclusion} \label{sec:conclusion}
We explored the unique topological structure that ReLU NNs create in the input space; identified the adjacency among the partitioned local polytopes; developed a traversing algorithm based on this adjacency; and proved the thoroughness of polytope traversing. Our polytope traversing algorithm could be extended to other piecewise linear networks such as those containing convolutional or maxpooling layers.

\section{Acknowledgments}
The authors would like to thank Lin Dong, Linwei Hu, Rahul Singh, and Han Wang from Wells Fargo, and Sihan Zeng from Georgia Institute of Technology for their valuable inputs and feedback on this project.

\bibliographystyle{IEEEbib}
\bibliography{references}

\section*{Appendix}
\subsection{Proof of Lemma \ref{def:redundant_ieq}}
\begin{lemma}
Given a set $\setR = \{ \vx | g_1(\vx) \leq 0,\ldots, g_M(\vx) \leq 0 \} \neq \emptyset$, then $g_m(\vx)$ is a redundant inequality if the new set formed by flipping this inequality is empty: $\hat{\setR} = \{ \vx | g_1(\vx) \leq 0, \ldots, g_{m}(\vx) \geq 0, \ldots, g_M(\vx) \leq 0 \} = \emptyset$.
\end{lemma}

\begin{proof}
Let $\tilde{\setR}$ be the set formed by removing inequality $g_m(\vx) \leq 0$: $\tilde{\setR} = \{ \vx | g_1(\vx) \leq 0, \ldots, g_{m-1}(\vx) \leq 0 ,g_{m+1}(\vx) \leq 0, \ldots, g_M(\vx) \leq 0 \}$. Then $\tilde{\setR} = \setR \cup \hat{\setR}$. If $ \hat{\setR}=\emptyset$, then $\setR = \tilde{\setR}$ and the inequality $g_m(\vx) \leq 0$ satisfies Definition \ref{def:redundant_ieq}.
\end{proof}

Note the other direction of Lemma \ref{them:redundant_ieq} may not hold. One example is when identical inequalities appear in the set: both inequalities in $\setR = \{ \vx | g_1(\vx)\leq0, g_2(\vx)\leq0 \}$ are redundant by definition if $g_1(\cdot)=g_2(\cdot)$. However, the procedure in Lemma\ref{them:redundant_ieq} will not identify them as redundant.

\subsection{Proof of Theorem \ref{them:traverseI}}
\begin{theorem}
Given a ReLU NN with one hidden layer of $M$ neurons as specified in (\ref{eq:relu_nn_I}), Algorithm \ref{algo:traverseI} covers all non-empty local polytopes created by the neural network. That is, for all $\vx \in \R^P$, there exists one $\setR_{\vc}$ as defined in (\ref{eq:polytope}) such that $\vx \in \setR_{\vc}$ and $\vc \in \setS_R$, where $\setS_R$ is the result returned by Algorithm \ref{algo:traverseI}.
\end{theorem}

\begin{proof}
Since each partitioning hyperplane divide $\R^P$ into two halfspaces, all $2^M$ activation patterns encoded by $\vc$ covers the entire input space. We construct a graph with $2^M$ nodes, each representing a possible polytope code. Some the nodes may correspond to an empty set due to conflicting inequalities. For each pair of non-empty polytope that are one-adjacent to each other, we add an edge to their corresponding nodes. What left to prove is that any pair of non-empty polytopes are connected.

W.l.o.g. assume two nodes with code $\vc$ and $\hat{\vc}$ that differ only in the first $K$ bits. Also assume the polytopes $\setR_{\vc}$ and $\setR_{\hat{\vc}}$ are both non-empty. We will show that there must exist a non-empty polytope $\setR_{\tilde{\vc}}$ that is one-adjacent to $\setR_{\vc}$ with code $\tilde{\vc}$ different from $\hat{\vc}$ in one of the first $K$ bits. As a result, $\tilde{\vc}$ is now one bit closer to $\hat{\vc}$.

We prove the claim above by contradiction. Assuming claim is not true, we flip any one of the first $K$ bits in $\setR_{\vc}$, and the corresponding polytope $\setR_{\tilde{\vc}^k}$ must be empty. By Definition \ref{def:redundant_ieq}, the inequality $(-1)^{c_m}\left(\vw_m^T\vx + b_m\right) \leq 0$, $m=1,2,\ldots,K$ must all be redundant, which means they can be removed from the set of constraints \cite{telgen1982minimal, telgen1983identifying}:
{\small
\begin{equation}
\begin{split}
\setR_{\vc} =& \{ \vx|(-1)^{c_m}\left(\vw_m^T\vx + b_m\right) \leq 0,\ m=1,2\ldots,M \} \\
=& \{ \vx|(-1)^{c_m}\left(\vw_m^T\vx + b_m\right) \leq 0,\ m=K+1,\ldots,M \}  \\
\supseteq &\{ \vx|(-1)^{c_m}\left(\vw_m^T\vx + b_m\right) \leq 0,\ m=1,2,\ldots,M \} \cup \\
&\{ \vx|(-1)^{c_m}\left(\vw_m^T\vx + b_m\right) \geq 0,\ m=1,\ldots,K, \\
&\quad\ \ (-1)^{c_m}\left(\vw_m^T\vx + b_m\right) \leq 0,\ m=K+1,\ldots,M \} \\
=& \setR_{\vc} \cup \setR_{\hat{\vc}} \ .
\end{split}
\label{eq:connected_proof}
\end{equation}
}%
The derived relationship in (\ref{eq:connected_proof}) plus the assumption that all $\setR_{\tilde{\vc}^k}$ must be empty lead to the conclusion that $\setR_{\hat{\vc}} = \emptyset$, which contradict with the non-empty assumption.

Therefore, for any two non-empty polytopes $\setR_{\vc}$ and $\setR_{\hat{\vc}}$, we can create a path from $\setR_{\vc}$ to $\setR_{\hat{\vc}}$ by iteratively finding an intermediate polytope whose code is one bit closer to $\hat{\vc}$. Since the polytope graph covers all input space and all non-empty polytopes are connected, BFS guarantees the thoroughness of traversing.
\end{proof}
\end{document}


Peer review:
A novel work for proposing polytope traversing algorithm
This paper points out that ReLU NNs will divide the input domain into many local polytopes. Based on this observation, this work develop a polytope traversing algorithm via BFS and apply in many aspects like local adversarial attacks. The idea is novel and exactly makes sense. 

However, I have some concerns about how this work can be applied to larger and more real tasks like image classification. First, the dimension of images is large, which may significantly increase the computational complexity of this algorithm. Second, since the sizes of different local polytopes may vary a lot, directly search polytopes via BFS may not find the optimal solution in some application like generating adversarial noises. Maybe the authors can attempt to test their methods in some image dataset like MNIST.
Rating: 7: Good paper, accept
Confidence: 3: The reviewer is fairly confident that the evaluation is correct