Provable Preimage Under-Approximation for Neural Networks (Full Version)
Neural network verification mainly focuses on local robustness properties,
which can be checked by bounding the image (set of outputs) of a given input
set. However, often it is important to know whether a given property holds
globally for the input domain, and if not then for what proportion of the input
the property is true. To analyze such properties requires computing preimage
abstractions of neural networks. In this work, we propose an efficient anytime
algorithm for generating symbolic under-approximations of the preimage of any
polyhedron output set for neural networks. Our algorithm combines a novel
technique for cheaply computing polytope preimage under-approximations using
linear relaxation, with a carefully-designed refinement procedure that
iteratively partitions the input region into subregions using input and ReLU
splitting in order to improve the approximation. Empirically, we validate the
efficacy of our method across a range of domains, including a high-dimensional
MNIST classification task beyond the reach of existing preimage computation
methods. Finally, as use cases, we showcase the application to quantitative
verification and robustness analysis. We present a sound and complete algorithm
for the former, which exploits our disjoint union of polytopes representation
to provide formal guarantees. For the latter, we find that our method can
provide useful quantitative information even when standard verifiers cannot
verify a robustness property.
Xiyue Zhang, Benjie Wang, Marta Kwiatkowska
# On Preimage Approximation for Neural Networks

###### Abstract

Neural network verification mainly focuses on local robustness properties. However, often it is important to know whether a given property holds globally for the whole input domain, and if not then for what proportion of the input the property is true. While exact preimage generation can construct an equivalent representation of neural networks that can aid such (quantitative) global robustness verification, it is intractable at scale. In this work, we propose an efficient and practical anytime algorithm for generating symbolic under-approximations of the preimage of neural networks based on linear relaxation. Our algorithm iteratively minimizes the volume approximation error by partitioning the input region into subregions, where the neural network relaxation bounds become tighter. We further employ sampling and differentiable approximations to the volume in order to prioritize regions to split and optimize the parameters of the relaxation, leading to faster improvement and more compact under-approximations. Evaluation results demonstrate that our approach is able to generate preimage approximations significantly faster than exact methods and scales to neural network controllers for which exact preimage generation is intractable. We also demonstrate an application of our approach to quantitative global verification.

 Neural networks, abstraction, verification, linear relaxation

## I Introduction

Despite the remarkable empirical success of neural networks, guaranteeing their correctness, especially when using them as decision-making components in safety-critical autonomous systems [1, 2, 3], is an important and challenging task. Towards this aim, various approaches have been developed for the verification of neural networks, with extensive effort devoted to local robustness verification [4, 5, 6, 7, 8, 9, 10, 11, 12]. While local robustness verification focuses on deciding the absence of adversarial examples within an \(\epsilon\)-perturbation neighbourhood, an alternative approach for neural network analysis is to construct the preimage of its predictions [13, 14]. By characterizing the preimage symbolically in an abstract representation, e.g., polyhedra, one can perform more complex analysis for a wider class of properties beyond local robustness.

However, the exact preimage generation method of [13] takes time exponential in the number of neurons in a network. Meanwhile, the approximate preimage generation method proposed in [14] bypasses the intractability of exact preimage generation by leveraging symbolic interpolants [15, 16] for abstraction of neural network layers. However, due to the complexity of interpolation, the time to compute the abstraction also scales exponentially with the number of neurons in hidden layers. Therefore, more efficient computation methods for (symbolic abstraction of) preimages of neural networks are needed.

This paper makes the following novel contributions. We propose an efficient and practical _anytime_ algorithm for generating symbolic under-approximations of the preimage of piecewise linear neural networks as a union of disjoint polytopes. The algorithm assumes a hyperrectangle input domain and relies on linear relaxation based perturbation analysis (LiRPA) algorithms [9, 10, 11], applied backward from a polyhedron output set. Our algorithm partitions the input region into disjoint subregions, which can be approximated independently in parallel in a divide-and-conquer approach. To assess and optimize the volume, our method bypasses the computational cost of exact volume computation by making use of statistical estimation and differentiable approximations. As an application, we show how to soundly and effectively verify _global quantitative_ properties of neural networks. We take advantage of the efficiency of our algorithm to iteratively generate a high-quality under-approximation to the property, while invoking expensive exact computation only at the end of the algorithm. Finally, we conduct an empirical analysis of our method on a range of control and reinforcement learning tasks, showing significant gains in efficiency compared to exact preimage generation, and demonstrating verification of quantitative properties of vehicle parking and aircraft collision avoidance systems.

## II Related Work

Our work is related to a series of works on local robustness verification of neural networks. To address the scalability issues with _complete_ verifiers [4, 5, 8] based on constraint solving, convex relaxation [17] has been used for developing highly efficient _incomplete_ verification methods [6, 9, 10, 18]. Later works employed the branch-and-bound (BaB) framework [7, 19] to achieve completeness, using convex relaxation for the bounding procedure [11, 12, 20]. We adapt convex relaxation for efficient preimage approximation. There are also works that have sought to define a weaker notion of local robustness known as _statistical robustness_[21, 22], which requires that a certain proportion of points under some perturbation distribution around an input point are classified in the same way. This has the advantage of providing quantitative information about when adversarial robustness does not hold, or cannot be proven. Verification of statistical robustness is typically achieved by sampling and statistical guarantees [23, 24, 21, 25]. In this work, we apply our symbolic approximation approach to quantitative analysis of neural networks, while providing exact quantitative rather than statistical guarantees [26].
While neural network verification has primarily focused on _local_ robustness properties in the vicinity of a specific input point, there has been recent interest in _global_ robustness properties [27, 28, 29], which measure the robustness of a neural network over the entire input space/data distribution, rather than the vicinity of a specific point. Our work is within the scope of global property analysis, aiming at analyzing neural network behaviour in a larger input subspace which may contain semantically distinct inputs.

Another line of related works considers neural network abstraction to derive exact or approximate representations of neural networks, which are more amenable to analysis. Abstraction techniques have been widely applied for explanations [30], verification [31, 32], reachability analysis [33], and preimage approximation [14, 34]. [14] leverages symbolic interpolants [16] to compute preimage approximations; however, checking interpolation condition suffers from exponential complexity w.r.t. the number of hidden neurons. [34] uses mixed integer linear programming (MILP) to compute a preimage over-approximation with the intersection of a finite set of (cutting-plane) constraints. In contrast, we construct an under-approximation abstraction of the preimage using a disjoint set of polytopes by employing efficient linear lower/upper bounds on nonlinear neurons of the neural network, and the derived symbolic abstraction can be further applied to provide provable guarantees on quantitative properties.

## III Preliminaries

We use \(f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}\) to denote a feedforward neural network. For layer \(i\), we use \(\textbf{W}^{(i)}\) to denote the weight matrix for the layer, \(\textbf{b}^{(i)}\) the bias, \(h^{(i)}\) the pre-activation neurons, and \(a^{(i)}\) the post-activation neurons, such that we have \(h^{(i)}=\textbf{W}^{(i)}a^{(i-1)}+\textbf{b}^{(i)}\). We focus on ReLU neural networks where \(\text{ReLU}(h):=\max(h,0)\).

**Linear Relaxation of Neural Networks.** Nonlinear activation functions lead to the NP-completeness of the neural network verification problem [5]. To address such intractability, linear relaxation is used to transform the nonconvex constraints into linear programs. As shown in Figure 1, given _concrete_ lower and upper bounds \(\textbf{l}^{(i)}\leq h^{(i)}(x)\leq\textbf{u}^{(i)}\) on the pre-activation values of layer \(i\), the post-activation neurons of layer \(i\), \(a^{(i)}_{j}(x)=ReLU(h^{(i)}_{j}(x))\) can be bounded by the following linear lower and upper bounding functions w.r.t. \(h^{(i)}_{j}(x)\):

\[\left\{\begin{array}{ll}0\cdot h^{(i)}_{j}(x)\leq a^{(i)}_{j}(x)\leq 0 \cdot h^{(i)}_{j}(x),&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad
where \(T\) consists of all values of \(x\) satisfying the first-order logic (FOL) formula \(\alpha(x):=\bigwedge_{i=1}^{K}\psi_{i}(x)\). We use the term polytope to refer to a bounded polyhedron, that is, a polyhedron \(T\) such that \(\exists R\in\mathbb{R}^{>0}:\forall x_{1},x_{2}\in T\), \(\left\|x_{1}-x_{2}\right\|_{2}\leq R\) holds. The abstract domain of polyhedra [35, 9, 36] has been widely used for the verification of neural networks and conventional programs. An important type of polytope is the hyperrectangle (box), which is a polytope defined by a closed and bounded interval \([\underline{x}_{i},\overline{x_{i}}]\) for each dimension, where \(\underline{x}_{i},\overline{x_{i}}\in\mathbb{Q}\). More formally, using the linear constraints \(\phi_{i}:=(\overline{x}_{i}\geq\underline{x}_{i})/\land(x_{i}\leq\overline{ x_{i}})\) for each dimension, the hyperrectangle takes the form \(\mathcal{C}=\{x\in\mathbb{R}^{d}|x\models\bigwedge_{i=1}^{d}\phi_{i}\}\).

## IV Problem Formulation

### _Approximate Preimage Generation_

We are interested in the problem of computing preimage abstraction for neural networks. Given a subset \(O\subset\mathbb{R}^{m}\) of the codomain, the preimage of a function \(f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}\) is defined to be the set of all inputs \(x\in\mathbb{R}^{d}\) that are mapped to an element of \(O\) by \(f\). For neural networks in particular, the input is typically restricted to some bounded input region \(\mathcal{C}\subset\mathbb{R}^{d}\). In this work, we restrict the output set \(O\) to be a polyhedron, and the input set \(\mathcal{C}\) to be an axis-aligned hyperrectangle region \(\mathcal{C}\subset\mathbb{R}^{d}\). We now define the notion of a restricted preimage:

**Definition 1** (Restricted Preimage).: _Given a neural network \(f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}\), and an input set \(\mathcal{C}\subset\mathbb{R}^{d}\), the restricted preimage of an output set \(O\subset\mathbb{R}^{m}\) is defined to be the set \(f_{\mathcal{C}}^{-1}(O):=\{x\in\mathbb{R}^{d}|f(x)\in O\wedge x\in\mathcal{C}\}\)._

**Example 1**.: _To illustrate our problem formulation and approach, we introduce a vehicle parking task [37] as a running example. In this task, there are four parking lots, located in each quadrant of a \(2\times 2\) grid \([0,2]^{2}\), and a neural network with two hidden layers of 10 ReLU neurons \(f:\mathbb{R}^{2}\rightarrow\mathbb{R}^{4}\) is trained to classify which parking lot an input point belongs to. To analyze the behaviour of the neural network in the input region \([0,1]\times[0,1]\) corresponding to parking lot 1, we set \(\mathcal{C}=\{x\in\mathbb{R}^{2}|(0\leq x_{1}\leq 1)\land(0\leq x_{2}\leq 1)\}\). Then the restricted preimage \(f_{\mathcal{C}}^{-1}(O)\) of the set \(O=\{\boldsymbol{y}\in\mathbb{R}^{4}|\bigwedge_{i\in\{2,3,4\}}y_{1}-y_{i}\geq 0\}\) is then the subspace of the region \([0,1]\times[0,1]\) that is labelled as parking lot \(1\) by the network._

Unfortunately, the (restricted) preimage of an output set of a neural network is expensive to compute or represent exactly. Thus, we resort to deriving approximations, or _abstractions_ of the preimage, that we can efficiently manipulate and analyze. In particular, we focus on _provable under-approximations_ of the preimage. Given a first-order formula \(A\), \(\alpha\) is an under-approximation of \(A\) if it holds that \(\forall x.\alpha(x)\implies A(x)\). In our context, the restricted preimage is defined by the formula \(A(x)=(f(x)\in O)\land(x\in\mathcal{C})\), and we restrict to under-approximations \(\alpha\) that take the form of a disjoint union of polytopes.

**Definition 2** (Disjoint Union of Polytopes).: _A disjoint union of polytopes (DUP) is a FOL formula \(\alpha\) of the form \(\alpha(x):=\bigcup_{i=1}^{D}\alpha_{i}(x)\) where each \(\alpha_{i}\) is a polytope formula (conjunction of linear half-space constraints), with the property that \(\alpha_{i}\land\alpha_{j}\) is unsatisfiable for any \(i\neq j\)._

We will also refer to the set of points \(\mathcal{T}=\{x\in\mathbb{R}^{d}|\alpha(x)\}\) satisfying \(\alpha\) as a disjoint union of polytopes. The goal of our method is to generate a DUP under-approximation \(\mathcal{T}\) that maximizes the _volume_\(\text{vol}(\mathcal{T})\) covered by the under-approximation. In our method, we use the volume both to optimize the under-approximation, and as a stopping criterion when the under-approximation reaches a sufficient coverage. Disjoint unions of polytopes have the advantage that we can often analyze the polytopes independently, allowing us to leverage parallel computation. For example, to compute the volume of a DUP \(\mathcal{T}\), it suffices to compute the sum of volumes of each individual polytope \(T_{i}=\{x\in\mathbb{R}^{d}|\alpha_{i}(x)\}\), which is easier to compute.

### _Quantitative Properties_

One of the most important verification problems for neural networks is that of proving guarantees on the output of a network for a given input set [38, 39, 40].

**Definition 3** (Network Properties).: _Given a neural network \(f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}\), a precondition (input set) \(I\subseteq\mathbb{R}^{d}\) and a postcondition (output set) \(O\subseteq\mathbb{R}^{m}\), we say that the neural network satisfies the property \((I,O)\) if \(x\in I\implies f(x)\in O\)._

This definition is commonly used to encode _safety_ properties, where the "unsafe region" is the complement of \(O\). A notable example of a safety problem is that of _local robustness_, where \(I\) is a perturbation region \(\left\|x^{\prime}-x\right\|_{p}\leq\epsilon\) around a fixed input \(x\), and \(O\) is the output region corresponding to a particular class label. The goal is then to prove whether the property holds, or else find a sufficiently small \(\epsilon\in\mathbb{R}^{>0}\) such that the property holds. This is usually achieved by _forward_ computation methods, such as abstract interpretation [38], which compute over-approximations of the _output_ set corresponding to the given input set \(I\); if the over-approximation is contained within \(O\), then the property is verified. Alternatively, we could verify such a property via _backward_ computation, by computing an _under-approximation_ of the pre-image of \(O\)[14]; if the input set \(I\) is contained within the under-approximation, then the property is verified.

When we consider _global_ properties that cover a larger input region and may consist of semantically distinct inputs, this property formulation often becomes inadequate. Firstly, if we cannot completely verify safety of the whole region, it is often preferable to obtain a quantitative guarantee of what proportion of the inputs satisfy the output condition, rather than restricting the size of the input region [23]. Secondly, many global properties are naturally expressed in a quantitative manner. For example, for an aircraft collision avoidance system, we might expect that the vast majority, but not all of the typical input space (say, \(>90\%\)) should result in a clear of conflict (COC) output.
**Definition 4** (Quantitative Property).: _Given a neural network \(f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}\), a measurable input set with non-zero measure (volume) \(I\subseteq\mathbb{R}^{d}\), a measurable output set \(O\subseteq\mathbb{R}^{m}\), and a rational proportion \(p\in[0,1]\) we say that the neural network satisfies the property \((I,O,p)\) if \(\frac{\text{vol}(f_{I}^{-1}(O))}{\text{vol}(I)}\geq p\). 1_

Footnote 1: In particular, the restricted preimage of a polyhedron under a neural network is Lebesgue measurable since polyhedra (intersection of half-spaces) are Borel measurable and NNs are continuous functions.

If the property \((I,O)\) holds, then the quantitative property \((I,O,1)\) holds with proportion \(1\), while quantitative properties provide more information when \((I,O)\) does not hold. Notice that forward over-approximations of the output cannot be used to directly verify quantitative properties, as the approximation is in the output space. On the other hand, backward under-approximations allow us to approximate \(\text{vol}(f_{I}^{-1}(O))\) using our preimage under-approximation. We will later show how to use our preimage approximation method to verify quantitative properties when the input set is a polytope and the output set is a polyhedron.

Borrowing from the definitions for non-quantitative neural network properties [41], we now define soundness and completeness of verification algorithms for quantitative properties:

**Definition 5** (Soundness).: _A verification algorithm \(QV\) is sound if, whenever \(QV\) outputs \(\text{True}\), the property \((I,O,p)\) holds._

**Definition 6** (Completeness).: _A verification algorithm \(QV\) is complete if (i) \(QV\) never returns \(\text{Unknown}\), and (ii) whenever \(QV\) outputs \(\text{False}\), the property \((I,O,p)\) does not hold._

## V Methodology

**Overview.** In this section, we present our _anytime_ algorithm for approximate preimage generation. The algorithm maintains and iteratively refines/grows a set of polytopes \(\mathcal{T}\) as a guaranteed under-approximation to the restricted preimage \(f_{\mathcal{C}}^{-1}(O)\), such that the volume of the approximation improves after each iteration. This is achieved by iteratively splitting the original input hyperrectangle \(\mathcal{C}\) via axis-aligned cuts into a set of hyperrectangle subregions, and approximating the preimage restricted to each subregion with a single polytope. As the subregions are disjoint, we obtain a DUP approximation after each iteration. The overall method is shown in Algorithm 1.

Our method consists of three main components. Firstly, in Section V-A, we show how to cheaply compute single polytope under-approximations to the restricted preimage, by adapting efficient linear relaxation methods. Next, in Section V-B, we introduce our refinement procedure that improves the approximation by partitioning a (sub)region into two subregions, which can then be approximated more accurately. We show how to choose the region to split and the input feature to cut on to maximize coverage (volume), leveraging parallel computation and Monte Carlo sampling for efficiency. In Section V-C, we further improve the method by optimizing the linear relaxation bounds to maximize volume coverage. Finally, in Section V-D, we show how our preimage generation method can be used to verify quantitative properties.

```
Input: Input region \(\mathcal{C}\), Output region \(O\), Volume threshold \(v\), Maximum iterations \(R\) Output: Disjoint union of polytopes \(\mathcal{T}\)
1\(T\leftarrow\text{GentUnderApprox}(\mathcal{C},O)\); \(\widetilde{\text{vol}_{T}},\widetilde{\text{vol}_{f_{\mathcal{C}}^{-1}(O)}} \leftarrow\text{EstimateVol}(T)\), EstimateVol(\(f_{\mathcal{C}}^{-1}(O)\)) ; \(\widetilde{\text{Dom}}\leftarrow\{(\mathcal{C},T,\text{vol}_{f_{\mathcal{C}}^{-1}(O )}-\widetilde{\text{vol}_{T}})\}\) ; // Priority queue // \(\mathcal{T}\)\(\widetilde{\text{Dom}}\) is the union of polytopes in Dom
2whileEstimateVol\((\mathcal{T}_{\text{Dom}})<v\)and Iterations \(\leq R\)do
3\(\mathcal{C}_{\text{sub}},T,\text{Priority}\leftarrow\text{Pop}(\text{Dom})\) ; // Subregion with highest priority
4\([\mathcal{C}_{\text{sub}}^{1},\mathcal{C}_{\text{sub}}^{1,\text{u}},\)..., \(\mathcal{C}_{\text{sub}}^{d,\text{u}},\mathcal{C}_{\text{sub}}^{d,\text{u}}] \leftarrow\text{SplitOnFeature}(\mathcal{C}_{\text{sub}})\); / Split selected subregion into two subregions w.r.t. all features
5\([T^{1,l},T^{1,u},...,T^{d,l},T^{d,u}]\leftarrow\)
6GenUnderApprox(\([\mathcal{C}_{\text{sub}}^{1,l},\mathcal{C}_{\text{sub}}^{1,\text{u}},\)..., \(\mathcal{C}_{\text{sub}}^{d,l},\mathcal{C}_{\text{sub}}^{d,u}]\), \(O\)) ; // Generate preimage in parallel
7\([\text{vol}_{T^{1,l}},\widetilde{\text{vol}_{T^{1,u}},...,\widetilde{\text{ vol}_{T^{d,l}}},\widetilde{\text{vol}_{T^{d,l}}}}]\leftarrow\) EstimateVol(\([T^{1,l},T^{1,u},...,T^{d,l},T^{d,u}]\)); \(id\leftarrow\text{arg}\max_{i}\left(\widetilde{\text{vol}_{T^{i,l}}}+\widetilde{ \text{vol}_{T^{i,u}}}\right)\) ; // Select the splitting feature
8\(\text{vol}_{T^{i-1}_{\mathcal{C}_{\text{sub}}^{i,l}(O)}}^{-1},\text{vol}_{T^{i -1}_{\mathcal{C}_{\text{sub}}^{i,l}(O)}}^{-1}\leftarrow\) EstimateVol(\(f_{\mathcal{C}_{\text{sub}}^{i,l}(O)}^{-1}\)), EstimateVol(\(f_{\mathcal{C}_{\text{sub}}^{i,l}(O)}^{-1}\)) ; \(\text{Dom}\leftarrow\text{Dom}\cup\{(\mathcal{C}_{\text{sub}}^{id,l},T^{id,u}, \widetilde{\text{vol}_{T^{-1}_{\mathcal{C}_{\text{sub}}^{i,l}(O)}}-\widetilde{ \text{vol}_{T^{id,l}}}})\}\) ; \(\cup\{(\mathcal{C}_{\text{sub}}^{id,u},T^{id,u},\widetilde{\text{vol}_{T^{-1}_{ \mathcal{C}_{\text{sub}}^{id,u}}}}(O)-\widetilde{\text{vol}_{T^{id,u}}})\}\) ; // Disjoint partition of the input
9return\(\mathcal{T}_{\text{Dom}}\)
```

**Algorithm 1**Preimage Under-Approximation

### _Polytope Under-Approximation via Linear Relaxation_

We begin by adapting linear relaxation techniques in order to cheaply generate valid under-approximations to the restricted preimage for a given input region \(\mathcal{C}\). Recall that LiRPA methods enable us to obtain linear lower and upper bounds on the output of a neural network \(f\), that is, \(\underline{\boldsymbol{\Lambda}}x+\underline{\boldsymbol{\mathsf{b}}}\leq f(x) \leq\overline{\boldsymbol{\Lambda}}x+\overline{\boldsymbol{\mathsf{b}}}\), where the linear coefficients depend on the input region \(\mathcal{C}\).

Now, suppose that we are interested in computing an under-approximation to the restricted preimage, given the input hyperrectangle \(\mathcal{C}=\{x\in\mathbb{R}^{d}|x\models\bigwedge_{i=1}^{d}\phi_{i}\}\), and the output polytope specified using the half-space constraints \(\psi_{i}(y)=(c_{i}^{T}y+d_{i}\geq 0)\) for \(i=1,...,K\) over the output space. For each such constraint, we append an additional linear layer at the end of the network \(f\), which maps \(y\mapsto c_{i}^{T}y+d_{i}\), such that the function \(g_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}\) represented by the new network is \(g_{i}(x)=c_{i}^{T}f(x)+d_{i}\). Then, applying LiRPA bounding to each \(g_{i}\), we obtain lower bounds \(\underline{g_{i}}(x)=\underline{a_{i}^{T}}x+\underline{b_{i}}\) for each \(i\), such that \(\underline{g_{i}}(x)\geq 0\implies g_{i}(x)\geq 0\) for \(x\in\mathcal{C}\). Notice that, for each \(i=1,...,K\), \(\underline{a_{i}^{T}}x+\underline{b_{i}}\geq 0\) is a half-space constraint
in the input space. We conjoin these constraints, along with the restriction to the input region \(\mathcal{C}\), to obtain a polytope

\[T_{\mathcal{C}}(O):=\{x|\bigwedge_{i=1}^{K}(\underline{g_{i}}(x)\geq 0) \wedge\bigwedge_{i=1}^{d}\phi_{i}(x)\} \tag{5}\]

**Proposition 1**.: \(T_{\mathcal{C}}(O)\) _is an under-approximation to the restricted preimage \(f_{\mathcal{C}}^{-1}(O)\)._

In our algorithm, we use CROWN [6] to generate the LiRPA bounds in parallel over the output polytope constraints \(i=1,...,K\), and store the resulting input polytope \(T_{\mathcal{C}}(O)\) as a list of constraints. This procedure is highly efficient, enabling us to employ it as a sub-routine LinearLowerBound in our overall algorithm (Line 4 of Algorithm 2).

### _Global Branching and Refinement_

As LiRPA performs crude linear relaxation, the resulting bounds can be quite loose, meaning that the polytope approximation \(T_{\mathcal{C}}(O)\) is unlikely to constitute a tight under-approximation to the restricted preimage. To address this challenge, we employ a divide-and-conquer approach that iteratively refines our under-approximation of the preimage. Starting from the single initial region \(\mathcal{C}\) represented at the root, our method generates a tree by iteratively partitioning the subregion \(\mathcal{C}_{sub}\) represented at a leaf node into two smaller subregions \(\mathcal{C}_{sub}^{l},\mathcal{C}_{sub}^{u}\) (via an axis-aligned bisection), which are then attached as children to that leaf node. In this way, the subregions represented by all leaves of the tree are disjoint, such that their union is the initial region \(\mathcal{C}\).

After each iteration in Algorithm 1, each leaf subregion \(\mathcal{C}_{sub}\) has an associated polytope, computed using LiRPA bounds in Algorithm 2, that approximates the restricted preimage in \(\mathcal{C}_{sub}\). Thus, the union of the polytopes corresponding to all leaves forms our refined, anytime DUP under-approximation \(\mathcal{T}\) to the preimage in the original region \(\mathcal{C}\). The algorithm terminates after a desired volume threshold \(v\) for \(\text{vol}(\mathcal{T})\) has been reached, or after a fixed maximum number of iterations (polytopes) has been reached.

Unfortunately, even with a moderate number of input dimensions \(d\), naively splitting along all dimensions quickly becomes computationally infeasible. For example, splitting a \(d\)-dimensional hyperrectangle using bisections along each dimension results in \(2^{d}\) subdomains to approximate. It thus becomes crucial to identify the subregion splits that have the most impact on the quality of the under-approximation. At each iteration, given the subregion tree, our algorithm must first select which leaf subregion to split, and then which input feature/dimesion to bisect on, both of which will impact the refinement quality.

**Subregion Selection.** Searching through all leaf subregions at each iteration is too computationally expensive. Thus, we propose a subregion selection strategy that prioritizes splitting subregions according to the difference in volume between the exact restricted preimage \(f_{C_{sub}}^{-1}(O)\) and the (already computed) polytope approximation \(T_{\mathcal{C}_{sub}}(O)\) on that subdomain, that is:

\[\text{Priority}(\mathcal{C}_{sub})=\text{vol}(f_{C_{sub}}^{-1}(O))-\text{ vol}(T_{\mathcal{C}_{sub}}(O)) \tag{6}\]

This measures the gap between the polytope under-approximation and the optimal approximation, namely, the preimage itself. A larger value suggests that, by refining the approximation of the preimage in that subregion, we can improve the volume of the under-approximation in this subregion significantly. Crucially, since the polytopes in separate subregions are all disjoint, this corresponds directly to improving the volume of the overall DUP under-approximation.

Suppose that a particular leaf subdomain attains the maximum of this metric among all leaves, and we split it into two subregions \(\mathcal{C}_{sub}^{l},\mathcal{C}_{sub}^{u}\), which we approximate with polytopes \(T_{\mathcal{C}_{sub}^{l}}(O),T_{\mathcal{C}_{sub}^{u}}(O)\). As tighter intermediate concrete bounds, and thus linear bounding functions, can be computed on the partitioned subregions, the polytope approximation on each subregion will be refined compared with the single polytope restricted to that subregion.

**Proposition 2**.: _Given any subregion \(\mathcal{C}_{sub}\) with polytope approximation \(T_{\mathcal{C}_{sub}}(O)\), and its children \(\mathcal{C}_{sub}^{l},\mathcal{C}_{sub}^{u}\) with polytope approximations \(T_{\mathcal{C}_{sub}^{l}}(O),T_{\mathcal{C}_{sub}^{u}}(O)\) respectively, it holds that:_

\[\text{vol}(T_{\mathcal{C}_{sub}^{l}}(O))+\text{vol}(T_{\mathcal{ C}_{sub}^{u}}(O))\geq\text{vol}(T_{\mathcal{C}_{sub}}(O)) \tag{7}\] \[\text{Priority}(\mathcal{C}_{sub})\geq\text{Priority}(\mathcal{ C}_{sub}^{l})+\text{Priority}(\mathcal{C}_{sub}^{u}) \tag{8}\]

That is, we replace the leaf with maximal priority with two leaves with lower priority, which sum to at most the priority of the original leaf. In the next iteration, we may choose either of these leaves to split, or another leaf in the tree which now has maximal priority. Further, given the volume guarantee in Equation 7, we have the following Corollary:

**Corollary 1**.: _In each iteration of Algorithm 1, the volume of the polytope approximation \(\mathcal{T}_{Dom}\) does not decrease._

In practice, we may not be able to compute the volumes in Equation 6 exactly. On the one hand, computing the volume of a polytope is a computationally expensive task
requiring specialized tools [42], and we do not want to add such runtime cost during each iteration; more pertinently, we cannot compute (the volume of) the exact restricted preimage corresponding to the output specifications. Thus, we resort to statistical estimation of the volumes. Since subdomains \(\mathcal{C}_{sub}\) are hyperrectangles, it is easy to sample \(N\) points \(x_{1},...,x_{N}\) uniformly from the subdomain. Then we can employ Monte Carlo estimation for both the preimage and the generated polytope approximation:

\[\widehat{\text{vol}}(f_{\mathcal{C}_{sub}}^{-1}(O))=\frac{\sum_{i=1}^{N} \mathds{1}_{x_{i}\in\mathcal{I}\bar{c}_{sub}^{-1}(O)}}{N} \tag{9}\]

\[\widehat{\text{vol}}(T_{\mathcal{C}_{sub}}(O))=\frac{\sum_{i=1}^{N}\mathds{1} _{x_{i}\in\mathcal{I}\bar{c}_{sub}(O)}}{N} \tag{10}\]

**Splitting Feature.** Once we have selected a subregion to refine, we then decide how to split the subregion. Given a subregion (hyperrectangle) defined by lower and upper bounds \(x_{i}\in[\underline{z}_{i},\overline{z}_{i}]\) for all dimensions \(i=1,...,d\), we split it into two subregions by bisecting along some feature \(i\). This splitting procedure will produce two subregions which are similar to the original subregion, but have updated bounds \([\underline{z}_{i},\frac{\underline{z}_{i}+\overline{z}_{i}}{2}],[\frac{ \underline{z}_{i}+\overline{z}_{i}}{2},\overline{z}_{i}]\) for feature \(i\) instead. In order to determine the weight feature/dimension to split on, we propose a greedy strategy. Specifically, for each feature, we generate a pair of polytopes via LiRPA for the two subregions resulting from the split, and choose the feature that results in the greatest total volume of the polytope pair. In practice, another commonly-adopted splitting heuristic is to select the dimension with the longest edge; that is to select feature \(i\) with the largest range: \(\arg\max_{i}(\overline{z}_{i}-\underline{z}_{i})\). Performance evaluation of these two splitting methods is presented in Section VI-B.

**Example 2**.: _We revisit the vehicle parking problem in Example 1. Figure 1(b) shows the polytope under-approximation computed on the input region \(\mathcal{C}\) before refinement, where each solid line represents the cutting plane for each output specification (\(y_{1}-y_{i}\geq 0\)). Figure 1(c) depicts the refined approximation by splitting the input region along the vertical axis, where the solid and dashed lines represent the cutting planes for the two resulting subregions. It can be seen that the total volume of the under-approximation has improved significantly._

### _Local Optimization_

One of the key components behind the effectiveness of LiRPA-based bounds is the ability to efficiently improve the tightness of the bounding function by optimizing the relaxation parameters \(\mathbf{\alpha}\), via projected gradient descent. In the context of local robustness verification, the goal is to optimize the concrete lower or upper bounds over the (sub)region \(\mathcal{C}\)[10], i.e., \(\min_{x\in\mathcal{C}}\underline{\mathbf{\Lambda}}(\mathbf{\alpha})x+\underline{ \textbf{b}}(\mathbf{\alpha})\), where we explicitly note the dependence of the linear coefficients on \(\mathbf{\alpha}\). In our case, we are instead interested in optimizing \(\mathbf{\alpha}\) to refine the polytope under-approximation, that is, increase its volume. Unfortunately, as before, computing the volume of a polytope exactly is computationally expensive, and further does not allow for computation of the gradients with respect to \(\mathbf{\alpha}\).

To address this challenge, we propose a loss function to encode a differentiable relaxation of the polytope volume. We have seen that we can estimate the volume of a polytope using a set of samples \(x_{1},...x_{N}\) as \(\widehat{\text{vol}}(T_{\mathcal{C}_{sub},\mathbf{\alpha}}(O))=\frac{1}{N}\sum_{i =1}^{N}\mathds{1}_{x_{i}\in T_{\mathcal{C}_{sub},\mathbf{\alpha}}}(O)\), where we have highlighted the dependence of the polytope \(T_{\mathcal{C}_{sub}}(O)=\{x|\bigwedge_{i=1}^{K}\underline{g}_{i}(x,\mathbf{\alpha} _{i})\geq 0\wedge\bigwedge_{i=1}^{d}\phi_{i}(x)\}\) on \(\mathbf{\alpha}=(\mathbf{\alpha}_{1},...,\mathbf{\alpha}_{K})\), and \(\mathbf{\alpha}_{i}\) are the \(\alpha\)-parameters for the linear relaxation of the neural network \(g_{i}\) corresponding to the \(i^{\text{th}}\) half-space constraint in \(O\). However, this is still non-differentiable w.r.t. \(\mathbf{\alpha}\) due to the identity function. We now show how to derive a differentiable relaxation which is amenable to gradient-based optimization:

\[\widehat{\text{vol}}(T_{\mathcal{C}_{sub}}(O))=\frac{1}{N}\sum_{j=1}^{N} \mathds{1}_{x_{j}\in T_{\mathcal{C}_{sub},\mathbf{\alpha}}(O)} \tag{11}\]

\[=\frac{1}{N}\sum_{j=1}^{N}\mathds{1}_{\min_{i=1,...K}\underline{g}_{i}(x_{j}, \mathbf{\alpha}_{i})\geq 0} \tag{12}\]

\[\approx\frac{1}{N}\sum_{j=1}^{N}\sigma\left(\min_{i=1,...K}\underline{g}_{i}(x_ {j},\mathbf{\alpha}_{i})\right) \tag{13}\]

\[\approx\frac{1}{N}\sum_{j=1}^{N}\sigma\left(-\text{LSE}(-\underline{g}_{1}(x_ {j},\mathbf{\alpha}_{1}),...,\underline{g}_{K}(x_{j},\mathbf{\alpha}_{K}))\right) \tag{14}\]

The second equality follows from the definition of the polytope \(T_{\mathcal{C}_{sub},\mathbf{\alpha}}(O)\); namely that a point is in the polytope if it satisfies \(\underline{g}_{i}(x_{j},\mathbf{\alpha}_{i})\geq 0\) for all \(i=1,...,K\), or equivalently, \(\min_{i=1,...K}\underline{g}_{i}(x_{j},\mathbf{\alpha}_{i})\geq 0\). After this, we approximate the

Fig. 2: Refinement and optimization for preimage approximation.
identity function using a sigmoid relaxation, where \(\sigma(y):=\frac{1}{1+e^{-y}}\), as is commonly done in machine learning to define classification losses. Finally, we approximate the minimum over specifications using the log-sum-exp (LSE) function. The log-sum-exp function is defined by \(LSE(y_{1},...,y_{K}):=\log(\sum_{i=1,...,K}e^{y_{i}})\), and is a differentiable approximation to the maximum function; we employ it to approximate the minimization by adding the appropriate sign changes. The final expression is now a differentiable function of \(\mathbf{\alpha}\). We employ this as the loss function in Algorithm 2 (Line 6), and optimize using projected gradient descent.

It is worth noting that employing \(\mathbf{\alpha}\)-optimization breaks the guarantee in Proposition 2, for two reasons. Firstly, due to the sampling and differentiable relaxation, the optimization objective is not exact volume. Secondly, gradient-based optimization cannot guarantee improvement in its objective after each update. Nonetheless, \(\mathbf{\alpha}\)-optimization can significantly improve the approximation in practice, and so we employ it in our method.

**Example 3**.: _We revisit the vehicle parking problem in Example 1. Figure 1(a) and 1(b) show the computed under-approximations before and after local optimization. We can see that the cutting lines for all three specifications are optimized, which effectively improves the approximation quality._

```
Input: Property \((I,O,p)\), Maximum iterations \(R\) Output: Verification Result \(\in\) {True, False, Unknown} \(\mathsf{vol}(I)\leftarrow\text{ExactVolume}(I)\); \(\mathcal{C}\leftarrow\text{QuerBox}(I)\); // For general polytopes \(I\)\(\mathcal{T}\leftarrow\text{InitialRun}(\mathcal{C},O)\); while Iterations \(\leq R\)do
1\(\mathcal{T}\leftarrow\text{Refine}(\mathcal{T},\mathcal{C},O)\);
2if\(\text{EstimateVolume}(\mathcal{T})\geq p\times\text{vol}(I)\)then
3if\(\text{ExactVolume}(\mathcal{T})\geq p\times\text{vol}(I)\)thenreturn True return Unknown
```

**Algorithm 3**Quantitative Verification

### _Quantitative Verification_

Given a quantitative property \((I,O,p)\), where \(O\) is a polyhedron and \(I\) a polytope, we now show how to use our efficient preimage under-approximation method to verify the property. Assume for now that \(I\) is a hyperrectangle, so that we can take \(\mathcal{C}=I\) (the case of a general polytope is discussed in the supplementary material).

In order to verify the property, we can set the volume threshold for Algorithm 1 to be \(p\times\text{vol}(I)\), such that we have \(\frac{\overline{\text{vol}(\mathcal{T})}}{\text{vol}(I)}\geq p\) if the algorithm terminates before reaching the maximum iterations. However, our Monte Carlo estimates of volume (Line 4) cannot provide a sound guarantee that \(\frac{\text{vol}(\mathcal{T})}{\text{vol}(I)}\geq p\). To resolve this problem, we propose to run exact volume computation [43] only when the Monte Carlo estimate reaches the threshold. If the exact volume \(\text{vol}(\mathcal{T})\geq p\times\text{vol}(I)\), then the property is verified. Otherwise, we continue running the preimage refinement. This procedure is shown in Algorithm 3, where InitialRun generates an initial approximation to the restricted preimage as in Lines 1-3 of Algorithm 1, and Refine performs one iteration of approximation refinement as in Lines 5-11 of Algorithm 1. Termination occurs when we have verified the quantitative property, or when the maximum number of iterations has been exceeded.

We see that the algorithm outputs _True_ only if the _exact_ volume of the DUP approximation \(\mathcal{T}\) exceeds the threshold, even though sampling is used during the approximation refinement.

**Theorem 1**.: _Algorithm 3 is sound for quantitative verification._

Algorithm 3 is not complete as it outputs Unknown when the iterations budget is exceeded. A possible modification to the algorithm is to remove the threshold, and terminate when the subregions are "sufficiently small". By sufficiently small, we mean that, for a subregion, the derived concrete bounds for each neuron in the neural network are such that the neuron is stable (Equation 1); such "sufficiently small" subregions exist by continuity. If this is the case, the linear lower (and upper) bound \(\underline{\mathbf{\Delta}x}+\underline{\mathbf{b}}\) is exact over that subregion. This means the polytope approximation to the restricted preimage will also be exact. Thus, we can split our original region until all subregions satisfy the condition.

Unfortunately, this is still not enough to guarantee completeness, because our choice of subregion to split is probabilistic (uses sampling). That is, there exists a run of the algorithm where it splits within the interval \([0,0.5]\) (say) in every iteration, but never splits on \([0.5,1]\). In practice, for even moderately sized neural networks, it is computationally infeasible to obtain sufficiently small subregions, and so terminating according to a time/iterations budget is more appropriate. Complete proofs for propositions and theorems are provided in the supplementary material.

## VI Experiments

We perform experimental evaluation of the proposed approach on a set of benchmark tasks and demonstrate its effectiveness in approximation precision and efficiency, as well as an application to quantitative analysis of global properties. We aim to answer the following research questions (RQs): \(\mathbf{\Theta}\) How effective is our approach in preimage approximation? \(\mathbf{\Theta}\) How effective are our global refinement and local optimization methods and how do the parameter configurations affect the performance of our approach? \(\mathbf{\Theta}\) Can our approach be applied to quantitative verification of neural networks?

### _Benchmark_

We focus on preimage analysis for neural networks on a benchmark of control and reinforcement learning (RL) tasks, where global analysis on the entire input space or a critical input subspace is especially needed. Besides the vehicle parking task [37] shown in the example, we use the following benchmarks:
* Aircraft collision avoidance system (VCAS) [44]: The VCAS system is used to provide advisory for collision avoidance between the ownship aircraft and the intruder. VCAS uses four input features \((h,h_{A},h_{B},t)\) representing the relative altitude of the aircrafts, vertical climbing rates of the ownship and intruder aircrafts, respectively, and time to the loss of horizontal separation. VCAS is implemented by nine feed-forward neural networks built in with a hidden layer of 21 neurons and an output layer of nine neurons corresponding to different advisories.
* Neural network controllers: We use the trained neural network controllers from VNN-COMP 2022 [45] for three RL tasks: Cartpole, Lunarlander, and Dubinsrejoin [46]. The neural networks for Cartpole and Lunarlander tasks have two hidden layers with 64 neurons, and the neural network for Dubinsrejoin has two hidden layers with 256 neurons. (1) The neural network controller for Cartpole aims to balance a pole atop a cart by controlling the movement of the cart, which has four input variables representing the position and velocity of the cart and the angle and angular velocity of the pole, and has two output actions. (2) The controller for Lunarlander tackles the task of correct landing of a moon lander on a pad, which has eight input features and four output actions. (3) The controller for Dubinsrejoin targets guiding a wingman craft to a certain radius around a lead aircraft, where the input and output spaces are both eight dimensional.

### _Evaluation_

#### Iv-B1 Evaluation Metric

To evaluate the quality of the preimage approximation generated from our approach, we define the _coverage ratio_ to be the ratio of volume covered to the volume of the restricted preimage, i.e. \(\text{cov}(\mathcal{T},f_{\mathcal{C}}^{-1}(O)):=\frac{\text{vol}(\mathcal{ T})}{\text{vol}(f_{\mathcal{C}}^{-1}(O))}\). Since \(\mathcal{T}\) is an under-approximation to \(f_{\mathcal{C}}^{-1}(O)\), we have \(\text{cov}(\mathcal{T},f_{\mathcal{C}}^{-1}(O))\in[0,1]\). This gives us a normalized measure of the quality of the approximation, which can be compared across preimages for different output sets (e.g., the preimage of different labels), where \(\text{vol}(f_{\mathcal{C}}^{-1}(O))\) may vary. In practice, before running preimage generation, we estimate \(\text{vol}(f_{\mathcal{C}}^{-1}(O))\) as \(\text{vol}(\widehat{f_{\mathcal{C}}^{-1}(O)})=\frac{1}{N}\sum_{i=1}^{N}\mathds {1}_{f(x_{i})\in O}\), where \(x_{1},...x_{N}\) are samples from \(\mathcal{C}\), while the estimate of \(\text{vol}(\mathcal{T})\) is updated in each iteration in Algorithm 1. We can set the target volume threshold (stopping criterion) to be \(v=r\times\text{vol}(\widehat{f_{\mathcal{C}}^{-1}}(O)\), where \(r\) is a _target coverage ratio_.

#### Iv-B2 RQ1. Effectiveness in Preimage Approximation

To investigate the efficiency and scalability of our approach, we perform comparison experiments with the exact preimage generation method [13]. We also evaluate the scalability of our approach w.r.t. network layer widths and depths.

Aircraft Collision AvoidanceIn our experiment, we use the following input region for the ownship and intruder aircraft as in [13]: \(h\in[-8000,8000]\), \(h_{A}\in[-100,100]\), \(h_{B}=30\), and \(t\in[0,40]\). We consider the output property \(O=\{y\in\mathbb{R}^{9}\mid\wedge_{i\in[2,9]}\;y_{1}\geq y_{i}\}\) and generate the preimage approximation for the nine neural networks VCAS 1 to 9. We set the target coverage ratio as 90% and the termination criterion - the iteration limit of our algorithm - as 500.

The experimental results are summarized in Table I. We compare our method with exact preimage generation, showing the number of polytopes (#Poly) in the under-approximation and exact preimage, respectively, and time in seconds (Time(s)). Column "PolyCov" shows the approximate coverage ratio of our approach when the algorithm terminates. For all VCAS networks, our approach effectively generates the preimage under-approximations with the polytope number varying from 6 to 18. Compared with the exact method, our approach realizes an average reduction of 91.1% (131 vs 12). Further, the computation time of our approach for all neural networks is less than 20s, demonstrating orders-of-magnitude improvement in efficiency (\(564\times\) faster on average).

Neural Network ControllersIn this experiment, we consider neural network controllers for the reinforcement learning tasks. The exact preimage generation method fails to scale to neural networks of these sizes due to exponential complexity. Table II summarizes the experimental results.

\begin{table}
\begin{tabular}{c|c c|c c} \hline \hline \multirow{2}{*}{**Models**} & \multicolumn{2}{c|}{**Exact**} & \multicolumn{2}{c}{**Our**} \\ \cline{2-5}  & \#**Poly** & **Time(s)** & **\#Poly** & **Time(s)** & **PolyCov(\%)** \\ \hline VCAS 1 & 49 & 6348.937 & 13 & 13.619 & 90.5 \\ VCAS 2 & 120 & 6325.712 & 11 & 10.431 & 90.4 \\ VCAS 3 & 253 & 6327.981 & 18 & 18.05 & 90.0 \\ VCAS 4 & 165 & 6435.46 & 11 & 10.384 & 90.4 \\ VCAS 5 & 122 & 6366.877 & 11 & 10.855 & 91.1 \\ VCAS 6 & 162 & 6382.198 & 10 & 9.202 & 92.0 \\ VCAS 7 & 62 & 6374.165 & 6 & 4.526 & 92.8 \\ VCAS 8 & 120 & 6341.173 & 14 & 13.677 & 91.3 \\ VCAS 9 & 125 & 6366.941 & 11 & 10.782 & 90.3 \\ \hline
**Average** & 131 & 6363.272 & 12 & 11.281 & 91.0 \\ \hline \hline \end{tabular}
\end{table} TABLE I: Performance of preimage generation for VCAS.

Fig. 3: Effect of different configurations on preimage approximation.
Firstly, we note that our approach can still effectively generate preimage under-approximations for the neural network controllers. Second, we see that, even with the same input region under analysis, there can be large differences in generation time and polytope number for different output sets to reach the target coverage ratio (e.g., 154 and 32 for Cartpole). This is because different output sets lead to differences in linear relaxation errors, as well as differences in the polytope tree/subregion splitting procedure.

Network SizeWe investigate the scalability of our approach by training neural networks of varying sizes for VCAS and evaluating the preimage approximation quality and runtime cost. Figure 2(a) and 2(b) depict the evaluation results w.r.t. layer width and layer depth. Regarding layer width, as the number of hidden nodes increases, more refinement iterations are needed to reach the target approximation quality, resulting in a greater number of polytopes and higher runtime cost. We observe that our approach is more sensitive to depth than width due to the accumulated convex relaxation errors. The approximation coverage decreases to 62.4% for 6 hidden layers when the algorithm reaches an iteration limit of 5000. Still, our method offers an advantage over the approximate method in [14], for which layer width (e.g., over 32) is computationally hard due to the exponential complexity.

Our approach can effectively generate approximations with orders-of-magnitude improvement in efficiency and scale beyond existing exact and approximate methods to neural network controllers for RL tasks.

#### Iv-B3 RQ2. Refinement Methods and Parameters

We examine the effectiveness of our global refinement and local optimization in improving the approximation quality, as well as the impact of parameter configurations.

Subregion SelectionAs a baseline for comparison, we conduct experiments with random region selection. The random strategy differs from our approach in selecting the next region to split randomly without prioritization. We perform comparison experiments on the benchmark tasks. Table III summarizes the evaluation results of random strategy (Column "Rand") and our method (Column "Our") for under-approximation refinement. We set the same target coverage ratio and iteration limit for both strategies. Note that, for _Dubinsrejoin_, random selection method hits the iteration limit and fails to reach the target coverage ratio. The results confirm the effectiveness of our region selection method in that _fewer_ iterations of the approximation refinement are required to reach the target coverage ratio, leading to (i) a smaller number of polytopes (#Polytope), reduced by 79.1% on average, and (ii) a 76.2% average runtime reduction.

Splitting FeatureWe compare our (greedy) splitting method with the heuristic splitting method, which chooses to split the selected subregion along the input index with the largest value range. We present the comparison results with the heuristic method (Column "Heuristic") in Table III. Our method requires splitting on all input features, computing the preimage approximations for all splits, and then choosing the dimension that refines the under-approximation the most. We find that, even with parallelization of the computation over input features, our approach leads to larger runtime overhead per-iteration compared with the heuristic method. Despite this, we find that our strategy actually requires _fewer_ refinement iterations to reach the target coverage, leading to a smaller number of polytopes (43.3% reduction on average) for the same approximation quality, demonstrating the per-iteration improvement in volume of the greedy vs heuristic strategy.

\begin{table}
\begin{tabular}{c|c c|c c|c c c} \hline \hline \multirow{2}{*}{**Task**} & \multirow{2}{*}{**Property**} & \multicolumn{3}{c|}{**Optim.**} & \multicolumn{3}{c}{**No Optim.**} \\ \cline{3-8}  & & **\#Poly** & **PolyCov(\%)** & **Time(s)** & **\#Poly** & **PolyCov(\%)** & **Time(s)** \\ \hline \multirow{2}{*}{Cartpole} & \(\{y\in\mathbb{R}^{2}|\ y_{1}\geq y_{2}\}\) & 154 & 75.2 & 98.839 & 5002 & 59.7 & 1139.853 \\  & \(\{y\in\mathbb{R}^{2}|\ y_{2}\geq y_{1}\}\) & 32 & 75.4 & 17.156 & 40 & 75.1 & 13.210 \\ \hline \multirow{2}{*}{Lunarlander} & \(\{y\in\mathbb{R}^{4}|\wedge_{i\in\{1,3,4\}}y_{2}\geq y_{i}\}\) & 130 & 75.0 & 154.882 & 5002 & 63.5 & 5023.098 \\ \hline \multirow{2}{*}{Dubinsrejoin} & \(\{y\in\mathbb{R}^{8}|\ \wedge_{i\in[2,4]}\ y_{1}\geq y_{i}\ \bigwedge_{i\in[6,8]}\ y_{5}\geq y_{i}\}\) & 78 & 75.2 & 102.673 & 5002 & 51.4 & 3985.861 \\ \hline \multicolumn{8}{c|}{**Average**} & 99 & 75.2 & 93.388 & 3762 & 62.4 & 2540.506 \\ \hline \hline \end{tabular}
\end{table} TABLE II: Performance of preimage generation for reinforcement learning tasks.

\begin{table}
\begin{tabular}{c c c|c c c|c c|c c} \hline \hline \multirow{2}{*}{**Model**} & \multicolumn{3}{c}{**\#Polytope**} & \multicolumn{3}{c}{**PolyCov(\%)**} & \multicolumn{3}{c}{**Time(s)**} \\ \cline{2-10}  & Rand & Heuristic & Our & Rand & Heuristic & Our & Rand & Heuristic & Our \\ \hline Vehicle Parking & 41 & 4 & 4 & 90.5 & 95.2 & 91.3 & 8.699 & 0.967 & 0.661 \\ \hline VCAS & 56 & 72 & 13 & 92.3 & 90.0 & 90.5 & 29.15 & 30.829 & 13.619 \\ \hline Cartpole & 151 & 34 & 32 & 75.3 & 75.3 & 75.4 & 84.639 & 12.313 & 17.156 \\ \hline Lunarlander & 481 & 238 & 130 & 75.1 & 75.1 & 75.0 & 505.747 & 119.744 & 154.882 \\ \hline Dubinsrejoin & 502 & 105 & 78 & 61.8 & 75.3 & 75.2 & 587.908 & 65.869 & 102.673 \\ \hline
**Average** & 246 & 91 & 51 & 79.0 & 82.2 & 81.5 & 243.229 & 45.944 & 57.798 \\ \hline \hline \end{tabular}
\end{table} TABLE III: Effectiveness of the refinement method.
Local OptimizationWe conduct an ablation analysis to assess the effectiveness of our local optimization method. We perform experiments on the RL tasks. We set the target coverage ratio as 75% and the iteration limit of our algorithm as 5000. Table II (Column "No Optim") summarizes the evaluation results when removing local optimization from the algorithm. The results reveal that removing local optimization leads to a substantial increase in polytope number (3670.9% on average) and runtime cost (2554.3%), while failing to reach the approximation quality when the algorithm terminates.

Parameter ConfigurationsWe investigate the impact of two parameters: the sample size \(N\) used in the Monte-Carlo volume estimates, and the targeted coverage for the tradeoff between approximation quality and efficiency.

Figure 2(c) shows the impact of sample size on the number of polytopes and time required to reach a fixed target coverage on the Cartpole task. It can be seen that both the number of polytopes and time required decrease as the sample size increases, up to a sample size of 10000. This reflects that a larger sample size leads to more accurate estimation of the volumes, improving the refinement procedure and resulting in fewer iterations to reach the target coverage. On the other hand, increasing the sample size further leads to diminishing returns in terms of estimation error, and an increase to per-iteration time when the volume estimation becomes the bottleneck. In our experiments, we use a sample size of 10000 to balance between the estimation accuracy and runtime efficiency.

The target coverage controls the preimage approximation quality of our approach. Figure 2(d) shows the relation between target coverage and the number of polytopes required (similar relation between target coverage and runtime). Overall, computing preimage approximation with a higher target coverage (better approximation quality) requires more refinement iterations on the input region, leading to more polytopes in the DUP approximation and higher runtime cost.

The proposed global refinement and local optimization methods are essential to improving preimage approximation quality and runtime efficiency. The sample size is a trade-off between the estimation accuracy and runtime overhead. The target coverage allows flexibility to control approximation precision based on available budget.

#### Vi-B4 RQ3. Quantitative Verification

Given a neural network, input set \(I\), and output set \(O\), we use our framework to perform quantitative verification of the property \((I,O,p)\); that is, to check whether \(f(x)\in O\) is true for a specified proportion \(p\) of input values \(x\in I\).

Vehicle Parking.We first perform quantitative verification on neural networks in the vehicle parking task. Consider the quantitative property with input set \(I=\{x\in\mathbb{R}^{2}\mid x\in[0,1]^{2}\}\), output set \(O=\{y\in\mathbb{R}^{4}|\bigwedge_{i=2}^{4}y_{1}-\ y_{i}\geq 0\}\), and quantitative proportion \(p=0.95\). We use Algorithm 3 to verify this property, with iteration limit 500. The computed under-approximation is a union of two polytopes, which takes 0.942s to reach the target coverage. We then compute the exact volume ratio of the under-approximation against the input region. The final quantitative proportion reached by our under-approximation is 95.2%, and thus verifies the quantitative property.

Aircraft Collision Avoidance.In this example, we consider the VCAS system and a scenario where the two aircrafts have negative relative altitude from intruder to ownship (\(h\in[-8000,0]\)), the ownship aircraft has a positive climbing rate \(h_{A}\in[0,100]\) and the intruder has a stable negative climbing rate \(h_{B}=-30\), and time to the loss of horizontal separation \(t\in[0,40]\), which defines the input region \(I\). For this scenario, it is clear that "COC" should be the advisory. Then we apply Algorithm 3 to verify the quantitative property where \(O=\{y\in\mathbb{R}^{9}|\bigwedge_{i=2}^{9}y_{1}-y_{i}\geq 0\}\) and the proportion \(p=0.9\), and we set an iteration limit of 500. The under-approximation computed is a union of 6 polytopes, which takes 5.620s to reach the target coverage. The exact quantitative proportion reached by the generated under-approximation is 90.8%, which thus verifies the quantitative property.

The proposed method can be effectively applied to provide provable guarantees for quantitative properties by enabling arithmetic computation over the input space.

## VII Threats to Validity

We discuss potential threats to the validity of our approach and possible ways to address them.

The _external_ threat mainly comes from the application domain. Our approach is designed for neural networks in the context of safety-critical control tasks, which have low input dimensionality. High-dimensional input spaces such as those found in computer vision tasks are challenging for the proposed method (and any preimage approximation method in general), as one can easily construct preimages that cannot be well-approximated by any set of polytopes. However, we believe that our approach can be adapted to investigate semantically meaningful input subspaces, e.g., a patch containing some visual attribute. To mitigate potential _internal_ threats, we perform a thorough analysis of our methodological choices and parameter configurations to evaluate their impact on the results.

## VIII Conclusion

We present an efficient and practical algorithm for preimage approximation of neural networks. Our algorithm is the first _anytime_ approach, with the output after each iteration being a valid under-approximation that is guaranteed to continually improve. This allows the user to freely configure the approximation according to a time budget or target coverage. As an application, we demonstrate verification of global quantitative properties of neural networks. Our evaluation on a range of benchmark tasks shows significant advantage in runtime efficiency compared with exact preimage generation. One promising direction for future work is to explore multi-neuron relaxation techniques to reduce accumulated errors through deep layers, further improving scalability.
## Acknowledgments

This project received funding from the ERC under the European Union's Horizon 2020 research and innovation programme (FUN2MODEL, grant agreement No. 834115) and ELSA: European Lighthouse on Secure and Safe AI project (grant agreement No. 101070617 under UK guarantee).
