Provable Preimage Under-Approximation for Neural Networks

. 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 show-case 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.


Introduction
Despite the remarkable empirical success of neural networks, guaranteeing their correctness, especially when using them as decision-making components in safetycritical autonomous systems [13,19,50], 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 [27,29,51,17,42,39,47,48,43].While local robustness verification focuses on deciding the absence of adversarial examples within an ϵ-perturbation neighbourhood, an alternative approach for neural network analysis is to construct the preimage of its predictions [34,21].Given a set of outputs, the preimage is defined as the set of all inputs mapped by the neural network to that output set.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, such as computing the proportion of inputs satisfying a property (quantitative verification) even if standard robustness verification fails.
Exact preimage generation [34] is intractable, taking time exponential in the number of neurons in a network; thus approximations are necessary.Unfortunately, existing methods are limited in their applicability.The inverse abstraction method in [21] bypasses the intractability of exact preimage generation by leveraging symbolic interpolants [20,6] 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.A concurrent work [30] proposed an input bounding algorithm targeting backward reachability analysis for control policies and out-of-distribution (OOD) detection in low-dimensional domains.Their method produces a preimage overapproximation, which cannot be used for quantitative verification.Therefore, more efficient and flexible computation methods for (symbolic abstraction of) preimages of neural networks are needed.
The main contribution of this paper is a scalable method for preimage approximation, which can be used for a variety of robustness analysis tasks.More specifically, we propose an efficient anytime algorithm for generating symbolic under-approximations of the preimage of piecewise linear neural networks as a union of disjoint polytopes.The algorithm computes a sound preimage underapproximation leveraging linear relaxation based perturbation analysis (LiRPA) [47,48,39], applied backwards from a polyhedron output set.It iteratively refines the preimage approximation by adding input and/or intermediate (ReLU) splitting (hyper)planes to partition the input region into disjoint subregions, which can be approximated independently in parallel in a divide-and-conquer approach.The refinement scheme uses a novel differential objective to optimize the quality (volume) of the polytope subregions.We also show that our method can be generalized to generate preimage over-approximations.We illustrate the application of our method to quantitative verification, input bounding for control tasks, and robustness analysis against adversarial and patch attacks.Finally, we conduct an empirical analysis on a range of control and computer vision tasks, showing significant gains in efficiency compared to exact preimage generation methods and scalability to high-input-dimensional tasks compared to existing preimage approximation methods.
Proofs and additional technical details are presented in the Appendix.

Preliminaries
We use f : R d → R m to denote a feedforward neural network.For layer i, we use W (i) to denote the weight matrix, b (i) the bias, h (i) the pre-activation neurons, and a (i) the post-activation neurons, such that we have h (i) = W (i) a (i−1) + b (i) .
In this paper, we focus on ReLU neural networks with a (i) (x) = ReLU (h (i) (x)), where ReLU(h) := max(h, 0) is applied element-wise.However, our method can be generalized to other activation functions bounded by linear relaxation [51].Linear Relaxation of Neural Networks.Nonlinear activation functions lead to the NP-completeness of the neural network verification problem [29].To address such intractability, linear relaxation is often used to transform the nonconvex constraints into linear programs.As shown in Figure 1, given concrete lower and upper bounds l (i) ≤ h (i) (x) ≤ u (i) on the pre-activation values of layer i, there are three cases to consider.In the inactive (u (i) j ≤ 0) and active (l (i) j ≥ 0) cases, the post-activation neurons a (i) j (x) are linear functions a (i) j (x) = 0 and a (i) j (x) = h (i) j (x) respectively.In the unstable case, a (i) j (x) can be bounded by α , where α (i) j is a configurable parameter that produces a valid lower bound for any value in [0, 1].Linear bounds can also be obtained for other non-piecewise linear activation functions [51].
Linear relaxation can be used to compute linear upper and lower bounds of the form Ax + b ≤ f (x) ≤ Ax + b on the output of a neural network, for a given bounded input region C.These methods are known as linear relaxation based perturbation analysis (LiRPA) algorithms [47,48,39].In particular, backward-mode LiRPA computes linear bounds on f by propagating linear bounding functions backward from the output, layer-by-layer, to the input layer.
Polytope Representations.Given an Euclidean space R d , a polyhedron T is defined to be the intersection of a set of half spaces.More formally, suppose we have a set of linear constraints defined by , where T consists of all values of x satisfying the first-order logic (FOL) formula α(x) := K i=1 ψ i (x).We use the term polytope to refer to a bounded polyhedron, that is, a polyhedron T such that ∃R ∈ R >0 : ∀x 1 , x 2 ∈ T ,∥x 1 − x 2 ∥ 2 ≤ R holds.The abstract domain of polyhedra [39,12,14] has been widely used for the verification of neural networks and computer programs.An important type of polytope is the hyperrectangle (box), which is a polytope defined by a closed and bounded interval [x i , x i ] for each dimension, where x i , x i ∈ Q.More formally, using the linear constraints 3 Problem Formulation

Preimage Approximation
In this work, we are interested in the problem of computing preimages for neural networks.Given a subset O ⊂ R m of the codomain, the preimage of a function f : R d → R m is defined to be the set of all inputs x ∈ 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 C ⊂ R d .In this work, we restrict the output set O to be a polyhedron, and the input set C to be an axis-aligned hyperrectangle region C ⊂ R d , as these are commonly used in neural network verification.We now define the notion of a restricted preimage: Example 1.To illustrate our problem formulation and approach, we introduce a vehicle parking task [8] as a running example.In this task, there are four parking lots, located in each quadrant of a 2 × 2 grid [0, 2] 2 , and a neural network with two hidden layers of 10 ReLU neurons f : R 2 → 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] × [0, 1] corresponding to parking lot 1, we set that is labelled as parking lot 1 by the network.
We focus on provable approximations of the preimage.Given a first-order formula A, α is an under-approximation (resp.over-approximation) of A if it holds that ∀x.α(x) =⇒ A(x) (resp.∀x.A(x) =⇒ α(x)).In our context, the restricted preimage is defined by the formula A(x) = (f (x) ∈ O) ∧ (x ∈ C), and we restrict to approximations α that take the form of a disjoint union of polytopes (DUP).The goal of our method is to generate a DUP approximation T that is as tight as possible; that is, to maximize the volume vol(T ) of an under-approximation, or minimize the volume vol(T ) of an over-approximation.

Definition 2 (Disjoint Union of Polytopes
, where each α i is a polytope formula (conjunction of a finite set of linear half-space constraints), with the property that α i ∧ α j is unsatisfiable for any i ̸ = j.

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 [24,25,37].This is often expressed as a property of the form (I, O) such that ∀x ∈ I =⇒ f (x) ∈ O.We can generalize this to quantitative properties: Definition 3 (Quantitative Property).Given a neural network f : R d → R m , a measurable input set with non-zero measure (volume) I ⊆ R d , a measurable output set O ⊆ R m , and a rational proportion p ∈ [0, 1] we say that the neural network satisfies the property (I, O, p) if Neural network verification algorithms [32] can be divided into two categories: sound, which always return correct results, and complete, guaranteed to reach a conclusion on any verification query.We now define soundness and completness of verification algorithms for quantitative properties.

Definition 4 (Soundness).
A verification algorithm QV is sound if, whenever QV outputs True, the property (I, O, p) holds.

Definition 5 (Completeness).
A verification algorithm QV is complete if (i) QV never returns Unknown, and (ii) whenever QV outputs False, the property (I, O, p) does not hold.
If the property (I, O) holds, then the quantitative property (I, O, 1) holds, while quantitative properties for 0 ≤ p < 1 provide more information when (I, O) does not hold.Most neural network verification methods produce approximations of the image of I in the output space, which cannot be used to verify quantitative properties.Preimage over-approximations include false regions, thereby not applicable for quantitative verification.In contrast, preimage under-approximations provide a lower bound on the volume of the preimage, allowing us to soundly verify quantitative properties.

Methodology
Overview.In this section we present the main components of our methodology.Firstly, in Section 4.1, we show how to cheaply and soundly under-approximate the (restricted) preimage with a single polytope, using linear relaxation methods (Algorithm 2).Secondly, in Section 4.2, we propose a novel differentiable objective to optimize the quality (volume) of the polytope under-approximation.Thirdly, in Section 4.3, we propose a refinement scheme that improves the approximation by partitioning a (sub)region into subregions with splitting planes, with each subregion then being under-approximated more accurately.The main contribution of this paper (Algorithm 1) integrates these three components and is described in Section 4.4.Finally, in Section 4.5, we apply our method to quantitative verification (Algorithm 3) and prove its soundness and completeness.

Polytope Under-Approximation via Linear Relaxation
We first show how to adapt linear relaxation techniques to efficiently generate valid under-approximations to the restricted preimage for a given input region C.

Algorithm 1: Preimage Approximation
Recall that LiRPA methods enable us to obtain linear lower and upper bounds on the output of a neural network f , that is, Ax + b ≤ f (x) ≤ Ax + b, where the linear coefficients depend on the input region C. Now, suppose that we are interested in computing an under-approximation to the restricted preimage, given the input hyperrectangle C = {x ∈ R d |x |= d i=1 ϕ i }, and the output polytope specified using the half-space constraints ψ i (y) = (c T i y + d i ≥ 0) for i = 1, ..., K over the output space.Given a constraint ψ i , we append an additional linear layer at the end of the network f , which maps y → c T i y +d i , such that the function Then, applying LiRPA bounding to each g i , we obtain lower bounds g i (x) = a T i x + b i for each i, such that g i (x) ≥ 0 =⇒ g i (x) ≥ 0 for x ∈ C. Notice that, for each i = 1, ..., K, a T i x+b i ≥ 0 is a half-space constraint in the input space.We conjoin these constraints, along with the restriction to the input region C, to obtain a polytope , where c i := e 1 − e i (where e i is the i th standard basis vector) and d i := 0. Applying LiRPA bounding, we obtain the linear lower bounds The intersection of these constraints, shown in Figure 2a, represents the region where any input is guaranteed to satisfy the output constraints.
We generate the linear bounds in parallel over the output polyhedron constraints i = 1, ..., K using the backward mode LiRPA [51], and store the resulting input polytope T C (O) as a list of constraints.This highly efficient procedure is used as a sub-routine LinearLowerBound when generating a preimage underapproximation as a polytope union using Algorithm 2 (Line 4).

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 α, 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 C [47], i.e., min x∈C A(α)x+b(α), where we explicitly note the dependence of the linear coefficients on α.In our case, we are instead interested in optimizing α to refine the polytope under-approximation, that is, increase its volume.Unfortunately, computing the volume of a polytope exactly is a computationally expensive task, and requires specialized tools [18] that do not permit easy optimization with respect to the α parameters.
To address this challenge, we propose to use statistical estimation.In particular, we sample N points x 1 , ..., x N uniformly from the input domain C then employ Monte Carlo estimation for the volume of the polytope approximation: where we highlight the dependence of , and α i are the α-parameters for the linear relaxation of the neural network g i corresponding to the i th half-space constraint in O.However, this is still non-differentiable w.r.t.α due to the identity function.We now show how to derive a differentiable relaxation which is amenable to gradient-based optimization: The second equality follows from the definition of the polytope T C,α (O); namely that a point is in the polytope if it satisfies g i (x j , α i ) ≥ 0 for all i = 1, ..., K, or equivalently, min i=1,...K g i (x j , α i ) ≥ 0. After this, we approximate the identity function using a sigmoid relaxation, where σ(y) := 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 logsum-exp function is defined by LSE(y 1 , ..., y K ) := log( i=1,...,K e yi ), 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 α.We employ this as the loss function in Algorithm 2 (Line 6) for generating a polytope approximation, and optimize volume using projected gradient descent.
Example 3. We revisit the vehicle parking problem in Example 1. Figure 2a and 2b show the computed under-approximations before and after local optimization.We can see that the bounding planes for all three specifications are optimized, which effectively improves the approximation quality.

Global Branching and Refinement
As LiRPA performs crude linear relaxation, the resulting bounds can be quite loose even with α-optimization, meaning that the polytope approximation T C (O) is unlikely to constitute a tight under-approximation to the preimage.To address this challenge, we employ a divide-and-conquer approach that iteratively refines our under-approximation of the preimage.Starting from the initial region C represented at the root, our method generates a tree by iteratively partitioning a subregion C sub represented at a leaf node into two smaller subregions C l sub , C u sub , 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 C.
For each leaf subregion C sub we compute, using LiRPA bounds (Line 4, Algorithm 2), an associated polytope that under-approximates the preimage in C sub .Thus, irrespective of the number of refinements performed, the union of the polytopes corresponding to all leaves forms an anytime DUP under-approximation T to the preimage in the original region C.The process of refining the subregions continues until an appropriate termination criterion is met.
Unfortunately, even with a moderate number of input dimensions or unstable ReLU nodes, naïvely splitting along all input-or ReLU-planes 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.Another important aspect is how to prioritize which leaf subregion to split.We describe these in turn.
Subregion Selection.Searching through all leaf subregions at each iteration is computationally too expensive.Thus, we propose a subregion selection strategy that prioritizes splitting subregions according to (an estimate of) the difference in volume between the exact preimage f −1 C sub (O) and the (already computed) polytope approximation T C sub (O) on that subdomain, that is: which measures the gap between the polytope under-approximation and the optimal approximation, namely, the preimage itself.Suppose that a particular leaf subdomain attains the maximum of this metric among all leaves, and we partition it into two subregions C l sub , C u sub , which we approximate with polytopes T C l sub (O), T C u sub (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 C sub with polytope approximation T C sub (O), and its children C l sub , C u sub with polytope approximations T C l sub (O), T C u sub (O) respectively, it holds that: Corollary 1.In each refinement iteration, the volume of the polytope approximation T Dom does not decrease.
Since computing the volumes in Equation 2 is intractable, we sample N points x 1 , ..., x N uniformly from the subdomain C sub and employ Monte Carlo estimation to estimate the volume for both the preimage and the polytope approximation using the same set of samples, i.e., vol(f  stress that volume estimation is only used to prioritize subregion selection, and does not affect the soundness of our method. Input Splitting.Given a subregion (hyperrectangle) defined by lower and upper bounds x i ∈ [x i , x i ] for all dimensions i = 1, ..., d, input splitting partitions it into two subregions by cutting along some feature i.This splitting procedure will produce two subregions which are similar to the original subregion, but have updated bounds [x i , for feature i instead.In order to determine which feature/dimension to split on, we propose a greedy strategy.Specifically, for each feature, we generate a pair of polytopes 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 [16], that is, to select feature i with the largest range: arg max i (x i −x i ).However, this method falls short in periteration approximation volume improvement compared to our greedy strategy.
Example 4. We revisit the vehicle parking problem in Example 1. Figure 2b shows the polytope under-approximation computed on the input region C before refinement, where each solid line represents the bounding plane for each output specification (y 1 − y i ≥ 0). Figure 2c depicts the refined approximation by splitting the input region along the vertical axis, where the solid and dashed lines represent the bounding planes for the two resulting subregions.It can be seen that the total volume of the under-approximation has improved significantly.
Intermediate ReLU Splitting.Refinement through splitting on input features is adequate for low-dimensional input problems such as reinforcement learning agents.However, it may be infeasible to generate sufficiently fine subregions for high-dimensional domains.We thus propose an algorithm for ReLU neural networks that uses intermediate ReLU splitting for preimage refinement.After determining a subregion for refinement, we partition the subregion based upon the pre-activation value of an intermediate unstable neuron z In this procedure, the order of splitting unstable ReLU neurons can greatly influence the refinement quality and efficiency.Existing heuristic methods of ReLU prioritization select ReLU nodes that lead to greater improvement in the final bound (maximum or minimum value) of the neuron network on the input domain [16], i.e. min x∈C f (x).However, these ReLU prioritization methods are not effective for preimage analysis, because our objective is instead to refine the overall preimage approximation.We thus propose a heuristic method to prioritize unstable ReLU nodes for preimage refinement.Specifically, we compute (an estimate of) the volume difference between the split subregions |vol(C + )|, using a single forward pass for a set of sampled datapoints from the input domain; note that this is bounded above by the total subregion volume vol(C sub ).We then propose to select the ReLU node that minimizes this difference.Intuitively, this choice results in balanced subdomains after splitting.
Another advantage of ReLU splitting is that we can replace the unstable neuron bound ch and a (i) j (x) = 0, respectively, as shown in Figure 1 (unstable to stable).This can then tighten the linear bounds for the other neurons, thus tightening the under-approximation on each subdomain.
Example 5. We now apply our algorithm with ReLU splitting to the vehicle parking problem in Example 1. Figure 2d shows the refined preimage polytope by adding the splitting plane (black solid line) along the direction of a selected unstable ReLU node.Compared with Figure 2b, we can see that the volume of the approximation is improved.
Remark 1 (Preimage Over-approximation).While Algorithms 1 and 2 focus on preimage under-approximations, they can be easily configured to generate overapproximations with two key modifications.Firstly, we generate polytope overapproximations by using LiRPA to propagate a linear upper bound g Secondly, the refinement and optimization objective is to minimize the volume of the over-approximation instead of maximizing the volume as in the case of under-approximation.

Overall Algorithm
Our overall preimage approximation method is summarized in Algorithm 1.It takes as input a neural network f , input region C, output region O, target polytope volume threshold v (a proxy for approximation precision), termination iteration number R, and a Boolean indicating whether to use input or ReLU splitting, and returns a disjoint polytope union T representing an underapproximation to the preimage.The algorithm initiates and maintains a priority queue of (sub)regions according to Equation 2. The initialization step (Lines 1-3) generates an initial polytope approximation on the whole region using Algorithm 2 (Sections 4.1, 2).Then, the preimage refinement loop (Lines 4-14) partitions a subregion in each iteration, with the preimage restricted to the child subregions then being re-approximated (Line 10-11).In each iteration, we choose the region to split (Line 5) and the splitting plane to cut on (Line 7 for input split and Line 9 for ReLU split), as detailed in Section 4.3.The preimage under-approximation is then updated by computing the priorities for each subregion (by approximating volumes) (Lines 12-14).The loop terminates and the approximation returned when the target volume threshold v or maximum iteration limit R is reached.

Quantitative Verification
We now show how to use our efficient preimage under-approximation method (Algorithm 1) to verify a given quantitative property (I, O, p), where O is a polyhedron, I a polytope and p the desired proportion value, summarized in Algorithm 3. To simplify assume that I is a hyperrectangle, so that we can take C = I (the case of general polytopes is discussed in the Appendix).We utilize Algorithm 1 by setting the volume threshold to p × vol(I), such that we have vol(T ) vol(I) ≥ p if the algorithm terminates before reaching the maximum number of iterations.However, the Monte Carlo estimates of volume cannot provide a sound guarantee that vol(T ) vol(I) ≥ p.To resolve this problem, we propose to run exact volume computation [11] only when the Monte Carlo estimate reaches the threshold.If the exact volume vol(T ) ≥ p × vol(I), then the property is verified.Otherwise, we continue running the preimage refinement.
In Algorithm 3, InitialRun generates an initial approximation to the preimage as in Lines 1-3 of Algorithm 1, and Refine performs one iteration of approximation refinement (Lines 5-14).Termination occurs when we have verified or falsified the quantitative property, or when the maximum number of iterations has been exceeded.Proposition 3. Algorithm 3 is sound for quantitative verification with input splitting.Proposition 4. Algorithm 3 is sound and complete for quantitative verification on piecewise linear neural networks with ReLU splitting.

Experiments
We have implemented our approach as a prototype tool3 for preimage approximation for polyhedron-type output sets/specifications.In this section, we perform experimental evaluation of the proposed approach on a set of benchmark tasks and demonstrate its effectiveness in approximation generation and its application to quantitative analysis of neural networks.

Benchmark and Evaluation Metric
We evaluate our preimage analysis approach on a benchmark of reinforcement learning and image classification tasks.Besides the vehicle parking task [8] shown in the running example, we use the following (trained) benchmarks: (1) aircraft collision avoidance system (VCAS) [28] with 9 feed-forward neural networks (FNNs); (2) neural network controllers from VNN-COMP 2022 [5] for three reinforcement learning tasks (Cartpole, Lunarlander, and Dubinsrejoin) [15]; and (3) the neural network from VNN-COMP 2022 for MNIST classification.Details of the models and additional experiments can be found in the Appendix.
Evaluation Metric To evaluate the quality of the preimage approximation, we define the coverage ratio to be the ratio of volume covered to the volume of the exact preimage, i.e., cov(T , f . Note that this is a normalized measure for assessing the quality of the approximation, as shown in Algorithm 3 when comparing with target coverage proportion p for termination of the refinement loop, and not as a measure for formal verification guarantees.In practice, we estimate vol(f , where r is the target coverage ratio.

Evaluation
Effectiveness in Preimage Approximation with Input Split We apply Algorithm 1 with input splitting to the input bounding problem for low-dimensional reinforcement learning tasks to evaluate its effectiveness.For comparison, we also run the exact preimage (Exact) [34] and preimage over-approximation (Invprop) [30,31]  in the preimage, computation time (Time(s)), and the approximate coverage ratio (Cov(%)) when the preimage approximation algorithm terminates with target coverage 90%.Compared with the exact method, our approach yields orders-ofmagnitude improvement in efficiency.It can also characterize the preimage with much fewer (and also disjoint) polytopes (average reduction of 91.1% for VCAS).The Invprop method [30] cannot be directly applied as it computes preimage over-approximations.We adapt it to produce an under-approximation by computing over-approximations for the complement of each output constraint; the resulting approximation is then the complement of a union of polytopes, rather than a DUP.On the 2D vehicle parking task, we find that the results (see Table 1) are comparable with ours in time and approximation coverage.Their implementation currently only supports two-dimensional input tasks [31].While their algorithm, which employs input splitting, can in theory be extended to higher-dimensional tasks, a significant unaddressed technical challenge is in how to choose the input splits effectively in high dimensions.This is confounded by the fact that, to generate an under-approximation, we need separate runs of their algorithm for each output constraint.In contrast, our method naturally incorporates a principled splitting and refinement strategy, and can also effectively employ ReLU splitting for further scalability, as we will show below.Our method can also be configured to generate over-approximations (Section 4.3, Remark 1).
Neural Network Controllers.In this experiment, we consider preimage underapproximation for neural network controllers in reinforcement learning tasks.Note that [34] (Exact) is unable to deal with neural networks of these sizes and To our knowledge, this is the first time preimage computation has been attempted for this challenging, high-dimensional task.Table 3 summarizes the evaluation results for two types of image attacks: l ∞ and patch attack.For L ∞ attacks, bounded perturbation noise is applied to all image pixels.The patch attack applies only to a smaller patch area but allows arbitrary perturbations covering the whole valid range [0, 1].The task is then to produce a DUP under-approximation of the perturbation region that is guaranteed to be classified correctly.For L ∞ attack, our approach generates a preimage approximation that achieves the targeted coverage of 75% for noise up to 0.08.Notice that, from e.g.0.05 to 0.07, the volume of the input region increases by tens of orders of magnitude due to the high dimensionality.The fact that the number of polytopes and computation time remains manageable is due to the effectiveness of ReLU splitting.Interestingly, for the patch attack, we observe that the number of polytopes required increases sharply when increasing the patch size at the center of the image, while this is not the case for patches in the corners of the image.We hypothesize this is due to the greater influence of central pixels on the neural network output, and correspondingly a greater number of unstable neurons over the input perturbation space.Comparison with Robustness Verifiers We now illustrate empirically the utility of preimage computation in robustness analysis compared to robustness verifiers.Table 4 shows comparison results with α, β-CROWN, winner of the VNN competition [5].We set the tasks according to the problem instances from VNN-COMP 2022 for local robustness verification (localized perturbation regions).For Cartpole, α, β-CROWN can provide a verification guarantee (yes/no or safe/unsafe) for both of the problem instances.However, in the case where the robustness property does not hold, our method explicitly generates a preimage approximation in the form of a disjoint polytope union (where correct classi- fication is guaranteed), and covers 94.9% of the exact preimage.For MNIST, while the smaller perturbation region is successfully verified, α, β-CROWN with tightened intermediate bounds by MIP solvers returns unknown with a timeout of 300s for the larger region.In comparison, our algorithm provides a concrete union of polytopes where the input is guaranteed to be correctly classified, which we find covers 100% of the input region (up to sampling error).Note also (Table 3) that our algorithm can produce non-trivial under-approximations for input regions far larger than α, β-CROWN can verify.
Quantitative Verification We now demonstrate the application of our preimage generation framework to quantitative verification of the property (I, O, p); that is, to check whether f (x) ∈ O for at least proportion p of input values x ∈ I.This leverages the disjointness of our approximation, such that we can exactly compute the volume covered by exactly computing the volume of each polytope.Vehicle Parking.We consider the quantitative property with input set , and quantitative proportion p = 0.95.We use Algorithm 3 to verify this property, with iteration limit 1000.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%, verifying the quantitative property.
Aircraft Collision Avoidance.In this example, we consider the VCAS system and a scenario where the two aircraft have negative relative altitude from intruder to ownship (h ∈ [−8000, 0]), the ownship aircraft has a positive climbing rate ḣA ∈ [0, 100] and the intruder has a stable negative climbing rate ḣ B = −30, and time to the loss of horizontal separation is t ∈ [0, 40], which defines the input region I.For this scenario, the correct advisory is "Clear Of Conflict" (COC).We apply Algorithm 3 to verify the quantitative property where O = {y ∈ R 9 | 9 i=2 y 1 − y i ≥ 0} and the proportion p = 0.9, with an iteration limit of 1000.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 verifies the quantitative property.

Related Work
Our paper is related to a series of works on robustness verification of neural networks.To address the scalability issues with complete verifiers [27,29,42] based on constraint solving, convex relaxation [38] has been used for developing highly efficient incomplete verification methods [51,46,39,47].Later works employed the branch-and-bound (BaB) framework [17,16] to achieve completeness, using incomplete methods for the bounding procedure [48,43,23].In this work, we adapt convex relaxation for efficient preimage approximation.Further, our divide-andconquer procedure is analogous to BaB, but focuses on maximizing covered volume rather than maximizing a function value.There are also works that have sought to define a weaker notion of local robustness known as statistical robustness [44,33], which requires that a proportion of points under some perturbation distribution around an input point are classified in the same way.Verification of statistical robustness is typically achieved by sampling and statistical guarantees [44,10,41,49].In this paper, we apply our symbolic approximation approach to quantitative analysis of neural networks, while providing exact quantitative rather than statistical guarantees [45].
Another line of related works considers deriving exact or approximate abstractions of neural networks, which are applied for explanation [40], verification [22,36], reachability analysis [35], and preimage approximation [21,30].[21] leverages symbolic interpolants [6] for preimage approximations, facing exponential complexity in the number of hidden neurons.Concurrently, [30] employs Lagrangian dual optimization for preimage over-approximations.Our anytime algorithm, which combines convex relaxation with principled splitting strategies for refinement, is applicable for both under-and over-approximations.Their work may benefit from our splitting strategies to scale to higher dimensions.

Conclusion
We present an efficient and flexible algorithm for preimage under-approximation of neural networks.Our anytime method derives from the observation that linear relaxation can be used to efficiently produce under-approximations, in conjunction with custom-designed strategies for iteratively decomposing the problem to rapidly improve the approximation quality.Unlike previous approaches, it is designed for, and scales to, both low and high-dimensional problems.Experimental evaluation on a range of benchmark tasks shows significant advantage in runtime efficiency and scalability, and the utility of our method for important applications in quantitative verification and robustness analysis.time to the loss of horizontal separation.VCAS is implemented by nine feedforward neural networks built with a hidden layer of 21 neurons and an output layer of nine neurons corresponding to different advisories.
The experimental results are summarized in Table 6.We compare our method with exact preimage generation, showing the number of polytopes (#Polytope) 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, and 564× faster than the exact method on average.

B.3 Neural Network Controllers & MNIST Classification
In this subsection, we summarize the detailed configuration for neural network controllers in reinforcement learning tasks and MNIST classification tasks.
Cartpole.The cartpole control problem considers balancing a pole atop a cart by controlling the movement of the cart.The neural network controller has two hidden layers with 64 neurons, and uses four input variables representing the position and velocity of the cart, the angle and angular velocity of the pole.The controller outputs are pushing the cart left or right.In the experiments, we set the following input region for the Cartpole task: (1) cart position [0, 1], (2) cart velocity [0, 2], (3) angle of the pole [−0.2, 0], and (4) angular velocity of the pole [−2, 0] (with varied feature length in the evaluation).We consider the output property for the action pushing left.
Dubinsrejoin.The Dubinsrejoin problem considers guiding a wingman craft to a certain radius around a lead aircraft.The neural network controller has two hidden layers with 256 neurons.The input space of the neural network controller is eight dimensional, with the input variables capturing the position, heading, velocity of the lead and wingman crafts, respectively.The outputs are also eight dimensional representing controlling actions of the wingman.Note that the eight network outputs are processed further as tuples of actuators (rudder, throttle) for controlling the wingman where each actuator has 4 options.The control action tuple is decided by taking the action with the maximum output value among the first four network outputs (the first actuator options) and the action with the maximum value among the second four network outputs (the second actuator options).In the experiments, we set the following input region: (1) horizontal and vertical position [−0.2, 0] × [0, 0.5], (2) heading and velocity [−1, 0] × [0, 0.2] for the lead aircraft, and (3) horizontal and vertical position [0.4,0.6] × [−0.3, 0.3] (with varied feature length for evaluation), (4) heading and velocity [0.2, 0.5] × [−1.5, 0] for the wingman aircraft.We consider the output property that both actuators (rudder, throttle) take the first option, i.e., {y ∈ R 8 | ∧ i∈{2,3,4} y 1 ≥ y i ∧ i∈{6,7,8} y 5 ≥ y i }.
MNIST Classification.We use the trained neural network from VNN-COMP 2022 [5] for digit image classification.The neural network has six layers with a hidden neuron size of 100 for each hidden layer.We consider two types of image attacks: l ∞ and patch attack.For L ∞ attack, a perturbation is applied to all pixels of the image.For the patch attack, it applies arbitrary perturbations to the patch area, i.e., the perturbation noise covers the whole valid range [0, 1], for which we set the patch area at the center and (upper-left) corner of the image with different sizes.

C Effect of Methodological Components and Parameter Configurations
In this section, we perform an ablation study on two methodological components, subregion prioritization and input feature selection, and investigate the impact of two parameters in our approach, the sample suite size N for approximation quality evaluation and the targeted coverage ratio for the tradeoff between approximation quality and efficiency.

C.1 Ablation Study
In this subsection, we examine the effectiveness of our subregion and input feature selection methods on the approximation quality.
Subregion Selection.To evaluate the impact of our region selection method, we conduct comparison 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 7 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 (1) a smaller number of polytopes (#Polytope), reduced by 78.7% on average, and (2) a 75.5% average runtime reduction.
Splitting Feature Selection.We 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 of our method with the heuristic method (Column "Heuristic") in Table 7.
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, and 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 (42.2% reduction on average) for the same approximation quality, demonstrating the per-iteration improvement in volume of the greedy vs heuristic strategy.

C.2 Effect of Parameter Configurations
Sample Suite Size.In our approach, the sample suite size N is used to estimate the volume of the polytope approximations and restricted preimage, which is then used to select the subregion to split and the feature to split on, as well as in the α-optimization procedure.Figure 3 shows the evaluation results of the impact of sample suite size on preimage approximation for Cartpole.As demonstrated in Figure 3, the number of polytopes (#Polytope) that comprises the under-approximation decreases as the sample suite size increases until it reaches a fairly stable number (with minor fluctuation) with a sample size of 10000, with diminishing returns beyond that point.This reflects that a smaller sample suite size leads to greater estimation error of the volumes which negatively impacts both the refinement procedure and the α-optimization, resulting in more iterations needed to reach the target coverage.On the other hand, a larger sample size provides more accurate estimation of the approximation quality, but requires larger runtime overhead, as reflected in the increase of runtime cost when the sample suite size continues to increase.In our experiments, we use a sample suite size of 10000 to balance between the estimation accuracy and runtime efficiency.
Target Coverage.The target coverage controls the preimage approximation quality of our approach.Figure 4 shows the evaluation results of the impact of target coverage ratio on the polytope number and runtime cost.Overall, computing preimage approximation with larger target coverage ratio (better approximation quality) requires more refinement iterations on the input region, thus leading to more polytopes in the DUP approximation and larger runtime cost.
Our approach tries to minimize the polytope number to reach the target coverage ratio by prioritizing refinement on the input subregion with the lowest polytope coverage and selecting the splitting feature that improves the approximation quality the most.

D Quantitative Verification for General Input Polytopes
In Section 4.5, we detailed how to use our preimage under-approximation method to verify quantitative properties (I, O, p), where I is a hyperrectangle.We now show how to extend this to a general polytope I. Firstly, in Line 2 of Algorithm 3, we derive a hyperrectangle C such that I ⊆ C, by converting the polytope I into its V-representation [26], which is a list of the vertices (extreme points) of the polytope, which can be computed as in [7,11].The computational expense of this operation is no greater than that of polytope exact volume computation.Once we have a V-representation, obtaining a bounding box can be achieved simply by computing the minimum and maximum value x i , x i of each dimension among all vertices.
Once we have the input region C, we can then run the preimage refinement as usual, but with the modification that, when defining the polytopes and restricted preimages, we must additionally add the polytope constraints from I. Practically, this means that during every call to EstimateVolume and ExactVolume in Algorithm 3, we add these polytope constraints, and in Line 8 of Algorithm 2, we add the polytope constraints from I, in addition to those derived from the output O and the box constraints from C sub .

E Proofs
Proposition 1. T C (O) is an under-approximation to the restricted preimage f −1 C (O). Proof.The approximating polytope is defined as: The LiRPA bound g i (x) ≤ g i (x) holds for any x ∈ C, and so we have Proof.We define T C sub (O)| l , T C sub (O)| r to be the restrictions of T C sub (O) to C l sub and C r sub respectively, that is: where we have replaced the constraint x ∈ C sub with x ∈ C l sub (resp.x ∈ C r sub ), and g i (x) is the LiRPA lower bound for the i th specification on the input region C sub .
On the other hand, we also have: where g l,i (x) (resp.g r,i (x)) is the LiRPA lower bound for the i th specification on the input region C l sub (resp.C r sub ).Now, it is sufficient to show that | r is entirely similar).Before proving this result in full, we outline the approach and a sketch proof.It suffices to prove (for all i) that g l,i (x) is a tighter bound than g i (x) on C l sub .That is, to show that g l,i (x) ≥ g i (x) for inputs x in C l sub , as then g i (x) ≥ 0 =⇒ g l,i (x) ≥ 0 for inputs x in C l sub , and so ) is tighter than g i (x) because the input region for LiRPA is smaller for g l,i (x), leading to tighter concrete neuron bounds, and thus tighter bound propagation through each layer of the neural network g i .We present the formal proof of greater bound tightness for input and ReLU splitting in the following.
Input split: We show g l,i (x) ≥ g i (x) for all x ∈ C l sub by induction (dropping the index i in the following as it is not important).Recall that LiRPA generates symbolic upper and lower bounds on the pre-activation values of each layer in terms of the input (i.e.treating that layer as output), which can then be converted into concrete bounds.
where h (j) (x) are the pre-activation values for the j th layer of the network g i , and A (j) , b (j) , A (j) , b (j) (resp.A (l,j) , b (l,j) , A (l,j) , b (l,j) ) are the linear bound coefficients, for input regions C sub (resp.C l sub ).Inductive Hypothesis For all layers j = 1, ..., L in the network, and for all x ∈ C l sub , it holds that: Base Case For the input layer, we have the trivial bounds Ix ≤ x ≤ Ix for both regions.
Inductive Step Suppose that the inductive hypothesis is true for layer j − 1 < L. Using the symbolic bounds in Equations 9, 10, we can derive concrete bounds l (j−1) ≤ h (j−1) (x) ≤ u (j−1) and l (l,j−1) ≤ h (l,j−1) (x) ≤ u (l,j−1) on the values of the pre-activation layer.By the inductive hypothesis, the bounds for region C l sub will be tighter, i.e. l (j−1) ≤ l (l,j−1) ≤ u (l,j−1) ≤ u (j−1) .Now, consider the backward bounding procedure for layer j as output.We begin by encoding the linear layer from post-activation layer j − 1 to pre-activation layer j as: Then, we bound a (j−1) (x) in terms of h (j−1) (x) using linear relaxation.Consider the three cases in Figure 5 (reproduced from main paper), where we have a bound ch for some scalars c, d, c, d.If the concrete bounds (horizontal axis) are tightened, then an unstable neuron may become inactive or active, but not vice versa.It can thus be seen that the new linear upper and lower bounds on h (j−1) k (x) will also be tighter.Substituting the linear relaxation bounds in Equation 12as in [48], we obtain bounds of the form A (j) j h (j−1) (x) + b (j) j ≤ h (j) (x) ≤ A (j) j h (j−1) (x) + b (j) j (13) A (l,j) j h (j−1) (x) + b (l,j) j ≤ h (j) (x) ≤ A (l,j) j h (j−1) (x) + b (l,j) j (14) such that A (j) j h (j−1) (x) + b (j) j ≤ A (l,j) j h (j−1) (x) + b (l,j) j ≤ A (l,j) j h (j−1) (x) + b (l,j) j ≤ A (j) j h (j−1) (x) + b (j) j for all l (l,j−1) ≤ h (j−1) (x) ≤ l (l,j−1) , by the fact that the concrete bounds are tighter for C l sub .Finally, substituting the bounds in Equations 9 and 10 (for h (j−1) ), and using the tightness result in the inductive hypothesis for j −1, we obtain linear bounds for h (j) (x) in terms of of the input x, such that the inductive hypothesis for j holds.
ReLU split: We use C l sub and C r sub to denote the input subregions when fixing unstable ReLU neuron z (x) < 0}.In the following, we prove that g l,i (x) ≥ g i (x) for all x ∈ C l sub .Assume we fix one unstable ReLU neuron of layer j − 1, then for all layers 1 ≤ m ≤ j − 1, for all x ∈ C l sub , it holds that: x + b (m) ≤ A (l,m) x + b (l,m) ≤ A (l,m) x + b (l,m) ≤ A (m) x + b (m) (15) where A (m) = A (l,m) , b (m) = b (l,m) and same for the upper bounding parameters.Now consider the bounding procedure for layer j.The linear layer from postactivation layer j − 1 to pre-activation layer j can be encoded as: W (j) a (j−1) (x) + b (j) ≤ h (j) (x) ≤ W (j) a (j−1) (x) + b (j) (16) Consider the post activation function a (j−1) (x) of the unstable neuron z where c = c = 1, d = d = 0, since the unstable neuron is fixed to be active.By substituting the linear relaxation bounds before and after splitting in Equation 12, we obtain the bounding functions with regard to h (j−1) (x) in the following form: A By the fact the relaxation is fixed to be exact for a j for C l sub .Finally, for the bound propagation procedure of layer L, substituting the tightened bounding for h (j−1) (x), we obtain that g l,i (x) = A (l,L) x + b (l,L) ≥ A (L) x + b (L) = g i (x).
Corollary 1.In each refinement iteration, the volume of the polytope approximation T Dom does not decrease.
Proof.In each iteration of Algorithm 1, we replace the polytope T C sub (O) in a leaf subregion with two polytopes T C l sub (O), T C r sub (O) in the DUP under-approximation.By Proposition 2, the total volume of the two new polytopes is at least that of the removed polytope.Thus the volume of the DUP approximation does not decrease.
Similarly, for ReLU splitting, we replace the polytope T C sub (O) in a leaf subregion with two polytopes T C l sub (O), T C r sub (O) where the relaxed bounding functions for one unstable neuron are replaced with exact linear functions, i.e., ch   Secondly, instead of maximizing the volume of the DUP approximation, when over-approximating we wish to minimize the volume.This affects the following steps of our method: -For subregion selection (Section 4.3), we define the priority (Equation 2) to instead be the volume of the over-approximation minus the true volume: The subregion with maximal priority then corresponds to the loosest overapproximation.-For local optimization (Section 4.2), we replace the lower bounds g i with the upper bounds g i constituting the over-approximation, and minimize rather than maximize the estimated volume vol(T C sub (O)) with respect to α.

jProposition 3 .
(x) + d is replaced with the exact linear function a(i) j (x) = h (i) j (x)and a (i) j (x) = 0, respectively, as shown in Figure5(from unstable to stable).By Proposition 2, the total volume of the two new polytopes is at least that of the removed polytope.Thus the volume of the DUP approximation does not decrease.Algorithm 3 is sound for quantitative verification with input splitting.Proof.Algorithm 3 outputs True only if, at some iteration, we have that the exact volume vol(T ) ≥ p × vol(I).Since T is an under-approximation to the restricted preimage f −1 I (O), we have thatvol(f −1 I (O)) vol(I) ≥ vol(T )vol(I) ≥ p, i.e. the quantitative property (I, O, p) holds.

Proposition 4 .
Algorithm 3 is sound and complete for quantitative verification on piecewise linear neural networks with ReLU splitting.
. We

Table 1 :
methods.Vehicle Parking & VCAS.Table 1 presents experimental results on the vehicle parking and VCAS tasks.In the table, we show the number of polytopes (#Poly) Performance comparison in preimage generation.

Table 2 :
Performance of preimage approximation for reinforcement learning tasks.

Table 3 :
Refinement with ReLU split for MNIST (FNN 6 × 100) Invprop) does not support these higher-dimensional input domains.Table 2 summarizes the experimental results.We evaluate Algorithm 1 with input split on a range of tasks/properties and configurations of the input region (e.g., angular velocity θ for Cartpole).Empirically, for the same coverage ratio, our method requires a number of polytopes and time roughly linear in the input region size, with the exception of Dubinsrejoin, where the larger number of output constraints and larger network size contribute to greater relaxation error.MNIST Preimage Approximation with ReLU Split Next, we evaluate the scalability of Algorithm 1 with ReLU splitting by applying it to MNIST image classifiers.

Table 4 :
Comparison with a robustness verifier.

Table 5 :
Performance of preimage generation for vehicle parking.

Table 6 :
Performance of preimage generation for VCAS.

Table 7 :
Ablation study on methodological components.
the polytope is an under-approximation.Given any subregion C sub with polytope approximation T C sub (O),