Hybrid Zonotopes Exactly Represent ReLU Neural Networks
We show that hybrid zonotopes offer an equivalent representation of
feed-forward fully connected neural networks with ReLU activation functions.
Our approach demonstrates that the complexity of binary variables is equal to
the total number of neurons in the network and hence grows linearly in the size
of the network. We demonstrate the utility of the hybrid zonotope formulation
through three case studies including nonlinear function approximation, MPC
closed-loop reachability and verification, and robustness of classification on
the MNIST dataset.
Joshua Ortiz, Alyssa Vellucci, Justin Koeln, Justin Ruths
# Hybrid Zonotopes Exactly Represent ReLU Neural Networks

###### Abstract

We show that hybrid zonotopes offer an equivalent representation of feed-forward fully connected neural networks with ReLU activation functions. Our approach demonstrates that the complexity of binary variables is equal to the total number of neurons in the network and hence grows linearly in the size of the network. We demonstrate the utility of the hybrid zonotope formulation through three case studies including nonlinear function approximation, MPC closed-loop reachability and verification, and robustness of classification on the MNIST dataset.

## I Introduction

We leverage the recently introduced concept of a hybrid zonotope to develop an equivalent representation of feed-forward fully connected neural networks whose activation functions are rectified linear activation functions (ReLU) [1]. We show that this exact and analytic expression of the input-output mapping of the neural network enables verification and robustness quantification for general neural networks and, when neural networks are embedded in closed-loop control applications, facilitates closed-loop reachability analysis and safety guarantees. This work was developed separately from, and in parallel to, a recently posted article [2]. The overall goal of our paper and use of hybrid zonotopes is similar, however, we show our approach leads to linear growth in complexity, whereas the complexity in [2] grows exponentially.

To verify neural networks, a wide range of methods have been proposed to upper and lower bound the output of neural networks. In [3], reactivation upper and lower bounds are calculated and propagated layer-by-layer to calculate an overall Lipschitz constant bound for the entire network. Due to their speed these iterative approaches continue to be used (e.g., [4]). However, it requires that bounds be recomputed if the input set changes and the iterative bounding erases any memory that exists between the layers. Semi-definite programming relaxations have been used to regain the layer-to-layer history in a computationally efficient manner [5] and mixed-integer linear programming (MILP) has been used to compute exact Lipschitz bounds [6]. In the context of control systems, closed-loop reachability has been investigated with similar methods [7, 8].

Hybrid zonotopes are a generalization of constrained zonotopes, which were introduced several years earlier [9]. Like our paper here, authors have leveraged constrained zonotopes [10, 11] or polytopes [12] to propagate sets through approximated ReLU activation functions to compute bounds on the output reachability of neural networks. [13] demonstrates how these analytic expression of relaxations of ReLU neural networks can be used to, for example, train networks to be more robust to adversarial attack. Star sets offer similar expressiveness as hybrid zonotopes and have also been used to provide approximate and exact reachability of feed-forward ReLU neural networks [14]. However, like [2], the computational complexity of the mixed-integer program (number of binary variables) grows exponentially.

At the heart of our approach is the ability of hybrid zonotopes to maintain an exact analytic expression from input to output. Crucially, this expression preserves how previous layers impact the propagation of the sets without iterative calculations or overapproximation. In contrast to MILP approaches (e.g., [15]), which also provide exact input-output relationships for ReLU neural network certification, the hybrid zonotope representation, due to the memory feature of zonotopes, enables it to be used in a variety of ways, including for closed loop analysis.

**Notation:** Matrices are denoted by uppercase letters, e.g., \(G\in\mathbb{R}^{n\times n_{g}}\), and sets by uppercase calligraphic letters, e.g., \(\mathcal{Z}\subset\mathbb{R}^{n}\). Vectors and scalars are denoted by lowercase letters. The \(n\)-dimensional unit hypercube is denoted by \(\mathcal{B}_{\infty}^{n}=\left\{x\in\mathbb{R}^{n}\mid\left\|x\right\|_{ \infty}\leq 1\right\}\). The set of all \(n\)-dimensional binary vectors is denoted by \(\{-1,1\}^{n}\). Matrices and vectors of all \(0\) and \(1\) elements are denoted by \(\mathbf{0}\) and \(\mathbf{1}\), respectively, of appropriate dimension. The Kronecker product of matrices \(A\in\mathbb{R}^{m\times n_{g}}\) and \(B\in\mathbb{R}^{p\times q}\) is given as \(A\otimes B\in\mathbb{R}^{pm\times qn}\).

Given the sets \(\mathcal{Z},\mathcal{W}\subset\mathbb{R}^{n},\mathcal{Y}\subset\mathbb{R}^{m}\), and matrix \(R\in\mathbb{R}^{m\times n}\), the linear mapping of \(\mathcal{Z}\) by \(R\) is \(R\mathcal{Z}=\left\{Rz\mid z\in\mathcal{Z}\right\}\), the Minkowski sum of \(\mathcal{Z}\) and \(\mathcal{W}\) is \(\mathcal{Z}\oplus\mathcal{W}=\left\{z+w\mid z\in\mathcal{Z},w\in\mathcal{W}\right\}\), the generalized intersection of \(\mathcal{Z}\) and \(\mathcal{Y}\) under \(R\) is \(\mathcal{Z}\cap_{R}\mathcal{Y}=\left\{z\in\mathcal{Z}\mid Rz\in\mathcal{Y}\right\}\), and the Cartesian product of \(\mathcal{Z}\) and \(\mathcal{Y}\) is \(\mathcal{Z}\times\mathcal{Y}=\left\{(z,y)\mid z\in\mathcal{Z},\,y\in\mathcal{Y}\right\}\).

## II Hybrid Zonotopes

**Definition 1**: _[_1_, Def. 3]_ _The set \(\mathcal{Z}_{h}\subset\mathbb{R}^{n}\) is a hybrid zonotope if there exists \(G^{c}\in\mathbb{R}^{n\times n_{g}}\), \(G^{b}\in\mathbb{R}^{n\times n_{b}}\), \(c\in\mathbb{R}^{n}\), \(A^{c}\in\mathbb{R}^{n_{c}n_{g}}\), \(A^{b}\in\mathbb{R}^{n_{c}n_{b}}\), and \(b\in\mathbb{R}^{n_{c}}\) such that_

\[\mathcal{Z}_{h}=\left\{\begin{bmatrix}G^{c}\,G^{b}\Big{]}\begin{bmatrix}\xi^{c} \\ \xi^{b}\end{bmatrix}+c\begin{bmatrix}\xi^{c}\\ \xi^{b}\end{bmatrix}\in\mathcal{B}_{\infty}^{n_{g}}\times\{-1,1\}^{n_{b}},\\ \begin{bmatrix}A^{c}\,A^{b}\end{bmatrix}\begin{bmatrix}\xi^{c}\\ \xi^{c}\end{bmatrix}=b\end{bmatrix}\right\}. \tag{1}\]

A hybrid zonotope is the union of \(2^{\mathbb{R}}\) constrained zonotopes corresponding to the possible combinations of binary factors, \(\xi^{b}\). The hybrid zonotope is given in _Hybrid Constrained Generator-representation_ and the shorthand notation of \(\mathcal{Z}_{h}=(G^{c},G^{b},c,A^{c},A^{b},b)\in\mathbb{R}^{n}\) is used to denote the set given by (1). Continuous and binary _generators_ refer to the columns of \(G^{c}\) and \(G^{b}\), respectively. A hybrid zonotope
with no binary generators is a constrained zonotope, \(\mathcal{Z}_{c}=(G,c,A,b)\subset\mathbb{R}^{n}\), and a hybrid zonotope with no binary generators and no constraints is a zonotope, \(\mathcal{Z}=(G,c)\subset\mathbb{R}^{n}\). Identities and time complexity of linear mappings, Minkowski sums, generalized intersections, and generalized half-space intersections are reported in [1, Section 3.2]. An identity and time complexity for Cartesian products is given in [16]. Preliminary methods for removing redundant generators and constraints of a hybrid zonotope were reported in [1] and further developed in [16].

_Example 1 (Zonotope Memory (see also [17])):_ If we consider a dynamic system \(x_{k+1}=Ax_{k}+Bu_{k}\), with \(x_{0}\in\mathcal{X}_{0}=(G_{x}^{c},G_{x}^{b},c_{x},A_{x}^{c},A_{x}^{b},b_{x})\) and \(u_{0}\in\mathcal{U}=(G_{u}^{c},G_{u}^{b},c_{u},A_{u}^{b},A_{u}^{b},b_{u})\), then by the linear mapping and Minkowski sum identities,

\[\begin{split} x_{1}\in\mathcal{X}_{1}=\Bigg{\{}& [AG_{x}^{c}\;BG_{u}^{c}],[AG_{x}^{b}\;BG_{u}^{b}],Ac_{x}+Bc_{u},\\ &\begin{bmatrix}A_{x}^{c}\\ &A_{u}^{c}\end{bmatrix},\begin{bmatrix}A_{x}^{b}\\ &A_{u}^{b}\end{bmatrix},\begin{bmatrix}b_{x}\\ b_{u}\end{bmatrix}\Bigg{\}}.\end{split} \tag{2}\]

The effect of the linear mapping is encoded in the transformation of the continuous and binary generators. The memory is, however, captured in the fact that the continuous and binary factors that define the sets \(\mathcal{X}_{0}\) and \(\mathcal{U}\) are embedded in the definition of \(\mathcal{X}_{1}\). To ignore this connection - what happens with iterative methods - is to erase this memory by treating the factors of \(\mathcal{X}_{1}\) as new unrelated factors. This feature of memory is especially visible when creating extended vectors, e.g.,

\[\begin{split}\begin{bmatrix}x_{0}\\ x_{1}\end{bmatrix}\in\Bigg{\{}&\begin{bmatrix}G_{x}^{c}\\ AG_{x}^{c}\;\;BG_{u}^{c}\end{bmatrix},&\begin{bmatrix}G_{x}^{b}\\ AG_{x}^{b}\;\;BG_{u}^{b}\end{bmatrix},&\begin{bmatrix}c_{x}\\ Ac_{x}+Bc_{u}\end{bmatrix},\\ &\begin{bmatrix}A_{x}^{c}\\ A_{x}^{c}\\ &A_{u}^{c}\end{bmatrix},&\begin{bmatrix}A_{x}^{b}\\ A_{x}^{b}\\ A_{u}^{b}\end{bmatrix},&\begin{bmatrix}b_{x}\\ b_{x}\\ b_{u}\end{bmatrix}\Bigg{\}}.\end{split}\]

Specifically, the continuous and binary factors that specify the hybrid zonotope \(\mathcal{X}_{1}\) are \([(\xi_{x}^{c})^{\top}\;(\xi_{u}^{b})^{\top}]^{\top}\) and \([(\xi_{x}^{b})^{\top}\;(\xi_{u}^{b})^{\top}]^{\top}\), respectively, where \(\xi_{x}^{c}\) and \(\xi_{x}^{b}\) specify \(\mathcal{X}_{0}\) and \(\xi_{u}^{c}\) and \(\xi_{u}^{b}\) specify \(\mathcal{U}\). There are no new factors introduced specific to \(\mathcal{X}_{1}\) - all are inherited from \(\mathcal{X}_{0}\) or \(\mathcal{U}\). In this case the redundant constraints can be dropped.

## III Neural Network Representation as a Hybrid Zonotope

Consider an \(L\)-layered feed-forward fully-connected neural network \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) mapping inputs \(x\in\mathbb{R}^{n}\) to outputs \(y=f(x)\in\mathbb{R}^{m}\) such that

\[\begin{split} x^{0}&=x,\\ x^{\ell+1}&=\phi(W^{\ell}x^{\ell}+b^{\ell}), \qquad\ell\in\{0,\cdots,L-1\},\\ y&=f(x)=W^{L}x^{L}+b^{L},\end{split} \tag{3}\]

where \(W^{\ell}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}}\) and \(b^{\ell}\in\mathbb{R}^{n_{\ell+1}}\) are the weight matrix and bias vector between layers \(\ell\) and \(\ell+1\), with \(n_{0}=n\) and \(n_{L+1}=m\). For this paper, all activation functions \(\phi\) are ReLU functions that operate element-wise, i.e., for the pre-activation vector \(v^{\ell+1}=W^{\ell}x^{\ell}+b^{\ell}\in\mathbb{R}^{n_{\ell+1}}\),

\[\phi(v^{\ell})=[\varphi(v_{1}^{\ell})\cdots\varphi(v_{n_{\ell}}^{\ell})]^{\top}, \tag{4}\]

where the ReLU function is defined as \(\varphi(v_{i})=\max\{0,v_{i}\}\).

The goal of this paper is to show that such a feed-forward neural network with only ReLU activation units can be exactly - and analytically - represented by a hybrid zonotope. Core to our approach is the idea that _sets_ (including hybrid zonotopes) can be used to describe _functions_. We can accomplish this by expressing an extended set that expresses the input-to-output relationship. Given the input set \(\mathcal{X}\subset\mathbb{R}^{n}\) and output set \(\mathcal{Y}=\{f(x):x\in\mathcal{X}\}\subset\mathbb{R}^{m}\), the function \(f\) over the input set \(\mathcal{X}\) can be expressed by the set \(\mathcal{F}\subset\mathbb{R}^{n+m}\) of extended vectors \([x^{\top}\;y^{\top}]^{\top}\in\mathbb{R}^{n+m}\) such that \(\mathcal{F}=\{[x^{\top}\;y^{\top}]^{\top}:y=f(x),x\in\mathcal{X}\}\). Zonotope-based sets especially enable function representations through the "memory" engendered by the factors. The crucial connection is that some of the factors that define the output zonotope set are factors derived from the input zonotope.

To begin, we represent a single ReLU activation function \(x_{i}=\varphi(v_{i})\) as a set of points \([v_{i}\;x_{i}]^{\top}\in\mathbb{R}^{2}\) over a predetermined domain \(v_{i}\in[-a,a]\). We assume \(a>0\) is chosen large enough to capture the largest anticipated absolute value of the input \(v_{i}\) such that \(|v_{i}|\leq a\).

_Lemma 1:_ The set \(\Phi\subset\mathbb{R}^{2}\) of points satisfying the ReLU activation function over the domain \(\mathcal{D}_{i}=[-a,a]\) can be exactly represented as a hybrid zonotope with \(n_{g}=4\) continuous generators, \(n_{b}=1\) binary generator, and \(n_{c}=2\) constraints.

To satisfy the ReLU function for \(v_{i}\in[-a,a]\), \(\Phi=\{[v_{i}\;x_{i}]^{\top}:x_{i}=\max\{0,v_{i}\},|v_{i}|\leq a\}\). Note that \(\Phi=\Phi_{1}\cup\Phi_{2}\), where \(\Phi_{1}=\{[v_{i}\;x_{i}]^{\top}:-a\leq v_{i}\leq 0,x_{i}=0\}\) and \(\Phi_{2}=\{[v_{i}\;x_{i}]^{\top}:0\leq v_{i}\leq a,x_{i}=v_{i}\}\). To represent this set, consider the zonotopes

\[\begin{split}\mathcal{Z}_{1}&=\left(\begin{bmatrix} \nicefrac{{a}}{{a}}&\nicefrac{{a}}{{2}}\\ \nicefrac{{a}}{{2}}&0\end{bmatrix},\begin{bmatrix}0\\ \nicefrac{{a}}{{2}}&0\end{bmatrix},\begin{bmatrix}-a\\ -a/2\end{bmatrix}\right),\\ \mathcal{Z}_{3}&=\left(\begin{bmatrix}\nicefrac{{a}}{{a}}&\nicefrac{{a}}{{2}}\\ \nicefrac{{a}}{{2}}&0\end{bmatrix},\begin{bmatrix}a\\ a/2\end{bmatrix}\right).\end{split} \tag{5}\]

As shown in Figure 1, \(\mathcal{Z}_{1}\), \(\mathcal{Z}_{2}\), and \(\mathcal{Z}_{3}\) are the same parallelogram with shifted centers. Additionally, Figure 1

Fig. 1: The intersection of the blue zonotope, \(\mathcal{Z}_{1}\), and red hybrid zonotope, \(\mathcal{Z}_{2}\cup\mathcal{Z}_{3}\), is equivalent to the ReLU function (purple) over the interval \([-a,a]\) (here shown for \(a=10\)).
shows that \(\mathcal{Z}_{1}\cap\mathcal{Z}_{2}=\Phi_{1}\) and \(\mathcal{Z}_{1}\cap\mathcal{Z}_{3}=\Phi_{2}\). Therefore, \(\Phi=\{\mathcal{Z}_{1}\cap\mathcal{Z}_{2}\}\cup(\mathcal{Z}_{1}\cap\mathcal{Z}_{ 3})=\mathcal{Z}_{1}\cap(\mathcal{Z}_{2}\cup\mathcal{Z}_{3})\). Note that \(\mathcal{Z}_{2}\cup\mathcal{Z}_{3}\) can be directly expressed as the unconstrained hybrid zonotope

\[\mathcal{Z}_{2}\cup\mathcal{Z}_{3}=\left\{\begin{bmatrix}a/2&a/2&0\\ a/2&0&0\end{bmatrix},\begin{bmatrix}a\\ a/2\end{bmatrix},\begin{bmatrix}0\\ 0\end{bmatrix},[\,],[\,],[\,]\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \,\
Proof:: The proof is similar to that of **Corollary 1**. First, we can write (10) as \(x_{k+1}=[A\;B][x_{k}^{\top}\;u^{\top}]^{\top}\). Then we note that \([x_{k}^{\top}\;u^{\top}]^{\top}\in\mathcal{F}\cap_{[\mathbf{I}_{n}\;\mathbf{0} _{n\times m}]}\mathcal{R}_{k}\) such that \(x_{k+1}\in[A\;B](\mathcal{F}\cap_{[\mathbf{I}_{n}\;\mathbf{0}_{n\times m}]} \mathcal{R}_{k})\).

The result of **Theorem 2** can be repeated recursively starting from \(\mathcal{R}_{0}=\mathcal{X}_{0}\) to compute the reachable set \(\mathcal{R}_{k}\) at any time step \(k\). We omit the simple proof to satisfy space constraints.

**Corollary 2**: _Assume the initial conditions set \(\mathcal{R}_{0}\) is represented as a hybrid zonotope with \(n_{g,0}\) continuous generators, \(n_{b,0}\) binary generators, and \(n_{c,0}\) constraints and that the neural network with \(n_{N}\) ReLU activation functions is represented as a hybrid zonotope \(\mathcal{F}\) with \(4n_{N}\) continuous generators, \(n_{N}\) binary generators, and \(3n_{N}\) constraints. Then the hybrid zonotope set representation complexity of the reachable set \(\mathcal{R}_{k}\) at time step \(k\) grows linearly with \(k\), resulting in \(n_{g,k}=n_{g,0}+4kn_{N}\) continuous generators, \(n_{b,k}=n_{b,0}+kn_{N}\) binary generators, and \(n_{c,k}=n_{c,0}+k(n+3n_{N})\) constraints._

In the following result, we demonstrate how the reachable sets over time can be stacked to provide additional interpretability with regard to the history of the control policy used. This result will be used in the context of the MPC case study in Section V-B.

**Corollary 3**: _Assume the initial conditions set \(\mathcal{R}_{0}\) is represented as a hybrid zonotope with \(n_{g,0}\) continuous generators, \(n_{b,0}\) binary generators, and \(n_{c,0}\) constraints and that the neural network with \(n_{N}\) ReLU activation functions is represented as a hybrid zonotope \(\mathcal{F}\) with \(4n_{N}\) continuous generators, \(n_{N}\) binary generators, and \(3n_{N}\) constraints. Then the hybrid zonotope set representation of the aggregate reachable set \([x_{0}^{\top}\;x_{1}^{\top}\;...\;x_{k}^{\top}]^{\top}\in\mathcal{R}_{0:k}= \mathcal{R}_{0}\times\mathcal{R}_{1}\times\cdots\times\mathcal{R}_{k}\) has complexity \(n_{g,k}=n_{g,0}+4kn_{N}\) continuous generators, \(n_{b,k}=n_{b,0}+kn_{N}\) binary generators, and \(n_{c,k}=n_{c,0}+(n+3n_{N})\sum_{k=0}^{k}\kappa\) constraints._

Proof:: Following from Corollary 2, the hybrid zonotope \(\mathcal{R}_{0:k}\) is the Cartesian product of sets \(\mathcal{R}_{\kappa}\) for \(\kappa=1,2,\ldots,k\) each with \(n_{g,\kappa}=n_{g,0}+4\kappa n_{N}\) continuous generators, \(n_{b,\kappa}=n_{b,0}+\kappa n_{N}\) binary generators, and \(n_{c,\kappa}=n+n_{c,0}+3\kappa n_{N}\) constraints. If these sets were unrelated, the resulting hybrid zonotope would have \(\sum_{\kappa=0}^{k}n_{g,\kappa}\) continuous generators, \(\sum_{\kappa=0}^{k}n_{b,\kappa}\) binary generators, and \(\sum_{\kappa=0}^{k}n_{c,\kappa}\).

However, from the proof of Theorem 2, \(\mathcal{R}_{\kappa+1}=[A\;B](\mathcal{F}\cap_{[\mathbf{I}_{n}\;\mathbf{0}_{ n\times m}]}\mathcal{R}_{\kappa})\), hence \(\mathcal{R}_{\kappa+1}\) and \(\mathcal{R}_{\kappa}\) share common factors. Namely, all continuous and binary factors that characterize \(\mathcal{R}_{\kappa}\) also contribute to characterize \(\mathcal{R}_{\kappa+1}\). Thus the extended vector \([x_{\kappa}^{\top}\;x_{\kappa+1}^{\top}]^{\top}\) belongs to a hybrid zonotope with \(\max(n_{g,\kappa},n_{g,\kappa+1})\) continuous and \(\max(n_{b,\kappa},n_{b,\kappa+1})\) binary generators. The constraints do stack and become \(n_{c,\kappa}+n_{c,\kappa+1}\).

Extending this pairwise relationship forward, \(\mathcal{R}_{0:k}\) has \(\max(n_{g,0},\ldots,n_{g,k})=n_{g,k}=n_{g,0}+4kn_{N}\) continuous generators, \(\max(n_{b,0},\ldots,n_{b,k})=n_{b,k}=n_{b,0}+4kn_{N}\) binary generators, and \(n_{c,0}+(n+3n_{N})\sum_{\kappa=0}^{k}\kappa\) constraints.

## V Applications & Demonstrations

The neural networks in this paper are feed-forward fully-connected neural networks with only ReLU activation units, and specified by the number of layers and width of each layer, e.g., [4,10,7,2] is a network with inputs in \(\mathbb{R}^{4}\), outputs in \(\mathbb{R}^{2}\), and with two hidden layers of 10 and 7 ReLU units, respectively. The networks are trained in MATLAB using the stochastic gradient descent optimizer with momentum (0.95) and 100 epochs. The hybrid zonotopes have been coded in MATLAB. Optimization problems have been solved using GUROBI [18]. These examples are conducted on a laptop computer using one core of an 1.9 GHz Intel i7 processor and 16GB of RAM.

### _Nonlinear Function Approximation_

Over the input domain \([-5,5]\times[-5,5]\) we train a ReLU feed-forward neural network with layer sizes [2,20,10,10,1] to approximate the function \(\cos(x_{1})+\sin(x_{2})\) by grid sampling 400 points.

There are \(20+10+10=40=n_{N}\) total ReLU functions in the network and we build a hybrid zonotope \(\mathcal{F}\) capturing the input-output mapping using the extended vector \([x^{\top}y^{\top}]^{\top}\). The hybrid zonotope has \(n+4n_{N}=2+4\times 40=162\) continuous generators, \(n_{N}=40\) binary generators, and \(3n_{N}=120\) constraints. Of the \(2^{40}\) possible combinations of binary factors, 829 satisfy the constraints, leading to 829 constrained zonotopes which appear as facets of the surface

Fig. 2: (Left) Hybrid zonotope of the neural network that approximates the function \(\cos(x_{1})+\sin(x_{2})\). (Right) Top-down view of the surface reveals the intricate faceting that corresponds to binary factors (equivalent to ReLU activations) being on/off. The colormap of the original function is overlaid to show how the faceting follows the contours.
plotted in Figure 2. There is a one-to-one correspondence between these 829 feasible binary factor combinations and the combinations of the 40 ReLU units being active (\(>0\)) or inactive (\(\leq 0\)).

### _Model Predictive Control Policy & Closed-Loop Reachability Analysis_

We consider a double integrator system discretized with a sampling time of 1 second,

\[x_{k+1}=\begin{bmatrix}1&1\\ 0&1\end{bmatrix}x_{k}+\begin{bmatrix}0.5\\ 1\end{bmatrix}u_{k}. \tag{11}\]

As in [2, 19, 8], an MPC control policy is designed to stablize the system to the origin while respecting the state and input constraints, \(x_{k}\in[-5,5]\times[-1,1]\) and \(u_{k}\in[-1,1]\). A prediction horizon of 10 steps is used, with state and input weighting matrices \(Q=I_{2}\) and \(R=I\), the terminal region \(\mathcal{O}_{\infty}^{\text{LOR}}\), and with terminal penalty matrix \(P_{\infty}\), the solution of the discrete time algebraic Riccati equation. We grid the state space over \([-0.5,3]\times[-1,1]\) to produce 10,201 state-input pairs to train a neural network with layer sizes \([2,8,4,1]\).

The mapping of state-to-control provided by the neural network is visualized as a hybrid zonotope surface in three dimensions in Figure 3 (right). This hybrid zonotope, with \(n=2\) and \(n_{N}=4+8=12\), has 50 continuous factors, 12 binary factors, and 36 constraints. Each facet of the surface corresponds to a different combination of binary factors. Although there are \(2^{12}\) different factor combinations, there are only 32 facets in the surface. The remaining binary factor combinations do not satisfy the linear constraints.

We now use the hybrid zonotope representation of the trained neural network to analyze the closed-loop reachability of the system, similar to the approach taken for an explicit MPC controller [20]. We consider an initial state (the gray square in Figure 3), \(x(0)\in[2.5,3]\times[-0.25,0.25]\). Following Corollary 3, we construct the hybrid zonotope for the extended vector of states across time, \([x_{0}^{\top}\cdots x_{5}^{\top}]^{\top}\). Despite the 60 binary generators in the extended state hybrid zonotope, there are only 8 binary factor combinations (8 different sequences of linear control policies) that satisfy the constraints. Part of the utility of the extended state vector is that the initial set (gray square) is originally specified as an unconstrained zonotope with a center and two continuous generators. As seen in Figure 3, the initial set is faceted by constraints at later time steps to reveal the sets of initial conditions that correspond to different sequences of linear feedback control policies.

Suppose that we wish to check that the system reaches the goal set \([-0.25,0.25]\times[-0.25,0.25]\) (green square in Figure 3) for any initial condition in the gray region. Although this can be confirmed by visual inspection, it can also be rigorously ensured by evaluating set containment. The containment check of two hybrid zonotopes can be posed as a feasibility evaluation of a mixed-integer linear program.

### _Classification Robustness on MNIST_

The MNIST dataset is a canonical classification problem in which \(28\times 28\) pixel images of handwritten numbers are classified into the digits 0-9. For clarity of presentation, we consider here the classification task to classify the digits "1" and "7". We train a ReLU feedforward network with layer sizes [784,5,5,1] on a corpus of 13,007 images of the digits 1 and 7, in which the input is the vectorized image (stacking columns) and the output value \(+1\) denotes a "1" and the output value \(-1\) denotes a "7". The trained network achieves an accuracy of 99.5% on a test bank of 2,163 images of the digits 1 and 7. Figure 4 (left) shows the neural network output over the test image set, with outliers denoting images that are classified poorly and - in a handful of cases - incorrectly. Inset are example images of "the 1 that looks most like a 7" (worst 1) and "the 7 that looks most like a 1" (worst 7). These examples show that the output range over the images that are 1s spans nearly the entire interval [-1,1] (same is true for 7s) despite the overall good performance of 99.5%.

The hybrid zonotope representation of the network allows us to explore this output range more rigorously as a function of the input space. We define the _1s input space_ as a zonotope

Fig. 3: (Right) The hybrid zonotope representation of the neural network trained to learn an MPC policy to reach the origin while avoiding constraints, e.g., avoid \(x\leq-1\). (Left) The evolution of the system state following the inputs produced by the neural network model. The initial state set at \(k=0\) is shown in gray, and steps \(k=1,\ldots,5\) alternate color between blue and red. The inset plot shows (green arrow) the break between the \(k=4\) and \(k=5\) reachable set to show that the final state set sits within the goal set (green box).
whose center is the mean value of the vectorized images of the digit 1 (both training and testing). Defining \(\sigma\in\mathbb{R}^{784}\) as the standard deviation of the same vectorized images, the (continuous) generator matrix of the 1s input set is \(G_{c}=\text{diag}(\sigma)\). It is important to note that while the input set has been defined by the images, the input set is now the interval box "containing" one standard deviation away from the mean image in 784-dimensional space. Thus this 1s input set contains the entire continuous region around the nominal image of the digit 1. We follow the same procedure to produce the _7s input space_. With these input spaces, the hybrid zonotope representing the neural network has \(784+4\times 10=824\) continuous generators, \(10\) binary generators, and \(30\) constraints.

We now consider an adjustment to these input sets such that \(G_{c}=\text{diag}(\alpha\sigma)\), where \(0<\alpha\leq 1\) specifies a fractional value of the standard deviation. When \(\alpha\) approaches zero, the input set approaches a single point, located at the center (mean) value. Given the classification accuracy of the neural network, this mean point will be reliably classified correctly, hence we expect the classification output range of the 1s input set with small \(\alpha\) to be very near \(+1\) and the output range of the 7s input set with small \(\alpha\) to be very near \(-1\). As \(\alpha\) grows, we expect that the output range to grow since we observed that misclassifications occur in the test dataset. Figure 4 (right) demonstrates how the hybrid zonotope representation of the neural network allows exact characterization of the upper and lower bounds of the output hybrid zonotope subject to the 1s input set and the 7s input set. The calculation of the upper and lower bounds along each dimension of a hybrid zonotope can be cast as a mixed-integer linear optimization problem.

Figure 4 (right) demonstrates that although the classification accuracy is good on the test images, relatively small perturbations (in this case perturbations on the order of \(0.3\sigma\)) are sufficient to compromise the classification of some images. The sample images from the 7s input set shows that sample images from input sets corresponding to \(\alpha=0.3\) are still human-identifiable as their correct digit. Being able to quantify the output bounds on a network over a given input domain is a powerful tool to probe the robustness of a neural network. This is especially true in the context of adversarial attacks on neural networks, in which small adjustments can lead to large changes in output value (e.g., classification outcome) [21].

## VI Conclusions

In this paper, we have shown that it is possible to exactly represent a ReLU feed-forward neural network as a hybrid zonotope. This approach provides the ability to characterize variation and robustness of the mapping given by the neural network and offers rigorous ways to certify safety and reachability in closed-loop control application.

A direction for future work is to formulate the Lipschitz constant for the neural network directly from the hybrid zonotope representation. Also, although the binary factors scale linearly with the number of ReLU neurons in the network, focusing on ways to best leverage advancements in mixed-integer linear optimization solvers will unlock the ability to analyze larger and larger networks.

## VII Acknowledgments

We thank Neera Jain and Trevor Bird for sharing their MATLAB code for hybrid zonotopes.
