Certified Invertibility in Neural Networks via Mixed-Integer Programming
Neural networks are known to be vulnerable to adversarial attacks, which are
small, imperceptible perturbations that can significantly alter the network's
output. Conversely, there may exist large, meaningful perturbations that do not
affect the network's decision (excessive invariance). In our research, we
investigate this latter phenomenon in two contexts: (a) discrete-time dynamical
system identification, and (b) the calibration of a neural network's output to
that of another network. We examine noninvertibility through the lens of
mathematical optimization, where the global solution measures the ``safety" of
the network predictions by their distance from the non-invertibility boundary.
We formulate mixed-integer programs (MIPs) for ReLU networks and $L_p$ norms
($p=1,2,\infty$) that apply to neural network approximators of dynamical
systems. We also discuss how our findings can be useful for invertibility
certification in transformations between neural networks, e.g. between
different levels of network pruning.
Tianqi Cui, Thomas Bertalan, George J. Pappas, Manfred Morari, Ioannis G. Kevrekidis, Mahyar Fazlyab
# Certified Invertibility in Neural Networks

###### Abstract

Neural networks are notoriously vulnerable to adversarial attacks - small imperceptible perturbations that can change the network's output drastically. In the reverse direction, there may exist large, meaningful perturbations that leave the network's decision unchanged (excessive invariance, noninvertibility). We study the latter phenomenon in two contexts: (a) discrete-time dynamical system identification, as well as (b) calibration of the output of one neural network to the output of another (neural network matching). We characterize noninvertibility through the lens of mathematical optimization, in which the global solution quantifies the "safety" of the network predictions: their distance from the noninvertibility boundary. For ReLU networks and \(L_{p}\) norms (\(p=1,2,\infty\)), we formulate these optimization problems as mixed-integer programs (MIPs) that apply to neural network approximators of dynamical systems. We also discuss the applicability of our results to invertibility certification in transformations between neural networks (e.g. at different levels of pruning).

ArXiv 1-24

## 1 Introduction

Despite achieving high performance in a variety of classification and regression tasks, neural networks are not always guaranteed to satisfy certain desired properties after training. A prominent example is adversarial robustness. Neural networks can be overly sensitive to carefully designed input perturbations (Szegedy et al. (2013)). This intriguing property holds in the reverse direction too. In classification problems, neural networks can also be excessively insensitive to large perturbations, causing two semantically different inputs (e.g., images) to be classified in the same category (Jacobsen et al. (2018)). Indeed, a fundamental trade-off has been shown between adversarial robustness and excessive invariance (Tramer et al. (2020)), which is mathematically related to the noninvertibility of the map defined by the neural network.
To mitigate noninvertibility, and hence excessive invariance, one can consider invertible-by-design architectures. Invertible neural networks (INNs) have been used to design generative models (Donahue and Simonyan (2019)), implement memory-saving gradient computation (Gomez et al. (2017)), and solve inverse problems (Ardizzone et al. (2018)). However, commonly-used INN architectures suffer from exploding inverses; in this paper, we therefore consider the problem of certifying the (possible) noninvertibility of conventional neural networks after training. Specifically, we study two relevant invertibility problems: _(i) local invertibility of neural networks:_ given a dynamical system whose time-\(\tau\) map is parameterized by a neural network, we verify whether it is locally invertible around a certain input (or trajectory) and compute the largest region of local invertibility; and _(ii) local invertibility of transformations between neural networks:_ we certify whether two (assumed "equivalent") neural networks (e.g., related through pruning) can be transformed (i.e. calibrated) to each other locally via an invertible transformation. We develop mathematical tools based on mixed-integer linear/quadratic programming for the characterization of noninvertibility that are applicable to both (a) neural network approximators of dynamics, as well as to (b) transformations between neural networks.

Related workNoninvertibility in neural networks was studied in the 1990s (Gicquel et al. (1998); Rico-Martinez et al. (1993)); more recently, several papers focus on the global invertibility property in neural networks (see Chang et al. (2018); Teshima et al. (2020); Chen et al. (2018); MacKay et al. (2018); Jaeger (2014)). Analyzing invertibility of neural networks (Behrmann et al. (2018)) and constructing invertible architectures arises in many contexts, such as generative modeling (Chen et al. (2019)), inverse problems (Ardizzone et al. (2019)) or probabilistic inference (Radev et al. (2020)). Neural networks invertible by design have been developed for these applications. Some of the these networks (e.g. RevNet (Gomez et al. (2017)), NICE (Dinh et al. (2015)), real NVP (Dinh et al. (2017))) partition the input domains and use affine or coupling transformations as the forward pass, keeping the Jacobians (block-)triangular with nonzero diagonal elements, resulting in nonzero determinants; others, like i-ResNet (Behrmann et al. (2019)) have no analytical forms for the inverse dynamics, yet their finite bi-Lipschitz constants can be derived: both methods can guarantee global invertibility. A comprehensive analysis is found in (Behrmann et al. (2021); Song et al. (2019)). However, a theoretical understanding of the expressiveness of these architectures, as well as of their universal approximation properties, is still incomplete. Compared to standard networks like multi-layer perceptrons (MLPs) or convolutional neural networks (CNNs), the novel invertible neural networks (INNs) become computationally demanding. Neural ODE (Chen et al. (2018)) use an alternative method to compute gradients for backward propagation; i-ResNet (Behrmann et al. (2019)) has restrictions on the norm of every weight matrix to be enforced during the training process. In most cases, the input domain of interest is a small subset of the whole space. For example, the grey-scale image domain in computer vision problems is \([0,1]^{H\times W}\) (where \(H\) and \(W\) are height and width of images), and it is unnecessary to consider the whole space \(\mathbb{R}^{H\times W}\). We thus focus on _local invertibility_: how do we know if our network is invertible on a given finite domain, and if not, how do we quantify noninvertibility?

Beyond classification problems, noninvertibility can also lead to catastrophic consequences in regression, and more specifically in dynamical systems prediction. The flow of smooth differential equations is invertible when it exists; yet traditional numerical integrators used to approximate them can be noninvertible. Neural network approximations of the corresponding time-\(\tau\) map also suffer
from this potential pathology. In this paper, we initially study noninvertibility in the context of dynamical systems predictions.

## 2 Local invertibility of dynamical systems and neural networks

Continuous-time dynamical systems, in particular autonomous ordinary differential equations (ODEs) have the form \(dX(t)/dt=f(X(t)),X(t=t_{0})=X_{0}\), where \(X(t)\in\mathbb{R}^{m}\) are the state variables of interest; \(f:\mathbb{R}^{m}\mapsto\mathbb{R}^{m}\) relates the states to their time derivatives and \(X_{0}\in\mathbb{R}^{m}\) is the initial condition at \(t_{0}\). If \(f\) is uniformly Lipschitz continuous in \(X\) and continuous in \(t\), the Cauchy-Lipschitz theorem guarantees the existence and uniqueness of the solution.

In practice, we observe the states \(X(t)\) at discrete points in time, starting at \(t_{0}=0\). For a fixed timestep \(\tau\in\mathbb{R}^{+}\), and \(\forall n\in\mathbb{N}\), \(t_{n}=n\tau\) denotes the \(n\)-th time stamp, and \(X_{n}=X(t=t_{n})\) the corresponding state values. Now we will have:

\[X_{n+1}:=F(X_{n})=X_{n}+\int_{t_{n}}^{t_{n+1}}f(X(t))dt;\;X_{n}=F^{-1}(X_{n+1}). \tag{1}\]

This equation also works as the starting point of many numerical ODE solvers.

For the time-\(\tau\) map in (1), the inverse function theorem provides a sufficient condition for its invertibility: If \(F\) is a continuously differentiable function from an open set \(\mathcal{B}\) of \(\mathbb{R}^{m}\) into \(\mathbb{R}^{m}\), and the Jacobian determinant of \(F\) at \(p\) is non-zero, then \(F\) is invertible near \(p\). Thus, if we define the _noninvertibility locus_ as the set \(J_{0}(F)=\{p\in\mathcal{B}:\det(\mathbf{J}_{F}(p))=0\}\), then the condition \(J_{0}(F)=\emptyset\) guarantees global invertibility of \(F\) (notice that this condition is not _necessary_: the scalar function \(F(X)=X^{3}\) provides a counterexample). If \(F\) is continuous over \(\mathcal{B}\) but not everywhere differentiable, then the definition of \(J_{0}\) set should be altered to:

\[J_{0}(F)=\{p\in\mathcal{B}:\forall N_{0}(p),\exists\,p_{1},p_{2}\in N_{0}(p), p_{1}\neq p_{2},\;\text{s.t.}\;\det(\mathbf{J}_{F}(p_{1}))\det(\mathbf{J}_{F}(p_{2})) \leq 0\}\,., \tag{2}\]

the set of points where the determinant discontinuously changes sign.

Numerical integrators are (often) noninvertibleNumerically approximating the finite integral in (1) can introduce noninvertibility in the transformation. Here is a simple one-dimensional illustrative ODE example: \(dX/dt=f(X)=X^{2}+bX+c,\quad X(t=0)=X_{0}\), where \(b,c\in\mathbb{R}\) are two fixed parameters. The analytical solution (1) is invertible; however a forward-Euler discretization with step \(\tau\) gives

\[X_{n+1}=F(X_{n})=X_{n}+\tau(X_{n}^{2}+bX_{n}+c)\Rightarrow\tau X_{n}^{2}+(\tau b +1)X_{n}+(\tau c-X_{n+1})=0. \tag{3}\]

Given a fixed \(X_{n+1}\), Equation (3) is quadratic w.r.t. \(X_{n}\); this determines the local invertibility of \(F\) based on \(\Delta=(\tau b+1)^{2}-4\tau(\tau c-X_{n+1})\): no real root if \(\Delta<0\); one real root with multiplicity 2 if \(\Delta=0\); and two distinct real roots if \(\Delta>0\). In practice, one uses small timesteps \(\tau\ll 1\) for accuracy/stability, leading to the last case: there will always exist a solution \(X_{n}\) close to \(X_{n+1}\), and a second preimage, far away from the region of our interest, and arguably physically irrelevant (this second \(X_{n}\rightarrow-\infty\) as \(\tau\to 0\)). On the other hand, as \(\tau\) grows, the two roots move closer to each other, \(J_{0}(F)\) moves close to the regime of our simulations, and noninvertibility can have visible implications on the predicted dynamics. Thus, choosing a small timestep in explicit integrators guarantees desirable accuracy, and simultaneously _practically_ mitigates noninvertibility pathologies in the dynamics.
Invertibility in transformations between neural networksTraining two neural networks for the same regression or classification task practically never gives identical network parameters. Numerous criteria exist for comparing the performance of different models (e.g. accuracy in classification, or mean-squared loss in regression). Here we explore whether two different models _can be calibrated to each other_ (leading to a _de facto_ implicit function problem). Extending our analysis provides invertibility guarantees for the transformation from the output of network 1 to the output of network 2 (and vice versa).

## 3 Invertibility certification of neural networks and of transformations between them

Here we pose the verification of local invertibility of continuous functions as an optimization problem. We then show that for ReLU networks, this leads to a mixed-integer linear/quadratic program. For an integer \(q\geq 1\), we denote the \(L_{q}\)-ball centered at \(x_{c}\) by \(\mathcal{B}_{q}(x_{c},r)=\{x\in\mathbb{R}^{n}\mid\|x-x_{c}\|_{q}\leq r\}\) (the notation also holds when \(q\to+\infty\)).

**Problem 1** (Local Invertibility of NNs): _Given a neural network \(f:\mathbb{R}^{m}\mapsto\mathbb{R}^{m}\) and a point \(x_{c}\in\mathbb{R}^{m}\) in the input space, we want to find the largest radius \(r>0\) such that \(f\) is invertible on \(\mathcal{B}_{q}(x_{c},r)\), i.e., \(f(x_{1})\neq f(x_{2})\) for all \(x_{1},x_{2}\in\mathcal{B}_{q}(x_{c},r)\), \(x_{1}\neq x_{2}\)._

Another relevant problem is to verify whether, for a particular point, a nearby point exists with the same forward image. This is of particular interest in assessing invertibility of discrete-time dynamical systems around a given trajectory. We formally state the problem as follows:

**Problem 2** (Pseudo Local Invertibility of NNs): _Given a neural network \(f:\mathbb{R}^{m}\mapsto\mathbb{R}^{m}\) and a point \(x_{c}\in\mathbb{R}^{m}\) in the input space, we want to find the largest radius \(R>0\) such that \(f(x)\neq f(x_{c})\) for all \(x\in\mathcal{B}_{q}(x_{c},R)\), \(x\neq x_{c}\)._

If \(r\) and \(R\) are the optimal radii in problems 1 and 2 respectively, we must have \(r\leq R\). For Problem 1, the ball \(\mathcal{B}_{q}(x_{c},r)\) just "touches" the \(J_{0}\) set; for Problem 2, the ball \(\mathcal{B}_{q}(x_{c},R)\) extends to the "other" closest preimage of \(f(x_{c})\). Figure 1 illustrates both concepts in the one-dimensional case. For the scalar function \(y=f(x)\) and around a particular input \(x_{c}\), we show the nearest bounds of local invertibility and pseudo invertibility. The points \(Q_{1}=(x_{Q_{1}},y_{Q_{1}})\) and \(Q_{2}=(x_{Q_{2}},y_{Q_{2}})\) are the two closest turning points (elements of the \(J_{0}\) set) to the point \(C=(x_{c},y_{c})\); \(f\) is uniquely invertible (bi-Lipschitz) on the open interval \((x_{Q_{1}},x_{Q_{2}})\), so that the optimal solution to Problem 1 is: \(r=\min\{|x_{Q_{1}}-x_{c}|,|x_{Q_{2}}-x_{c}|\}=|x_{Q_{1}}-x_{c}|\). Noting that \(M_{1}=(x_{M_{1}},y_{M_{1}})\) and \(M_{2}=(x_{M_{2}},y_{M_{2}})\) are the two closest points that have the same \(y\)-coordinate as the point \(C=(x_{c},y_{c})\), the optimal solution to Problem 2 is \(R=\min\{|x_{M_{1}}-x_{c}|,|x_{M_{2}}-x_{c}|\}=|x_{M_{1}}-x_{c}|\).
We now state our first result, posing the local invertibility of a function (such as a neural network) as a constrained optimization problem.

**Theorem 1** (Local Invertibility of Continuous Functions): _Let \(f\colon\mathbb{R}^{m}\to\mathbb{R}^{m}\) be a continuous function and \(\mathcal{B}\subset\mathbb{R}^{m}\) be a compact set. Consider the following optimization problem,_

\[p^{\star}\leftarrow\max\quad\|x_{1}-x_{2}\|\quad\text{subject to }x_{1},x_{2} \in\mathcal{B},\quad f(x_{1})=f(x_{2}). \tag{4}\]

_Then \(f\) is invertible on \(\mathcal{B}\) if and only if \(p^{\star}=0\)._

**Theorem 2** (Pseudo Local Invertibility): _Let \(f\colon\mathbb{R}^{m}\to\mathbb{R}^{m}\) be a continuous function and \(\mathcal{B}\subset\mathbb{R}^{m}\) be a compact set. Suppose \(x_{c}\in\mathcal{B}\). Consider the following optimization problem,_

\[P^{\star}\leftarrow\max\quad\|x-x_{c}\|\quad\text{subject to }x\in\mathcal{B}, \quad f(x)=f(x_{c}). \tag{5}\]

_Then we have \(f(x)\neq f(x_{c})\) for all \(x\in\mathcal{B}\setminus\{x_{c}\}\) if and only if \(P^{\star}=0\)._

Note that by adding the equality constraints \(x=x_{1},x_{c}=x_{2}\) to the optimization problem (4), we obtain the optimization problem (5). Hence, we will only focus on (4) in what follows.

Invertibility certification of ReLU networks via mixed-integer programmingWe now show that for a given ball \(\mathcal{B}_{\infty}(x_{c},r)\) in the input space, and piecewise linear networks with ReLU activations, the optimization problem in (4) can be cast as an MILP.

A single ReLU constraint \(y=\max(0,x)\) with pre-activation bounds \(\underline{x}\leq x\leq\bar{x}\) can be equivalently described by the following mixed-integer linear constraints (Tjeng et al. (2017)):

\[y=\max(0,x),\ \underline{x}\leq x\leq\bar{x}\iff\{y\geq 0,\ y\geq x,y\leq x- \underline{x}(1-t),\ y\leq\bar{x}t,\ t\in\{1,0\}\}, \tag{6}\]

where the binary variable \(t\in\{1,0\}\) is an indicator of the activation function being active (\(y=x\)) or inactive (\(y=0\)). Now consider an \(\ell\)-layer feed-forward fully-connected ReLU network with input \(x\) given by the following recursions,

\[x^{(k+1)}=\max(W^{(k)}x^{(k)}+b^{(k)},0)\text{ for }k=0,\cdots,\ell-1;\ f(x^{(0)} )=W^{(\ell)}x^{(\ell)}+b^{(\ell)}, \tag{7}\]

Figure 1: Illustration of problems 1 (distance to invertibility boundary, red) and 2 (distance to pseudo invertibility boundary, blue).
where \(x^{(k)}\in\mathbb{R}^{n_{k}}\) gives the input to the \((k+1)\)-th layer (specifically, we have \(x=x^{(0)}\) and \(n_{0}=m\)), \(W^{(k)}\in\mathbb{R}^{n_{k+1}\times n_{k}},b^{(k)}\in\mathbb{R}^{n_{k+1}}\) are the weight matrices and bias vectors of the affine layers. We denote \(n=\sum_{k=1}^{\ell}n_{k}\) the total number of neurons. Suppose \(l^{(k)}\) and \(u^{(k)}\) are known elementwise lower and upper bounds on the input to the \((k+1)\)-th activation layer, i.e., \(l^{(k)}\leq W^{(k)}x^{(k)}+b^{(k)}\leq u^{(k)}\). Then the neural network equations are equivalent to a set of mixed-integer constraints as follows:

\[x^{(k+1)}=\max(W^{(k)}x^{(k)}+b^{(k)},0)\Leftrightarrow\begin{cases}x^{(k+1)} \geq W^{(k)}x^{(k)}+b^{(k)}\\ x^{(k+1)}\leq W^{(k)}x^{(k)}+b^{(k)}-l^{(k)}\odot(1_{n_{k+1}}-t^{(k)})\\ x^{(k+1)}\geq 0,\quad x^{(k+1)}\leq u^{(k)}\odot t^{(k)},\end{cases} \tag{8}\]

where \(t^{(k)}\in\{1,0\}^{n_{k+1}}\) is a vector of binary variables for the \((k+1)\)-th activation layer and \(1_{n_{k+1}}\) denotes vector of all \(1\)'s in \(\mathbb{R}^{n_{k+1}}\). We note that the element-wise pre-activation bounds \(\{l^{(k)},u^{(k)}\}\) can be precomputed by, for example, interval bound propagation or linear programming, assuming known bounds on the input of the neural network (Weng et al. (2018); Zhang et al. (2018); Hein and Andriushchenko (2017); Wang et al. (2018); Wong and Kolter (2018)). Since the state-of-the-art solvers for mixed-integer programming are based on branch \(\&\) bound algorithms (Land and Doig (1960); Beasley (1996)), tight pre-activation bounds will allow the algorithm to prune branches more efficiently and reduce the total running time.

Local invertibility certificates via mixed-integer programmingHaving represented the neural network equations by mixed-integer constraints, it remains to encode the objective function \(\|x_{1}-x_{2}\|\) of (4) as well as the set \(\mathcal{B}\). We assume that \(\mathcal{B}\) is an \(L_{\infty}\) ball around a given point \(x_{c}\), i.e., \(\mathcal{B}=\mathcal{B}_{\infty}(x_{c},r)\). Furthermore, for the sake of space, we only consider \(L_{\infty}\) norms for the objective function. Specifically, consider the equality \(w=\|x_{1}-x_{2}\|_{\infty}\). This equality can be encoded as mixed-integer linear constraints by introducing \(2n_{0}\) mutually exclusive indicator variables, which leads to the following MILP:

\[p^{\star}\leftarrow\max w\text{ subject to }\|x_{1}-x_{c}\|_{ \infty}\leq r,\ \ \|x_{2}-x_{c}\|_{\infty}\leq r\] \[\text{(I)}: \begin{cases}(x_{1}-x_{2})\leq w1_{n_{0}}\leq(x_{1}-x_{2})+4r(1_ {n_{0}}-f)\\ -(x_{1}-x_{2})\leq w1_{n_{0}}\leq-(x_{1}-x_{2})+4r(1_{n_{0}}-f^{\prime})\\ f+f^{\prime}\leq 1_{n_{0}},1_{n_{0}}^{\top}(f+f^{\prime})=1,f,f^{\prime}\in\{0,1\}^ {n_{0}}\end{cases}\] \[\text{(II)}: W^{(\ell)}x_{1}^{(\ell)}=W^{(\ell)}x_{2}^{(\ell)}\] \[\text{for }k=0,\cdots,\ell-1:\] \[\text{(III)}: \begin{cases}x_{1}^{(k+1)}\geq W^{(k)}x_{1}^{(k)}+b^{(k)},x_{2}^{ (k+1)}\geq W^{(k)}x_{2}^{(k)}+b^{(k)}\\ x_{1}^{(k+1)}\leq W^{(k)}x_{1}^{(k)}+b^{(k)}-l^{(k)}\odot(1-t^{(k)}),x_{2}^{(k +1)}\leq W^{(k)}x_{2}^{(k)}+b^{(k)}-l^{(k)}\odot(1-t^{(k)})\\ x_{1}^{(k+1)}\geq 0,x_{2}^{(k+1)}\geq 0,x_{1}^{(k+1)}\leq u^{(k)}\odot t^{(k)},x_{2}^{ (k+1)}\leq u^{(k)}\odot t^{(k)};t^{(k)},s^{(k)}\in\{0,1\}^{n_{k}+1},\end{cases}\]

where the set of constraints in (I) model the objective function \(\|x_{1}-x_{2}\|_{\infty}\), and the set of constraints \(\text{(III)}\) encode the network \(x_{1}^{(k+1)}=\max(W^{(k)}x_{1}^{(k)}+b^{(k)},0)\) and \(x_{2}^{(k+1)}=\max(W^{(k)}x_{2}^{(k)}+b^{(k)},0)\). The constraint \(\text{(II)}\) enforces that \(f(x_{1})=f(x_{2})\). This optimization problem (4) has \(2(n_{0}+n)\) integer variables.
**Remark 3**: _If we instead use the \(\ell_{2}\) norm both for the objective function and the ball \(\mathcal{B}_{2}(x_{c},r)\), we will arrive at a mixed-integer quadratic program (MIQP). However, (9) remains an MILP if we change them to \(\ell_{1}\) norms._

_Largest region of invertibility_For a fixed radius \(r\geq 0\), the optimization problem (9) either verifies whether \(f\) is invertible on \(\mathcal{B}_{\infty}(x_{c},r)\) or it finds counterexamples \(x_{1}\neq x_{2}\) such that \(f(x_{1})=f(x_{2})\). Thus, we can find the maximal \(r\) by performing a bisection search on \(r\) (Problem 1).

To close this section, we consider the problem of invertibility certification in transformations between two functions (and in particular two neural networks).

**Problem 3** (Transformation Invertibility): _Given two functions \(f_{1},f_{2}\colon\mathbb{R}^{m}\to\mathbb{R}^{m}\) and a particular point \(x_{c}\in\mathbb{R}^{m}\) in the input space, we would like to find the largest ball \(\mathcal{B}_{q}(x_{c},r)\) over which the output of \(f_{2}\) is a function of the output of \(f_{1}\) (and vice versa)._

**Theorem 4**: _Let \(f_{1}\colon\mathbb{R}^{m}\to\mathbb{R}^{n}\), \(f_{2}\colon\mathbb{R}^{m}\to\mathbb{R}^{n}\) be two continuous functions and \(\mathcal{B}\subset\mathbb{R}^{m}\) be a compact set. Consider the following optimization problem,_

\[p_{12}^{\star}\leftarrow\max\quad\|f_{2}(x_{1})-f_{2}(x_{2})\|\quad\text{ subject to }x_{1},x_{2}\in\mathcal{B},\quad f_{1}(x_{1})=f_{1}(x_{2}). \tag{10}\]

_Then the output of \(f_{2}\) is a function of the output of \(f_{1}\) on \(\mathcal{B}\) if and only if \(p_{12}^{\star}=0\)._

Similar to Problem 1, we can pose Problem 3 as a mixed-integer program. Furthermore, we can also define \(p_{21}^{\star}\), whose zero value determines whether output of \(f_{1}\) is a function of output of \(f_{2}\) over \(\mathcal{B}\). It is straightforward to see that \(p_{12}^{\star}=p_{21}^{\star}=0\) if and only if output of \(f_{2}\) is an invertible function of output of \(f_{1}\).

## 4 Numerical Experiments

We now present experiments with \(\mathrm{ReLU}\) multi-layer perceptrons (MLPs) in both (a) regression problems, and also in (b) transformations between two \(\mathrm{ReLU}\) networks.

1D ExampleWe use a 1-10-10-1 randomly generated fully-connected neural network \(f(x)\) with \(\mathrm{ReLU}\) activations. We find the largest interval around the points \(x=-1.8;-1;-0.3\) on which \(f\) is invertible (Problem 1); we also find the largest interval around the point \(x=-1\) for which no other interior points map to \(f(-1)\) (Problem 2). The results are plotted in Figure 2, where intervals in red and blue respectively represent the optimal solutions for the two problems. The largest certified radii are 0.157, 0.322 and 0.214 for Problem 1 and 0.553 for Problem 2.
2D Example: a discrete-time integrator.The Brusselator (Tyson (1973)) is a system of two ODEs for the two variables \((x,y)\), depending on the parameters \((a,b)\); it describes oscillatory dynamics in a theoretical chemical reaction scheme. We use its forward-Euler discretization with step \(\tau\),

\[x_{n+1}=x_{n}+\tau(a+x_{n}^{2}y_{n}-(b+1)x_{n}),\;y_{n+1}=y_{n}+\tau(bx_{n}-x_{n }^{2}y_{n}). \tag{11}\]

Rearranging and eliminating \(y_{n+1}\) in (11) we obtain:

\[\tau(1-\tau)x_{n}^{3}+\tau(\tau a-x_{n+1}-y_{n+1})x_{n}^{2}+(\tau b+\tau-1)x_{n }+(x_{n+1}-\tau a)=0. \tag{12}\]

Equation (12) is a cubic for \(x_{n}\) given \((x_{n+1},y_{n+1})\) when \(\tau\neq 1\). By varying the parameters \(a\), \(b\) and \(\tau\), we see the past states \((x_{n},y_{n})^{T}\) of point \((x_{n+1},y_{n+1})^{T}\) (also called "inverses" or "preimages") may be multi-valued, so that this discrete-time system is, in general, noninvertible. We fix \(a=1\) and consider how inverses will be changing (a) with \(b\) for fixed \(\tau=0.15\); and (b) with \(\tau\), for fixed \(b=2\).

We are interested in training a neural network that learns this time-\(\tau\) mapping; for a fixed set of parameter values, this is a network from 3D to 2D: \((x_{n+1},y_{n+1})^{T}\approx\mathcal{N}(x_{n},y_{n};p)^{T}\), where \(p\in\mathbb{R}\) is the parameter. The network dynamics will be parameter-dependent if we set \(p\equiv b\), or timestep-dependent if \(p\equiv\tau\). The first layer of such an MLP reads

\[W^{(0)}\begin{bmatrix}x_{n}\\ y_{n}\\ p\end{bmatrix}+b^{(0)}=(W^{(0)}(e_{1}+e_{2}))\begin{bmatrix}x_{n}\\ y_{n}\end{bmatrix}+(pW^{(0)}e_{3}+b^{(0)}), \tag{13}\]

where \(e_{1,2,3}\in\mathbb{R}^{3}\) are indicator vectors. Here we trained two separate MLPs, tone with \(b\) and one with \(\tau\) dependence. For fixed \(p\) (either \(b\) or \(\tau\)) each of these two networks \(\mathcal{N}\) can be thought of as a MLP mapping from \(\mathbb{R}^{2}\) to \(\mathbb{R}^{2}\), by slightly modifying the weights and biases in the first linear layer.

Parameter-dependent InversesIt is useful to start with a brief discussion of the dynamics and noninvertibilities in the ground-truth system (see Figure 3). Consider a state located on the invariant circle (IC, shown in orange), for we therefore know there exists at least one preimage _also on this IC_. In Figure 3 we indeed see that every point on the IC has three preimages: one still on the IC, and two extra inverses (in green and purple) after one iteration, all three loops map to the orange one,

Figure 2: Solutions to Problem 1 (left, red) and Problem 2 (right, blue) for the MLP corresponding to a randomly-generated ReLU network (see text).
and then remain forward invariant. The phase space, upon iteration, _folds_ along the two branches of the \(J_{0}\) curve (sets of red points). For lower values of \(b\), these three closed loops _do not intersect each other_. As \(b\) increases the (orange) attractor will become tangent to, and subsequently intersect \(J_{0}\), leading to an interaction with the other (green) preimage branch. At this point the dynamics predicted by the network become unphysical (beyond just inaccurate).

After convergence of training, we employ our algorithm to obtain noninvertibility certificates for the resulting MLP, and plot results for \(b=2.1\) in Figure 4. In Figure 4, we arbitrarily select one representative point, marked by triangle (\(\triangle\)), on the attractor (the orange invariant circle); we know there exists one inverse _also_ located on the attractor, see the nearby cross (\(+\)); we call this the _primal_ inverse. Our algorithm will produce two regions for this point, one for each of our problems (squares of constant \(L_{\infty}\) distance in 2D). As a sanity check, we also compute the \(J_{0}\) sets (the red point), as well as a few additional inverses, beyond the primal ones with the help of a numerical root solver and automatic differentiation (Baydin et al. (2017)). Clearly, the smaller square neighborhood just hits the \(J_{0}\) curve, while the larger one extends to the closest non-primal inverse of the attractor.

Timestep-dependent InversesIn the right two subfigures of Figure 4, we explore the effect of varying the time horizon \(\tau\). We compare a single Euler step of the ground truth ODE to the MLP approximating the same time \(\tau\) map, and find that, for both of them, smaller time horizons lead to larger regions of invertibility.

Figure 3: Attractors (and their multiple inverses) for several parameter values of the discrete Brusselator neural network for \(\tau=0.15\). Notice the relative positions of the \(J_{0}\) curves (red), the “main” preimage locus (yellow), and the “extra” preimages (green, purple). When the attractor starts interacting with the \(J_{0}\) curve and, therefore, with these extra preimages, the dynamic behavior degenerates quantitatively and qualitatively (see also Rico-Martinez et al. (1993)).
Network Transformation Example: Learning the Van der Pol EquationHere, to test our algorithm on the problem of transformations between networks 3, we trained two networks on the same regression task. Our data comes from the 2D Van der Pol equation \(dx_{1}/dt=x_{2},dx_{2}/dt=\mu(1-x_{1}^{2})x_{2}-x_{1}\), where the input and output are the initial and final states of 1000 short solution trajectories of duration 0.2 for \(\mu=1\), when a stable limit cycle exists. The initial states are uniformly sampled in the region \([-3,3]\times[-3,3]\). The neural network A used to learn the time-\(\tau=0.2\) map is a 2-32-32-2 MLP, while the neural network B is a retrained sparse version of A, where half of the weight entries are pruned (set to zero) based on Zhu and Gupta (2018). To visualize the performances of the two networks, two trajectories, generated by respectively iterating each network function for a fixed number of times starting from a common given initial state have been plotted in the left subplot of Figure 5. The ODE solution trajectory starting at the same initial state with same overall time duration is also shown. We see that both network functions A and B exhibit long term oscillations; the shapes of both attractors appear to only have small visual differences from the true ODE solution (the red curve).

These two network functions were then used to illustrate the algorithm for Problem 3. Here we chose a center point \(x_{c}=(0,0)^{T}\), computed and plotted the mappable regions (the regions over which there is a one-to-one mapping between the output of one network and the output of the other, i.e. where one network can be _calibrated_ to the other). This was done for two subcases (see the right subfigure of Figure 5): (a) where the output of network \(B\) is a function of the output of network \(A\) (the square with white bounds centered at the red point, radius 3.0820), and vice versa, where the output of network \(A\) is a function of the output of the network \(B\) (the square with black bounds centered at the red point, radius 3.6484). This also gives us the "common" region (the interior of the white square) where both networks can be calibrated _to each other_. For validation we also computed the Jacobian values of network \(A\) and network \(B\) on every grid point of the input domain, and shown that the white square touches the \(J_{0}\) curve of network \(A\), while the black square touches the \(J_{0}\) curve of network \(B\). Inside the black square the Jacobian of network \(B\) remains positive, so that network \(B\) is invertible (i.e. there exists a mapping from \(f_{B}(x)\) to \(x\), or equivalently, \(f_{B}^{-1}(x)\));

Figure 4: Left: illustration of our solution to Problems 1 and 2 for the Brusselator network with \((a,b)=(1,2.1)\). For a particular reference point on the attractor, we show the neighborhoods found by our algorithms. They clearly locate the closest point on the \(J_{0}\) curve / the closest “extra preimage” of the point of interest. Last two: plots of \(J_{0}\) curves at different \(\tau\) with \((a,b)=(1,2)\), for both the Euler integrator and our Brusselator ReLU network. Small timesteps lead to progressively more remote \(J_{0}\) curves. Notice also the piecewise linear nature of the \(J_{0}\) curve for the ReLU network; its accurate computation constitutes an interesting challenge by itself.
therefore we can find the mapping from \(f_{B}(x)\) to \(f_{A}(x)\) by composing the mapping from \(f_{B}(x)\) to \(x\) with the mapping from \(x\) to \(f_{A}(x)\) (the function \(f_{A}(x)\) itself). The size of the white square can be similarly rationalized, validating our computation.

As a sanity check, we consructed eight more pruned networks; two of them have \(50\%\) sparsity (networks \(B_{5}\) and \(B_{6}\)), three have \(40\%\) sparsity (networks \(B_{1},B_{2}\) and \(B_{3}\)) and the others have \(60\%\) sparsity (networks \(B_{7},B_{8}\) and \(B_{9}\)). Above we discussed network \(B_{4}\) For each pruned network, we computed the radii of the regions of interest (aka \(r_{AB}\) and \(r_{BA}\)). The results are listed in Table 1. All pruned networks \(\{B_{i}\}\) share the same radii \(r_{AB}\), consistent with the invertibility of \(A\) itself. Since \(r_{A}=3.0820\), \(A\) is invertible in the ball we computed, and the existence of the mapping \(f_{A}(x)\mapsto f_{B}(x)\) follows by composition of \(f_{A}(x)\mapsto x\) and \(x\mapsto f_{B}(x)\). Based on these few computational experiments one might very tentatively surmise a trend: the higher the pruning (e.g. \(60\%\)) the larger the invertibility guarantee for the pruned network. In our work the input and output

\begin{table}
\begin{tabular}{c|c c c|c c c|c c c} \hline \hline Sparsity & \multicolumn{4}{c|}{40 \%} & \multicolumn{4}{c|}{50 \%} & \multicolumn{4}{c}{60 \%} \\ \hline Network \(B\) & \(B_{1}\) & \(B_{2}\) & \(B_{3}\) & \(B_{4}\) & \(B_{5}\) & \(B_{6}\) & \(B_{7}\) & \(B_{8}\) & \(B_{9}\) \\ \hline \(r_{AB}\) & 3.0820 & 3.0820 & 3.0820 & 3.0820 & 3.0820 & 3.0820 & 3.0820 & 3.0820 & 3.0820 & 3.0820 \\ \hline \(r_{BA}\) & 3.4609 & 3.1055 & 3.8555 & 3.6484 & 2.6523 & 3.8203 & 3.6328 & 3.9727 & 4.5547 \\ \hline \hline \end{tabular}
\end{table}
Table 1: The radii of the mappable regions between the original network \(A\) and its pruned versions \(B\). \(r_{AB}\) relates to the region within which \(f_{B}(x)\) is a function of \(f_{A}(x)\).

Figure 5: Left: Trajectories of the ODE solution for the Van der Pol system (red), and their discrete-time neural network approximations (blue and green). All three trajectories begin at the same initial state. While the ODE solution curve is smooth due to its continuous-time nature, the others are just straight line segments connecting consecutive states (discrete-time dynamics). However, it is clear that all three systems have visually nearby long-time dynamic attractors, corroborating the good performance of the network and its pruned version. Right: visualization of MILP computation results, along with signs of Jacobian values of networks on the grid points of the input domain. Here, the center of the region is shown in red, while the white and black boundaries quantify the mappable region between outputs of network A and network B.
dimensions are the same (e.g. \(m=n\) in Problem 3). However, this condition is not necessary, and our algorithm can be conceptually extended to classification problems, where in general \(m\gg n\).

## 5 Conclusions

In this paper, we revisited noninvertibility issues that arise in discrete-time dynamical systems integrators) as well as in neural networks that perform approximations of the same (time-series related) task. We argued that such noninvertibility may have dramatic pathological consequences, going beyond just inaccuracies, in the dynamics predicted by the networks. We also extended the analysis to transformations between different neural networks. We formulated three problems that provide a quantifiable assessment of "local" invertibility for any given, arbitrarily selected input. Specifically, for functions like MLPs with ReLU activations, these problems were formulated as mixed-integer programs. We then performed experiments on regression tasks. An extension of our algorithm to ResNets. can be found in the Appendix.

Future directions include developing structure-exploiting methods to globally solve these MIPs more efficiently, and for larger networks. On the other hand, given that convolution and average pooling are linear operations, while max pooling is piecewise linear, it is natural to adapt our algorithms to convolutional neural networks like AlexNet (Krizhevsky et al. (2017)) or VGG (Simonyan and Zisserman (2015)). The successful application of our algorithm to ResNet architectures (He et al. (2016)) holds promise for applicability also to recursive architectures (Lu et al. (2018); E (2017)), such as fractal networks (Larsson et al. (2017)), poly-inception networks (Zhang et al. (2016)), and RevNet (Gomez et al. (2017)). We are working on making the algorithm practical for continuous differentiable activations like tanh or Swish (Ramachandran et al. (2017)), and for other piecewise activations like gaussian error linear units (GELUs, Hendrycks and Gimpel (2016)). We are particularly interested in the case when the input and output domains are of different dimension (e.g., classifiers).
