Partition-Based Formulations for Mixed-Integer Optimization of Trained ReLU Neural Networks
Year: 2021

Abstract:
This paper introduces a class of mixed-integer formulations for trained ReLU neural networks. The approach balances model size and tightness by partitioning node inputs into a number of groups and forming the convex hull over the partitions via disjunctive programming. At one extreme, one partition per input recovers the convex hull of a node, i.e., the tightest possible formulation for each node. For fewer partitions, we develop smaller relaxations that approximate the convex hull, and show that they outperform existing formulations. Specifically, we propose strategies for partitioning variables based on theoretical motivations and validate these strategies using extensive computational experiments. Furthermore, the proposed scheme complements known algorithmic approaches, e.g., optimizationbased bound tightening captures dependencies within a partition.

heading: 1 Introduction
text: Many applications use mixed-integer linear programming (MILP) to optimize over trained feedforward ReLU neural networks (NNs) [5, 14, 18, 23, 32, 36]. A MILP encoding of a ReLU-NN enables network properties to be rigorously analyzed, e.g., maximizing a neural acquisition function [33] or verifying robustness of an output (often classification) within a restricted input domain [6]. MILP encodings of ReLU-NNS have also been used to determine robust perturbation bounds [9], compress NNs [28], count linear regions [27], and find adversarial examples [17]. The so-called bigM formulation is the main approach for encoding NNs as MILPs in the above references. Optimizing the resulting MILPs remains challenging for large networks, even with state-of-the-art software.
Effectively solving a MILP hinges on the strength of its continuous relaxation [10]; weak relaxations can render MILPs computationally intractable. For NNs, Anderson et al. [3] showed that the big-M formulation is not tight and presented formulations for the convex hull (i.e., the tightest possible formulation) of individual nodes. However, these formulations require either an exponential (w.r.t. inputs) number of constraints or many additional/auxiliary variables. So, despite its weaker continuous relaxation, the big-M formulation can be computationally advantageous owing to its smaller size.
Given these challenges, we present a novel class of MILP formulations for ReLU-NNs. The formulations are hierarchical: their relaxations start at a big-M equivalent and converge to the convex hull. Intermediate formulations can closely approximate the convex hull with many fewer variables/constraints. The formulations are constructed by viewing each ReLU node as a two-part
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
disjunction. Kronqvist et al. [21] proposed hierarchical relaxations for general disjunctive programs. This work develops a similar hierarchy to construct strong and efficient MILP formulations specific to ReLU-NNs. In particular, we partition the inputs of each node into groups and formulate the convex hull over the resulting groups. With fewer groups than inputs, this approach results in MILPs that are smaller than convex-hull formulations, yet have stronger relaxations than big-M.
Three optimization tasks evaluate the new formulations: optimal adversarial examples, robust verification, and `1-minimally distorted adversaries. Extensive computation, including with convexhull-based cuts, shows that our formulations outperform the standard big-M approach with 25% more problems solved within a 1h time limit (average 2.2X speedups for solved problems).
Related work. Techniques for obtaining strong relaxations of ReLU-NNs include linear programming [16, 34, 35], semidefinite programming [12, 25], Lagrangian decomposition [7], combined relaxations [30], and relaxations over multiple nodes [29]. Recently, Tjandraatmadja et al. [31] derived the tightest possible convex relaxation for a single ReLU node by considering its multivariate input space. These relaxation techniques do not exactly represent ReLU-NNs, but rather derive valid bounds for the network in general. These techniques might fail to verify some properties, due to their non-exactness, but they can be much faster than MILP-based techniques. Strong MILP encoding of ReLU-NNs was also studied by Anderson et al. [3], and a related dual algorithm was later presented by De Palma et al. [13]. Our approach is fundamentally different, as it constructs computationally cheaper formulations that approximate the convex hull, and we start by deriving a stronger relaxation instead of strengthening the relaxation via cutting planes. Furthermore, our formulation enables input node dependencies to easily be incorporated.
Contributions of this paper. We present a new class of strong, yet compact, MILP formulations for feed-forward ReLU-NNs in Section 3. Section 3.1 observes how, in conjunction with optimizationbased bound tightening, partitioning input variables can efficiently incorporate dependencies into MILP formulations. Section 3.2 builds on the observations of Kronqvist et al. [21] to prove the hierarchical nature of the proposed ReLU-specific formulations, with relaxations spanning between big-M and convex-hull formulations. Sections 3.3–3.4 show that formulation tightness depends on the specific choice of variable partitions, and we present efficient partitioning strategies based on both theoretical and computational motivations. The advantages of the new formulations are demonstrated via extensive computational experiments in Section 4.

heading: 2 Background

heading: 2.1 Feed-forward Neural Networks
text: A feed-forward neural network (NN) is a directed acyclic graph with nodes structured into k layers. Layer k receives the outputs of nodes in the preceding layer k − 1 as its inputs (layer 0 represents inputs to the NN). Each node in a layer computes a weighted sum of its inputs (known as the preactivation), and applies an activation function. This work considers the ReLU activation function:
y = max(0,wTx + b) (1)
where x ∈ Rη and y ∈ [0,∞) are, respectively, the inputs and output of a node (wTx + b is termed the preactivation). Parameters w ∈ Rη and b ∈ R are the node weights and bias, respectively.

heading: 2.2 ReLU Optimization Formulations
text: In contrast to the training of NNs (where parameters w and b are optimized), optimization over a NN seeks extreme cases for a trained model. Therefore, model parameters (w, b) are fixed, and the inputs/outputs of nodes in the network (x, y) are optimization variables instead.
Big-M Formulation. The ReLU function (1) is commonly modeled [9, 23]:
y ≥ (wTx + b) (2) y ≤ (wTx + b)− (1− σ)LB0 (3) y ≤ σUB0 (4)
where y ∈ [0,∞) is the node output and σ ∈ {0, 1} is a binary variable corresponding to the on-off state of the neuron. The formulation requires the bounds (big-M coefficients) LB0,UB0 ∈ R, which should be as tight as possible, such that (wTx + b) ∈ [LB0,UB0]. Disjunctive Programming [4]. We observe that (1) can be modeled as a disjunctive program:[
y = 0 wTx + b ≤ 0 ] ∨ [ y = wTx + b wTx + b ≥ 0 ] (5)
The extended formulation for disjunctive programs introduces auxiliary variables for each disjunction. Instead of directly applying the standard extended formulation, we first define z := wTx and assume z is bounded. The auxiliary variables za ∈ R and zb ∈ R can then be introduced to model (5):
wTx = za + zb (6) za + σb ≤ 0 (7) zb + (1− σ)b ≥ 0 (8) y = zb + (1− σ)b (9) σLBa ≤ za ≤ σUBa (10) (1− σ)LBb ≤ zb ≤ (1− σ)UBb (11)
where again y ∈ [0,∞) and σ ∈ {0, 1}. Bounds LBa,UBa,LBb,UBb ∈ R are required for this formulation, such that za ∈ σ[LBa,UBa] and zb ∈ (1−σ)[LBb,UBb]. Note that the inequalities in (5) may cause the domains of za and zb to differ. The summation in (6) can be used to eliminate either za or zb; therefore, in practice, only one auxiliary variable is introduced by formulation (6)–(11).
Importance of relaxation strength. MILP is often solved with branch-and-bound, a strategy that bounds the objective function between the best feasible point and its tightest optimal relaxation. The integral search space is explored by “branching” until the gap between bounds falls below a given tolerance. A tighter, or stronger, relaxation can reduce this search tree considerably.

heading: 3 Disaggregated Disjunctions: Between Big-M and the Convex Hull
text: Our proposed formulations split the sum z = wTx into partitions: we will show these formulations have tighter continuous relaxations than (6)–(11). In particular, we partition the set {1, ..., η} into subsets S1 ∪ S2 ∪ ... ∪ SN = {1, ..., η}; Sn ∩ Sn′ = ∅ ∀n 6= n′. An auxiliary variable is then introduced for each partition, i.e., zn = ∑ i∈Sn wixi. Replacing z = w Tx with ∑N n=1 zn, the disjunction (5) becomes:[ y = 0∑N
n=1 zn + b ≤ 0
] ∨ [ y = ∑N n=1 zn + b y ≥ 0 ] (12)
Assuming now that each zn is bounded, the extended formulation can be expressed using auxiliary variables zan and z b n for each zn. Eliminating z a n via the summation zn = z a n + z b n results in our proposed formulation (complete derivation in appendix A.1):∑ n (∑ i∈Sn wixi − zbn ) + σb ≤ 0 (13)
∑ n zbn + (1− σ)b ≥ 0 (14)
y = ∑ n zbn + (1− σ)b (15)
σLBan ≤ ∑ i∈Sn wixi − zbn ≤ σUB a n, ∀n = 1, ..., N (16)
(1− σ)LBbn ≤ zbn ≤ (1− σ)UB b n, ∀n = 1, ..., N (17)
where y ∈ [0,∞) and σ ∈ {0, 1}. We observe that (13)–(17) exactly represents the ReLU node (given that the partitions zn are bounded in the original disjunction): the sum zn = ∑ i∈Sn wixi is
substituted in the extended convex hull formulation [4] of disjunction (12). Note that, compared to the general case presented by Kronqvist et al. [21], both sides of disjunction (12) can be modeled using the same partitions, resulting in fewer auxiliary variables. We further note that domains [LBan,UB a n] and [LBbn,UB b n] may not be equivalent, owing to inequalities in (12).

heading: 3.1 Obtaining and Tightening Bounds
text: The big-M formulation (2)–(4) requires valid bounds (wTx + b) ∈ [LB0,UB0]. Given bounds for each input variable, xi ∈ [x¯ i , x̄i], interval arithmetic gives valid bounds:
LB0 = ∑ i (x ¯ i max(0, wi) + x̄imin(0, wi)) + b (18)
UB0 = ∑ i (x̄imax(0, wi) + x¯ i min(0, wi)) + b (19)
But (18)–(19) do not provide the tightest valid bounds in general, as dependencies between the input nodes are ignored. Propagating the resulting over-approximated bounds from layer to layer leads to increasingly large over-approximations, i.e., propagating weak bounds through layers results in a significantly weakened model. These bounds remain in the proposed formulation (13)–(17) in the form of bounds on both the output y and the original variables x (i.e., outputs of the previous layer).
Optimization-Based Bound Tightening (OBBT) or progressive bounds tightening [32], tightens variable bounds and constraints [15]. For example, solving the optimization problem with the objective set to minimize/maximize (wTx + b) gives bounds: min(wTx + b) ≤ wTx + b ≤ max(wTx + b). To simplify these problems, OBBT can be performed using the relaxed model (i.e., σ ∈ [0, 1] rather than σ ∈ {0, 1}), resulting in a linear program (LP). In contrast to (18)–(19), bounds from OBBT incorporate variable dependencies. We apply OBBT by solving one LP per bound.
The partitioned formulation (13)–(17) requires bounds such that zan ∈ σ[LB a n,UB a n], z b n ∈ (1 − σ)[LBbn,UB b n]. In other words, [LB a n,UB a n] is a valid domain for z a n when the node is inactive
(σ = 1) and vice versa. These bounds can also be obtained via OBBT: min( ∑ i∈Sn wixi) ≤ z a n ≤
max( ∑ i∈Sn wixi);w
Tx + b ≤ 0. The constraint on the right-hand side of disjunction (12) can be similarly enforced in OBBT problems for zbn. In our framework, OBBT additionally captures dependencies within each partition Sn. Specifically, we observe that the partitioned OBBT problems effectively form the convex hull over a given polyhedron of the input variables, in contrast to the convex hull formulation, which only considers the box domain defined by the min/max of each input node [3]. Since ∑ i∈Sn wixi = z a n + z b n and z a nz b n = 0, the bounds [min( ∑ i∈Sn wixi),
max( ∑ i∈Sn wixi)] are valid for both z a n and z b n. These bounds can be from (18)–(19), or by solving two OBBT problems for each partition (2N LPs total). This simplification uses equivalent bounds for zan and z b n, but tighter bounds can potentially be found by performing OBBT for z a n and z b n separately (4N LPs total).
Figure 1 compares the proposed formulation with bounds from interval arithmetic (top row) vs OBBT (bottom row). The true model outputs (blue), and minimum (orange) and maximum (green) outputs of the relaxed model are shown over the input domain. The NNs comprise two inputs, three hidden layers with 20 nodes each, and one output. As expected, OBBT greatly improves relaxation tightness.

heading: 3.2 Tightness of the Proposed Formulation
text: Proposition 1. (13)–(17) has the equivalent non-lifted (i.e., without auxiliary variables) formulation: y ≤ ∑ i∈Ij wixi + σ(b+ ∑ i∈I\Ij UB i) + (σ − 1)( ∑ i∈Ij LB i), ∀j = 1, ..., 2N (20)
y ≥ wTx + b (21) y ≤ σUB0 (22)
where UB i and LB i denote the upper and lower bounds of wixi. The set I denotes the input indices {1, ..., η}, and the subset Ij denotes the indices contained by the union of the j-th combination of partitions in {S1, ...,SN}.
Proof. Formulation (13)–(17) introduces N auxiliary variables zbn, which can be projected out using Fourier-Motzkin elimination (appendix A.2), resulting in combinations I1, ..., IJ , J = 2N .
Remark. When N < η, the family of constraints in (20) represent a subset of the inequalities used to define the convex hull by Anderson et al. [3], where UB i and LB i would be obtained using interval arithmetic. Therefore, though a lifted formulation is proposed in this work, the proposed formulations have non-lifted relaxations equivalent to pre-selecting a subset of the convex hull inequalities.
Proposition 2. Formulation (13)–(17) is equivalent to the big-M formulation (2)–(4) when N = 1 . Proof. When N = 1, it follows that S1 = {1, .., η}. Therefore, z1 = ∑η i=1 wixi = w
Tx, and∑ n zn = z1 = z. Conversely, big-M can be seen as the convex hull over a single aggregated “variable,” z = wTx.
Proposition 3. Formulation (13)–(17) represents the convex hull of (1), given the inputs x are bounded, for the case of N = η.
Proof. When N = η, it follows that Sn = {n},∀n = 1, .., η, and, consequently, zn = wnxn. Since zn are linear transformations of each xn (as are their respective bounds), forming the convex hull over zn recovers the convex hull over x. An extended proof is provided in appendix A.3.
Proposition 4. A formulation with N partitions is strictly tighter than any formulation with (N − 1) partitions that is derived by combining two partitions in the former.
Proof. When combining two partitions, i.e., S′N−1 := SN−1∪SN , constraints in (20) where S′N−1 ⊆ Ij are also obtained by {SN−1,SN} ⊆ Ij . In contrast, those obtained by SN−1 ∨ SN ⊆ Ij cannot be modeled by S′N−1. Since each constraint in (20) is facet-defining [3] and distinct from each other, omissions result in a less tight formulation.
Remark. The convex hull can be formulated with η auxiliary variables using (13)–(17) for N = η, or 2η constraints using (20). While these formulations have tighter relaxations, they can perform worse than big-M due to having more difficult branch-and-bound subproblems. Our proposed formulation balances this tradeoff by introducing a hierarchy of relaxations with increasing tightness and size. The convex hull is created over partitions zn, n = 1, ..., N , rather than the input variables xn, n = 1, ..., η. Therefore, onlyN auxiliary variables are introduced, withN ≤ η. The results in this work are obtained using this lifted formulation (13)–(17); Appendix A.5 compares the computational performance of equivalent non-lifted formulations, i.e., involving 2N constraints.
Figure 1 shows a hierarchy of increasingly tight formulations from N = 1 (equiv. big-M ) to N = 20 (convex hull of each node over a box input domain). The intermediate formulations approximate the convex-hull (N = 20) formulation well, but need fewer auxiliary variables/constraints.

heading: 3.3 Effect of Input Partitioning Choice
text: Formulation (13)–(17) creates the convex hull over zn = ∑ i∈Sn wixi, n = 1, ..., N . The hyperparameter N dictates model size, and the choice of subsets, S1, ...,SN , can strongly impact the strength of its relaxation. By Proposition 4, (13)–(17) can in fact give multiple hierarchies of formulations. While all hierarchies eventually converge to the convex hull, we are primarily interested in those with tight relaxations for small N .
Bounds and Bounds Tightening. Bounds on the partitions play a key role in the proposed formulation. For example, consider when a node is inactive: σ = 1, zbn = 0, and (16) gives the bounds σLBan ≤ ∑ i∈Sn wixi ≤ σUB a n. Intuitively, the proposed formulation represents the
convex hull over the auxiliary variables, zn = ∑ i∈Sn wixi, and their bounds play a key role in model tightness. We hypothesize these bounds are most effective when partitions Sn are selected such that wixi,∀i ∈ Sn are of similar orders of magnitude. Consider for instance the case of w = [1, 1, 100, 100] and xi ∈ [0, 1], i = 1, ..., 4. As all weights are positive, interval arithmetic gives 0 ≤ ∑ xiwi ≤
∑ x̄iwi. With two partitions, the choices of S1 = {1, 2} vs S1 = {1, 3} give:[
x1 + x2 ≤ σ2 x3 + x4 ≤ σ2
] vs. [ x1 + 100x3 ≤ σ101 x2 + 100x4 ≤ σ101 ] (23)
where σ is a binary variable. The constraints on the right closely approximate the η-partition (i.e., convex hull) bounds: x3 ≤ σ and x4 ≤ σ. But x1 and x2 are relatively unaffected by a perturbation σ = 1 − δ (when σ is relaxed). Whereas the formulation on the left constrains the four variables equivalently. If the behavior of the node is dominated by a few inputs, the formulation on the right is strong, as it approximates the convex hull over those inputs (z1 ≈ x3 and z2 ≈ x4 in this case). For the practical case of N << η, there are likely fewer partitions than dominating variables.
The size of the partitions can also be selected to be uneven:[ x1 ≤ σ1
x2 + 100x3 + 100x4 ≤ σ201
] vs. [ x3 ≤ σ1
x1 + x2 + 100x4 ≤ σ102
] (24)
Similar tradeoffs are seen here: the first formulation provides the tightest bound for x1, but x2 is effectively unconstrained and x3, x4 approach “equal treatment.” The second formulation provides the tightest bound for x3, and a tight bound for x4, but x1, x2 are effectively unbounded for fractional σ. The above analyses also apply to the case of OBBT. For the above example, solving a relaxed model for max(x1+x2) obtains a bound that affects the two variables equivalently, while the same procedure for max(x1 + 100x3) obtains a bound that is much stronger for x3 than for x1. Similarly, max(x1 + x2) captures dependency between the two variables, while max(x1 + 100x3) ≈ max(100x3). Relaxation Tightness. The partitions (and their bounds) also directly affect the tightness of the relaxation for the output variable y. The equivalent non-lifted realization (20) reveals the tightness of the above simple example. With two partitions, the choices of S1 = {1, 2} vs S1 = {1, 3} result in the equivalent non-lifted constraints:[
y ≤ x1 + 100x3 + σ(b+ 101) y ≤ x2 + 100x4 + σ(b+ 101)
] vs. [ y ≤ x1 + x2 + σ(b+ 200)
y ≤ 100x3 + 100x4 + σ(b+ 2)
] (25)
Note that combinations Ij = ∅ and Ij = {1, 2, 3, 4} in (20) are not analyzed here, as they correspond to the big-M/1-partition model and are unaffected by choice of partitions. The 4-partition model is the tightest formulation and (in addition to all possible 2-partition constraints) includes:
y ≤ xi + σ(b+ 201), i = 1, 2 (26) y ≤ 100xi + σ(b+ 102), i = 3, 4 (27)
The unbalanced 2-partition, i.e., the left of (25), closely approximates two of the 4-partition (i.e., convex hull) constraints (27). Analogous to the tradeoffs in terms of bound tightness, we see that this formulation essentially neglects the dependence of y on x1, x2 and instead creates the convex hull over z1 = x1 + 100x3 ≈ x3 and z2 ≈ x4. For this simple example, the behavior of y is dominated by x3 and x4, and this turns out to be a relatively strong formulation. However, when N << η, we expect neglecting the dependence of y on some input variables to weaken the model.
The alternative formulation, i.e., the right of (25), handles the four variables similarly and creates the convex hull in two new directions: z1 = x1 + x2 and z2 = 100(x3 + x4). While this does not
model the individual effect of either x3 or x4 on y as closely, it includes dependency between x3 and x4. Furthermore, x1 and x2 are modeled equivalently (i.e., less tightly than individual constraints). Analyzing partitions with unequal size reveals similar tradeoffs. This section shows that the proposed formulation has a relaxation equivalent to selecting a subset of constraints from (20), with a tradeoff between modeling the effect of individual variables well vs the effect of many variables weakly.
Deriving Cutting Planes from Convex Hull. The convex-hull constraints (20) not implicitly modeled by a partition-based formulation can be viewed as potential tightening constraints, or cuts. Anderson et al. [3] provide a linear-time method for selecting the most violated of the exponentially many constraints (20), which is naturally compatible with our proposed formulation (some constraints will be always satisfied). Note that the computational expense can still be significant for practical instances; Botoeva et al. [5] found adding cuts to <0.025% of branch-and-bound nodes to balance their expense and tightening, while De Palma et al. [13] proposed adding cuts at the root node only. Remark. While the above 4-D case may seem to favor “unbalanced” partitions, it is difficult to illustrate the case where N << η. Our experiments confirm “balanced” partitions perform better.

heading: 3.4 Strategies for Selecting Partitions
text: The above rationale motivates selecting partitions that result in a model that treats input variables (approximately) equivalently for the practical case of N << η. Specifically, we seek to evenly distribute tradeoffs in model tightness among inputs. Section 3.3 suggests a reasonable approach is to select partitions such that weights in each are approximately equal (weights are fixed during optimization). We propose two such strategies below, as well as two strategies in red for comparison.
Equal-Size Partitions. One strategy is to create partitions of equal size, i.e., |S1| = |S2| = ... = |SN | (note that they may differ by up to one if η is not a multiple of N ). Indices are then assigned to partitions to keep the weights in each as close as possible. This is accomplished by sorting the weights w and distributing them evenly among partitions (array_split(argsort(w), N ) in Numpy).
Equal-Range Partitions. A second strategy is to partition with equal range, i.e., range i∈S1 (wi) = ... = range i∈SN (wi). We define thresholds v ∈ RN+1 such that v1 and vN+1 are min(w) and max(w). To reduce the effect of outliers, we define v2 and vN as the 0.05 and 0.95 quantiles of w, respectively. The remaining elements of v are distributed evenly in (v2, vN ). Indices i ∈ {1, ..., η} are then assigned to Sn : wi ∈ [vn, vn+1). This strategy requires N ≥ 3, but, for a symmetrically distributed weight vector, w, two partitions of equal size are also of equal range.
We compare our proposed strategies against the following:
Random Partitions. Input indices {1, ..., η} are assigned randomly to partitions S1, ...,SN . Uneven Magnitudes. Weights are sorted and “dealt” in reverse to partitions in snake-draft order.

heading: 4 Experiments
text: Computational experiments were performed using Gurobi v 9.1 [19], and models were implemented and trained using PyTorch [24]. The computational set-up, solver settings, and models for MNIST [22] and CIFAR-10 [20] are described in appendix A.4.

heading: 4.1 Optimal Adversary Results
text: The optimal adversary problem takes a target image x̄, its correct label j, and an adversarial label k, and finds the image within a range of perturbations maximizing the difference in predictions. Mathematically, this is formulated as max
x (fk(x)−fj(x));x ∈ X , where fk and fj correspond to the
k- and j-th elements of the NN output layer, respectively, and X defines the domain of perturbations. We consider perturbations defined by the `1-norm (||x − x̄||1 ≤ 1 ∈ R), which promotes sparse perturbations [8]. For each dataset, we use the first 100 test-dataset images and randomly selected adversarial labels (the same 100 instances are used for models trained on the same dataset).
Table 1 gives the optimal adversary results. Perturbations 1 were selected such that some big-M problems were solvable within 3600s (problems become more difficult as 1 increases). While our
formulations consistently outperform big-M in terms of problems solved and solution times, the best choice of N is problem-dependent. For instance, the 4-partition formulation performs best for the 2× 100 network trained on CIFAR-10; Figure 2 (right) shows the number of problems solved for this case. The 2-partition formulation is best for easy problems, but is soon overtaken by larger values of N . Intuitively, simpler problems are solved with fewer branch-and-bound nodes and benefit more from smaller subproblems. Performance declines again near N ≥ 7, supporting observations that the convex-hull formulation (N = η) is not always best [3].
Convex-Hull-Based Cuts. The cut-generation strategy for big-M by Anderson et al. [3] is compatible with our formulations (Section 3.3), i.e., we begin with a partition-based formulation and apply cuts during optimization. Figure 2 (left) gives results with these cuts added (via Gurobi callbacks) at various cut frequencies. Specifically, a cut frequency of 1/k corresponds to adding the most violated cut (if any) to each NN node at every k branch-and-bound nodes. Performance is improved at certain frequencies; however, our formulations consistently outperform even the best cases using big-M with cuts. Moreover, cut generation is not always straightforward to integrate with off-the-shelf MILP solvers, and our formulations often perform just as well (if not better) without added cuts. Appendix A.5 gives results with most-violated cuts directly added to the model, as in De Palma et al. [13].
Partitioning Strategies. Figure 3 shows the result of the input partitioning strategies from Section 3.4 on the MNIST 2×100 model for varyingN . Both proposed strategies (blue) outperform formulations with random and uneven partitions (red). With OBBT, big-M (N = 1) outperforms partition-based formulations when partitions are selected poorly. Figure 3 again shows that our formulations perform best for some intermediate N ; random/uneven partitions worsen with increasing N .
Optimization-Based Bounds Tightening. OBBT was implemented by tightening bounds for all zbn. We found that adding the bounds from the 1-partition model (i.e., bounds on ∑ zbn) improved all models, as they account for dependencies among all inputs. Therefore these bounds were used in all models, resulting in 2N + 2 LPs per node (N ≥ 2). We limited the OBBT LPs to 5s; interval bounds were used if an LP was not solved. Figure 3 shows that OBBT greatly improves the optimization performance of all formulations. OBBT problems for each layer are independent and could, in practice, be solved in parallel. Therefore, at a minimum, OBBT requires the sum of max solution times found in each layer (5s × # layers in this case). This represents an avenue to significantly improve MILP optimization of NNs via parallelization. In contrast, parallelizing branch-and-bound is known to be challenging and may have limited benefits [1, 26].

heading: 4.2 Verification Results
text: The verification problem is similar to the optimal adversary problem, but terminates when the sign of the objective function is known (i.e., the lower/upper bounds of the MILP have the same sign). This problem is typically solved for perturbations defined by the `∞-norm (||x− x̄||∞ ≤ ∞ ∈ R). Here, problems are difficult for moderate ∞: at large ∞ a mis-classified example (positive objective) is easily found. Several verification tools rely on an underlying big-M formulation, e.g., MIPVerify [32], NSVerify [2], making big-M an especially relevant point of comparison. Owing to the early termination, larger NNs can be tested compared to the optimal adversary problems. We turned off cuts (cuts= 0) for the partition-based formulations, as the models are relatively tight over the box-domain perturbations and do not seem to benefit from additional cuts. On the other hand, removing cuts improved some problems using big-M and worsened others. Results are presented in Table 2. Our formulations again generally outperform big-M (N=1), except for a few of the 4-partition problems.

heading: 4.3 Minimally Distorted Adversary Results
text: In a similar vein as Croce and Hein [11], we define the `1-minimally distorted adversary problem: given a target image x̄ and its correct label j, find the smallest perturbation over which the NN predicts an adversarial label k. We formulate this as min
1,x 1; ||x− x̄||1 ≤ 1; fk(x) ≥ fj(x). The
adversarial label k is selected as the second-likeliest class of the target image. Figure 4 provides examples illustrating that the `1-norm promotes sparse perturbations, unlike the `∞-norm. Note that the sizes of perturbations 1 are dependent on input dimensionality; if it were distributed evenly, 1 = 5 would correspond to an `∞-norm perturbation ∞ ≈ 0.006 for MNIST models.
Table 3 presents results for the minimally distorted adversary problems. As input domains are unbounded, these problems are considerably more difficult than the above optimal adversary problems. Therefore, only smaller MNIST networks were manageable (with OBBT) for all formulations. Again, partition-based formulations consistently outperform big-M, solving more problems and in less time.

heading: 5 Conclusions
text: This work presented MILP formulations for ReLU NNs that balance having a tight relaxation and manageable size. The approach forms the convex hull over partitions of node inputs; we presented theoretical and computational motivations for obtaining good partitions for ReLU nodes. Furthermore, our approach expands the benefits of OBBT, which, unlike conventional MILP tools, can easily be parallelized. Results on three classes of optimization tasks show that the proposed formulations consistently outperform standard MILP encodings, allowing us to solve 25% more of the problems (average >2X speedup for solved problems). Implementations of the proposed formulations and partitioning strategies are available at https://github.com/cog-imperial/ PartitionedFormulations_NN.

heading: Acknowledgments
text: This work was supported by Engineering & Physical Sciences Research Council (EPSRC) Fellowships to CT and RM (grants EP/T001577/1 and EP/P016871/1), an Imperial College Research Fellowship to CT, a Royal Society Newton International Fellowship (NIF\R1\182194) to JK, a grant by the Swedish Cultural Foundation in Finland to JK, and a PhD studentship funded by BASF to AT.

Review summaries:
This paper proposes a mixed-integer programming formulation for trained ReLU neural networks, which will enable people to solve many related problems like optimal adversary, verification and minimally distorted adversary. The main advantage of this formulation comes from partitioning node inputs into a number of groups and use the disjunctive programming to form the convex hull. By selecting different partition size, one can recover the convex hull or the big-M formulation. For the purpose of balancing model size and tightness of the formulation, numerical experiments suggest that moderate size of partitions performs the best. 

This paper presents a family of mixed-integer optimization formulations for trained neural networks (or, rather, individual neurons in the network), with primary application to verification and adversarial analysis/training. The family traces out a spectrum between a small, weak formulation (the big-M) and a larger, stronger formulation (the "convex hull"). The family is parameterized by a partition of the multivariate input to a single ReLU neuron. A handful of strategies for doing this partition are presented, along with computational experiments.

The paper introduces a class of novel mixed-integer programming formulations for optimizing the output of a trained ReLU network. Its approach is to model a ReLU as a disjunction where the preactivation function is partitioned in its indices. This results a parameterized formulation that can range from the traditional big-M formulation to an ideal formulation at its extremes. Since the former is weak but light and the latter is strong but large, this allows us to find a trade-off where the formulation is both reasonably strong and not too large so that its performance is better than either extreme. Moreover, the paper presents an optimization-based bounds tightening approach that can be applied to these formulation. The authors show that this formulation has good computational performance compared to its extremes for three classes of problems in the area of neural network verification.