Provably Bounding Neural Network Preimages
Most work on the formal verification of neural networks has focused on
bounding the set of outputs that correspond to a given set of inputs (for
example, bounded perturbations of a nominal input). However, many use cases of
neural network verification require solving the inverse problem, or
over-approximating the set of inputs that lead to certain outputs. We present
the INVPROP algorithm for verifying properties over the preimage of a linearly
constrained output set, which can be combined with branch-and-bound to increase
precision. Contrary to other approaches, our efficient algorithm is
GPU-accelerated and does not require a linear programming solver. We
demonstrate our algorithm for identifying safe control regions for a dynamical
system via backward reachability analysis, verifying adversarial robustness,
and detecting out-of-distribution inputs to a neural network. Our results show
that in certain settings, we find over-approximations over 2500x tighter than
prior work while being 2.5x faster. By strengthening robustness verification
with output constraints, we consistently verify more properties than the
previous state-of-the-art on multiple benchmarks, including a large model with
167k neurons in VNN-COMP 2023. Our algorithm has been incorporated into the
$\alpha,\!\beta$-CROWN verifier, available at https://abcrown.org.
Suhas Kotha, Christopher Brix, Zico Kolter, Krishnamurthy Dvijotham, Huan Zhang
# Provably Bounding Neural Network Preimages

###### Abstract

Most work on the formal verification of neural networks has focused on bounding forward images of neural networks, i.e., the set of outputs of a neural network that correspond to a given set of inputs (for example, bounded perturbations of a nominal input). However, many use cases of neural network verification require solving the inverse problem, i.e, over-approximating the set of inputs that lead to certain outputs. In this work, we present the first efficient bound propagation algorithm, INVPROP, for verifying properties over the preimage of a linearly constrained output set of a neural network, which can be combined with branch-and-bound to achieve completeness. Our efficient algorithm allows multiple passes of intermediate bound refinements, which are crucial for tight inverse verification because the bounds of an intermediate layer depend on relaxations both before and after this layer. We demonstrate our algorithm on applications related to quantifying safe control regions for a dynamical system and detecting out-of-distribution inputs to a neural network. Our results show that in certain settings, we can find over-approximations that are over \(2500\times\) tighter than prior work while being \(2.5\times\) faster on the same hardware.

## 1 Introduction

Applications of neural networks to safety-critical settings often require reasoning about the set of inputs that will produce a particular set of outputs. For example, for a physical system controlled by a policy parameterized as a neural network, it is of interest to understand the states that would cause the policy to choose control actions that are outside acceptable limits (for example torques that are beyond the safe operating range of a motor), or that would result in a future state that violates safety constraints. Another example is quantifying the out-of-distribution (OOD) detection behavior of a neural network: here, it is of interest to bound the set of inputs that leads to a classifier making a confident prediction, i.e, where the highest scoring logit exceeds the second highest scoring logit by a large margin. Given this set, one can assess whether there are inputs far from the training data that are assigned high confidence predictions, which would be an undesirable property.

Formal verification of neural networks seeks to provide provable guarantees demonstrating that networks satisfy formal specifications on their inputs and outputs. Practical algorithms often rephrase the above as finding provably correct bounds on the inputs or outputs of a neural network subject to specific constraints. Most work to date has focused on developing algorithms that can bound the outputs of a neural network given constraints on the inputs, which is challenging due to high
dimensionality and non-convexity. For example, to analyze the robustness of a neural network to perturbations of a given input, these algorithms compute bounds on the outputs of the network given that the inputs are bounded perturbations of a nominal input (Wong and Kolter, 2018; Dvijotham et al., 2018; Zhang et al., 2018; Raghunathan et al., 2018; Gehr et al., 2018).

In this work, we address the inverse problem, motivated by the use-cases outlined above. Our goal is to find bounds on the pre-image of a neural network: given a set of outputs \(\mathcal{S}_{\text{out}}\) (described by linear constraints on the output of a neural network), we seek to find a set that provably contains all inputs that lead to such outputs. The verification problem is already challenging due to non-convexity and high dimensionality. This new problem is increasingly difficult since neural networks are not invertible.

Specifically, representative efficient verifiers (such as state of the art bound-propagation-based methods (Zhang et al., 2022)) can only compute bounds in the forward direction, and the tightness of the linear programming relaxations they rely on to perform formal verification critically depends on having tight bounds on intermediate activations. In the setting of this paper, however, the bounds available on the intermediate activations derived from the input constraints are weak, since the only constraints on the input are that it should be from the valid input domain (for example, normalized pixels from an image lie between \(0\) and \(1\)). In fact, in some applications like control, the input domain may be unbounded. Given this, the bounds on intermediate activations derived from the input constraints may be vacuous or very loose, impacting the tightness of the intermediate activation relaxations. We efficiently solve this problem by significantly generalizing the existing bound propagation-based verification framework. Our contributions are as follows:

\(\bullet\) We formulate the _inverse verification problem_ for neural networks, i.e, the problem of over-approximating the set of inputs that leads to a given set of outputs, with a mixed integer linear programming formulation. This motivates our development of an efficient bound propagation method.

\(\bullet\) We develop an effective bound propagation framework, Inverse Propagation for Neural Network Verification (INVPROP), which computes provable over-approximations of neural network preimages. We study the convex outer approximation of this problem which leads to a novel inverse bound propagation method for verification. We also unify INVPROP and forward verification bound propagation into a more general verification framework, allowing us to connect our method to standard tools and benefit from prior progress in the verification community, such as efficient GPU acceleration.

\(\bullet\) We demonstrate that tight inverse verification requires multiple refinement passes of intermediate bounds. This is unique to the inverse verification problem because bounds for an intermediate layer depend on not only intermediate bounds before this layer, but also constraints imposed on downstream layers or the output of the network. In other words, bounds cannot be simply propagated in a single forward pass first tightening all neurons in a given layer and then tightening the next layer based on the previous layer bounds. Instead, inverse verification requires an iterative process where, in each pass, all bounds are tightened based on previous bounds (including output constraints), and the iterative process continues until bounds converge.

\(\bullet\) We improve the state-of-the-art on a control benchmark Rober et al. (2022a,b) by providing \(2500\times\) tighter bounds with \(2.5\times\) faster computation. We also demonstrate applicability in OOD detection.

## 2 Setup

### Notation

Throughout this paper, we will use \([L]\) for \(L\in\mathbb{N}\) to refer to the set \(\{1,2,\ldots,L\}\), \(\mathbf{W}_{:,j}^{(i)}\) to refer to column \(j\) of the matrix \(\mathbf{W}^{(i)}\), \([\cdot]_{+}\) to refer to \(\max(0,\cdot)\), and \([\cdot]_{-}\) to refer to \(-\min(0,\cdot)\). We will also boldface symbols for vectors and matrices (such as \(\mathbf{x}^{(i)}\) and \(\mathbf{W}^{(i)}\)) and use regular symbols for scalars (such as \(x_{j}^{(i)}\)). We use \(\mathbf{x}\odot\mathbf{y}\) to denote element-wise multiplication of vectors \(\mathbf{x},\mathbf{y}\).

We define an \(L\) layer ReLU neural network by its weight matrices \(\mathbf{W}^{(i)}\) and bias vectors \(\mathbf{b}^{(i)}\) for \(i\in[L]\). The output of the neural network for the input \(\hat{\mathbf{x}}^{(0)}\) from a bounded input domain \(\mathcal{X}\) is computed by alternately applying linear layers \(\mathbf{x}^{(i)}=\mathbf{W}^{(i)}\hat{\mathbf{x}}^{(i-1)}+\mathbf{b}^{(i)}\) and ReLU layers
\(\hat{\mathbf{x}}^{(i)}=\max(\mathbf{0},\mathbf{x}^{(i)})\) until we receive the output \(\mathbf{x}^{(L)}\) (which we refer to as the logits). Note that we treat softmax as a component of the loss function, not the neural network.

The inputs to our network will lie in the Euclidean space \(\mathbb{R}^{\text{in}}\) and the outputs will lie in \(\mathbb{R}^{\text{out}}\). For any subset \(S\) of \(\mathbb{R}^{d}\) (\(d\) dimensional Euclidean space), \(\operatorname{CONV}(S)\) denotes the convex hull of \(S\), i.e, the smallest convex set containing \(S\).

### Problem Statement

Given a neural network \(f:\mathcal{X}\subseteq\mathbb{R}^{\text{in}}\to\mathbb{R}^{\text{out}}\) and an output constraint \(\mathcal{S}_{\text{out}}\subseteq\mathbb{R}^{\text{out}}\), we want to compute \(f^{-1}(\mathcal{S}_{\text{out}})\subseteq\mathcal{X}\). Precisely computing or expressing \(f^{-1}(\mathcal{S}_{\text{out}})\) is an intractable problem in general. Therefore, we strive to compute a tight over-approximation \(\mathcal{S}_{\text{over}}\) such that \(f^{-1}(\mathcal{S}_{\text{out}})\subseteq\mathcal{S}_{\text{over}}\). In particular, we will develop an algorithm to compute \(\mathcal{S}_{\text{over}}=\operatorname{CONV}\left(f^{-1}(\mathcal{S}_{\text {out}})\right)\), the convex hull of the pre-image, via a cutting-plane representation. In this paper, \(\mathcal{S}_{\text{out}}\) will always be defined by a set of linear constraints parameterized by \(\mathbf{H}f\left(\mathbf{x}\right)+\mathbf{d}\leq\mathbf{0}\).

### Applications

We outline two concrete applications of our approach.

Backward Reachability Analysis for Neural Feedback Loops.For control problems involving linear dynamical systems with quadratic cost functions, the optimal feedback policies can be proven to be linear Bertsekas (2000). However, many tasks require non-quadratic cost functions (for example, tasks requiring driving a system towards a desired safe set). In this case, linear feedback policies are no longer sufficient. Neural networks constitute a flexible class of nonlinear feedback policies. However, establishing safety guarantees for neural network policies is a challenging task. In particular, one problem of interest is to find a set of initial states that does not reach a particular set of future states under the neural network policy. This can be helpful in collision avoidance or control with safety constraints.

Under linear time-invariant dynamics, the state of a system is determined by the state transition function

\[\mathbf{x}_{t+1}=f(\mathbf{x}_{t})=\mathbf{A}\mathbf{x}_{t}+\mathbf{B}\pi(\mathbf{x}_{t})+ \mathbf{c}\]

where for \(x_{n},x_{u}\in\mathbb{N}\), \(\mathbf{x}_{t}\in\mathcal{X}\subseteq\mathbb{R}^{x_{n}}\) is the state at discrete timestep \(t\), \(\mathbf{A}\in\mathbb{R}^{x_{n}\times x_{n}}\) and \(\mathbf{B}\in\mathbb{R}^{x_{n}\times x_{n}}\) are known system matrices, \(\mathbf{c}\in\mathbb{R}^{x_{n}}\) is a known exogenous input, and \(\pi:\mathcal{X}\to\mathbb{R}^{n_{u}}\) is the state-feedback control policy (the decision the model takes given the state).

In the controls benchmark we use for our evaluation Hu et al. (2020); Everett et al. (2021, 2022), a robot is moving on the floor of a room with coordinates spanning \(\mathcal{X}=[-5,5]\times[-5,5]\). The position of the robot at time \(t+1\) can be directly computed based on the position at time \(t\), and follows the equation

\[\mathbf{x}_{t+1}=f(\mathbf{x}_{t})=\begin{bmatrix}1&1\\ 0&1\end{bmatrix}\mathbf{x}_{t}+\begin{bmatrix}0.5\\ 1\end{bmatrix}\pi(\mathbf{x}_{t})\]

with \(\pi:\mathbb{R}^{2}\to\mathbb{R}\). There is an obstacle in the room covering the region \([4.5,5.0]\times[-0.25,0.25]\) that the robot needs to avoid. It is of interest to understand whether, starting in a given state, the system will enter the unsafe region in the next time-step. We can represent this obstacle set with the linear constraints

\[\mathcal{S}_{\text{out}}=\left\{\mathbf{x}_{t+1}:\begin{bmatrix}-1&0\\ 0&-1\\ 1&0\\ 0&1\end{bmatrix}\mathbf{x}_{t+1}+\begin{bmatrix}4.5\\ -0.25\\ -5.0\\ -0.25\end{bmatrix}\leq\begin{bmatrix}0\\ 0\\ 0\\ 0\end{bmatrix}\right\}\]

\(f^{-1}(\mathcal{S}_{\text{out}})\) denotes the set of states \(\mathbf{x_{t}}\) such that \(\mathbf{x_{t+1}}\) lies in the unsafe region given the control policy \(\pi\). Overapproximating this set allows us to define the set of states to avoid one timestep in advance. We can compose \(f\) with itself \(t\) times to obtain the set of initial states where \(\pi\) would drive the system to the unsafe set after \(t\) steps.
Computing the set of inputs that leads to confident predictions.Training a classifier that knows when it doesn't know is a challenging problem. A simple and effective approach to doing this is to train a classifier with an additional logit representing (out of distribution) OOD-ness. Labeled data for this additional class is generated via outlier exposure, i.e., adding synthetic training data with features drawn from an OOD distribution (known as outlier exposure (Hendrycks et al., 2018)) and labels corresponding to the additional OOD logit (Chen et al., 2021).

Consider the logits \(y\) produced by a classifier trained in this manner for a binary classification task. The classifier can classify in-distribution data by comparing \(y_{0}\) and \(y_{1}\), the logits corresponding to the two in-distribution labels. The quantity \(\max(y_{0},y_{1})-y_{2}\), known as the logit gap, can be used to identify out-of-distribution data (Fig. 1).

A natural question to ask is: Can it correctly identify datapoints far from the training data as being OOD? In order to do so, one needs to compute the preimage of the set

\[\mathcal{S}_{\text{out}}=\{\mathbf{y}:\max(y_{0},y_{1})\geq y_{2}\}\]

We detail how to model this set with linear constraints in Section D.1. The algorithm INVPROP we develop in this paper enables us to compute an overapproximation of \(f^{-1}(\mathcal{S}_{\text{out}})\), thereby answering the question of whether the classifier learned to correctly identify OOD data.

## 3 Approach

We first provide an algorithm to find convex over-approximations of feed-forward ReLU networks. We then extend our method to non-convex over-approximations of general computational graphs.

### Convex over-approximation

Suppose we want to find a convex over-approximation of \(f^{-1}(\mathcal{S}_{\text{out}})\). The tightest such set is its convex hull, which is the intersection of all half-spaces that contain \(f^{-1}(\mathcal{S}_{\text{out}})\)(Boyd et al., 2004). This intersection is equivalent to \(\bigcap_{\mathbf{c}\in\mathbb{R}^{n}}\{\mathbf{x}:\mathbf{c}^{\top}\mathbf{x}\geq\min_{f(\mathbf{ x}^{\prime})\in\mathcal{S}_{\text{out}}}\mathbf{c}^{\top}\mathbf{x}^{\prime}\}\), which means that we can build an over-approximation by taking this intersection for finitely many \(\mathbf{c}\). Furthermore, replacing the minimization problem with a lower bound to its value still yields a valid half-space, where a tighter bound yields a tighter approximation.

We focus on convex over-approximations for two reasons: First, checking that the pre-image satisfies a linear constraint \(\mathbf{c}^{\top}\mathbf{x}+d\geq 0\) is equivalent to minimizing the linear function \(\mathbf{x}\mapsto\mathbf{c}^{\top}\mathbf{x}\) over the pre-image, which is in turn equivalent to minimizing the linear function over \(\operatorname{CONV}\left(f^{-1}(\mathcal{S}_{\text{out}})\right)\)(Boyd et al., 2004). Second, the convex hull and its tractable over-approximations are conveniently represented as intersections of linear constraints on the pre-image, each of which can be quickly computed after a single precomputation phase to tighten bounds on intermediate neurons using the INVPROP algorithm outlined in the next section.

Figure 1: To improve our model’s ability to detect data outside the training distribution (of green and red points), we randomly sample points that are far from every training point. In the right contour, we see that the model’s confidence of in-distribution increases as we approach the centers of the distribution, demonstrating our model is more calibrated.
For the above reasons, we will focus on solving the following constrained optimization problem.

\[\min_{\mathbf{x}} \mathbf{c}^{\top}\mathbf{x}\] (1a) s.t. \[\mathbf{x}\in\mathcal{X};\quad f\left(\mathbf{x}\right)\in\mathcal{S}_{ \text{out}} \tag{1b}\]

Note that this differs from the forward verification problem widely studied in the literature, which can be phrased as

\[\min_{\mathbf{x}\in\mathcal{S}_{\text{in}}}\quad\mathbf{c}^{\top}f\left(\mathbf{x}\right)\]

where \(\mathcal{S}_{\text{in}}\) is a set representing constraints on the input, and the goal of the verification is to establish that the output \(f\left(\mathbf{x}\right)\) satisfies a linear constraint of the form \(\mathbf{c}^{\top}f\left(\mathbf{x}\right)\geq d\) for all \(x\in\mathcal{S}_{\text{in}}\).

### The INVPROP Algorithm

As a brief overview of our method, we first construct the Mixed Integer Linear Program (MILP) for solving this optimization problem, which generates the over-approximation \(\mathcal{S}_{\text{MILP}}\) when solved for finitely many \(\mathbf{c}\). Then, we relax the MILP to a Linear Program (LP), which can construct the over-approximation \(\mathcal{S}_{\text{LP}}\). Finally, we relax the LP via its Lagrangian dual, which will be used to construct our over-approximation \(\mathcal{S}_{\text{over}}\). This chain of relaxations is visualized in Figure 2. Though the diagram displays looseness, in the coming sections we provide a comprehensive methodology to arbitrarily reduce error in all three relaxations.

The Mixed Integer Linear Programming (MILP) FormulationFor feed-forward neural networks with ReLU non-linearities, this problem admits a MILP encoding similar to prior work in adversarial robustness Tjeng et al. (2017). The weights and biases at layer \(i\) are denoted by \(\mathbf{W}^{(i)},\mathbf{b}^{(i)}\). Let \(\hat{\mathbf{x}}^{(i)}\) denote the post-activations (following the application of the ReLU nonlinearity) at layer \(i\) and \(\mathbf{x}^{(i)}\) denote the pre-activations. Then, the problem (1) is equivalent to:

\[\min_{\mathbf{x},\hat{\mathbf{x}},\mathbf{x}} \mathbf{c}^{\top}\mathbf{x}\] (2a) s.t. \[\mathbf{x}\in\mathcal{X};\hat{\mathbf{x}}^{(0)}=\mathbf{x};\quad\mathbf{H} \mathbf{x}^{(L)}+\mathbf{d}\leq\mathbf{0} \tag{2b}\] \[\mathbf{x}^{(i)}=\mathbf{W}^{(i)}\hat{\mathbf{x}}^{(i-1)}+\mathbf{b}^{(i )};\quad i\in[L],\] (2c) \[\hat{\mathbf{x}}^{(i)}=\max(0,\mathbf{x}^{(i)});\quad i\in[L-1] \tag{2d}\]

We note that, as opposed to prior work on verifying forward properties of neural networks (Zhang et al., 2022), our approach seeks to bound a linear function of the _inputs_ to the neural network, subject to constraints on the _output_ of the neural network. This makes it particularly challenging to directly apply bound propagation algorithms derived in prior work, as these algorithms rely on propagating bounds on the inputs to obtain tight bounds on intermediate activations. In our problem, we need to modify these algorithms to leverage constraints on the outputs of downstream layers to obtain bounds on intermediate or input neurons.

Figure 2: Visualization of relaxations. The inner green region depicts the true \(f^{-1}(\mathcal{S}_{\text{out}})\), the blue relaxation depicts the intersection of finite half-spaces solved via MILP, the red relaxation displays the same via LP, and the purple relaxation displays the same via INVPROP.
The Linear Programming (LP) FormulationUnfortunately, finding an exact solution to this MILP is NP-complete Katz et al. (2017). To sidestep the intractability of exact verification, we can compute lower bounds for this program via its convex relaxation. In order to do so, we begin by finding bounds on the outputs of intermediate layers:

\[l_{j}^{(i)}\leq x_{j}^{(i)}\leq u_{j}^{(i)}\]

Based on these bounds, we can rewrite the last constraint of (2) as follows:

\[0\leq\hat{x}_{j}^{(i)} \leq u_{j}^{(i)}z_{j}^{(i)}\] \[x_{j}^{(i)}\leq\hat{x}_{j}^{(i)} \leq x_{j}^{(i)}-l_{j}^{(i)}(1-z_{j}^{(i)})\]

Relaxing the constraints on \(z\) to \(z_{j}^{(i)}\in[0,1]\), we obtain the following LP relaxation of (2):

\[\min_{\mathbf{x},\hat{\mathbf{x}},\mathbf{x}} \mathbf{c}^{\top}\mathbf{x}\] (3a) s.t. \[\mathbf{x}\in\mathcal{X};\mathbf{l}^{(0)}\leq\mathbf{x}\leq\mathbf{u}^{(0)};\hat {\mathbf{x}}^{(0)}=\mathbf{x} \tag{3b}\] \[\mathbf{H}\mathbf{x}^{(L)}+\mathbf{d}\leq\mathbf{0}\] (3c) \[\mathbf{x}^{(i)}=\mathbf{W}^{(i)}\hat{\mathbf{x}}^{(i-1)}+\mathbf{b}^{( i)};\quad i\in[L],\] (3d) \[\mathbf{0}\leq\hat{\mathbf{x}}^{(i)}\leq\mathbf{u}^{(i)}\odot\mathbf{x}^{(i)} \quad i\in[L]\] (3e) \[\mathbf{x}^{(i)}\leq\hat{\mathbf{x}}_{j}^{(i)}\leq\mathbf{x}^{(i)}-\mathbf{l}^{ (i)}\odot\left(1-\mathbf{z}^{(i)}\right)\quad i\in[L],\] (3f) \[\mathbf{0}\leq\mathbf{z}^{(i)}\leq\mathbf{1}\quad i\in[L] \tag{3g}\]

where the bounds \(\mathbf{l}^{(0)}\leq\mathbf{x}\leq\mathbf{u}^{(0)}\) are either the bounds already implicit in \(\mathcal{X}\), or a refinement of these obtained via previous rounds of bound propagation or input branching (see Algorithm 1).

Most efficient neural network verifiers do not solve an LP formulation of verification directly because LP solvers are often slow for problem instances involving neural networks. Instead, the bound propagation framework (Zhang et al., 2018; Wang et al., 2021; Zhang et al., 2022) is a practical and efficient way to lower bound the LP relaxation of the forward verification problem. However, there are _two major roadblocks_ to applying existing methods here: First, typical methods cannot directly handle the output constraints (Eq. 3c) and the objective involving input layer variables (Eq. 3a); second, the intermediate bounds \(\mathbf{l}\) and \(\mathbf{u}\) must be computed considering the output constraint (Eq. 3c), which is impossible with existing bound propagation methods, as discussed below.

**Example 1**.: Consider the toy neural network in Figure 3 under \(\mathcal{X}=[-2,2]\) and \(\mathcal{S}_{\text{out}}=\{y:1\leq y\leq 1.02\}\). Suppose we wanted to find the tightest \([l_{0}^{(1)},u_{0}^{(1)}]\) (bounds on \(x_{0}^{(1)}\)). If we only enforce the constraints preceding the ReLU (standard bound propagation), we get \([-2,2]\). If we only enforce the constraints following the ReLU, we get \((-\infty,1.02]\). However, when we utilize all of the constraints, we find that the true intermediate bounds are \([0,0.01]\), which is over \(100\) times tighter than the intersection of the two previous methods (all derivations provided in Appendix C). Therefore, new techniques are necessary for optimizing intermediate bounds in this setting.

The Inverse Propagation (INVPROP) FormulationBy changing the LP in Section 3.2 to optimize the quantity \(\hat{x}_{j}^{(0)}\) or \(x_{j}^{(i)}\) for \(i\in[L-1]\), the bounds for the \(j\)th neuron of layer \(i\) can be iteratively tightened in separate LP calls. However, this highly constrained program is too expensive to be run multiple times for each neuron in a neural network.

Inspired by the success of CROWN-family robustness verifiers Zhang et al. (2018); Xu et al. (2021); Wang et al. (2021); Zhang et al. (2022), we efficiently provide lower bounds to the program by optimizing its Lagrangian dual. This dual is highly structured, allowing us to formulate quantities such as the input cutting plane and intermediate bounds as closed-form expressions of the dual variables. Our main generalization of the existing method is the ability to optimize variables with

Figure 3: A simple neural network from \(\mathbb{R}\) to \(\mathbb{R}\). Every node is the sum of its incoming nodes unless the edge is labeled ReLU.
respect to constraints that are sequentially **after** it in the neural network. We present this contribution in the following theorem.

**Theorem 1** (Bounding cutting plane).: _Given an output set \(\mathcal{S}_{\text{out}}=\{\mathbf{y}:\mathbf{H}\mathbf{y}+\mathbf{d}\leq\mathbf{0}\}\) and vector \(\mathbf{c}\), \(g_{\mathbf{c}}(\mathbf{\alpha},\mathbf{\gamma})\) is a lower bound to the linear program in (3) for \(\mathbf{0}\leq\mathbf{\alpha}\leq\mathbf{1}\), \(\mathbf{\gamma}\geq\mathbf{0}\), and \(g_{\mathbf{c}}\) defined via_

\[g_{\mathbf{c}}(\mathbf{\alpha},\mathbf{\gamma}) =\mathsf{ReLU}\left(\mathbf{c}-\mathbf{\nu}^{(1)\top}\mathbf{W}^{(1)} \right)\mathbf{I}^{(0)}\] \[-\mathsf{ReLU}\left(-\mathbf{c}+\mathbf{\nu}^{(1)\top}\mathbf{W}^{(1)} \right)\mathbf{u}^{(0)}\] \[-\sum_{i=1}^{L}\mathbf{\nu}^{(i)\top}\mathbf{b}^{(i)}+\sum_{i=1}^{L-1 }\sum_{j\in\mathcal{I}^{\pm(i)}}\left[\frac{u_{j}^{(i)}l_{j}^{(i)}[\hat{\nu}_{ j}^{(i)}]_{+}}{u_{i}^{(j)}-l_{i}^{(j)}}\right]\]

_where every term can be directly recursively computed via_

\[\mathcal{I}^{-(i)} =\{j:u_{j}^{(i)}\leq 0\}\text{ (Provably Off ReLUs)}\] \[\mathcal{I}^{+(i)} =\{j:l_{j}^{(i)}\geq 0\}\text{ (Provably On ReLUs)}\] \[\mathcal{I}^{\pm(i)} =\{j:l_{j}^{(i)}<0<u_{j}^{(i)}\}\text{ (Ambiguous ReLUs)}\] \[\mathbf{\nu}^{(L)} =-\mathbf{\gamma}\] \[\hat{\mathbf{\nu}}^{(i)} =\left(\mathbf{W}^{(i+1)\top}\right)\mathbf{\nu}^{(i+1)}\quad\forall i \in[L-1]\] \[\nu_{j}^{(i)} =\left\{\begin{array}{ll}\hat{\nu}_{j}^{(i)},&j\in\mathcal{I}^{ +(i)}\\ 0,&j\in\mathcal{I}^{-(i)}\\ \frac{u_{j}^{(i)}}{u_{j}^{(j)}-l_{i}^{(j)}}[\hat{\nu}_{j}^{(i)}]_{+}-\alpha_{ j}^{(i)}[\hat{\nu}_{j}^{(i)}]_{-},&j\in\mathcal{I}^{\pm(i)}\end{array}\right.\]

Proof.: Full proof is presented in Appendix B. 

In Appendix A, we show how to also compute bounds on the values of intermediate neurons in the network using a similar approach. Since the intermediate bounds might depend on bounds on neurons in subsequent layers (due to the output constraint), we cannot simply optimize bounds in a single forward pass layer by layer, unlike prior work (Zhang et al., 2022). Instead, we must iteratively tighten these bounds, where in each pass all bounds are tightened given the tightest bounds on all neurons computed thus far. This iterative approach is able to tighten the bounds by several orders of magnitude, as shown in Figure 5. After performing this procedure once, the intermediate bounds can be used to tightly lower bound \(\mathbf{c}^{\top}\mathbf{x}\) for any \(\mathbf{c}\) via Theorem 1. Therefore, this computation can be shared across all the constraints \(\mathbf{c}\) we use to describe \(\mathcal{S}_{\text{over}}\). Our algorithm can be expressed in terms of forward/backward passes through layers of the neural network and implemented via standard deep learning modules in libraries like PyTorch (Paszke et al., 2019). Since all of the operations are auto-differentiable, we can tighten our lower bound using standard gradient ascent (projected by the dual variable constraints).

Connection to forward verificationOur bound in Theorem 1 introduces the dual variable \(\mathbf{\gamma}\), which enforces the output constraint during optimization. In fact, we can use this variable to get a better conceptual interpretation of our result. Consider the optimization problem in (1). Dualizing the output constraint and taking the dual wrt the constraint \(\mathbf{H}f\left(\mathbf{x}\right)+\mathbf{d}\leq\mathbf{0}\) yields the following lower bound:

\[\max_{\mathbf{\gamma}}\min_{\mathbf{x}} \mathbf{c}^{\top}\mathbf{x}+\mathbf{\gamma}^{\top}\left(\mathbf{H}f\left(\mathbf{ x}\right)+\mathbf{d}\right)\] s.t. \[\mathbf{x}\in\mathcal{X};\quad\mathbf{\gamma}\geq 0\]

The objective for the inner minimization can be represented as minimizing a linear function of the output of a residual neural network with a skip connection from the input to the output layer, subject to constraints on the input. This is precisely the standard forward verification problem.1 This unifies the inverse verification problem with the forward verification problem, connecting our method to standard tools Brix et al. (2023); Xu et al. (2021); Ferrari et al. (2022). For example, as long as \(f\) is a
general computation graph with ReLU activations, Auto-LiRPA Xu et al. (2020) can be used to lower bound the objective function with a closed form expression parameterized by relaxation variables \(\mathbf{\alpha}\). With this, we can directly optimize our lower bound with respect to \(\mathbf{\alpha},\mathbf{\gamma}\).

When \(f\) is a feedforward ReLU network, the lower bound described in this section is precisely the same as Theorem 1 (since both are solving the dual of the same linear program). This shows that the introduction of \(\mathbf{\gamma}\) constitutes our generalization to the bound propagation framework.

### Branch and Bound

Our current formulation still faces two sources of looseness: the gap in \(\mathcal{S}_{\text{LP}}\supseteq\mathcal{S}_{\text{MILP}}\) and the gap in \(\mathcal{S}_{\text{MILP}}\supseteq f^{-1}(\mathcal{S}_{\text{out}})\). To overcome both of these issues, we can make use of branching. While there are several possibilities here, we focus on input branching, which gave the biggest empirical gains in our experiments described in Section 4. More concretely, we divide the input domain \(\mathcal{X}=[\mathbf{l}^{(0)},\mathbf{u}^{(0)}]\) into several regions by choosing a coordinate \(i\) and computing

\[\mathcal{X}_{a}=\mathcal{X}\cap\{\mathbf{x}:\mathbf{x}_{i}\geq s_{i}\},\mathcal{X}_{b} =\mathcal{X}\cap\{\mathbf{x}:\mathbf{x}_{i}\leq s_{i}\}\]

so that \(\mathcal{X}=\mathcal{X}_{a}\cup\mathcal{X}_{b}\). Doing this recursively, we obtain several regions and can compute an overapproximation of \(f^{-1}(\mathcal{S}_{\text{out}})\) when the input is in each of those regions, and take a union of the resulting sets. The finer the input domain we compute the over-approximation over, the tighter the approximation is.

## 4 Results

We demonstrate the impact of our approach on a controls benchmark Rober et al. (2022a,b) and OOD detection. For the established controls benchmark, we significantly outperform the state-of-the-art.

We measure the tightness of an over-approximation using its Approximation Ratio, defined as \(\frac{\text{vol}(\mathcal{S}_{\text{out}})}{\text{vol}(f^{-1}(\mathcal{S}_{ \text{out}}))}\). Both of these volumes are heuristically estimated with 1 million evenly distributed samples of the input space.

While INVPROP allows the optimization of arbitrary cutting planes, all experiments were performed using 40 planes with slopes of equally distributed angles. All implementation details are described in Section D.
### Backward Reachability Analysis for Neural Feedback Loops

After combining (see Section D.2) the state-feedback control policy \(\pi\) of the benchmark with the state transition function (see Section 2.3), we get a three layer MLP with \(12\), \(7\), and \(2\) neurons which is typical for applications in this setting. We model multiple time steps via composing copies of this function. For example, the \(10\) time step transition function can be represented as a \(21\) layer MLP (after fusing consecutive linear layers in the composition). Moreover, we leverage the fact that intermediate bounds for the \(t\) time step transition can be used for the \(t+1\) time step transition.

We find that INVPROP is significantly tighter and faster than ReBreach-LP, the SOTA available method for this problem Rober et al. (2022a,b). As evident in Figure 4, ReBreach-LP suffers from increasingly weak bounds as \(t\) increases, whereas our approach is able to compute tight bounds for all \(t\). Different to their implementation, we are able to factor in the output constraint to improve the intermediate neuron bounds (see Figure 6). Furthermore, we can iteratively tighten these bounds with respect to each other. Even though we optimize many more quantities, our efficient bound propagation allows us to solve the problem faster overall. Figure 5 visualizes the improvement of the approximation ratio score and intermediate bounds over time. We also support the partitioning over the input space introduced by Rober et al. 2022a and the polytope bounds introduced by Everett et al. 2022. Both our increased tightness and speed are quantified in Table 1.

### OOD Detection

Consider the calibrated OOD detector presented in Figure 1, encoded by a four layer MLP with \(200\), \(200\), \(3\), and \(2\) neurons. We compute the set of inputs which induce a sufficiently high ID confidence (measured by \(\max\{y_{0},y_{1}\}>y_{2}\)), pictured in green in Figure 7. This set is non-convex, making the convex hull a poor over-approximation. With some simple input space branching, we can get a much tighter over-approximation, as shown in the right plot. We compare the performance of our approach with and without branching over the input space with the MILP baseline (see Table 2).

\begin{table}
\begin{tabular}{l c c c} \hline \hline Method & Input & Approx & Time (sec) \\ \cline{2-3} ReBreach-LP & yes & 4021.47 & \(42.86\pm 0.04\) \\ INVPROP & no & 1.46 & \(17.89\pm 0.03\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Comparison of our method vs SOTA in bounding the inputs that will collide with an obstacle in 10 time steps (as shown in Figure 4). We observe that we are significantly tighter and faster than their method.

Figure 4: In both of the above plots, we over-approximate the set of states that would collide with the obstacle (black square) at \([4.5,5.0]\times[-0.25,0.25]\). We carry out this task for up to \(10\) steps back in time. Each green region is an approximation of the preimages sets via \(1\) million samples. The error of ReBreach-LP exponentially blows up further back in time while our bounding polytopes of 40 edges continue to stay at near perfect tightness.
## 5 Related Work

Formally Verified Neural Networks.There has been a large body of work on the formal verification of neural networks, tackling the problem from the angles of Interval Bound Propagation Gowal et al. (2018); Gehr et al. (2018); Convex Relaxations Wong and Kolter (2018); Shiqi et al. (2018); Salman et al. (2019); Dvijotham et al. (2018); Abstract Interpretation Singh et al. (2018); Linear Programming Ehlers (2017); Semidefinite Programming Raghunathan et al. (2018); Satisfiability Modulo Theories Katz et al. (2017), and Mixed Integer Linear Programming Tjeng et al. (2017). However, most of this work is in the different setting of forward verification. Our work strictly generalizes the bound propagation framework presented in Xu et al. 2021 since setting \(\mathbf{\gamma}=\mathbf{0}\) in Theorem 1 recovers its results.

Formal Verification for Neural Feedback Loops.Our work was motivated by the growing body of work on backward reachability analysis and over-approximating states that result in a target set Everett et al. (2021). The original method solely utilized input constraints for deriving intermediate bounds. The later development of Everett et al. 2022 improved upon this by optimizing \(\mathbf{t}^{(0)}\) and \(\mathbf{u}^{(0)}\)

Figure 5: Our approach enables us to iteratively tighten the bounds of all layers, with each iteration allowing for a smaller approximation ratio with respect to the true preimage. Green: Sum of bound intervals for all neurons in the third layer (second hidden layer); Blue: resulting approximation ratio. Measured for step \(t=10\) of the control benchmark. One iteration updates all bounds of layers not previously optimized in step \(t=9\) as well as all cutting hyperplanes. Note that the improvement from iterative tightening is two orders of magnitude for intermediate bounds, and four orders of magnitude for the approximation ratio.

Figure 6: Comparison of constraints enforced in forward propagation, SOTA controls methods, and our approach. Forward propagation (blue arrows) cannot tighten intermediate bounds with respect to an output constraint. ReBreach-LP and DRIP (blue and red arrows) use the output constraint to tighten the input bounds, and use these input bounds to tighten the intermediate bounds. We (blue, red, and green arrows) tighten all intermediate bounds with respect to input and output constraints simultaneously.
with respect to output bounds. The work of Vincent and Schwager 2021 explores utilizing an LP with complete neuron branching to verify control safety, which can be viewed as a domain-specific implementation of our MILP forumalation of the inverse verification problem. Both approaches share the ability to be terminated at any time for a looser approximation.

Certified OOD detection.There is a wide variety of OOD detection systems Yang et al. (2021); Salehi et al. (2022). Commonly, they are evaluated empirically based on known OOD data Wang et al. (2022); Liang et al. (2018); Chen et al. (2021). Therefore, they fail to provide verified guarantees for their effectiveness. In fact, many OOD detection systems are susceptible to adversarial attacks Sehwag et al. (2019); Chen et al. (2022). Meinke et al. (2021) show how to verify robustness to adversarial examples around given input samples. Berrada et al. (2021) develop a general framework for certifying properties of output distributions of neural networks given constraints on the input distribution. They study the problem of distributionally robust OOD detection, where OOD data can be perturbed by noise drawn from a collection of distributions. However, this work is still constrained to verifying properties of the outputs given constraints (albeit probabilistic) on the inputs. In contrast, INVPROP is able to certify arbitrary regions of the input space that lead to confident predictions.

## 6 Discussion

In this work we present the challenge of over-approximating neural network preimages and provide an efficient algorithm to solve it. By doing so, we demonstrate strong performance on multiple application areas. This work provides a first attempt at this problem and we believe there is a large scope for future investigation and new applications of our contribution.

\begin{table}
\begin{tabular}{l c c c} \hline \hline \multirow{2}{*}{Method} & Input & Approx & \multirow{2}{*}{Time (sec)} \\  & Branching & Ratio & \\ \hline MILP & no & 1.47 & 1562.26 \\ INVPROP & no & 4.39 & \(3.77\pm 0.02\) \\ INVPROP & yes & 1.14 & \(12.02\pm 0.06\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Comparison of methods for over-approximating the preimage of the OOD region shown in Figure 1 using 40 half-spaces.

Figure 7: In both plots, the green regions approximate inputs which the model believes are in-distribution (generated from 1 million samples). On the left, the blue border traces out the intersection of the half-spaces we tightened without branching. On the right, we do the same with input space branching and display the union of multiple sets. We observe that branching produced a tighter and more expressive over-approximation.
LimitationsOur current implementation does not implement neuron and output branching. The input branching strategy we implemented will be less effective for high dimensional inputs.

Future WorkTo scale this verification to higher-dimensional tasks such as standard image datasets would require addressing additional challenges. The preimage may not be well represented by the intersection of half-spaces. Moreover, the difficulty of uniformly sampling cutting planes \(\mathbf{c}\) will complicate measuring the efficacy of our methods.

Potential Negative Social ImpactOur work improves reliable ML through facilitating the deployment of systems that provably align with practitioner expectations. As such, we expect INVPROP to have positive societal impact. We acknowledge that our improvements to verification may be repurposed as improvements to attacks given access to a model though we believe the positive use cases of our technique greatly outweigh current speculation of misusage.

## 7 Acknowledgements

We thank Michael Everett and Nicholas Rober for helpful discussion and feedback on the paper.