ReLU Neural Networks, Polyhedral Decompositions, and Persistent Homolog
A ReLU neural network leads to a finite polyhedral decomposition of input
space and a corresponding finite dual graph. We show that while this dual graph
is a coarse quantization of input space, it is sufficiently robust that it can
be combined with persistent homology to detect homological signals of manifolds
in the input space from samples. This property holds for a variety of networks
trained for a wide range of purposes that have nothing to do with this
topological application. We found this feature to be surprising and
interesting; we hope it will also be useful.
Yajing Liu, Christina M Cole, Chris Peterson, Michael Kirby
# ReLU Neural Networks, Polyhedral Decompositions, and Persistent Homology

###### Abstract

A ReLU neural network leads to a finite polyhedral decomposition of input space and a corresponding finite dual graph. We show that while this dual graph is a coarse quantization of input space, it is sufficiently robust that it can be combined with persistent homology to detect homological signals of manifolds in the input space from samples. This property holds for a variety of networks trained for a wide range of purposes that have nothing to do with this topological application. We found this feature to be surprising and interesting; we hope it will also be useful.

Machine Learning, Polyhedral Decompositions, Persistent Homology, Persistent Homology, Persistent Homology

## 1 Introduction

The rectified linear unit (ReLU) function is the default activation function for many well known, deep, feed forward neural networks (AlexNet, ResNet, Inception, SqueezeNet, etc). These networks frequently use some convolutional layers to decrease the number of parameters being trained and often utilize skips to decrease issues from overfitting. The importance of ReLU feed forward neural networks (FFNNs) has motivated an extensive investigation of their properties from a variety of aspects to demystify their "black-box nature". In this paper, we start from the observation that, regardless of their use of skips and/or convolutions, a ReLU FFNN \(F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}\) decomposes the input space (\(\mathbb{R}^{m}\)) into convex polyhedra and assigns to each polyhedron a unique binary vector that encodes the ReLU activation pattern of all nodes in the ReLU layers of the network. More precisely, the neural network assigns a binary vector to each input point, assigning identical binary vectors to input points lying in the same polyhedra. Despite this granularity imposed on input space, we found that the homology of manifolds embedded in the input space of the network can still be detected via persistent homology applied solely to distance measures built from the binary vectors of the points in the sampling. After explaining the connection between polyhedra and binary vectors in a ReLU FFNN, we describe methods for finding the binary vectors associated with the nearby neighbors of a given polyhedron in the decomposition. For small networks, one can (in principle) use these methods to determine all the polyhedra in the decomposition and determine their proximity within the decomposition. Toward this goal, we exploit the following observation: _two polyhedra share a facet if and only if their associated binary vectors differ in exactly one bit_. Due to this property, we propose the Hamming distance between binary vectors as a proxy for proximity between polyhedra in the decomposition. More precisely, if \(\mathcal{X}\) is a collection of data points in the input space, each \(x\in\mathcal{X}\) lies in one of the polyhedra and can be assigned the binary vector, \(s(x)\), that labels the polyhedron. Using the Hamming metric, we build a network driven distance matrix between pairs of binary vectors associated to data points in \(\mathcal{X}\). Using a software package such as Ripser (see [https://live.ripser.org](https://live.ripser.org)), we illustrate how the coarse distance measure corresponding to a polyhedral decomposition of input space can be used to extract homological information about a manifold lying in the input space from sample points on the manifold. It is worth noting that the Hamming distance between polyhedral bit vectors is an approximation for the smallest number of polyhedral steps between two polyhedra (i.e. for the shortest path between vertices on the dual graph of the polyhedral decomposition). (Jamil et al., 2023) has demonstrated the efficacy of the Hamming distance in illuminating the mechanisms underlying neural network functionality.

The paper is organized as follows: Section 2 introduces the basic background on ReLU FFNNs, the polyhedral decomposition associated to a network, the binary vectors assigned to each polyhedron, the linear model that determines each polyhedron, the linear program for filtering out redundant inequalities in the linear model, and the affine linear function that is attached to each polyhedron. Section 3 outlines algorithms for determining the binary vectors that occur for polyhedra in the entire input space and for polyhedra that occur within a bounded region. Section 4 describes how the polyhedral decomposition, realized through binary vectors, can be combined with persistent homology to uncover topological features in data. Section 5 concludes the work and
sketches directions for future research.

## 2 Neural Networks and Polyhedral Decompositions

This section gives the basic definitions of a ReLU neural network and describes the connection to polyhedral decompositions and binary vectors.

### Notation for Neural Networks

Consider an \((L+1)-\)layer ReLU feedforward neural network:

\[\mathbb{R}^{m}\ \xrightarrow[]{(W_{1},b_{1})}\ \mathbb{R}^{h_{1}}\ \xrightarrow[]{(W_{2},b_{2})} \mathbb{R}^{h_{2}}\to\ldots\to\mathbb{R}^{h_{L-1}}\] \[\xrightarrow[]{(W_{L},b_{L})}\ \mathbb{R}^{h_{L}}\xrightarrow[]{(W_{L+1},b_{L+1})} \mathbb{R}^{n}. \tag{1}\]

In this model, \(\mathbb{R}^{m}\) is the input space, \(\mathbb{R}^{n}\) is the output space, and \(h_{i}\) corresponds to the number of nodes at layer \(i\). Layer 0 corresponds to the input space and layer \(L+1\) corresponds to the output space (so \(h_{0}=m\) and \(h_{L+1}=n\)). We let \(W_{i}\in\mathbb{R}^{h_{i}\times h_{i-1}}\) and \(b_{i}\in\mathbb{R}^{h_{i}}\) denote the weight matrix and bias vector of layer \(i\), respectively. The activation functions for the hidden layers (layers \(1,\ldots,L\)) are assumed to be ReLU functions (applied coordinate-wise) while the map to the last layer (the output layer) is assumed to be affine linear (without a ReLU function being applied to the image). Recall that the ReLU function is the map \(ReLU:\mathbb{R}\to\mathbb{R}\) given by

\[ReLU(a)=\begin{cases}a&\text{if }a>0\\ 0&\text{if }a\leq 0.\end{cases} \tag{2}\]

The ReLU map is a map on real numbers that is piecewise linear and continuous. It can be naturally extended to a piecewise linear continuous map on vector spaces (which we also denote as ReLU). More precisely, we define \(ReLU:\mathbb{R}^{h_{i}}\to\mathbb{R}^{h_{i}}\), by applying the ReLU function to each coordinate of \(x\in\mathbb{R}^{h_{i}}\). Let \(w_{i,j}\) denote the \(j^{th}\) row of \(W_{i}\) and let \(b_{i,j}\) denote the \(j^{th}\) entry of \(b_{i}\). Given an input data point \(x\in\mathbb{R}^{m}\), we denote the output of \(x\) in layer \(i\) as \(F_{i}(x)\). Thus, with this notation we have \(F_{i}(x)\in\mathbb{R}^{h_{i}}\), \(F_{0}(x)=x\), and

\[\begin{split} F_{i}(x)&=\text{ReLU}(W_{i}F_{i-1}(x )+b_{i})\\ &=\begin{bmatrix}\max\{0,w_{i,1}F_{i-1}(x)+b_{i,1}\}\\ \vdots\\ \max\{0,w_{i,h_{i}}F_{i-1}(x)+b_{i,h_{i}}\}\end{bmatrix}.\end{split} \tag{3}\]

### Definitions of Binary Vectors

Consider model (1). Given an input data point \(x\in\mathbb{R}^{m}\), for each hidden layer \(i\) (so \(1\leq i\leq L\)), we introduce a binary (bit) vector

\[s_{i}(x)=[s_{i,1}(x)\ \ldots\ s_{i,h_{i}}(x)]^{\top}\in\mathbb{R}^{h_{i}},\]

where \(s_{i,j}(x)\) (with \(1\leq j\leq h_{i}\)) is defined as follows:

\[s_{i,j}(x)=\begin{cases}1&\text{if }w_{i,j}F_{i-1}(x)+b_{i,j}>0\\ 0&\text{if }w_{i,j}F_{i-1}(x)+b_{i,j}\leq 0.\end{cases} \tag{4}\]

Thus, for each point \(x\in\mathbb{R}^{m}\), we have a sequence of binary vectors \(s_{1}(x),s_{2}(x),\ldots,s_{L}(x)\). We can stack the binary vectors associated to \(x\) to make a long column vector

\[s(x)=[s_{1}^{\top}(x)\ \ldots\ s_{L}^{\top}(x)]^{\top}\in\mathbb{R}^{h}, \tag{5}\]

where \(h=\sum_{i=1}^{L}h_{i}\) is the total number of nodes in the hidden layers. We call \(s(x)\) the binary vector of \(x\).

Different points from the input space \(\mathbb{R}^{m}\) can have the same binary vector. We next show that the set of points that have the same binary vector, \(\{x^{\prime}:s(x^{\prime})=s(x),x^{\prime}\in\mathbb{R}^{m}\}\), form a convex polyhedron in \(\mathbb{R}^{m}\).

### Linear Model for Binary Vectors

Let \(s_{1},s_{2},\ldots,s_{L}\) denote a given sequence of binary vectors for model (1). For each layer \(i\)\((1\leq i\leq L)\), we describe inequality constraints that must be satisfied for an input data point to have the same binary vector \(s_{i}\). To describe the inequalities in a consistent manner, we introduce a sign vector of \(1^{\prime}s\) and \(-1^{\prime}s\) for each hidden layer. For layer \(i\), define \(s^{\prime}_{i}=[s^{\prime}_{i,1}\ \ldots\ s^{\prime}_{i,h_{i}}]^{\top}\) with

\[s^{\prime}_{i,j}=\begin{cases}1&\text{if }s_{i,j}=0\\ -1&\text{if }s_{i,j}=1.\end{cases} \tag{6}\]

In layer 1, because \(F_{1}(x)=\text{ReLU}(W_{1}x+b_{1})\), any data point \(x\) that has the bit vector \(s_{1}\) satisfies the following linear inequality:

\[\text{diag}(s^{\prime}_{1})(W_{1}x+b_{1})\leq 0 \tag{7}\]

where \(\text{diag}(v)\) is a square diagonal matrix with the elements of vector \(v\) on the main diagonal. The inequality (7) is established by a fundamental principle: the affine output generated by the first hidden layer, in response to input \(x\), must satisfy a greater-than-zero condition for nodes where the corresponding binary vector is set to 1, and a less-than-or-equal-to-zero condition for nodes where the corresponding binary vector is set to 0.

Suppose that \(x\) has the bit vector sequence \(s_{1},s_{2},\ldots,s_{L}\). Let \(\hat{W}_{j}=W_{j}\text{diag}(s_{j-1})\hat{W}_{j-1}\) and \(\hat{b}_{j}=W_{j}\text{diag}(s_{j-1})\hat{b}_{j-1}+b_{j}\) for \(2\leq j\leq L\) with \(\hat{W}_{1}=W_{1},\hat{b}_{1}=b_{1}\). By model (1), the following equation holds for \(1\leq j\leq L\):

\[F_{j}(x)=\text{ReLU}(W_{j}F_{j-1}(x)+b_{j})=\text{diag}(s_{j})(\hat{W}_{j}x+ \hat{b}_{j})\]
where \(F_{0}(x)=x\).

More generally, any data point \(x\) that has \(s_{j}\) as its bit vector for layer \(j\) should satisfy the following linear inequalities:

\[\text{diag}(s^{\prime}_{j})\hat{W}_{j}x\leq\text{diag}(s^{\prime}_{j})(-\hat{b}_ {j}). \tag{8}\]

Let \(A_{j}=\text{diag}(s^{\prime}_{j})\hat{W}_{j}\) and \(c_{j}=\text{diag}(s^{\prime}_{j})(-\hat{b}_{j})\) for \(1\leq j\leq L\). Combining (7) and (8), we have

\[Ax\leq c, \tag{9}\]

where

\[A=[A_{1}^{\top}\ A_{2}^{\top}\ \dots\ A_{L}^{\top}]^{\top}\ \text{and}\ c=[c_{1}^{\top}\ c_{2}^{ \top}\ \dots\ c_{L}^{\top}]^{\top}\]

with \(A_{i}\in\mathbb{R}^{h_{i}\times m}\) and \(c_{i}\in\mathbb{R}^{h_{i}}\). Note that the set \(P=\{x\in\mathbb{R}^{m}:Ax\leq c\}\) is a convex polyhedron (potentially empty, potentially unbounded). Considering the totality of points in the input space, the weights of a ReLU neural network lead to a decomposition of the input space into a collection of bounded and unbounded polyhedra and each polyhedron has a corresponding bit vector. Note that a random bit vector may or may not correspond to a non-empty polyhedra.

It is not difficult to see that for a polyhedron \(P\) of full dimension (the same dimension as the ambient space), there is a unique smallest subset of inequalities, which we denote by \((A^{\prime},c^{\prime})\), that one can obtain from \((A,c)\) that leaves the polyhedron \(P\) unchanged. Thus, \(A^{\prime}\) is built from a subset of rows of \(A\) and \(c^{\prime}\) is the corresponding subset of rows of \(c\) such that the set of points that satisfy \(Ax\leq c\) is the same as the set of points that satisfy \(A^{\prime}x\leq c^{\prime}\).

Write

\[A=\begin{bmatrix}a_{1}\ a_{2}\ \dots\ a_{h}\end{bmatrix}^{\top}\ \text{ and }\ \ c= \begin{bmatrix}c^{1}\ c^{2}\ \dots\ c^{h}\end{bmatrix}^{\top}\]

with \(a_{i}\in\mathbb{R}^{m}\) and \(c^{i}\in\mathbb{R}\). To determine if \(a_{i}x\leq c^{i}\) is a redundant constraint, first define

\[\tilde{A}=\begin{bmatrix}a_{1}\ a_{2}\ \dots\ a_{i-1}\ a_{i+1}\dots\ a_{h} \end{bmatrix}^{\top}\]

and

\[\tilde{c}=\begin{bmatrix}c^{1}\ c^{2}\ \dots\ c^{i-1}\ c^{i+1}\ \dots\ c^{h} \end{bmatrix}^{\top}.\]

Next, consider the following linear program

\[\begin{array}{l}\operatorname*{maximize}_{x}\ \ a_{i}^{\top}x\\ \text{s. t.}\ \ \tilde{A}x\leq\tilde{c}.\end{array} \tag{10}\]

If the optimal objective value of (10) is less than or equal to \(c^{i}\), then the \(i\)th linear inequality is redundant, and we can remove \(a_{i}^{\top}\) and \(c^{i}\) from \(A\) and \(c\), respectively. We determine \((A^{\prime},c^{\prime})\) by iterating this process to remove all redundant constraints.

Given a bit vector \(s\), the \(i\)-th entry of \(s\) is called active if the \(i\)th row of \(A\) is in \(A^{\prime}\), and inactive otherwise. For any input \(x\) in the polyhedron determined by the bit vector \(s\), the output of \(x\) is determined by the single affine map: \(G(x)=W_{L+1}\text{diag}(s_{L})\hat{W}_{L}x+W_{L+1}\text{diag}(s_{L})\hat{b}_{ L}+b_{L+1}\).

We summarize some of the points of this section, together with implications (many of which we leave to the reader) in the bullet points that follow. We emphasize that the polyhedral decomposition of input space has an associated dual graph. This graph has vertices corresponding to polyhedra and edges corresponding to polyhedra that share an \(m-1\) dimensional face. It is not hard to show that this dual graph is bipartite (refer to Appendix B). We will be using this graph, implicitly, to detect homological signals via persistent homology. For the interested reader, we point out the following references related to topological features and combinatorial features of the type of neural networks we are considering (Grigsby et al., 2022; Masden, 2022). We would also like to acknowledge the quickly-growing library of work in the area of polyhedral theory, a survey of which can be found in (Huchette et al., 2023).

The following is a summary of the key ideas presented either explicitly or implicitly in this section:

* A ReLU FFNN determines a continuous piecewise linear map \(F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}\).
* The neural network leads to a decomposition of \(\mathbb{R}^{m}\) into a collection of bounded and unbounded convex polyhedra.
* A binary vector, supported on the hidden nodes of the network, can be attached to each point in the domain of the network. The value of the binary vector associated to a given input \(x\in\mathbb{R}^{m}\), at a given node, is \(1\) if ReLU was not applied at the node and \(0\) if ReLU was applied at the node (the ReLU activation pattern).
* If \(h\) denotes the number of hidden nodes in the network, then there are, a priori, \(2^{h}\) possible binary vectors and thus \(2^{h}\) possible polyhedra in the decomposition of input space. In reality, the number of realized convex polyhedra in the domain is much smaller.
* If two points lie in the interior of the same convex polyhedron then they determine the same binary vector.
* If two points are in the interior of distinct convex polyhedra then they determine distinct binary vectors.
* From the ReLU activation pattern for points in a convex polyhedron, one can determine an affine linear equation that represents the behavior of the neural network on points in the polyhedron. If the polyhedron is denoted by \(P\), then for points in \(P\), the function \(F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}\) can be expressed as \(F(x)=A_{P}x+b_{P}\) for some matrix \(A_{P}\) and some vector \(b_{P}\).
* Input values that determine a value of \(0\) on some node, before applying ReLU, are precisely the input values that lie on the boundary of distinct convex polyhedra. This is due to the fact that applying ReLU to \(0\) has the same effect as not applying ReLU to \(0\).
* If two convex polyhedra, \(P_{1},P_{2}\), share an \((m-1)\)-dimensional face, then the binary vectors associated to each of the polyhedra differ in one bit. Furthermore, the affine linear functions \(A_{P_{1}}+b_{P_{1}}\) and \(A_{P_{2}}+b_{P_{2}}\) agree on this \((m-1)\)-dimensional face.
* Any polyhedral decomposition of \(\mathbb{R}^{m}\) has a natural dual graph with vertices corresponding to \(m\)-dimensional polyhedra and edges corresponding to polyhedra sharing an \(m-1\) dimensional facet.
* The Hamming distance between binary vectors (that represent two polyhedra) can be used as an approximation for the smallest number of polyhedral steps between the two polyhedra (i.e. the length of a minimal geodesic on the dual graph).

## 3 Algorithms for Bit Vector Search and Examples

The number of polyhedra in the input space or within a bounded region of the input space provides a measure of the network's expressivity and complexity. Upper and lower bounds on the maximal number of polyhedra obtainable from a given ReLU FFNN architecture can be found in (Pascanu et al., 2014), (Montufar et al., 2014), (Raghu et al., 2017), (Arora et al., 2018), and (Serra et al., 2018). Several algorithms ((Xiang et al., 2018), (Yang et al., 2020), and (Xu et al., 2021)) have been developed to compute the exact polyhedra decomposition of the input space through layer-by-layer linear inequality solving. Larger decomposition examples can be computed using a method developed by (Vincent and Schwager, 2021), which enumerates all polyhedra in the input space as follows:

* Start with a random point \(x\in\mathbb{R}^{m}\) and determine its bit vector \(s(x)\). This bit vector labels a polyhedron \(P\).
* Find the active bits in \(s(x)\) for the polyhedron \(P\).
* Each active bit corresponds to a neighboring polyhedron. Each neighboring polyhedron has a bit vector that can be obtained by "flipping" one of the active bits for \(P\). Thus, one can find all of the neighboring polyhedra for \(P\), in terms of their binary vectors.
* Repeat the process to find the neighbors of each of these newly identified polyhedra.
* The number of active bits for a polyhedron is equal to the number of nearest neighbors of the polyhedron. The previous steps continue until we have a list of polyhedra \(\mathcal{P}\) that satisfies the property that for each \(P\in\mathcal{P}\), the set of nearest neighbors of \(P\) is itself a subset of \(\mathcal{P}\).

This process leads to a set \(\mathcal{P}\) of convex polyhedra that decompose input space. Such decompositions have been used to define network-imposed distance metrics for model inputs (Balestriero and Baraniuk, 2018). We coarsen these metrics and show they can still be used to detect homological signals. Equating neural networks of various types (convolutional neural networks, residual networks, skip-connected networks, recurrent neural networks) to max-affine spline operators has been carried out by (Balestriero and Baraniuk, 2018). The paper (Sattelberg et al., 2020) investigated the behavior of the ReLU FFNNs based on the structure of the decomposition, together with the affine map attached to each polyhedron.

While there are multiple methods that can be used to determine the polyhedral decomposition, and their associated binary vectors, imposed by the weights of a ReLU neural network, all of these methods are woefully inadequate for the large and deep networks that are commonplace today. This is not a fault of the algorithms, many modern networks contain many millions of hidden nodes and have input spaces with dimension well beyond 100,000. There are just too many polyhedra represented by such networks. However, for small networks, any of the methods will work. We present two straightforward methods. The first is a brute-force method formulated as a Linear programming problem. The second, which is the traversal-and-pruning method outlined above, computationally improves upon the first. We demonstrate our methods with examples. Pseudocode for these two methods are provided in Appendix A.

### Algorithms

A mixed-integer linear program to count the number of polyhedra in a bounded region can be found in (Serra et al., 2018). We present a linear program that can also count polyhedra, determine if they share a facet, and whose implementational simplicity is useful. Our proposed method yields not only the count of polyhedra in the entire input space \(\mathbb{R}^{m}\) but also their respective binary vectors. To count the number of polyhedra in a bounded region, the linear inequality constraints in (15) need to be augmented to include bounds for variables. For example, when the input is 2-dimensional and the bounded region is determined by \(a_{1}\leq x_{1}\leq b_{1},\ a_{2}\leq x_{2}\leq b_{2}\), the constraints become \(\tilde{A}^{j}x\leq\tilde{c}^{j}\) where \(\tilde{A}^{j}\) is the concatenation of \(A^{j}\) and \(B\) while \(\tilde{c}^{j}\) is the concatenation of \(c^{j}\) and \(d\), with

\[B=\begin{bmatrix}1&-1&0&0\\ 0&0&1&-1\end{bmatrix}^{T}\ \operatorname{and}\ d=\begin{bmatrix}b_{1}\ -a_{1}\ b_{2}\ -a_{2}\end{bmatrix}^{T}.\]

Algorithm 1 is fast for small neural network structure. However, (Hanin and Rolnick, 2019) proved that the number of bit vectors that correspond to polyhedra is approximately \(h^{m}\) which is much smaller than \(2^{h}\) for large \(h\), highlighting the computational savings that can be procured by only traversing active bit vectors instead of all possibilities as is done by
the brute-force method. These savings are what Algorithm 2 has to offer.

Algorithm 1 is easy to implement but slow. The idea for Algorithm 2 is based on the fact that each polyhedron corresponding to \(Ax\leq c\) is determined by \(A^{\prime}x\leq c^{\prime}\) and two adjacent polyhedra that share one facet differ by one active bit. This active bit corresponds to the hyperplane containing this facet. The key step for Algorithm 2 is to find the active bits for a given bit vector using the method mentioned in Section 2.2. Once identified, these active bits can then be used to find all neighboring polyhedra, starting a scheme that ripples through the desired domain until all polyhedra are found. The idea behind this algorithm also motivated the reachable polyhedral marching algorithm in (Vincent and Schwager, 2021). We make a small modification to enumerate the polyhedra in a bounded region rather than in the entire, unbounded domain. For a bounded region, the bit vector can be expanded by checking whether the corresponding polyhedron hits the boundary or not. For example, when \(m=2\) and the bounded region is defined by \(a_{1}\leq x_{1}\leq b_{1},a_{2}\leq x_{2}\leq b_{2}\), the bit vector \(s\) can be expanded by 1 more bit with 1's or 0's by checking if its \(Ax\leq c\) solutions hit any of the domain boundaries defined by \(x_{1}=a_{1}\), \(x_{1}=b_{1}\), \(x_{2}=a_{2}\), or \(x_{2}=b_{2}\) (1 if true, 0 otherwise). For the bit vectors that have \(1\) for their last bit, we update their labels with 1 after they are added to \(\mathcal{P}\). The reason is that only if the bit vector of the initial point does not hit the region boundaries, the bit vectors that hit the region boundaries will be derived by flipping the active bits of their neighbors, which means that all its adjacent polyhedra in the bounded region are already in \(\mathcal{P}\).

### Examples

In this section, we visualize the polyhedral decomposition, determined by a ReLU FFNN, on a bounded region for two basic models:

\[\mathbb{R}^{2}\xrightarrow[\text{ReLU}]{(W_{1},b_{1})}\mathbb{R}^{3} \xrightarrow[\text{ReLU}]{(W_{2},b_{2})}\mathbb{R}^{3}\xrightarrow[\text{ReLU}]{(W_{3},b_{3})}\mathbb{R} \tag{11}\]

\[\mathbb{R}^{3}\xrightarrow[\text{ReLU}]{(W_{1},b_{1})}\mathbb{R}^{10} \xrightarrow[\text{ReLU}]{(W_{2},b_{2})}\mathbb{R}^{10}\xrightarrow[\text{ReLU }]{(W_{3},b_{3})}\mathbb{R}^{10}\xrightarrow[\text{ReLU}]{(W_{4},b_{4})} \mathbb{R} \tag{12}\]

We first applied model (11) to fit \(f_{1}(x_{1},x_{2})=x_{1}^{2}+x_{2}^{2}-2/3\) using 10000 points uniformly sampled in \([-1,1]^{2}\) and model (12) to fit \(f_{2}(x_{1},x_{2},x_{3})=(x_{1}-1)^{2}+2x_{2}^{2}+x_{3}^{2}+1\) using 125000 points uniformly sampled in \([-1,1]^{3}\). We used TensorFlow to train the above models with batch size of 50 and early stopping criterion based on the convergence of validation loss.

We used Algorithm 1 to enumerate all binary vectors in \(\mathbb{R}^{2}\) and in \([-1,1]^{2}\). We used Algorithm 2 to enumerate all binary vectors in \(\mathbb{R}^{3}\) and in \([-1,1]^{3}\). We plotted the polyhedra by finding their extremal vertices.

Figure 1(a) and Figure 1(b) provide visualizations of the polyhedral decomposition derived from model (11) and model (12), respectively. For Figure 1(a), there are 25 polyhedra in the bounded region \([-1,1]^{2}\) and 27 polyhedra in \(\mathbb{R}^{2}\). Of these polyhedra, 21 are bounded and 6 are unbounded. The binary vectors are superimposed onto the polyhedra. Note that the binary vectors associated to different polyhedra differ in exactly one bit if and only if they share a facet. For Figure 1(b), there are 1858 polyhedra in the bounded region \([-1,1]^{3}\) and 3331 polyhedra in \(\mathbb{R}^{3}\).

## 4 Persistent Homology and Polyhedral Decompositions

It has been well documented that the Euclidean distance between sampled points on a manifold in \(\mathbb{R}^{n}\) can be employed to detect the topology of the manifold. In this section, we provide a description of Vietoris-Rips persistent homology and illustrate how it can be effectively combined with the non-Euclidean distance measure, associated to the polyhedral decomposition, to also identify homological features. Persistent homology has been a rapidly developing branch of topology largely due to its usefulness in data analysis and machine learning (Ghrist, 2008; Zomorodian and Carlsson, 2004; Carlsson, 2009; Edelsbrunner and Harer, 2008) (a collection of additional resources and videos can be found at [https://www.aatrn.net](https://www.aatrn.net)). Work linking persistent homology and neural networks has been appearing with increasing frequency, please see the following for a sampling of some recent works in this direction (Rieck et al., 2018; Zhao et al., 2020; Carriere et al., 2020; Birdal et al., 2021). The description of persistent homology given below is extremely condensed. The interested reader is encouraged to read the survey article (Ghrist, 2008) where additional details can be found.

Figure 1: Visualizations of low-dimensional polyhedral decompositions. A detailed discussion on (a) can be found in Appendix B along with the exact binary vectors that label all of the regions of the decomposition.
### Persistent Homology

Consider a metric space, \(C\), with distance function \(d:C\times C\to\mathbb{R}\). Let \(x^{(1)},x^{(2)},\ldots x^{(N)}\) be points in \(C\). We can utilize \(d\) to build an \(N\times N\) matrix \(D\) by setting \(D_{i,j}=d(x^{(i)},x^{(j)})\). \(D\) will be hollow (i.e. \(D_{i,i}=0\)), symmetric, and non-negative. From \(D\) we can build a family, \(A(t)\), of \(N\times N\)\(\{0,1\}\)-matrices parameterized by a real parameter \(t\) using the rule

\[A(t)_{i,j}=\begin{cases}0&\text{if $D(i,j)>t$ or if $i=j$}\\ 1&\text{else}.\end{cases} \tag{13}\]

For each \(t\), \(A(t)\) can be viewed as an adjacency matrix of a graph \(G(t)\). Let \(CL(t)\) denote the clique complex of \(G(t)\)(Hausmann, 1995). By construction, \(CL(t)\) is a simplicial complex and there is a natural inclusion map

\[CL(t_{1})\hookrightarrow CL(t_{2})\]

whenever \(t_{1}<t_{2}\). Let \(\mathbb{F}\) be a field. A simplicial complex \(S\) has an associated chain complex \(S_{\bullet}\) of \(\mathbb{F}\)-vector spaces. The failure of the chain complex to be exact at location \(i\) is measured by the \(i^{th}\) homology \(H_{i}(S_{\bullet},\mathbb{F})\) (which itself is an \(\mathbb{F}\)-vector space). The inclusion map \(CL(t_{1})\hookrightarrow CL(t_{2})\) induces a chain map

\[CL(t_{1})_{\bullet}\to CL(t_{2})_{\bullet}.\]

Whenever you have a chain map \(F:S_{\bullet}\to T_{\bullet}\) between chain complexes, \(S_{\bullet},T_{\bullet}\) you get associated linear maps between \(H_{i}(S_{\bullet},\mathbb{F})\) and \(H_{i}(T_{\bullet},\mathbb{F})\) for each \(i\). Thus, the chain map \(CL(t_{1})_{\bullet}\to CL(t_{2})_{\bullet}\) induces, for each \(i\), a linear map

\[H_{i}(CL(t_{1})_{\bullet},\mathbb{F})\to H_{i}(CL(t_{2})_{\bullet},\mathbb{F}).\]

If we pick values \(t_{1}<t_{2}<\cdots<t_{k}\), then we can build a nested sequence of simplicial complexes

\[CL(t_{1})\subset CL(t_{2})\subset\cdots\subset CL(t_{k})\]

which leads to

\[H_{i}(CL(t_{1})_{\bullet},\mathbb{F})\to H_{i}(CL(t_{2})_{\bullet},\mathbb{F}) \rightarrow\cdots\to H_{i}(CL(t_{k})_{\bullet},\mathbb{F}), \tag{14}\]

where the arrows denote linear maps.

Each of the \(H_{i}(CL(t_{j})_{\bullet},\mathbb{F})\) are finite dimensional \(\mathbb{F}\)-vector spaces and we can view (14) as a finitely generated graded \(\mathbb{F}[x]\) module, where \(\mathbb{F}[x]\) is the ring of polynomials. This module is frequently called the \(i^{th}\) persistence module associated to the nested sequence of simplicial complexes. Since \(\mathbb{F}[x]\) is a principal ideal domain, the \(F[x]\) module has a decomposition into its invariant factors by the well known structure theorem for finitely generated graded modules over a principal ideal domain. Each invariant factor can be viewed as a homology class that has a birth time \(t_{b}\) and a death time \(t_{d}\) (possibly infinite) and this invariant factor can be represented as an interval \([t_{b},t_{d}]\). This collection of intervals corresponds to the barcode representation of the invariant factors of a persistence module.

### Combining Persistence with a Polyhedral Decomposition

First we provide two examples based on the polyhedral decomposition of model (12).

_Example 4.1_.: We consider the circle in \(\mathbb{R}^{3}\) (the input space) with parameterization given by \((0,cos(t),sin(t))\). We sampled the circle at \(20\) evenly spaced points \(x^{(1)},x^{(2)},\ldots,x^{(20)}\). Each of these points lie in one of the polyhedra from Figure 1(b). A picture of the polyhedra encountered by the \(20\) sample points is found in Figure 2(a). We recorded the binary vectors for each of the encountered polyhedra. This gave a total of \(20\) binary vectors but only \(19\) were distinct. We labeled these distinct binary vectors \(s^{(1)},s^{(2)},\ldots s^{(19)}\). We built a \(19\times 19\) matrix \(E\) by setting \(E_{i,j}\) equal to the Hamming distance between \(s^{(i)}\) and \(s^{(j)}\). The resulting matrix is hollow, symmetric, and has positive integers in entries off the diagonal. We input this matrix into Ripser (see [https://live.ripser.org](https://live.ripser.org)) and asked for the \(H_{0}\) and \(H_{1}\) barcodes. The result of this experiment can be found in Figure 3(a).

_Example 4.2_.: We considered the same circle in \(\mathbb{R}^{3}\) as in the previous example but sampled at \(500\) evenly spaced points \(x^{(1)},x^{(2)},\ldots,x^{(500)}\). A picture of the polyhedra encountered by the \(500\) sample points is found in Figure 2(b). We recorded the \(500\) binary vectors for each of the polyhedra that was hit by a data point. Only \(41\) of the bit vectors were distinct (corresponding to \(41\) distinct polyhedra) and we labeled these \(s^{(1)},s^{(2)},\ldots s^{(41)}\). We built a \(41\times 41\) matrix \(F\) by setting \(F_{i,j}\) equal to the Hamming distance between \(s^{(i)}\) and \(s^{(j)}\). We input this matrix into Ripser (see [https://live.ripser.org](https://live.ripser.org)) and asked for the \(H_{0}\) and \(H_{1}\) barcodes. The result of this experiment can be found in Figure 3(b).

The \(H_{0}\) barcode from the first example indicates connectivity at Hamming distance 4. The \(H_{1}\) barcode indicates a spurious (i.e. short lived) closed loop occurring at Hamming distance 3. The homological signal of the circle appears (and is quite strong) at Hamming distance 4. The \(H_{0}\) barcode from the second example indicates connectivity at

Figure 2: Polyhedra from two sample sizes
Hamming distance 2. The \(H_{1}\) barcode indicates a long lived loop beginning at Hamming distance 2.

In the following two examples, we will consider a similar example but utilizing real images from the ImageNet validation dataset (which contains 50K images). We will calculate the Hamming distance between bit vectors of data points via a much deeper neural network. The network we use is known as ResNet-50, it is a 50-layer convolutional neural network and was pre-trained on ImageNet (Deng et al., 2009). The training images are 224x224x3 thus the input space has dimension 150528. It uses ReLU as an activation function on the outputs from the convolutional layers and contains more than 6,000,000 nodes in its many layers. Like in the previous examples, the activation pattern of each ReLU layer is stored as a bit vector (as defined in Section 2.2).

_Example 4.3_.: We started with 3 pictures, denoted by \(A_{1},A_{2},A_{3}\), chosen from the ImageNet validation dataset. They are photos of a miniature poodle, a Persian cat, and a Saluki (see Figure 4). Each photo is represented by a 224 x 224 x 3 array. We generated \(50\) data points using the formulas \(\sin\theta A_{1}+\cos\theta A_{2}\) and \(A_{3}+\sin\theta A_{1}+\cos\theta A_{2}\), respectively, with \(\theta\) consisting of \(50\) points uniformly sampled from \([0,2\pi]\). We calculated the bit vectors via ResNet-\(50\) for the \(50\) data points and label them as \(s^{(1)},s^{(2)},\ldots s^{(50)}\). We built a \(50\times 50\) distance matrix \(G\) by setting \(G_{i,j}\) equal to the Hamming distance between \(s^{(i)}\) and \(s^{(j)}\). We input this matrix into Ripser which returned the \(H_{0}\) and \(H_{1}\) barcodes. The result of this experiment can be found in Figure 5.

We note that there is homological noise (or homological dust) in Figure 5(a) but none in Figure 5(b). In both examples, there is a strong signal representing the circle. The sample points in the second example have all positive entries while the sample points in the first example do not have this property. We were unsure why (or if) this is related to the homological noise but we found it interesting nevertheless. It is a potentially useful feature that the input space to a neural network simultaneously has two kinds of distance measures. The first derives from the standard Euclidean geometry and the second derives from the coarse geometry implied by the Hamming distance. If one applies an isometry to a data set then its pairwise distance matrix will not change. However, if one applies an isometry with respect to one metric but measure distance via the second metric then one can definitely observe a change. It may be useful to combine information from multiple measurements. One way to carry this out is by using a monotone Boolean function (i.e Boolean functions built using only **and** and **or** operations). The next example illustrates this approach using a small rotation as the isometry in order to produce two distance matrices arising from essentially the same data set.

_Example 4.4_.: Using the same \(A_{1},A_{2}\) with Example 4.3, we uniformly generated 50 points using \(\sin\theta A_{1}+\cos\theta A_{2}\) from \([0,2\pi]\) and \([1,2\pi+1]\), respectively. We calculated the bit vectors via ResNet-50 for the \(50\) data points sampled from \([0,2\pi]\) and labelled them as \(s^{(1)},s^{(2)},\ldots s^{(50)}\). Similarly, we label the bit vectors of the \(50\) data points sampled from \([1,2\pi+1]\) as \(\bar{s}^{(1)},\bar{s}^{(2)},\ldots\bar{s}^{(50)}\). We built four distance matrices \(G^{1},G^{2},G^{3}\), and \(G^{4}\) with \(G^{1}_{i,j}\) equal to the Hamming distance between \(s^{i}\) and \(s^{j}\), \(G^{2}_{i,j}\) equal to the Hamming distance between \(\bar{s}^{i}\) and \(\bar{s}^{j}\), \(G^{3}_{i,j}\) equal to the maximum of the Hamming distance between \(s^{i}\) and \(s^{j}\) and between \(\bar{s}^{i}\) and \(\bar{s}^{j}\), and \(G^{4}_{i,j}\) equal to minimum of the Hamming distance between \(s^{i}\) and \(s^{j}\) and between \(\bar{s}^{i}\) and \(\bar{s}^{j}\). The \(H_{1}\) barcodes for the four cases are presented in Figure 6

The distance matrix corresponding to the max function corresponds to requiring that distance matrices \(G_{1}\) and \(G_{2}\) both have their corresponding entry below some threshold. In other words, this corresponds to applying the **and** Boolean function. Similarly, the min function corresponds to applying the **or** Boolean function. Each seems to improve some feature related to the homological noise but it is hard to

Figure 4: Images of \(A_{1},A_{2},\) and \(A_{3}\)

Figure 3: \(H_{0}\) and \(H_{1}\) barcode plots
say which of the two is better. In the next example, we see an overall strengthening of the length of the homological signal for the max function and we see an earlier start of the homological signal for the min function. Homological noise did not make an appearance in the next example.

_Example 4.5_.: Using the same \(A_{1},A_{2},A_{3}\) with Example 4.3, we uniformly generated 50 points using \(A_{3}+\sin\theta A_{1}+\cos\theta A_{2}\) from \([0,2\pi]\) and \([1,2\pi+1]\), respectively. As in Example 4.4, we generated 4 different distance matrices \(\hat{G}^{1}\), \(\hat{G}^{2}\), \(\hat{G}^{3}\), and \(\hat{G}^{4}\), respectively. The \(H_{1}\) barcodes for the four cases are presented in Figure 7.

## 5 Conclusion and Future Work

A ReLU feedforward neural network induces a finite polyhedral decomposition of the input space. The corresponding dual graph represents this decomposition, where vertices correspond to polyhedra and edges represent shared facets. Data in the input space gets mapped to vertices of this graph. The geometry of the graph can be exploited, via persistent homology, to detect homological signals of manifolds from points sampled from the manifold. Many techniques in data analysis build on the premise that data, sharing a collection of common features, tends to be well approximated by a low dimensional geometric object. It is conjectured that the coarseness of the polyhedral decomposition can be helpful in dealing with the noise present in many data sets (maybe by including monotone Boolean functions to combine multiple distance matrices). In future work, we hope to extend the examples in this paper with the goal of detecting manifolds, known to exist in data, such as \(SO(3)\), higher dimensional tori, and various fiber bundles involving these manifolds. A torus example can be found in Appendix C.

Neural networks of the kind considered in this paper have a finite number of "polyhedral resources". The training function, the distribution of training data, the network architecture, and the training method are some of the factors in how the network utilizes its polyhedra. In the networks we considered, we observed the general pattern that there was a dense collection of smaller volume polytopes near the training data, larger volume polytopes in the general vicinity but away from the training data, with a shell of unbounded polyhedra surrounding the polytopes. You can see hints of this tendency in Figure 1(a). In future work, we hope to make more precise characterizations in the direction of this observation. In particular, we are eager to extend this work to depictions that more accurately quantify the number of polytopes that emerge from a neural network and their distributions of volumes, the "polytopes landscape".

## 6 Acknowledgements

This work is partially supported by the United States Air Force under Contract No. FA865020C1121 and the DARPA Geometries of Learning Program under Award No. HR00112290074.