# Solution-Aware Vs Global ReLU Selection: Partial MILP Strikes Back for DNN Verification

Yuke Liao, Blaise Genest, Kuldeep Meel, Shaan Aryaman

► To cite this version:

Yuke Liao, Blaise Genest, Kuldeep Meel, Shaan Aryaman. Solution-Aware Vs Global ReLU Selection: Partial MILP Strikes Back for DNN Verification. 23rd International Symposium Automated Technology for Verification and Analysis (ATVA 2025), Oct 2025, Bangalore, India. pp.299-320, 10.1007/978-3-032-08707-2_14 . hal-05384407

HAL Id: hal-05384407

https://hal.science/hal-05384407v1

Submitted on 27 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Solution-aware vs global ReLU selection: partial MILP strikes back for DNN verification

Yuke Liao<sup>1</sup>[0009-0004-3763-686X], Blaise Genest<sup>2,3</sup>[0000-0002-5758-1876], Kuldeep Meel<sup>4</sup>[0000-0001-9423-5270], and Shaan Aryaman<sup>5</sup>[0000-0001-7576-0766]

<sup>1</sup> CNRS@CREATE, Singapore yuke.liao@cnrsatcreate.sg

<sup>2</sup> CNRS@CREATE & IPAL, Singapore blaise.genest@cnrsatcreate.sg

<sup>3</sup> CNRS, IPAL, France blaise.genest@cnrs.fr

<sup>4</sup> University of Toronto, Toronto, Canada
meel@cs.toronto.edu

<sup>5</sup> NYU Courant Institute of Mathematical Sciences, New York, USA
aryaman.shaan@gmail.com

**Abstract.** Branch and Bound (BaB) is considered as the most efficient technique for DNN verification: it can propagate bounds over numerous branches, to accurately approximate values a given neuron can take even in large DNNs, enabling formal verification of properties such as local robustness. Nevertheless, the number of branches grows exponentially with important variables, and there are complex instances for which the number of branches is too large to handle even using BaB. In these cases, providing more time to BaB is not efficient, as the number of branches treated is *linear* with the time-out. Such cases arise with verification-agnostic DNNs, non-local properties (e.g. global robustness, computing Lipschitz bound), etc.

To handle complex instances, we revisit a divide-and-conquer approach to break down the complexity: instead of few complex BaB calls, we rely on many small *partial* MILP calls. The crucial step is to select very few but very important ReLUs to treat using (costly) binary variables. The previous attempts were suboptimal in that respect. To select these important ReLU variables, we propose a novel *solution-aware* ReLU scoring (SAS), as well as adapt the BaB-SR and BaB-FSB branching functions as *global* ReLU scoring (GS) functions. We compare them theoretically as well as experimentally, and SAS is more efficient at selecting a set of variables to open using binary variables. Compared with previous attempts, SAS reduces the number of binary variables by around 6 times, while maintaining the same level of accuracy. Implemented in *Hybrid MILP*, calling first $\alpha$, $\beta$-CROWN with a short time-out to solve easier instances, and then partial MILP, produces a very accurate yet efficient verifier, reducing by up to 40% the number of undecided instances to low levels (8-15%), while keeping a reasonable runtime (46s-417s on average per instance), even for fairly large CNNs with 2 million parameters.

# 1 Introduction

Deep neural networks (DNNs for short) have demonstrated remarkable capabilities, achieving human-like or even superior performance across a wide range
<table><caption>Table 1: Accuracy of DNN (class predicted vs ground truth), upper bound on robustness (robustness attacks found on remaining images), and % of images verified by α, β-CROWN with different time-outs (TO) on 7 DNNs, and average runtime per image. The 6 first DNNs are complex instances. The last DNN (ResNet) is an easy instance (trained using Wong to be easy to verify, but with a very low accuracy level), provided for reference.</caption><thead><tr><th>Network<br>Perturbation</th><th>nbr<br>activ.</th><th>Accur.</th><th>Upper<br>Bound</th><th>α, β-CROWN<br>TO=10s</th><th>α, β-CROWN<br>TO=30s</th><th>α, β-CROWN<br>TO=2000s</th></tr></thead><tbody><tr><td>MNIST 5×100<br>ε = 0.026</td><td>500<br>ReLU</td><td>99%</td><td>90%</td><td>33%<br>6.9s</td><td>35%<br>18.9s</td><td>40%<br>1026s</td></tr><tr><td>MNIST 5×200<br>ε = 0.015</td><td>1000<br>ReLU</td><td>99%</td><td>96%</td><td>46%<br>6.5s</td><td>49%<br>16.6s</td><td>50%<br>930s</td></tr><tr><td>MNIST 8×100<br>ε = 0.026</td><td>800<br>ReLU</td><td>97%</td><td>86%</td><td>23%<br>7.2s</td><td>28%<br>20.1s</td><td>28%<br>930s</td></tr><tr><td>MNIST 8×200<br>ε = 0.015</td><td>1600<br>ReLU</td><td>97%</td><td>91%</td><td>35%<br>6.8s</td><td>36%<br>18.2s</td><td>37%<br>1083s</td></tr><tr><td>MNIST 6×500<br>ε = 0.035</td><td>3000<br>ReLU</td><td>100%</td><td>94%</td><td>41%<br>6.4s</td><td>43%<br>16.4s</td><td>44%<br>1003s</td></tr><tr><td>CIFAR CNN-B-Adv<br>ε = 2/255</td><td>16634<br>ReLU</td><td>78%</td><td>62%</td><td>34%<br>4.3s</td><td>40%<br>8.7s</td><td>42%<br>373s</td></tr><tr><td>CIFAR ResNet<br>ε = 8/255</td><td>107496<br>ReLU</td><td>29%</td><td>25%</td><td>25%<br>2s</td><td>25%<br>2s</td><td>25%<br>2s</td></tr></tbody></table>

of tasks. However, their robustness is often compromised by their susceptibility to input perturbations [24]. This vulnerability has catalyzed the verification community to develop various methodologies, each presenting a unique balance between completeness and computational efficiency [18,17,22]. This surge in innovation has also led to the inception of competitions such as VNNComp [6], which aim to systematically evaluate the performance of neural network verification tools. While the verification engines are generic, the benchmarks usually focus on local robustness, i.e. given a DNN, an image and a small neighbourhood around this image, is it the case that all the images in the neighbourhood are classified in the same way. For the past 5 years, VNNcomp has focused on rather easy instances, that can be solved within tens of seconds (the typical hard time-out is 300s). For this reason, DNN verifiers in the past years have mainly focused on optimizing for such easy instances. Among them, NNenum [2], Marabou [18,29], and PyRAT [12], respectively 4th, 3rd and 2nd of the last VNNcomp'24 [5] and 5th, 2nd and 3rd of the VNNcomp'23 [4]; MnBAB [15], 2nd in VNNcomp'22 [20], built upon ERAN [22] and PRIMA [19]; and importantly, α, β-CROWN [27,31], the winner of the last 4 VNNcomp, benefiting from branch-and-bound based methodology [32,8]. We will thus compare mainly with α, β-CROWN experiments as gold standard in the following<fnref n="6" />.

<footnote n="6">⁶ GCP-CROWN [32] is slightly more accurate than α, β-CROWN on the DNNs we tested, but necessitates IBM CPLEX solver, which is not available to us</footnote>
α, β-CROWN, as well as BaBSR [8] and MN-BaB [15], rely on Branch and Bound technique (BaB), which call BaB once per output neuron (few calls). In the worst case, this involves considering all possible ReLU configurations, though branch and bound typically circumvents most possibilities. For easy instances, BaB is highly efficient as all branches can be pruned early. However, BaB methods hit a complexity barrier when verifying more complex instances, due to an overwhelming number of branches (exponential in the height of branches that cannot be pruned as they need too many variables to branch over). This can be clearly witnessed on the verification-agnostic [9] DNNs of Table 1 (6 first DNNs), where vastly enlarging the time-out only enables to verify few more % of images, leaving a large proportion (20% - 50%) of images undecided despite the large runtime. As argued in [9], there are many situations (workflow, no access to the dataset...) where using specific trainers to learn easy to verify DNN is simply not possible, leading to verification-agnostic networks, and such cases should be treated as well as DNNs specifically trained to be easy to verify, e.g. using [30]. Verification-agnostic are the simplest instances to demonstrate the scaling behavior of BaB on complex instances using standard local robustness implementations. Other complex instances include solving non-local properties, e.g. global robustness computed through Lipschitz bound [28], etc. The bottom line is that one cannot expect to have only easy instances to verify. It is important to notice that the number of activation functions of the DNN is a poor indicator of the hardness of the instance, e.g. 5 × 100 with 500 ReLUs is far more complex to certify (50% undecided images) than 100 times bigger ResNet (0% undecided images), see Table 1.

Other standard non-BaB methods such as Marabou, NNenum or a Full MILP encoding, show similar poor performance on such complex instances as well, even with a large 10 000s Time-out: Table 2 reveals that only NNenum succeeds to verify images not verified by α, β-CROWN, limited to 9% more images out of the 50% undecided images on 5 × 100, and with a very large runtime of almost 5000s per image. It appeared pointless to test these verifiers on larger networks.

Our main contributions address the challenges to verify *complex* instances efficiently, as current methods are not appropriate to verify such instances:

1. We revisit the idea from [16] to consider small calls to a partial MILP (pMILP) solver, i.e. with few binary variables encoding few ReLU functions exactly (other being encoded with more efficient but less accurate linear variables), to compute bounds for each neuron inductively, hence many ($O(n)$,

<table><caption>Table 2: Result of non-BaB methods on the hard 5 × 100 with TO = 10 000s. Only NNenum verifies more instances (9% out of 50% undecided images) than α, β-CROWN (40%), at the cost of a much larger runtime (4995s vs 1026s).</caption><thead><tr><th>Network</th><th>Accuracy</th><th>Upper</th><th>Marabou 2.0</th><th>NNenum</th><th>Full MILP</th></tr></thead><tbody><tr><td>MNIST 5×100</td><td>99%</td><td>90%</td><td>28%</td><td>49%</td><td>40%</td></tr><tr><td>ε = 0.026</td><td></td><td></td><td>6200s</td><td>4995s</td><td>6100s</td></tr></tbody></table>
the number of neuron) small calls, with a complexity exponential only in
the few binary variables (Section 4). Compared to the few (one per output
neuron) complex call to BaB, each with a worst case complexity exponential
in the number of neurons of the DNN (which is far from the actual complex-
ity thanks to pruning branches in BaB - but which can be too large as we
shown in the 6 first DNNs of Table 1). Two questions arise: how to select
few very important ReLUs?, and is computing the bounds of intermediate
neurons a good trade-off compared with the theoretical loss of accuracy due
to selecting only some binary ReLUs? Answer to these questions were not
looking very promising judging by previous attempt [16], which was using a
simple selection heuristic of nodes in the previous layer only.

2. On the first question, we adapted from BaB-SR [8] and FSB [11], which choose branching nodes for BaB, *global scoring (GS)* functions to choose important ReLUs (Section 5). These revealed to be much more accurate than the simple heuristic in [16]. However, we also uncover that the *improvement function* that GS tries to approximate depends heavily upon the mode of the ReLU function, and as this mode is unavailable to GS, there are many cases in which GS is far from the improvement (with both under or over-approximation of the improvement function).

3. We thus designed a *novel solution-aware scoring (SAS)*, which uses the solution of a unique LP call, that provides the mode to consider (Section 6). Theoretically, we show that *SAS* is always an over-approximation of the improvement (Proposition 2), which implies that a small *SAS* value implies that the ReLU is unnecessary. Experimentally, we further show that *SAS* is very close to the actual improvement, closer than *GS*, and that overall, the accuracy from *SAS* is significantly better than *GS*. Compared with the heuristic in [16], *SAS* is much more efficient (≈ 6 times less binary variables for same accuracy (Fig. 4)).

4. Compared with many calls using full MILP, where all the ReLUs are encoded as binary variables, *SAS* (and *GS*) encode only a subset of ReLUs as binary and others as linear variables. The model in full MILP is fully accurate, while *SAS* (and *GS*) are abstractions, and thus much faster to solve. While full MILP is thus *asymptotically* more accurate than *SAS* and *GS*, *experimentally*, every reasonable time-out leads to much better practical accuracy of *SAS* (and *GS*) (see Fig. 5).

5. For the second question, we propose a new verifier, called *Hybrid MILP*, in-
voking first α, β-CROWN with short time-out to settle the easy instances. On
those (*hard*) instances which are neither certified nor falsified, we call pMILP
with few neurons encoded as binary variables. Experimental evaluation re-
veals that Hybrid MILP achieves a beneficial balance between accuracy and
completeness compared to prevailing methods. It reduces the proportion of
undecided inputs from 20–58% (α, β-CROWN with 2000s TO) to 8–15%,
while taking a reasonable average time per instance (46–420s), Table 4. It
scales to fairly large networks such as CIFAR-10 CNN-B-Adv [9], with more
than 2 million parameters.
Limitation: We consider DNNs employing the standard ReLU activation function, though our findings can be extended to other activation functions, following similar extension by [16], with updated MILP models e.g. for maxpool.

## 1.1 Related Work

We compare Hybrid MILP with major verification tools for DNNs to clarify our methodology and its distinction from the existing state-of-the-art. It scales while preserving good accuracy, through targeting a limited number of binary variables, stricking a good balance between exact encoding of a DNN using MILP [26] (too slow) and LP relaxation (too inaccurate). MIPplanet [14] opts for a different selection of binary variables, and execute one large MILP encoding instead of Hybrid MILP's many small encodings, which significantly reduce the number of binary variables necessary for each encoding. In [16], small encodings are also considered, however with a straightforward choice of binary nodes based on the weight of outgoing edges, which need much more binary variables (thus runtime) to reach the same accuracy.

Compared with [1], which uses pMILP in an abstraction refinement loop, they iteratively call pMILP to obtain bounds for the same output neuron, opening more and more ReLUs. This scales only to limited size DNN (500 neurons), because of the fact that many ReLUs need to be open (and then Gurobi takes a lot of time) and the iterative nature which cannot be parallelized, unlike our method which scales up to 20.000 neurons.

The fact that pure BaB is not that efficient for e.g. verification-agnostic (even very small) DNNs has been witnessed before [23]. The workaround, e.g. in *refined* $\alpha$, $\beta$-CROWN, was to precompute very accurate bounds for the first few neurons of the DNN using a complete full MILP encoding, and then rely on a BaB call from that refined bounds (more complex calls to full MILP and BaB). Non-surprisingly, this very slow technique does not scale but to small DNNs (max 2000 ReLU activation functions). Hybrid MILP on the other hand relies only on small calls: it is much more efficient on small DNNs, and it can scale to larger DNNs as well: we demonstrated strong performance with at least one order of magnitude larger networks (CNN-B-Adv).

Last, ERAN-DeepPoly [22] computes bounds on values very quickly, by abstracting the weight of every node using two functions: an upper function and a lower function. While the upper function is fixed, the lower function offers two choices. It relates to the LP encoding through Proposition 1 [10]: the LP relaxation precisely matches the intersection of these two choices. Consequently, LP is more accurate (but slower) than DeepPoly, and Hybrid MILP is considerably more precise. Regarding PRIMA [19], the approach involves explicitly maintaining dependencies between neurons.

Finally, methods such as Reluplex / Marabou [17,18] abstract the network: they diverge significantly from those abstracting values such as PRIMA, $\alpha$, $\beta$-CROWN) [19,27], Hybrid MILP. These network-abstraction algorithms are designed to be *intrinsically complete* (rather than asymptotically complete as BaB),
but this comes at the price of significant scalability challenges, and in practice they time-out on complex instances as shown in Table 2.

## 2 Notations and Preliminaries

In this paper, we will use lower case latin $a$ for scalars, bold $z$ for vectors, capitalized bold $W$ for matrices, similar to notations in [27]. To simplify the notations, we restrict the presentation to feed-forward, fully connected ReLU Deep Neural Networks (DNN for short), where the ReLU function is ReLU : $\mathbb{R} \to \mathbb{R}$ with ReLU($x$) = $x$ for $x \ge 0$ and ReLU($x$) = 0 for $x \le 0$, which we extend componentwise on vectors.

An $\ell$-layer DNN is provided by $\ell$ weight matrices $W^i \in \mathbb{R}^{d_i \times d_{i-1}}$ and $\ell$ bias vectors $b^i \in \mathbb{R}^{d_i}$, for $i = 1, \dots, \ell$. We call $d_i$ the number of neurons of hidden layer $i \in \{1, \dots, \ell - 1\}$, $d_0$ the input dimension, and $d_\ell$ the output dimension.

Given an input vector $z^0 \in \mathbb{R}^{d_0}$, denoting $\hat{z}^0 = z^0$, we define inductively the value vectors $z^i, \hat{z}^i$ at layer $1 \le i \le \ell$ with

$$z^i = W^i \cdot \hat{z}^{i-1} + b^i \qquad \hat{z}^i = \text{ReLU}(z^i). \qquad (1)$$

The vector $\hat{z}$ is called post-activation values, $z$ is called pre-activation values, and $z^i$ is used to call the $j$-th neuron in the $i$-th layer. For $x = z^0$ the (vector of) input, we denote by $f(x) = z^\ell$ the output. Finally, pre- and post-activation neurons are called nodes, and when we refer to a specific node/neuron, we use $a, b, c, d, n$ to denote them, and $W_{a,b} \in \mathbb{R}$ to denote the weight from neuron $a$ to $b$. Similarly, for input $x$, we denote by value$_x(a)$ the value of neuron $a$ when the input is $x$.

Concerning the verification problem, we focus on the well studied local robustness question. Local robustness asks to determine whether the output of a neural network will be affected under small perturbations to the input. Formally, for an input $x$ perturbed by $\varepsilon > 0$ under distance $d$, then the DNN is locally $\varepsilon$-robust in $x$ whenever:

$$\forall x' \text{ s.t. } d(x, x') \le \varepsilon, \text{ we have } \operatorname{argmax}(f(x')) = \operatorname{argmax}(f(x)) \qquad (2)$$

## 3 Value Abstraction for DNN verification

In this section, we describe different value (over-)abstractions on $z$ that are used by efficient algorithms to certify robustness around an input $x$. Over-abstractions of values include all values for $z$ in the neighbourhood of $x$, and thus a certificate for safety in the over-abstraction is a proof of safety for the original input $x$.
<figure><img src="image_2.png" /><figcaption>Fig. 1: A DNN. Every neuron is separated into 2 nodes, $n$ pre- and $\hat{n}$ post-ReLU activation.</figcaption></figure>

## 3.1 The Box Abstractions

The concept of value abstraction involves calculating upper and lower bounds for the values of certain neurons in a Deep Neural Network (DNN) when inputs fall within a specified range. This approach aims to assess the network's robustness without precisely computing the values for every input within that range.

Firstly, it's important to note that weighted sums represent a linear function, which can be explicitly expressed with relative ease. However, the ReLU (Recified Linear Unit) function presents a challenge in terms of accurate representation. Although ReLU is a relatively straightforward piecewise linear function with two modes (one for $x < 0$ and another for $x \ge 0$), it is not linear. The complexity arises when considering the compounded effects of the ReLU function across the various layers of a ReLU DNN. It's worth noting that representing ReLU($x$) precisely is feasible when $x$ is "stable", meaning it's consistently positive or consistently negative, as there's only one linear mode involved in each scenario. Consequently, the primary challenge lies in addressing "unstable" neurons, where the linearity of the function does not hold consistently.

Consider the simpler abstraction, termed “Box abstraction”, recalled e.g. in [22]: it inductively computes the bounds for each neuron in the subsequent layer independently. This is achieved by considering the weighted sum of the bounds from the previous layer, followed by clipping the lower bound at $\max(0, \text{lower bound})$ to represent the ReLU function, and so forth. For all $i \ge 3$, define $x_i = \text{value}_x(n_i)$, where $x = (x_1, x_2)$. Taking the DNN example from Fig 1, assume $x_1, x_2 \in [-1, 1]$. This implies that $x_3, x_4 \in [-2, 2]$. After applying the ReLU function, $\hat{x}_3, \hat{x}_4$ are constrained to $[1.5, 3.5]$, leading to $x_5 \in [0, 6]$ and $x_6 \in [0, 2]$. The bounds for $n_1, \dots, n_4$ are exact, meaning for every $\alpha$ within the range, an input $y$ can be found such that $\text{value}_y(n_i) = \alpha$. However, this precision is lost from the next layer (beginning with $n_5, n_6$) due to potential dependencies among preceding neurons. For example, it is impossible for $x_5 = \text{value}_x(n_5)$ to reach 6, as it would necessitate both $x_3 = 2$ and $x_4 = 2$, which is not possible at the same time as $x_3 = 2$ implies $x_1 = x_2 = 1$ and $x_4 = 2$ implies $x_2 = -1$ (and $x_1 = 1$), a contradiction.
In [13,22] and others, the triangular abstraction was proposed:

$$ \mathrm{ReLU}(x) = \max(0, x) \le \hat{x} \le \mathrm{UB} \frac{x - \mathrm{LB}}{\mathrm{UB} - \mathrm{LB}} \quad (3) $$

It has two lower bounds (the 0 and identity functions), and one upper bound. DeepPoly [22] chooses one of the two lower bounds for each neuron $x$, giving rise to a greedy quadratic-time algorithm to compute very fast an abstraction of the value of $\hat{x}$ (but not that accurately).

## 3.2 MILP and LP encodings for DNNs

At the other end of the spectrum, we find the Mixed Integer Linear Programming (MILP) value abstraction, which is a complete (but inefficient) method. Consider an unstable neuron $n$, whose value $x \in [\mathrm{LB}, \mathrm{UB}]$ with $\mathrm{LB} < 0 < \mathrm{UB}$. The value $\hat{x}$ of $ReLU(x)$ can be encoded exactly in an MILP formula with one integer (actually even binary) variable $a$ valued in $\{0, 1\}$, using constants $\mathrm{UB}, \mathrm{LB}$ with 4 constraints [26]:

$$ \hat{x} \ge x \wedge \hat{x} \ge 0 \wedge \hat{x} \le \mathrm{UB} \cdot a \wedge \hat{x} \le x - \mathrm{LB} \cdot (1 - a) \quad (4) $$

For all $x \in [\mathrm{LB}, \mathrm{UB}] \setminus 0$, there exists a unique solution $(a, \hat{x})$ that meets these constraints, with $\hat{x} = \mathrm{ReLU}(x)$ [26]. The value of $a$ is 0 if $x < 0$, and 1 if $x > 0$, and can be either if $x = 0$. This encoding approach can be applied to every (unstable) ReLU node, and optimizing its value can help getting more accurate bounds. However, for networks with hundreds of unstable nodes, the resulting MILP formulation will contain numerous binary variables and generally bounds obtained will not be accurate, even using powerful commercial solvers such as Gurobi.

MILP instances can be linearly relaxed into LP over-abstraction, where variables originally restricted to integers in $\{0, 1\}$ (binary) are relaxed to real numbers in the interval $[0, 1]$, while maintaining the same encoding. As solving LP instances is polynomial time, this optimization is significantly more efficient. However, this efficiency comes at the cost of precision, often resulting in less stringent bounds. This approach is termed the LP abstraction. We invoke a follore result on the LP relaxation of (4), for which we provide a direct and explicit proof:

**Proposition 1.** [10] The LP relaxation of (4) is equivalent with the triangular abstraction (3).

*Proof.* Consider an unstable neuron $n$, that is $\mathrm{LB} < 0 < \mathrm{UB}$. The lower bound on $\hat{x}$ is simple, as $\hat{x} \ge 0 \wedge \hat{x} \ge x$ is immediately equivalent with $\hat{x} \ge \mathrm{ReLU}(x)$.

We now show that the three constraints $\hat{x} \le \mathrm{UB} \cdot a \wedge \hat{x} \le x - \mathrm{LB} \cdot (1 - a) \wedge a \in [0, 1]$ translates into $\hat{x} \le \mathrm{UB} \frac{x - \mathrm{LB}}{\mathrm{UB} - \mathrm{LB}}$. We have $\hat{x}$ is upper bounded by $\max_{a \in [0, 1]}(\min(\mathrm{UB} \cdot a, x - \mathrm{LB}(1 - a)))$, and this bound can be reached. Furthermore, using standard function analysis tools (derivative...), we can show
<caption>**Algorithm 1: pMILP<sub>K</sub>**</caption>

**Input:** Bounds [LB(*m*), UB(*m*)] for nodes *m* at layer 0
**Output:** Bounds [LB, UB] for every node *n*
1 **for** layer *k* = 1, ..., *ℓ* **do**
2 &nbsp;&nbsp;&nbsp;&nbsp;**for** neuron *n* in layer *k* **do**
3 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute *X* a set of *K* nodes most important for target neuron *n*.
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Run pMILP<sub>*X*</sub>(*n*) to obtain [LB(*n*), UB(*n*)] using MILP<sub>*X*</sub> and
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;additional constraints $\bigvee_{m \text{ in layer } <k} \text{value}(m) \in [\text{LB}(m), \text{UB}(m)]$.

that the function $a \mapsto \min(\text{UB} \cdot a, x - \text{LB}(1 - a))$ attains its maximum when
$\text{UB} \cdot a = x - \text{LB}(1 - a)$, leading to the equation $(\text{UB} - \text{LB})a = x - \text{LB}$ and
consequently $a = \frac{x - \text{LB}}{\text{UB} - \text{LB}}$. This results in an upper bound $\hat{x} \le \text{UB} \frac{x - \text{LB}}{\text{UB} - \text{LB}}$, which
can be reached, hence the equivalence. □

4 Partial MILP

In this paper, we revisit the use of *partial MILP* (pMILP for short, see [16]), to
get interesting trade-offs between accuracy and runtime. Let *X* be a subset of the
set of unstable neurons, and *n* a neuron for which we want to compute upper and
lower bounds on values: pMILP<sub>*X*</sub>(*n*) simply calls Gurobi to minimize/maximize
the value of *n* with the MILP model encoding (equation (2) Section 3.2), where
variable *a* is:

- binary for neurons in *X* (exact encoding of the ReLU),

- linear for neurons not in X (linear relaxation).

Formally, we denote by MILP<sub>X</sub> the MILP encoding where the variable *a*
encoding ReLU(*y*) is binary for *y* ∈ *X*, and linear for *y* /∈ *X* (*y* being unstable
or stable). We say that nodes in *X* are *opened*. That is, if *X* is a strict subset of
the set of all unstable neurons, then MILP<sub>X</sub> is an abstraction of the constraints
in the full MILP model. If *X* covers all unstable neurons, then MILP<sub>X</sub> is exact
as there is only one linear mode for stable variables. The worst-case complexity
to solve MILP<sub>X</sub> is NP-complete in |*X*|, ie the number of binary variables.

To reduce the runtime, we will limit the size of subset X. This a priori hurts
accuracy. To recover some of this accuracy, we use an iterative approach similar
to the box abstraction or DeepPoly [22], computing lower and upper bounds
LB, UB for neurons n of each layer iteratively, that are used when computing
values of the next layer. So we trade-off the NP-complexity in |X| with a linear
complexity in the number of neurons. Notice that all the neurons of a layer can
be treated in parallel.

We provide the pseudo code for pMILP<sub>K</sub> in Algorithm 1. pMILP<sub>K</sub> has a
worst case complexity bounded by O(N · MILP(N, K)), where N is the number
of nodes of the DNN, and MILP(N, K) is the complexity of solving a MILP
program with $K$ binary variables and $N$ linear variables. We have MILP($N$, $K$) $\le 2^K$ LP($N$) where LP($N$) is the Polynomial time to solve a Linear Program with $N$ variables, $2^K$ being an upper bound. Solvers like Gurobi are quite adept and usually do not need to evaluate all $2^K$ ReLU configurations to deduce the bounds. It is worth mentioning that the "for" loop iterating over neurons $n$ in layer $k$ (line 2) can be executed in parallel, because the computation only depends on bounds from preceding layers, not the current layer $k$. This worst-case complexity compares favorably when $K \ll N$ with the worst-case complexity of BaB which would be $O(out \times 2^N)$, where out is the number of output neurons. Of course, actual complexity is always much better than the worst case, thanks to different heuristics, such as pruning, which is what happens for easy instances, for which the linear cost $N$ in pMILP (which cannot be avoided) is actually the limiting factor. Only experimental results can tell whether the accuracy lost using reduced variables in pMILP can be a good trade off vs efficiency, both depending on the number $K$ chosen.

If $K$ is sufficiently small, this approach is expected to be efficient. The crucial factor in such an approach is to select few opened ReLU nodes in $X$ which are the most important for the accuracy. An extreme strategy was adopted in [16], where only ReLU nodes of the immediate previous layer can be opened, and the measure to choose ReLU $a$ when computing the bounds for neuron $b$ was to consider $|W_{ab}|(UB(a) - LB(a))$.

# 5 Global Scoring Functions

To select nodes $X$ for pMILP$_X$, we first adapt from the BaB-SR [8] and FSB [11] functions, which originally iteratively select one node to branch on for BaB. Intuitively, we extract the scoring $s_{SR}$ and $s_{FSB}$ from both BaB-SR and FSB, as the BaB bounding step is not adapted for pMILP. We will call such functions global scoring (GS) functions.

They both work by backpropagating gradients vectors $\lambda$ from the neurons under consideration in layer $n$, back to neurons to be potentially selected. To do so, they consider the rate of ReLU($u_k$) to be $r(u_k) = \frac{\max(0,UB(u_k))}{\max(0,UB(u_k))-\min(0,LB(u_k))} \in [0, 1]$, with $r(b) = 0$ iff $UB(b) \le 0$ and $r(b) = 1$ iff $LB(b) \ge 0$.

$$ \lambda_{n-1} = -(W^n)^T 1, \quad \lambda_{k-1} = (W^k)^T (r(u_k) \odot \lambda_k) \quad k \in [n-1, 2] \quad (5) $$

Then, the scoring functions $s_{SR}$ and $s_{FSB}$ for ReLUs in layer $k$ are computed by approximating how each node would impact the neurons in layer $n$, using the computed $\lambda$, differing only in how they handle the bias $b_k$:

$$ s_{SR}(k) = \left| r(u_k) \odot \text{LB}(u_k) \odot \max(\lambda_k, 0) + \max\{0, \lambda_k \odot b_k\} - r(u_k) \odot \lambda_k \odot b_k \right| \\ s_{FSB}(k) = \left| r(u_k) \odot \text{LB}(u_k) \odot \max(\lambda_k, 0) + \min\{0, \lambda_k \odot b_k\} - r(u_k) \odot \lambda_k \odot b_k \right| $$
<figure><img src="image_3.png" /><figcaption>Fig. 2: A running example with parametric weights *c* and *d*</figcaption></figure>

Then, ReLU are ranked using these scores, to select the most important unstable ReLUs (with $LB(u_k) < 0 < UB(u_k)$).

*Running Example:* Consider Fig. 2. It has no bias, so $s_{SR} = s_{FSB}$. We have $\lambda(b) = d$, $\lambda(a) = \frac{cd}{2}$. The value of $s_{FSB}(a)$ depends on the signs of *c*, *d*.

We perform a novel comparison between $s_{FSB}(a)$ and $\Delta(z)$, the difference on the maximal bound computed by pMILP$_X(z)$ when opening the ReLU of node *a* ($X = \{a\}$), yielding value($\hat{a}$) = ReLU(value(*a*)), versus having $X = \emptyset$, for which value($\hat{a}$) can be any value in the triangle approximation (Prop. 1).

The most favorable cases are $c < 0 < d$ and $d < 0 < c$: as $cd < 0$, we have $s_{FSB}(a) = \max(0, \frac{cd}{4}) = 0$. Because $cd < 0$, both *a* and $\hat{a}$ need to be minimized by MILP$_X$ in order to maximize the value of *z*. For $X = \emptyset$, the LP approximation ($=\text{MILP}_{\emptyset}$) will thus set value($\hat{a}$) = ReLU(value(*a*)) as this is the minimum value in the triangle approximation. Notice that opening $X = \{a\}$ yields the same value($\hat{a}$) = ReLU(value(*a*)). That is, opening node *a* is not improving the bound $UB(z)$, as correctly predicted by the score $s_{FSB}(a) = \Delta(a) = 0$.

Case $c > 0, d > 0$: we have $s_{FSB}(a) = \frac{cd}{4}$. The value of $a, \hat{a}, b, \hat{b}$ should be maximized to maximize the value of *z*, because $W_{bz} = d > 0$ and $W_{ab} = c > 0$. Now, let us call *val(a)* the maximum value for *a* (same for MILP$_\emptyset$ and MILP$_{\{a\}$}). The maximum value for $\hat{a}$ under LP ($X = \emptyset$) is $\frac{1}{2} \cdot (\text{value}(a) - \text{LB}(a))$ following the triangle approximation. Now, if $X = \{a\}$, then the value of $\hat{a}$ is ReLU(value(*a*)): the difference $\Delta(\hat{a})$ is between 0 (for value(*a*) ∈ {LB, UB}) and $\frac{1}{2}$ (for value(*a*) = 0). The difference $\Delta(b)$ for the value of *b* is thus between 0 and $\frac{c}{2}$, which means $\Delta(\hat{b}) \in [0, \frac{cd}{4}]$, using the upper function of the triangle approximation of rate $r(b) = \frac{1}{2}$, as the value of $\hat{b}$ should be maximized. This means a difference $\Delta(z) \in [0, \frac{cd}{4}]$ (depending on the maximal value *val(a)*), to compare with the fixed $s_{FSB}(a) = \frac{cd}{4}$.

The last case is the most problematic: Case $c < 0, d < 0$, implying that $s_{FSB}(a) = \frac{cd}{4}$, because $cd > 0$. Opening ReLU(*a*) will have the same impact on the value of $\hat{a}, b$ than in case $c > 0, d > 0$, with $\Delta(b) \in [0, \frac{c}{2}]$. However, as $d < 0$, value of $\hat{b}$ needs to be minimized to maximize the value of *z*. That is, the value of $\hat{b}$ will be the value of ReLU(*sol(b)*), and the change $\Delta(\hat{b})$ will be either 0 in case value(*b*) < 0, or $\Delta(b) \in [0, \frac{c}{2}]$ for value(*b*) > 0. That is, *z* will be modified by either $\Delta(z) = 0$ or $\Delta(z) = d\Delta(b)$, to be compared with
the fix value $s_{FSB}(a) = \frac{cd}{4}$, which is not always an overapproximation: we have $\Delta(z) = \frac{cd}{2} > s_{FSB}(a)$, if value(a) = 0 and value(b) > 0.

We call global scoring (GS) these functions $s_{FSB}, s_{SR}$ because they score ReLUs as accurately as possible, considering that they do not have access to the values value(a), value(b) maximizing z. Following this analysis of $s_{FSB}, s_{SR}$, next Section presents a novel scoring function more accurate wrt $\Delta(z)$.

# 6 Solution-Aware Scoring.

In this section, we propose Solution-Aware Scoring (SAS), to evaluate accurately how opening a ReLU impacts the accuracy. To do so, SAS considers explicitly a solution to a unique LP call, which is reasonably fast to obtain as there is no binary variables (polynomial time). Assume that we want to compute an upper bound for neuron z on layer $\ell_z$. We write $n < z$ if neuron n is on a layer before $\ell_z$, and $n \le z$ if $n < z$ or $n = z$. We denote $(\text{Sol\_max}_X^z(n))_{n \le z}$ a solution of $M_X$ maximizing z: $\text{Sol\_max}_X^z(z)$ is the maximum of z under $M_X$.

Consider $(\text{sol}(n))_{n \le z} = (\text{Sol\_max}_\emptyset^z(n))_{n \le z}$, a solution maximizing the value for z when all ReLU use the LP relaxation. Function $\text{Improve\_max}^z(n) = \text{sol}(z) - \text{Sol\_max}_{\{n\}}^z(z)$, accurately represents how much opening neuron $n < z$ reduces the maximum computed for z compared with using only LP. We have $\text{Improve\_max}^z(n) \ge 0$ as $\text{Sol\_max}_{\{n\}}^z$ fulfills all the constraints of $M_\emptyset$, so $\text{Sol\_max}_{\{n\}}^z(z) \le \text{sol}(z)$. Computing exactly $\text{Improve\_max}^z(n)$ would need a MILP call on $M_{\{n\}}$ for every neuron $n \le z$, which would be very time consuming. Instead, the SAS function uses a (single) LP call to compute $(\text{sol}(n))_{n \le z}$, with negligible runtime wrt the forthcoming $\text{MILP}_X$ call, and yet accurately approximates $\text{Improve\_max}^z(n)$ (Fig. 4).

For a neuron b on the layer before layer $\ell_z$, we define:

$$ \text{SAS\_max}^z(b) = W_{bz} \times (\text{sol}(\hat{b}) - \text{ReLU}(\text{sol}(b))) \quad (6) $$

*Comparison:* consider b with $W_{bz} < 0$: to maximize z, the value of $\text{sol}(\hat{b})$ is minimized by LP, ie $\text{sol}(\hat{b}) = \text{ReLU}(\text{sol}(b))$ thanks to Proposition 1. Thus, we have $\text{SAS\_max}^z(b) = 0 = \text{Improve}_{\text{max}}^z(b)$. Notice that the original scoring function $|W_{bz}|(\text{UB}(b) - \text{LB}(b))$ [16] would be possibly very large in this case. However, GS scoring functions from BaB-SR and FSB would also accurately compute $s_{FSB}(b) = s_{SR} = 0$. Notice that SAS does not need to consider bias explicitly, unlike $s_{FSB}, s_{SR}$, as they are already accounted for in the solution considered.

Consider a neuron a two layers before $\ell_z$, b denoting neurons in the layer $\ell$ just before $\ell_z$. Recall the rate $r(b) = \frac{\text{max}(0, \text{UB}(b))}{\text{max}(0, \text{UB}(b)) - \min(0, \text{LB}(b))} \in [0, 1]$. We define:
<figure><img src="image_4.png" /><figcaption>Fig. 3: Different cases for ReLU(b)</figcaption></figure>

$$ \Delta(\hat{a}) = \text{ReLU}(\text{sol}(a)) - \text{sol}(\hat{a}) \qquad (7) $$

$$ \forall b \in \ell, \Delta(b) = W_{ab}\Delta(\hat{a}) \qquad (8) $$

$$ \forall b \in \ell, \Delta(\hat{b}) = \begin{cases} r(b)\Delta(b), & \text{for } W_{bz} > 0 \end{cases} \qquad (9a) $$

$$ \max(0, \Delta(b) + \text{sol}(b)), \quad \text{for } W_{bz} < 0 \text{ and } \text{sol}(b) \ge 0 \qquad (9b) $$

$$ \max(0, \Delta(b) + \text{sol}(b)), \quad \text{for } W_{bz} < 0 \text{ and } \text{sol}(b) < 0 \qquad (9c) $$

$$ \text{SAS\_max}^z(a) = \Delta(z) = - \sum_{b \in \ell} W_{bz} \Delta(\hat{b}) \qquad (10) $$

*Comparison:* First, the original [16] does not propose a formula for node a two layers before z. So we will compare SAS with GS. Consider again the running example of Fig. 2.

In the case $c \cdot d < 0$, we have $\text{SAS\_max}^z(a) = \text{Improve}_{\text{max}}^z(b) = s_{\text{FSB}}(a) = s_{\text{SR}}(a) = 0$, as $\Delta(\hat{a}) = 0$.

In the case $c > 0, d > 0$, $\Delta(\hat{a}) = \text{ReLU}(\text{sol}(a)) - \text{sol}(\hat{a})$ is precise, whereas the corresponding $\Delta_{\text{FSB}}(\hat{a}) = \text{LB}(a)r(a)$ is only an upperbound.

The last case $c < 0, d < 0$ is the most extreme: $\Delta(\hat{b})$ adapts to the case (b),(c) in Fig. 3 leveraging the value $\text{sol}(b)$, which yields very different values, whereas the corresponding $\Delta_{\text{FSB}}(\hat{b})$ is always $r(b)\Delta_{\text{FSB}}(b)$:

- For $\text{sol}(b) \ll 0$, we will have $\Delta(\hat{b}) = 0 < \Delta_{\text{FSB}}(\hat{b})$.

- For $\text{sol}(b) \gg 0$, we will have $\Delta(\hat{b}) = \Delta(b) > \frac{1}{2}\Delta_{\text{FSB}}(\hat{b})$ as $r(b) = \frac{1}{2}$.

Further, we can show that SAS is a safe overapproximation of $\text{Improve\_max}^z(a)$, which does not hold for $s_{\text{FSB}}, s_{\text{SR}}$ (because of the case $\text{sol}(b) \gg 0$):

**Proposition 2.** $0 \le \text{Improve\_max}^z(a) \le \text{SAS\_max}^z(a)$.

In particular, for all nodes a with $\text{SAS\_max}^z(a) = 0$, we are sure that this node is not having any impact on $\text{Sol\_max}^z_{\{a\}}(z)$.
Proof. Consider $sol'(n)_{n \le z}$ with $sol'(n) = sol(n)$ for all $n \notin \{z, \hat{a}\} \cup \{b, \hat{b} | b \in \ell\}$. In particular, $sol'(a) = sol(a)$. Now, define $sol'(\hat{a}) = ReLU sol(a)$. That is, $sol'(\hat{a})$ is the correct value for $\hat{a}$, obtained if we open neuron $a$, compared to the LP abstraction for $sol(\hat{a})$. We define $sol'(b) = sol(b) + \Delta(b)$ and $sol'(\hat{b}) = sol(\hat{b}) + \Delta(\hat{b})$. Last, $sol'(z) = sol(z) + \sum_{b \in \ell} W_{bz} \Delta(\hat{b})$. We will show:

$$(sol'(n))_{n \le z} \text{ satisfies the constraints in } \mathcal{M}_{\{a\}} \qquad (11)$$

This suffices to conclude: as $sol'(z)$ is a solution of $\mathcal{M}_{\{a\}}$, it is smaller or equal to the maximal solution: $sol'(z) \le Sol\_max_{\{a\}}(z)$. That is, $sol(z) - sol'(z) \ge sol(z) - Sol\_max_{\{a\}}(z)$, i.e. $SAS\_max^z(a) \ge Improve\_max^z(a)$. In particular, we have that $SAS\_max^z(a) \ge 0$, which was not obvious from the definition.

Finally, we show (11). First, opening $a$ changes the value of $\hat{a}$ from $sol(\hat{a})$ to ReLU($sol(a) = sol(\hat{a}) + \Delta(a)$), and from $sol(b)$ to $sol(b) + \Delta(b)$. The case of $\Delta(\hat{b})$ is the most interesting: If (a) $W_{bz} > 0$, to maximize $z$, the LP solver sets $sol(\hat{b})$ to the maximal possible value, which is $r(b)sol(b) + Cst$ according to Proposition 1. Changing $b$ by $\Delta(b)$ thus results in changing $sol(\hat{b})$ by $r(b)\Delta(b)$. If $W_{bz} \le 0$, then the LP solver sets $sol(\hat{b})$ to the lowest possible value to maximize $z$, which happens to be ReLU($b$) according to Proposition 1. If (b) $sol(b) > 0$, then $sol(\hat{b}) = ReLUsol(b) = sol(b)$, and the change to $\hat{b}$ will be the full $\Delta(b)$, unless $\Delta(b) < -sol(b) < 0$ in which case it is $-sol(b)$. If (c) $sol(b) < 0$, then we have $sol(\hat{b}) = ReLU(b) = 0$ and opening $a$ moves away from 0 only if $sol(b) + \Delta(b) > 0$. □

# 7 Experimental Evaluation

We implemented Hybrid MILP in Python 3.8: it first calls $\alpha, \beta$-CROWN with small time out (10s for small and 30s for larger DNNs), and then call pMILP on the undecided inputs. Gurobi 9.52 was used for solving LP and MILP problems. We conducted our evaluation on an AMD Threadripper 7970X (32 cores@4.0GHz) with 256 GB of main memory and 2 NVIDIA RTX 4090. The objectives of our evaluation was to answer the following questions:

1. How does the choice of the set $X$ impacts the accuracy of MILP$_X$?

2. How accurate is Hybrid MILP, and how efficient is it?

## 7.1 Comparison between different scoring functions (accuracy).

To measure the impact of the scoring function to select neurons to open, we considered the complex CNN-B-Adv, which has 16634 nodes. We tested over the $x$ = 85th image in the CIFAR-10 dataset. To measure the accuracy, we measure the uncertainty of all nodes in the 4th layer: the uncertainty of a node is the range between its computed lower and upper bound. We then average the uncertainty among all the nodes of the layer. Formally, the uncertainty of a node
<figure><img src="image_5.png" /><figcaption>Fig. 4: Average uncertainty of pMILP for neurons of the fourth layer of CNN-B-Adv, when selecting between 0 and 30 ReLUs with different scoring functions.</figcaption></figure>

a with bounds [LB, UB] is $uncert(a) = UB(a) - LB(a)$. The average uncertainty of layer $\ell$ is $\frac{\sum_{a \in \ell} \text{uncert}(a)}{\text{size}(\ell)}$.

We report in Figure 4 the average uncertainty of MILP$_X$ following the choice of the K heaviest neurons of SAS, compared with a random choice, with Huang [16], based on strength($n$) = $(UB(n) - LB(n)) \cdot |W_{nz}|$, and with GS (here, $s_{FSB}$).

Overall, SAS is more accurate than GS, more accurate than Huang, more accurate than random choice of variables. SAS significantly outperforms other solutions, with 2 times less binary variables for the same accuracy (10 vs 22) vs GS and 6 times vs Huang [16] (4 vs 24). Each node of the 4th layer displays a similar pattern, and the pattern is similar for different images.

Gurobi relies on a Branch and Bound procedure to compute bounds. We report in Table 3 experiments on changing the ordering on variables in Gurobi, when the selection of ReLUs is fixed by SAS on CNN-B-Adv Image 85. We compared the standard Gurobi ordering of variables, which is adaptive in each

<table><caption>Table 3: Comparison of different orderings with selection of ReLUs by SAS.</caption><thead><tr><th>Image 85 order</th><th>45 open ReLUs distance</th><th>45 open ReLUs time</th><th>50 open ReLUs distance</th><th>50 open ReLUs time</th><th>55 open ReLUs distance</th><th>55 open ReLUs time</th><th>60 open ReLUs distance</th><th>60 open ReLUs time</th></tr></thead><tbody><tr><td>Gurobi</td><td>0.384</td><td>252s</td><td>0.368</td><td>265s</td><td>0.342</td><td>287s</td><td>0.315</td><td>306s</td></tr><tr><td>SAS (static)</td><td>0.384</td><td>448s</td><td>0.368</td><td>759s</td><td>0.365</td><td>971s</td><td>0.355</td><td>971s</td></tr><tr><td>GS (static)</td><td>0.384</td><td><strong>251s</strong></td><td>0.368</td><td><strong>253s</strong></td><td>0.342</td><td><strong>270s</strong></td><td>0.315</td><td><strong>281s</strong></td></tr></tbody></table>
branch, with the static order provided by SAS, as well as the static order provided by GS, wrt runtime at a fixed accuracy (we report the *distance* to verify the image).

SAS ordering is particularly counter-productive because it is accurate only for one branch (the most complicated one), whereas Gurobi can adapt the order to each branch. However, the GS *ordering* (with the SAS *selection*) is better as it is general and not local to a solution, with better runtime than Gurobi, despite its staticness whereas Gurobi order is adaptative.

## 7.2 Comparison with MILP (time and accuracy)

Restricting the set $X$ of open ReLU nodes potentially hurts accuracy. Another strategy could be to still compute inductively accurate bounds for each nodes, replacing pMILP (SAS, GS, etc.) with a full MILP model, and stopping it early after some fixed time-out (100s, 250s, ..., 2000s).

We evaluate in Figure 5 the full MILP and different pMILP models on the output layer of CNN-B-Adv. The bounds for hidden layers have been computed with the same SAS method. We compare the *distance to verify* the same Image 85, specifically the lower bound of an output neuron, the furthest away from verification. The curve is not as smooth as in Fig. 4, as a unique neuron is considered rather than an average over neurons.

SAS (even $GS = s_{FSB}$) is much more accurate than full MILP for every reasonable time-outs considered, with > 10 times fastest runtime at better accuracy.

<figure><img src="image_6.png" /><figcaption>Fig. 5: Distance to verify vs runtime: comparison between SAS, GS, full MILP and Huang's method for different number of opened ReLUs / time-outs.</figcaption></figure>
The reason is that even the advanced MILP solver Gurobi struggles to sort out all ReLUs. It is thus particulary important to have accurate scoring functions as SAS, in order to optimize both accuracy and runtime. Huang [16] is limited in ReLUs in the previous layer, reason why it is stuck at relatively poor accuracy. The accuracy/runtime curve from GS is closer to SAS than the number of nodes / accuracy curve of Fig. 4. This is because many ReLU nodes deemed important by GS are not relevant for accuracy, but they are also not penalizing runtime. Still, SAS is faster than GS at every accuracy. Importantly, SAS is more deterministic in its runtime given a number of ReLU nodes, with half the variance in runtime over different output neurons compared with GS for similar accuracy, which helps setting a number of ReLU nodes to open.

7.3 Ablation study: on usefulness of computing previous layers
accurately

We explore the usefulness of computing accurately bounds for each neuron in-
ductively on the layers, even on small networks. For that, we consider MNIST
5 × 100, computing bounds for nodes of layer 3, comparing when bounds for
neurons of layer 2 have been computed inaccurately using LP rather than with
the more accurate (partial) MILP.

<figure><img src="image_7.png" /><figcaption>Fig. 6: Comparison of accuracy in layer 3 when bounds for neurons in layer 2 is computed inaccurately using LP vs when bounds of layer 2 are computed accurately using pMILP. Time is using logscale.</figcaption></figure>
This experiment explains the rationale to use divide and conquer methodology, using many calls (one for each neuron) of pMILP with relatively small number $|X|$ of open ReLUs rather than few calls (one per output neuron) of pMILP with larger number $|X|$ of open nodes. The benefit is clear already from layer 3, obtaining much tighter bounds (lower uncertainty) when bounds for neurons in layer 2 have been computed accurately using pMILP.

## 7.4 Comparison with α, β-CROWN

To assess the verification efficiency and runtime vs α, β-CROWN, we conducted evaluations on neural networks tested in [27] which are mostly *complex* to verify (as easy instances are already appropriately taken care of). Namely, 6 ReLU-DNNs: 5 MNIST DNN that can be found in the ERAN GitHub (the 4th to the 8th DNNs provided) as well as 1 CIFAR CNN from [3], see also [9], which can be downloaded from the α, β-CROWN GitHub. We commit to the same $\epsilon$ settings as in [27], that are recalled in Table 1. For reference, we also report an easy but very large ResNet Network for CIFAR10, already tested with α, β CROWN. We report in Table 4 the % of undecided images, that is the % of images than can be neither falsified (by α, β-CROWN) nor verified by the tested verifier, among the 100 first images for each MNIST or CIFAR10 benchmark. The exact same DNNs, image set and $\epsilon$ are used in Tables 1 and 4.

*Analysis:* overall, Hybrid MILP is very accurate, only leaving 8%-15% of images undecided, with runtime taking less than 500s in average per image, and even 10 times less on smaller DNNs. It can scale up to quite large hard DNNs, such as CNN-B-Adv with 2M parameters.

Compared with α, β-CROWN: on easy instances, Hybrid MILP is virtually similar to α, β-CROWN (e.g. even on the very large ResNet), since Hybrid MILP calls first α, β-CROWN as it is very efficient for easy instances.

On hard instances (the 6 first DNNs tested), compared with α, β-CROWN with a time-out of TO=2000s, Hybrid MILP is much more accurate, with a

<table><caption>Table 4: Undecided images (%, *lower is better*) as computed by α, β-CROWN, Refined β-CROWN and Hybrid MILP on 7 DNNs (average runtime per image). The 6 first DNNs are hard instances. The last DNN (ResNet) is an easy instance (trained using Wong to be easy to verify), provided for reference.</caption><thead><tr><th>Network</th><th>α, β-CROWN<br>TO=10s</th><th>α, β-CROWN<br>TO=30s</th><th>α, β-CROWN<br>TO=2000s</th><th>Refined<br>β-CROWN</th><th>Hybrid<br>MILP</th></tr></thead><tbody><tr><td>MNIST 5×100</td><td>57% (6.9s)</td><td>55% (18.9s)</td><td>50% (1026s)</td><td>13% (92s)</td><td>13% (46s)</td></tr><tr><td>MNIST 5×200</td><td>50% (6.5s)</td><td>47% (17s)</td><td>46% (930s)</td><td>9% (80s)</td><td>8% (71s)</td></tr><tr><td>MNIST 8×100</td><td>63% (7.2s)</td><td>58% (20s)</td><td>58% (1163s)</td><td>21% (102s)</td><td>15% (61s)</td></tr><tr><td>MNIST 8×200</td><td>56% (6.8s)</td><td>55% (18s)</td><td>54% (1083s)</td><td>16% (83s)</td><td>8% (78s)</td></tr><tr><td>MNIST 6×500</td><td>53% (6.4s)</td><td>51% (16s)</td><td>50% (1002s)</td><td>—</td><td>10% (402s)</td></tr><tr><td>CNN-B-Adv</td><td>28% (4.3s)</td><td>22% (8.7s)</td><td>20% (373s)</td><td>—</td><td>11% (417s)</td></tr><tr><td>ResNet</td><td>0% (2s)</td><td>0% (2s)</td><td>0% (2s)</td><td>—</td><td>0% (2s)</td></tr></tbody></table>
reduction of undecided images by 9% - 43%. It is also from 20x faster on smaller networks to similar time on the largest DNN. Compared with α, β-CROWN with a time-out of TO=30s, the accuracy gap is even larger (e.g. 11% for CNN-B-Adv, i.e. half the undecided images), although the average runtime is also obviously larger (solving hard instances takes longer than solving easy instances).

Last, compared with *Refined β-CROWN*, we can observe three patterns: on the shallowest DNNs (5×100, 5×200), *Refined β-CROWN* can run full MILP on almost all nodes, reaching almost the same accuracy than Hybrid MILP, but with longer runtime (up to 2 times on 5×100). As size of DNNs grows (8×100, 8×200), full MILP invoked by *Refined β-CROWN* can only be run on a fraction of the neurons, and the accuracy is not as good as Hybrid MILP, with 6% - 8% more undecided images (that is double on 8×200), while having longer runtime. Last but not least, *Refined β-CROWN* cannot scale to larger instances (6×500, CNN-B-Adv), while Hybrid MILP can.

## 7.5 Comparison with other Verifiers

We voluntarily limited the comparison so far to α, β-CROWN because it is one of the most efficient verifier to date, which allowed us to consider a spectrum of parameters to understand α, β-CROWN scaling without too much clutter.

Interestingly, GCP-CROWN [32] is slightly more accurate than α, β-CROWN on the DNNs we tested, but necessitates IBM CPLEX solver, which is not available to us. We provide in Table 5 results from the Table 3 page 9 of [32], experimenting with different verifiers on the most interesting CIFAR CNN-B-Adv. Notice that results are not exactly comparable with ours because the testing images are not the same (ours has 62% upper bound while [32] has 65% upper bound, so the image set in [32] are slightly easier). There, GCP-CROWN is 2% more accurate (with higher runtime) than α, β-CROWN with 90s TO, so we can deduce that Hybrid MILP, which is 9% more accurate than α, β-CROWN with 2000s TO, is significantly more accurate than GCP-CROWN.

Results on other networks from other verifiers (PRIMA [19], SDP-FO [9], etc) were already reported [27] on other tested DNNs, with unfavorable comparison vs α, β-CROWN. Further, we reported accuracy of NNenum [2], Marabou [18,29], respectively 4th, 3rd of the last VNNcomp'24 [5], as well as full MILP [26] in Table 2, showing that these verifiers are not competitive on complex (even small) instances. Concerning MnBAB [15], and it compares slightly unfavorably in time and accuracy towards α, β-CROWN on CNN-B-Adv and *complex* MNIST DNNs

<table><caption>Table 5: Images verified by different verifiers (%), higher is better) on CIFAR CNN-B-Adv. Results, from [32], are not fully comparable with Tables 1,4.</caption><thead><tr><th>Upper Bound</th><th>SDP-FO</th><th>PRIMA</th><th>refined β-CROWN</th><th>β-CROWN</th><th>GCP-CROWN</th></tr></thead><tbody><tr><td>65%</td><td>32.8%</td><td>38%</td><td>27%</td><td>46.5%</td><td>48.5%</td></tr><tr><td>ε = 2/255</td><td>&gt; 25h</td><td>344s</td><td>361s</td><td>32s</td><td>58s</td></tr></tbody></table>
at several time-out settings. Last, Pyrat [12] (2nd in the latest VNNComp) is
not open source, which made running it impossible.

8 Conclusion

In this paper, we developed a novel solution-aware scoring (SAS) function to select few ReLU nodes to consider with binary variables to compute accurately bounds in DNNs. The solution awareness allows SAS to compute an accurate score for each ReLU, which enables partial MILP to be very efficient, necessitating ≈ 6x less binary variables than previous proposals [16] for the same accuracy, and ≈ 2x less than GS scoring adapted from FSB [11]. As the worst-case complexity is exponential in the number of binary variables, this has large implication in terms of scalability to larger DNNs, making it possible to verify accurately quite large DNNs such as CNN-B-Adv with 2M parameters.

While α, β-CROWN is known to be extremely efficient to solve easier verifica-
tion instances, we exhibit many cases (complex instances) where its worst-case
exponential complexity in the number of ReLUs is tangible, with unfavorable
scaling (Table 1). Resorting to Hybrid MILP, a divide-and-conquer approach
[16], revisited thanks to the very efficient SAS, revealed to be a much better
trade-off than augmenting α, β-CROWN time-outs, with 8% to 40% less unde-
cided images at ISO runtime. Currently, for hard instances, there is no alternative
to partial MILP, other methods being > 10 times slower.

This opens up interesting future research directions, to verify global [28], [7]
rather than local (robustness) properties, which need very accurate methodology
and give rise to hard instances as the range of each neuron is no more local to a
narrow neighborhood (and thus almost all ReLUs are unstable, with both modes
possible).

**Acknowledgement:** This research was conducted as part of the DesCartes program and was supported by the National Research Foundation, Prime Minister's Office, Singapore, under the Campus for Research Excellence and Technological Enterprise (CREATE) program, and partially supported by ANR-23-PEIA-0006 SAIF.

References

1. Mohammad Afzal, Ashutosh Gupta, and S. Akshay. Using Counterexamples to Improve Robustness Verification in Neural Networks. In *Automated Technology for Verification and Analysis (ATVA'23)*, LNCS 14215, pages 422–443, 2023.

2. Stanley Bak. nnenum: Verification of relu neural networks with optimized abstraction refinement. In *NASA Formal Methods Symposium*, pages 19–36. Springer, 2021.

3. Mislav Balunovic and Martin Vechev. Adversarial training and provable defenses: Bridging the gap. In *International Conference on Learning Representations (ICLR'20)*, 2020.
4. Christopher Brix, Stanley Bak, Changliu Liu, and Taylor T. Johnson. The fourth international verification of neural networks competition (vnn-comp 2023): Summary and results, 2023.

5. Christopher Brix, Stanley Bak, Changliu Liu, Taylor T. Johnson, David Shriver, and Haoze (Andrew) Wu. 5th international verification of neural networks competition (vnn-comp'24), 2024.

6. Christopher Brix, Mark Niklas Müller, Stanley Bak, Taylor T. Johnson, and Changliu Liu. First three years of the international verification of neural networks competition (vnn-comp), 2023.

7. Rudy Bunel, Krishnamurthy Dvijotham, M. Pawan Kumar, Alessandro De Palma, and Robert Stanforth. Verified neural compressed sensing, 2024.

8. Rudy Bunel, Jingyue Lu, Ilker Turkaslan, Philip HS Torr, Pushmeet Kohli, and M Pawan Kumar. Branch and bound for piecewise linear neural network verification. *Journal of Machine Learning Research*, 21(42):1–39, 2020.

9. Sumanth Dathathri, Krishnamurthy Dvijotham, Alexey Kurakin, Aditi Raghunathan, Jonathan Uesato, Rudy R Bunel, Shreya Shankar, Jacob Steinhardt, Ian Goodfellow, Percy S Liang, et al. Enabling certification of verification-agnostic networks via memory-efficient semidefinite programming. *Advances in Neural Information Processing Systems*, 33:5318–5331, 2020.

10. Alessandro De Palma, Harkirat S Behl, Rudy Bunel, Philip Torr, and M Pawan Kumar. Scaling the convex barrier with active sets. In *International Conference on Learning Representations (ICLR’21)*. Open Review, 2021.

11. Alessandro De Palma, Rudy Bunel, Alban Desmaison, Krishnamurthy Dvijotham, Pushmeet Kohli, Philip HS Torr, and M Pawan Kumar. Improved branch and bound for neural network verification via lagrangian decomposition. *arXiv preprint arXiv:2104.06718*, 2021.

12. Serge Durand, Augustin Lemesle, Zakaria Chihani, Caterina Urban, and Francois Terrier. Reciph: Relational coefficients for input partitioning heuristic. In *1st Workshop on Formal Verification of Machine Learning (WFVML 2022)*, 2022.

13. Rüdiger Ehlers. Formal verification of piece-wise linear feed-forward neural networks. In Deepak D’Souza and K. Narayan Kumar, editors, *Automated Technology for Verification and Analysis - 15th International Symposium, ATVA 2017, Pune, India, October 3-6, 2017, Proceedings*, volume 10482 of Lecture Notes in Computer Science, pages 269–286. Springer, 2017.

14. Ruediger Ehlers. Formal verification of piece-wise linear feed-forward neural networks. In *Automated Technology for Verification and Analysis: 15th International Symposium, ATVA 2017, Pune, India, October 3-6, 2017, Proceedings 15*, pages 269–286. Springer, 2017.

15. Claudio Ferrari, Mark Niklas Mueller, Nikola Jovanović, and Martin Vechev. Complete verification via multi-neuron relaxation guided branch-and-bound. In *International Conference on Learning Representations (ICLR’22)*, 2022.

16. Chao Huang, Jiameng Fan, Xin Chen, Wenchao Li, and Qi Zhu. Divide and slide: Layer-wise refinement for output range analysis of deep neural networks. *IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems*, 39(11):3323–3335, 2020.

17. Guy Katz, Clark Barrett, David L. Dill, Kyle Julian, and Mykel J. Kochenderfer. Reluplex: An efficient smt solver for verifying deep neural networks. In Rupak Majumdar and Viktor Kunčak, editors, *Computer Aided Verification*, pages 97–117, Cham, 2017.
18. Guy Katz, Derek A. Huang, Duligur Ibeling, Kyle Julian, Christopher Lazarus, Rachel Lim, Parth Shah, Shantanu Thakoor, Haoze Wu, Aleksandar Zeljić, David L. Dill, Mykel J. Kochenderfer, and Clark Barrett. The marabou framework for verification and analysis of deep neural networks. In *Isil Dillig and Serdar Tasiran, editors, *Computer Aided Verification*, pages 443–452, Cham, 2019. Springer International Publishing.

19. Mark Niklas Muller, Gleb Makarchuk, Gagandeep Singh, Markus Püschel, and Martin Vechev. Prima: General and precise neural network certification via scalable convex hull approximations. volume 6, New York, NY, USA, jan 2022. Association for Computing Machinery.

20. Mark Niklas Müller, Christopher Brix, Stanley Bak, Changliu Liu, and Taylor T. Johnson. The third international verification of neural networks competition (vnn-comp 2022): Summary and results, 2022.

21. Gagandeep Singh, Rupanshu Ganvir, Markus Püschel, and Martin Vechev. Beyond the single neuron convex barrier for neural network certification. *Advances in Neural Information Processing Systems*, 32, 2019.

22. Gagandeep Singh, Timon Gehr, Markus Püschel, and Martin Vechev. An abstract domain for certifying neural networks. *Proc. ACM Program. Lang.*, 3(POPL), jan 2019.

23. Gagandeep Singh, Timon Gehr, Markus Püschel, and Martin Vechev. Robustness certification with refinement. In *International Conference on Learning Representations (ICLR'19)*, 2019.

24. Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus. Intriguing properties of neural networks. *International Conference on Learning Representations (ICLR'14)*, 2014.

25. Christian Tjandraatmadja, Ross Anderson, Joey Huchette, Will Ma, Krunal Kishor Patel, and Juan Pablo Vielma. The convex relaxation barrier, revisited: Tightened single-neuron relaxations for neural network verification. *Advances in Neural Information Processing Systems*, 33:21675–21686, 2020.

26. Vincent Tjeng, Kai Xiao, and Russ Tedrake. Evaluating robustness of neural networks with mixed integer programming. *International Conference on Learning Representations (ICLR'19)*, 2019.

27. Shiqi Wang, Huan Zhang, Kaidi Xu, Xue Lin, Suman Jana, Cho-Jui Hsieh, and J. Zico Kolter. Beta-crown: Efficient bound propagation with per-neuron split constraints for neural network robustness verification. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, *Advances in Neural Information Processing Systems*, volume 34, pages 29909–29921. Curran Associates, Inc., 2021.

28. Zhilu Wang, Chao Huang, and Qi Zhu. Efficient global robustness certification of neural networks via interleaving twin-network encoding. In *2022 Design, Automation & Test in Europe Conference & Exhibition (DATE)*, pages 1087–1092. IEEE, 2022.

29. Haoze Wu, Omri Isac, Aleksandar Zeljić, Teruhiro Tagomori, Matthew Daggitt, Wen Kokke, Idan Refaeli, Guy Amir, Kyle Julian, Shahaf Bassan, Pei Huang, Ori Lahav, Min Wu, Min Zhang, Ekaterina Komendantskaya, Guy Katz, and Clark Barrett. Marabou 2.0: A versatile formal analyzer of neural networks. In *Arie Gurfinkel and Vijay Ganesh, editors, *Proceedings of the 36th International Conference on Computer Aided Verification (CAV '24)*, volume 14681 of *Lecture Notes in Computer Science*, pages 249–264. Springer, July 2024. Montreal, Canada.
30. Dong Xu, Nusrat Jahan Mozumder, Hai Duong, and Matthew B. Dwyer. Training for verification: Increasing neuron stability to scale dnn verification. In Bernd Finkbeiner and Laura Kovács, editors, *Tools and Algorithms for the Construction and Analysis of Systems*, pages 24–44, Cham, 2024. Springer Nature Switzerland.

31. Kaidi Xu, Huan Zhang, Shiqi Wang, Yihan Wang, Suman Jana, Xue Lin, and Cho-Jui Hsieh. Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers. In *International Conference on Learning Representations (ICLR'21)*. OpenReview.net, 2021.

32. Huan Zhang, Shiqi Wang, Kaidi Xu, Linyi Li, Bo Li, Suman Jana, Cho-Jui Hsieh, and J. Zico Kolter. General cutting planes for bound-propagation-based neural network verification. In *Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022, New Orleans, LA, USA, November 28 - December 9, 2022, 2022*.
# Appendix

## A Parameter settings

### Setting for Hybrid MILP

Hybird MILP first call $\alpha$, $\beta$-CROWN with short time-out (TO), then call partial MILP on those inputs which was neither certified nor falsified by this run of $\alpha$, $\beta$-CROWN. We are using two settings of TO, for smaller DNNs we use TO= 10s, and for the two larger ones, we use TO= 30s.

Partial MILP uses 20 CPU-threads, while $\alpha$, $\beta$-CROWN uses massively parallel (>4096 threads) GPU,

The setting for partial MILP for fully-connected DNNs is about how many neurons need to be opened (once set, the selection is automatic). The runtime depending crucially upon the number of open ReLU neurons, we set it quite tightly, only allowing few neuron deviation to accommodate to a particularly accurate/inaccurate bound computation (measured by the weight of the remaining SAS score). As complexity increases with the layer considered, as the size of the MILP model grows, we lower this number with the depth, only committing to an intermediate number for the output neuron (the number of output neurons is smaller than hidden layer, and this is the most important computation). We experimentally set this number so that each computing the bounds in each hidden layer takes around the same time. Remember that in layer 1, partial MILP is not necessary and propagating bounds using interval arithmetic is already exact. We open [48,48] to compute bounds for hidden layer 2, [21,24] for layer 3, [11,14] for layer 4, [6,9] for layer 5, [3,6] for layer 6, [2,5] for layer 7, [1,4] for hidden layer 8 (if any), and we open [14,17] for the output layer. The exact number of open nodes in the range $[a,a+3]$ is decided automatically for each neuron being computed: ReLUs are ranked according to their value by SAS, and the top ReLUs are open. Then, ReLUs ranked $a+1, a+2, a+3$ are opened if their SAS value is larger than a small threshold. We set the threshold at 0.01. It should be seen as a way to save runtime when SAS knows that the next node by ranking ($a+i$) will not impact accuracy much (thanks to the upper bound from Proposition 2).

For convolutional CNNs, the strategy is adapted, as there is much more neurons, but in a shallower architecture and not fully connected. The second

<table><caption>Table 6: Settings of Hybrid MILP for the different *hard* instances</caption><thead><tr><td>Network</td><td>TO for α, β-CROWN</td><td>Minimum number of Open neurons</td></tr></thead><tbody><tr><td>MNIST 5 × 100</td><td>10s</td><td>48,21,11,6,14</td></tr><tr><td>MNIST 5 × 200</td><td>10s</td><td>48,21,11,6,14</td></tr><tr><td>MNIST 8 × 100</td><td>10s</td><td>48,21,11,6,3,2,1,14</td></tr><tr><td>MNIST 8 × 200</td><td>10s</td><td>48,21,11,6,3,2,1,14</td></tr><tr><td>MNIST 6 × 500</td><td>30s</td><td>48,21,11,6,3,14</td></tr><tr><td>CIFAR CNN-B-Adv</td><td>30s</td><td>200, 0, 45</td></tr></tbody></table>
layer is computed accurately, opening 200 neurons, which is manageable as there is only one ReLU layer to consider, and accuracy here is crucial. We do not open any nodes in the third layer (the first fully connected layer) if the output layer is the next one (which is the case for CNN-B-Adv), and instead rely on the choice of important nodes for the output layer. Otherwise, we open 20 neurons. In the output layer, we open at least 45 neurons (there is less output neurons than nodes in the previous layer), and enlarge the number of open neurons (up to 300) till we find an upper bound, that is a best current MILP solution, of around +0.1 (this 0.1 was experimentally set as target, a good balance between accuracy and efficiency), and compute a guaranteed lower bound (the goal is to guarantee the bound is > 0).

Table 6 sums up the TO and the minimum numbers of ReLU opened.

Last, for Gurobi, we use a custom MIP-Gap (from 0.001 to 0.1) and time-out parameters, depending on the seen improvement and the possibility to make a node stable. This is low level implementation details that will be available in the code once the paper is accepted.

Notice that a different balance between accuracy and runtime could be set. For instance, we set up the numbers of open neurons to have similar runtime as Refined $\beta$-CROWN for the first 4 DNNs (50s – 100s). We could easily target better accuracy (e.g. for 8 × 100 with a relatively high 15% undecided images) by increasing the number of open neurons, with a trade-off on runtime (current runtime is at 61s). By comparison, the sweet spot for $\alpha, \beta$-CROWN seems to be around TO= 30s, enlarging the time-out having very little impact on accuracy but large impact on runtime (Table 1).

## Setting for $\alpha, \beta$-CROWN

The networks were already tested by $\alpha, \beta$-CROWN [27]. We thus simply reused the parameter files from their Github, except for time-out which we explicitly mention: e.g. for CNN-B-Adv: "solver: batch size: 512 beta-crown: iteration: 20"; and for MNIST 5x100: "solver: batch size: 1024 beta-crown: iteration: 20".

## B Additional experiments including Ablation studies

First, we provide Fig. 7, similar to Fig. 4 but on a different image, namely 37. It displays similar pattern as image 85, with an even larger difference between SAS and GS (around 3x more binary variables are necessary for GS to match the accuracy of SAS).

Then, we provide Table 7 to show the runtime variability of GS and SAS. The number of open nodes for SAS is 14 and the number of open nodes for GS is 30, to account for the lack of accuracy of SAS, making the average distance to verify similar between both (−3.16 vs −3.13). The runtime variability between different output nodes is twice as high for GS vs SAS, because the number of irrelevant nodes chosen by GS is unpredictable. If this number is high, then the runtime can be particularly short (node 8), with also a worse accuracy; while if
<figure><img src="image_8.png" /><figcaption>Fig. 7: Comparison of different methods with image 37.</figcaption></figure>

<table><caption>Table 7: Distance to verify and runtime (lower is better) for different output nodes for image 85. Here, negative distance means the node is verified.</caption><thead><tr><th>Output Node</th><th>0</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>Avg</th><th>Var</th></tr></thead><tbody><tr><td>SAS distance</td><td>-1.47</td><td>-4.56</td><td>-4.48</td><td>-3.46</td><td>-5.6</td><td>-3.06</td><td>-4.17</td><td>0.639</td><td>-2.33</td><td>-3.16</td><td></td></tr><tr><td>Time (s)</td><td>66.9</td><td>54.8</td><td>69.9</td><td>88.8</td><td>94.8</td><td>63.3</td><td>73.4</td><td>56.7</td><td>64.6</td><td>70.3</td><td>184</td></tr><tr><td>GS distance</td><td>-1.38</td><td>-4.54</td><td>-4.45</td><td>-3.48</td><td>-5.7</td><td>-2.98</td><td>-4.08</td><td>0.726</td><td>-2.28</td><td>-3.13</td><td></td></tr><tr><td>Time (s)</td><td>71.6</td><td>84.9</td><td>71.2</td><td>113</td><td>93.7</td><td>68.5</td><td>88.6</td><td>46.5</td><td>69.4</td><td>78.6</td><td>359</td></tr></tbody></table>

it is low (e.g. 4), then the runtime is particularly long. This does not necessarily translates into better accuracy, e.g. node 2 where the runtime is much longer for GS vs SAS, while the accuracy is also lower. SAS displays much less variability in time between different output nodes as most nodes are relevant.

Further, we consider ablation studies to understand how each feature enables the efficiency of pMILP.

## Time scaling with open nodes

First, we explore the time scaling with different number of open nodes, for SAS using nodes in the last two layers (Layer 1 and 2) wrt nodes of layer 3 of 5 × 100 on image 59 of MNIST, presented in Table 8 and Fig. 8.

The exponential complexity with the number of nodes can be seen on Figure 8, where time is represented using logarithmic scale. The flat area in the middle
<table><caption>Table 8: Time and uncertainty scaling of pMILP with number of nodes.</caption><thead><tr><th>X</th><th>Time</th><th>Uncertainty</th></tr></thead><tbody><tr><td>0</td><td>2.6</td><td>1.760946128</td></tr><tr><td>1</td><td>7.3</td><td>1.702986873</td></tr><tr><td>2</td><td>11.1</td><td>1.65469034</td></tr><tr><td>3</td><td>16.3</td><td>1.612137282</td></tr><tr><td>4</td><td>15.5</td><td>1.571001109</td></tr><tr><td>5</td><td>15.7</td><td>1.531925404</td></tr><tr><td>6</td><td>15.8</td><td>1.49535638</td></tr><tr><td>7</td><td>16.4</td><td>1.46189314</td></tr><tr><td>8</td><td>15.8</td><td>1.4299535</td></tr><tr><td>9</td><td>17.2</td><td>1.4006364</td></tr><tr><td>10</td><td>22.5</td><td>1.3711203</td></tr><tr><td>11</td><td>27.2</td><td>1.3438245</td></tr><tr><td>12</td><td>21.6</td><td>1.3183356</td></tr><tr><td>13</td><td>28.7</td><td>1.2938690</td></tr><tr><td>14</td><td>29.6</td><td>1.2690507</td></tr><tr><td>15</td><td>24.5</td><td>1.2475106</td></tr></tbody></table><table><thead><tr><th>X</th><th>Time</th><th>Uncertainty</th></tr></thead><tbody><tr><td>16</td><td>31.9</td><td>1.2243065</td></tr><tr><td>17</td><td>28.6</td><td>1.2031791</td></tr><tr><td>18</td><td>30.4</td><td>1.1839474</td></tr><tr><td>19</td><td>34.0</td><td>1.1644653</td></tr><tr><td>20</td><td>42.1</td><td>1.1456181</td></tr><tr><td>21</td><td>47.6</td><td>1.1261252</td></tr><tr><td>22</td><td>62.7</td><td>1.1089745</td></tr><tr><td>23</td><td>70.0</td><td>1.0931242</td></tr><tr><td>24</td><td>70.8</td><td>1.0773088</td></tr><tr><td>25</td><td>139.9</td><td>1.060928</td></tr><tr><td>26</td><td>154.2</td><td>1.045715</td></tr><tr><td>27</td><td>213.1</td><td>1.030605</td></tr><tr><td>28</td><td>211.3</td><td>1.016058</td></tr><tr><td>29</td><td>373.1</td><td>1.001374</td></tr><tr><td>max=116</td><td>3300</td><td>0.895</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Time and uncertainty scaling of pMILP with number of nodes. Time is using logscale.</figcaption></figure>

is Gurobi having good heuristic to avoid considering all $2^K$ cases when $K < 21$
is not too large, but not working so well for $K > 25$. Notice that when certifying, pMILP uses $|X| \in 21-24$, which is a good trade off between time and accuracy.

**restricting number of open nodes (pMILP) vs setting time-outs (full MILP)** Running full MILP till a small MIP-Gap (typically 0.001) is reached is extremely time inefficient.

Instead, the standard strategy is to set a reasonable time-out and use what-ever bound has been generated. We compare this standard strategy with the pMILP strategy of setting a priori a number of open nodes.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Comparison of uncertainty at layer 7 for full MILP with different time-outs vs pMILP with different number of open nodes. Time is using logscale.</figcaption></figure>

pMILP obtains 2.8 accuracy in < 21 seconds (with 7 open nodes), while full MILP needs 400 seconds to obtain it, a 19x speed up. For 2.7 accuracy, the speedup is >> 22.

Figure 9 shows that choosing nodes is much more efficient for time/accuracy trade-off than setting time outs and use full MILP. And this is for the smallest DNN we considered (500 hidden neurons, far from the biggest 20k neuron DNN we experimented with)
<figure><img src="image_11.png" /><figcaption>Table 9: Comparison of bounding the number of nodes for pMILP and using different time outs for full MILP. In both settings, lower and upper bounds of previous layers are the same (computed by pMILP).</figcaption></figure>

## C Comparison with other DNN verifiers

In the following, we provide results comparing α, β-CROWN to other verifiers, to justify our use of α, β-CROWN as state of the art for efficient verifiers as main source of comparison to Hybrid MILP for hard DNN instance.

### Comparison α, β-CROWN vs PRIMA

<table><caption>Table 10: Undecided images (%, lower is better), as computed by α, β-CROWN, Refined β-CROWN, and PRIMA, as reported in [27], except for 6 × 500 that we run ourselves. N/A means that [27] did not report the numbers, while — means that Refined β-CROWN cannot be run on these DNNs.</caption><thead><tr><th>Network</th><th>α, β-CROWN</th><th>Refined β-CROWN</th><th>PRIMA</th></tr></thead><tbody><tr><td>MNIST 5 × 100</td><td>N/A</td><td>14.3% (102s)</td><td>33.2% (159s)</td></tr><tr><td>MNIST 5 × 200</td><td>N/A</td><td>13.7% (86s)</td><td>21.1% (224s)</td></tr><tr><td>MNIST 8 × 100</td><td>N/A</td><td>20.0% (103s)</td><td>39.2% (301s)</td></tr><tr><td>MNIST 8 × 200</td><td>N/A</td><td>17.6% (95s)</td><td>28.7% (395s)</td></tr><tr><td>MNIST 6 × 500</td><td>51% (16s)</td><td>—</td><td>64% (117s)</td></tr><tr><td>CIFAR CNN-B-Adv</td><td>18.5% (32s)</td><td>—</td><td>27% (344s)</td></tr><tr><td>CIFAR ResNet</td><td>0% (2s)</td><td>—</td><td>0% (2s)</td></tr></tbody></table>

Most data is directly from [27]. N/A means no data either in [27] or by our running.

<footnote>* The data in this row is from our own running on first 100 images of the MNIST dataset.</footnote>

<footnote>** The data is from [27] on first 200 images of the CIFAR10 dataset.</footnote>
PRIMA [19] is a major verifier in the ERAN toolkit. In Table 10, we report the comparison between PRIMA and α, β-CROWN, mainly from [27]. The setting is mainly similar from ours, but numbers are not perfectly comparable as the images tested are not exactly the same (1000 first or 200 first images for CNN-B-Adv), vs 100 first in Tables 4, 1. Also, time-out settings and hardware are slightly different. The overall picture is anyway the same.

Analysis: On the 4 smallest MNIST networks, PRIMA uses a refined path comparable with Refined β-CROWN. However, it is slower and less accurate than Refined β-CROWN. On larger *hard* networks, PRIMA has also more undecided images than α, β-CROWN, while the runtime is > 5 times larger. Hence, Hybrid MILP is more accurate than PRIMA with similar runtime or faster.

Notice that kPoly [21], OptC2V [25], SDP-FO [9] numbers were also reported in [27] on these networks, with even more unfavorable results.

## Comparison α, β-CROWN vs MN-BaB

MN-BaB [15] is an improvement built over PRIMA, using a similar Branch and Bound technique as used in α, β-CROWN. Results in [15] are close to those of α, β-CROWN. However, none of the *hard* networks from [27] that we consider have been tested. We thus tested three representative *hard* DNNs (first 100 images) to understand how MN-BaB fairs on such hard instances, and report the numbers in Table 11. Results are directly comparable with Table 4.

<table><caption>Table 11: Undecided images (%, lower is better), as computed by α, β-CROWN, and MN-BaB</caption><thead><tr><th rowspan="2">Network</th><th colspan="2">α, β-CROWN</th><th colspan="2">MN-BaB</th></tr><tr><th>TO=30s</th><th>TO=2000s</th><th>TO=30s</th><th>TO=2000s</th></tr></thead><tbody><tr><td>MNIST 5 × 100</td><td>55% (19s)</td><td>50% (1026s)</td><td>60% (19s)</td><td>50% (1027s)</td></tr><tr><td>MNIST 6 × 500</td><td>51% (16s)</td><td>50% (1002s)</td><td>58% (18s)</td><td>55% (1036s)</td></tr><tr><td>CIFAR CNN-B-Adv</td><td>22% (8.7s)</td><td>20% (373s)</td><td>43% (14s)</td><td>24% (576s)</td></tr></tbody></table>

Analysis: results reveal that MN-BaB is slightly slower and slightly less accurate than α, β-CROWN. Notice the specially high number of undecided images for CNN-B-Adv with TO=30s, probably meaning that 30s is too small for MN-BaB on this large DNN. Hence, Hybrid MILP is more accurate than MN-BaB with similar runtime or faster.

## Comparison α, β-CROWN vs NNenum

NNenum [2] is a complete verifier with good performance according to VN-Ncomp. It was the only complete verifier tested in Table 2 to verify more images than α, β-CROWN. The experiments section in [2] does not report the *hard* DNNs we are considering. We tried to experiment it on the same MNIST 6 × 500
and CIFAR CNN-B-Adv as we did in Table 11 for MN-BaB. Unfortunately, on 6 × 500, buffer overflow were reported. We report in Table 12 experiments with the same 2000s Time-out (it was 10000s in Table 2) for a fair comparison with α, β-CROWN, on both MNIST 5 × 100 and CIFAR CNN-B-Adv. On MNIST 5 × 100, NNenum is slightly more accurate than α, β-CROWN, but far from the accuracy Hybrid MILP. On CIFAR CNN-B-Adv, NNenum was much less accurate than α, β-CROWN, and thus of Hybrid MILP. In both test, the runtime of NNenum was also much longer than for Hybrid MILP.

<table><caption>Table 12: Undecided images (%, lower is better), as computed by α, β-CROWN and NNenum with 2000s time-out, and Hybrid MILP</caption><thead><tr><th>Network</th><th>α, β-CROWN<br>TO=2000s</th><th>NNenum<br>TO=2000s</th><th>Hybrid<br>MILP</th></tr></thead><tbody><tr><td>MNIST 5 × 100<br>ε = 0.026</td><td>50%(1026s)</td><td>44% (1046s)</td><td><b>13% (46s)</b></td></tr><tr><td>CIFAR CNN-B-Adv<br>ε = 0.015</td><td>20% (373s)</td><td>40% (1020s)</td><td><b>11% (417s)</b></td></tr></tbody></table>

.

## D Average vs max time per pMILP call

We provide in Table 13 the average as well as maximum time to perform MILP<sub>X</sub> calls as called by pMILP, on a given input: image 3 for MNIST, and image 76 for CIFAR10. For 6x500, we provide results for two different ε.

<table><caption>Table 13: average and maximum time per MILP<sub>X</sub> calls for image 3 (MNIST) and image 76 (CIFAR10).</caption><thead><tr><th>Network</th><th>average time</th><th>maximum time</th></tr></thead><tbody><tr><td>MNIST 5×100<br>ε = 0.026</td><td>0.41s</td><td>1.87</td></tr><tr><td>MNIST 5×200<br>ε = 0.015</td><td>0.75s</td><td>5.31s</td></tr><tr><td>MNIST 8×100<br>ε = 0.026</td><td>0.39s</td><td>1.41s</td></tr><tr><td>MNIST 8×200<br>ε = 0.015</td><td>0.49s</td><td>1.63s</td></tr><tr><td>MNIST 6×500<br>ε = 0.035</td><td>1.4s</td><td>3.5s</td></tr><tr><td>ε = 0.1</td><td>44.6s</td><td>310s</td></tr><tr><td>CIFAR CNN-B-Adv<br>ε = 2/255</td><td>1s</td><td>609s</td></tr></tbody></table>
Notice that DNN 6× 500 and $\epsilon = 0.1$ is a very hard instance as being very close to the falsification $\epsilon \approx 0.11$. This is thus not representative of the average case. Also, on this image 3, pMILP succeeds to verify $\epsilon = 1.054$, while $\alpha, \beta$-CROWN can only certify $\epsilon = 0.0467$ within the 10 000s Time-out.

For CNN-B-Adv, the very long maximum time for a MILP call is an outlier: it happens only for one output layer, for which the number $K$ of open nodes is particularly large (around 200 out of 20000 neurons) to certify this hard image 76. Indeed, the average time is at 1s. Notice that this does not lead to a runtime of 20.000s, as 20 threads are used by pMILP in parallel (similar to competing solutions, except $\alpha, \beta$-CROWN which uses > 4096 GPU cores).