DeepSplit: Scalable Verification of Deep Neural Networks via Operator Splitting
Year: 2021

Abstract:
Analyzing the worst-case performance of deep neural networks against input perturbations amounts to solving a large-scale non-convex optimization problem, for which several past works have proposed convex relaxations as a promising alternative. However, even for reasonably-sized neural networks, these relaxations are not tractable, and so must be replaced by even weaker relaxations in practice. In this work, we propose a novel operator splitting method that can directly solve a convex relaxation of the problem to high accuracy, by splitting it into smaller sub-problems that often have analytical solutions. The method is modular, scales to very large problem instances, and compromises of operations that are amenable to fast parallelization with GPU acceleration. We demonstrate our method in obtaining tighter bounds on the worst-case performance of large convolutional networks in image classification and reinforcement learning settings.

heading: 1 Introduction
text: Despite their superior performance, neural networks lack formal guarantees, raising serious concerns about their adoption in safety-critical applications such as autonomous vehicles (Cao et al., 2019) and medical machine learning (Finlayson et al., 2019). Motivated by this drawback, there has been an increasing interest in developing tools to verify desirable properties for neural networks, such as robustness to adversarial attacks. Neural network verification refers to the problem of verifying whether the output of a neural network satisfies certain properties for a bounded set of input perturbations. This problem can be framed as optimization problems of the form
minimize J(f(x)) subject to x ∈ X , (1) where f is given by a deep neural network, J is a real-valued function representing a performance measure (or a specification), and X is a set of inputs to be verified. In this formulation, verifying the neural network amounts to certifying whether the optimal value of (1) is bounded below by a certain threshold. This problem is large-scale and non-convex, making it extremely difficult to solve efficiently–both in terms of time and memory. For ReLU activation functions and linear objectives, the problem in (1) can be cast as a MixedInteger Linear Program (MILP) (Lomuscio & Maganti, 2017; Cheng et al., 2017; Dutta et al., 2018; Fischetti & Jo, 2018), which can be solved for the global solution via, for example, Branch-and-Bound (BaB) methods. While we do not expect these approaches to scale to large problems, for small neural networks they can still be practical.
Instead of solving (1) for its global minimum, one can instead find guaranteed lower bounds on the optimal value via convex relaxations, such as Linear Programming (LP) (Wong & Kolter, 2017) and Semidefinite Programming (SDP) (Raghunathan et al., 2018; Fazlyab et al., 2019; 2020). Verification methods based on convex relaxations are sound but incomplete, i.e., they are guaranteed to detect all false negatives but also produce false positives, whose rate depends on the tightness of the relaxation. Although convex relaxations are polynomial-time solvable (in terms of number of decision variables), in practice they are not computationally tractable for large-scale neural networks. To improve scalability, these relaxations must typically be further relaxed (Wong & Kolter, 2017; Weng et al., 2018; Zhang et al., 2018; Dvijotham et al., 2018; Bunel et al., 2020a).
Contributions In this work, we propose an algorithm for solving convex relaxations exactly for large-scale neural networks. Our starting point is to express (1) as a constrained optimization problem whose constraints are imposed by the forward passes in the network. We then introduce additional decision variables and consensus constraints that naturally split the corresponding problem into independent subproblems, which often have closed-form solutions. Finally, we employ an operator splitting technique based on the Alternating Direction Method of Multipliers (ADMM) (Boyd et al., 2011), to solve the corresponding Lagrangian relaxation of the problem. This approach has several favorable properties. First, the method requires minimal parameter tuning and relies on simple operations, which scale to very large problems and can achieve a good trade-off between runtime and solution accuracy. Second, all the solver operations are amenable to fast parallelization with GPU acceleration. Third, our method is fully modular and applies to standard network architectures.
We employ our method to compute exact solutions to LP relaxations on the worst-case performance of adversarially trained deep networks, with a focus on networks whose convex relaxations were previously impossible to solve exactly due to their size. Specifically, we perform extensive experiments in the `∞ perturbation setting, where we verify robustness properties of image classifiers for CIFAR10 and deep Q-networks (DQNs) in Atari games (Zhang et al., 2020). Our method is able to quickly solve LP relaxations at scales that are too large for exact BaB verifiers, SDP relaxations, or commercial LP solvers such as Gurobi.

heading: 1.1 Related Work
text: Convex relaxations LP relaxations are relatively the most scalable form of convex relaxations (Ehlers, 2017). However, even solving LPs can become computationally prohibitive for small convolutional networks (Salman et al., 2019). One line of work studies computationally cheaper but looser bounds of the LP relaxation (Wong & Kolter, 2017), which have been extended to larger and more general networks and settings (Wong et al., 2018; Weng et al., 2018; Zhang et al., 2018; Xu et al., 2020). These bounds tend to be loose unless optimized during training, which typically comes at a significant cost to standard performance. Further work has aimed to tighten these bounds (Singh et al., 2019; Tjandraatmadja et al., 2020; de2), however these works focus primarily on small convolutional networks and struggle to scale to more typical deep networks. Other work has studied the limits of these convex relaxations on these small networks using vast amounts of CPU-compute (Salman et al., 2019). Recent SDP-based approaches (Dathathri et al., 2020) can produce much tighter bounds on these small networks.
Lagrangian-based bounds Related to our work is that which solves the Lagrangian of the LP relaxation (Dvijotham et al., 2018; Bunel et al., 2020a), which can tighten the bound but do not aim to solve the LP exactly due to relatively slow convergence. However, these works primarily study small networks whose LP relaxation can still be solved exactly with Gurobi. Although these works could in theory be used on larger networks, only the faster, linear-based bounds (Xu et al., 2020) have demonstrated applicability to standard deep learning architectures. In our work, we solve the LP relaxation exactly in large network settings that previously have only been studied with loose bounds of the LP relaxation such as LiRPA (Xu et al., 2020).
Complete verification methods These methods verify properties of deep networks exactly (without false positives or false negatives) using methods such as SMT solvers (Scheibler et al., 2015; Ehlers, 2017; Katz et al., 2017) and MILP solvers (Lomuscio & Maganti, 2017; Cheng et al., 2017; Dutta et al., 2018; Fischetti & Jo, 2018). Complete verification methods typically rely on BaB algorithms (Bunel et al., 2017), in which the verification problem is divided into subproblems (branching) that can be verified using incomplete verification methods (bounding) (Bunel et al., 2020b). However, these methods have a worst-case exponential runtime and have difficulty scaling beyond relatively small convolutional networks. Motivated by this, several recent works have improved the practical running time of BaB methods with custom solvers for medium-sized networks (Bunel et al., 2020a; De Palma et al., 2021; Xu et al., 2021; Wang et al., 2021), but have yet to scale to larger deep networks.
Operator splitting methods Operator splitting, and in particular the ADMM method, is a powerful technique in solving structured convex optimization problems and has applications in numerous settings ranging from optimal control (O’Donoghue et al., 2013) to training neural networks (Taylor et al., 2016). These methods scale well with the problem dimensions, can exploit sparsity in the problem data efficiently (Zheng et al., 2017), are amenable to parallelization on GPU (Schubiger et al., 2020), and have well-understood convergence properties under minimal regularity assumptions (Boyd et al., 2011). The benefit of ADMM as an alternative to interior-point solvers has been shown in various classes of optimization problems (O’donoghue et al., 2016). Our operator splitting method is specifically tailored for neural network verification in order to fully exploit the problem structure.
Notation. We denote the set of real numbers by R, the set of nonnegative real numbers by R+, the set of real n-dimensional vectors by Rn, the set of m× n-dimensional matrices by Rm×n, and the n-dimensional identity matrix by In. The p-norm (p ≥ 1) is denoted by ‖ · ‖p : Rn → R+. For a set S, we define the indicator function IS(x) of S as IS(x) = 0 if x ∈ S and IS(x) = +∞ otherwise. Given a function f : X → Y, the graph of f is the set Gf = {(x, f(x)) | x ∈ X}.

heading: 2 Scalable Neural Network Verification via Operator Splitting
text: We consider an `-layer feed-forward neural network f(x) : Rn0 → Rn` described by the following recursive equations,
x0 = x, xk+1 = φk(xk) k = 0, · · · , `− 1, f(x) = x` (2) where x0 ∈ Rn0 is the input to the neural network, xk ∈ Rnk is the input to the k-th layer, and φk : Rnk → Rnk+1 is the operator of the k-th layer, which can represent any commonly-used operator in deep networks, such as linear (convolutional) layers, MaxPooling units, and activation functions.
Given the neural network f , a specification function J : Rn` 7→ R, and an input set X ⊂ Rn0 , we say that f satisfies the specification J if J(f(x)) ≥ 0 for all x ∈ X . This is equivalent to verifying that the optimal value of (1) is non-negative. We assume X ⊂ Rn0 is a closed convex set and J : Rn` → R∪{+∞} is a convex function, defining a performance measure on the output of the network. Note that our formulation generalizes to arbitrary computational graphs (e.g. residual blocks) and general objective functions (see Appendix B).
Using the sequential representation of the neural network in (2), we may rewrite the optimization problem in (1) as the following constrained optimization problem, J? ← minimize J(x`) subject to xk+1 = φk(xk), k = 0, · · · , `− 1, x0 ∈ X , (3)
with n := ∑` k=0 nk decision variables. We can rewrite (3) equivalently as
J? ← minimize J(x`) subject to (xk, xk+1) ∈ Gφk , k = 0, · · · , `− 1, x0 ∈ X , (4) where Gφk = {(xk, xk+1) | xk+1 = φk(xk), xk ≤ xk ≤ x̄k} is the graph of φk. Here xk and x̄k are a priori known bounds on xk when x0 ∈ X . The problem in (4) is non-convex due to presence of nonlinear operators in the network, such as activation layers. By overapproximating Gφk by a convex set (or ideally by its convex hull), we arrive at a direct layer-wise convex relaxation of the problem. However, solving this relaxation directly cannot scale to even medium-sized neural networks Salman et al. (2019). In this section, by exploiting the sequential structure of the constraints, we propose a reformulation of (3) whose convex relaxation can be solved efficiently and in a scalable manner.

heading: 2.1 Variable Splitting
text: By introducing the intermediate variables yk and zk, we can rewrite (3) as J? ←minimize J(x`) (5)
subject to yk = xk, zk = φk(yk), xk+1 = zk, k = 0, · · · , `− 1, x0 ∈ X
which has now 3n− n0 − n` decision variables. Intuitively, we have introduced additional “identity layers” between consecutive layers (see Figure 1). By overapproximating Gφk by a convex set Sφk , we obtain the convex relaxation
J?relaxed ←minimize J(x`) (6) subject to yk = xk, (yk, zk) ∈ Sφk , xk+1 = zk, k = 0, · · · , `− 1,
x0 ∈ X ,
for which J?relaxed ≤ J?. This form is known as consensus as yk and zk−1 are just copies of the variable xk. As shown below, this “overparameterization” allows us to split the optimization problem into smaller sub-problems that can be solved in parallel and often in closed form.

heading: 2.2 Lagrangian Relaxation and Operator Splitting
text: We use x = (x0, · · · , x`), y = (y0, · · · , y`−1) and z = (z0, · · · , z`−1) to denote the concatenated variables. By relaxing the equality constraints with Lagrangian multipliers, we define the augmented Lagrangian for (6) as follows,
L(x,y, z,λ,µ)=J(x`)+ `−1∑ k=0 ISφk (yk, zk)+IX (x0) + (ρ/2) `−1∑ k=0 ( ‖xk − yk + λk‖22 − ‖λk‖22 ) + (ρ/2)
`−1∑ k=0 ( ‖xk+1 − zk + µk‖22 − ‖µk‖22 ) . (7)
where λ = (λ0, · · · , λ`−1) and µ = (µ0, · · · , µ`−1) are the scaled dual variables (by 1/ρ) and ρ > 0 is the augmentation constant. Note that we have only relaxed the equality constraints in (6), and the constraints describing the sets Sφk as well as X are kept intact. Furthermore, the inclusion of augmentation will render the dual function differentiable (see the Appendix for more details), and hence, easier to optimize. For the Lagrangian in (7), the dual function, which provides a lower bound to J?relaxed, is given by g(λ,µ) = inf(x,y,z) L(x,y, z,λ,µ). The best lower bound can then be found by maximizing the dual function. 1
Solving the inner problem jointly over (x,y, z) to find the dual function is as difficult as solving a direct convex relaxation of (3). Instead, we split the primal variables (x,y, z) into x and (y, z) and apply the classical ADMM algorithm to obtain the following iterations
1If strong duality holds, then this best lower bound would match J?relaxed. As shown in Salman et al. (2019), strong duality holds under mild conditions.
Algorithm 1: DeepSplit Algorithm Data: neural network f (Eq. 2), bounded convex input set X , convex function J . Result: lower bound J?relaxed on Problem (1). Initialization: x0 ∈ X , xk+1 = φk(xk), yk = xk, zk = xk+1, λk = 0, µk = 0, k = 0, · · · , `− 1, augmentation constant ρ > 0. repeat
Step 1: x-update (9) Step 2: (y, z)-update (10) Step 3: dual update (13)
until stopping criterion is met ; Output: J(x`)
(shown in Figure 1) for updating the primal and dual variables,
x+ ∈ argminx L(x,y, z,λ,µ) (8a) (y+, z+) ∈ argmin(y,z) L(x+,y, z,λ,µ) (8b) (λ+,µ+)=(λ,µ)+∇(λ,µ) L(x+,y+, z+,λ,µ). (8c)
As we show below, the Lagrangian has a separable structure by construction that can be exploited in order to efficiently implement each step of (8).
2.3 The x-update
The Lagrangian in (7) is separable across the xk variables; hence, the minimization in (8a) can be done independently for each xk. Specifically, for k = 0, we obtain the following update rule for x0,
x+0 = ProjX (y0 − λ0). (9a) Projections onto the `∞ and `2 balls can be done in closed-form. For the `1 ball, we can use the efficient projection scheme proposed in (Duchi et al., 2008), which has O(n0) complexity in expectation. For subsequent layers k = 1, · · · , `, we obtain the updates
x+k = 1
2 (yk − λk + zk−1 − µk−1) (9b)
x+` =arg min x`
J(x`)+ ρ
2 ‖x`−z`−1+µ`−1‖22. (9c)
For convex J and ρ > 0, the optimization problem for updating x` is strongly convex with a unique optimal solution. Indeed, its solution is the proximal operator of J/ρ evaluated at z`−1−µ`−1. For the special case of linear objectives, J(x`) = c>x`, we obtain the closed-form solution x+` = −c/ρ+ (z`−1 − µ`−1).
2.4 The (y, z)-update
Similarly, the Lagrangian is also separable across the (yk, zk) variables. Updating these variables in (8b) corresponds to the following projection operations per layer,
(y+k , z + k ) = ProjSφk (x+k + λk, x + k+1 + µk), (10)
for k = 0, · · · , ` − 1. Depending on the type of the layer (linear, activation, convolution, etc.), we obtain different projections, which we describe below.
Affine Layers. Suppose φk(yk) = Wkyk+bk is an affine layer representing a fully-connected, convolutional, or an average pooling layer. Then the graph of φk is already a convex set given by Gφk = {(yk, zk) | zk = Wkyk + bk}. Choosing Sφk = Gφk , the projection in (10) takes the form
y+k = (Ink+W > k Wk) −1(x+k + λk+W > k (x + k+1+µk−bk))
z+k = Wky + k + bk.
(11)
The matrix (Ink + W>k Wk) −1 can be pre-computed and cached for subsequent iterations. We can do this efficiently for convolutional layers using the fast Fourier transform, which we discuss later in Section D.3.
Activation Layers. For an activation layer of the form φ(x) := [ϕ1(x1) · · · ϕn(xn)]>, the convex relaxation of Gφ is given by the Cartesian product of individual convex relaxations i.e., Sφ = Sϕ1 × · · · × Sϕn . For a generic activation function ϕ : R→ R, suppose we have a concave upper bound ϕ̄ and a convex lower bound ϕ on ϕ over an interval I = [x, x̄], i.e., ϕ(x) ≤ ϕ(x) ≤ ϕ̄(x) ∀ x ∈ [x, x̄]. A convex overapproximation of Gϕ is
Sϕ={(x, y) | ϕ(x) ≤ y ≤ ϕ̄(x), x ≤ x ≤ x̄}, (12) which turns out to be the convex hull of G(ϕ) when ϕ̄ and ϕ are concave and convex envelopes of ϕ, respectively. As an example, when x < 0 < x̄, the ReLU activation function y = max(0, x) admits the envelopes ϕ(x) = max(0, x), ϕ̄(x) = y + ȳ−yx̄−x (x − x) on [x , x̄ ], where y = max(0, x) and ȳ = max(0, x̄) (Ehlers, 2017; Wong & Kolter, 2017) (see Figure 1). The assumed pre-activation bounds x and x̄ used to relax the activation functions can be obtained a priori via a variety of existing techniques (Wong & Kolter, 2017; Zhang et al., 2018; 2019). We present projection operators of additional layers in Appendix C.

heading: 2.5 The (λ,µ)-update
text: Finally, we update the scaled dual variables as follows, λ+k = λk + (x + k − y + k ), k = 0, · · · , `− 1,
µ+k = µk + (x + k+1 − z + k ), k = 0, · · · , `− 1.
(13)
The DeepSplit Algorithm is summarized in Algorithm 1. We provide further details such as parameter selection and efficient implementations for convolutional layers in Appendix D.

heading: 3 Connection to Lagrangian-based methods
text: In this section, we draw connections to two related methods relying on Lagrangian relaxation. Specifically, we relate our approach to an earlier dual method (Dvijotham et al., 2018) (Section 3.1), as well as a recent Lagrangian decomposition method (Bunel et al., 2020a) (Section 3.2). Overall, these approaches use a similar Lagrangian formulation, but the specific choices in splitting and augmentation of the Lagrangian result in poorer theoretical convergence guarantees when solving the convex relaxation to optimality.

heading: 3.1 Dual method via Lagrangian Relaxation of the nonconvex problem
text: Instead of splitting the neural network equations (3) with auxiliary variables, an alternative strategy is to directly relax (3) with Lagrangian multipliers (Dvijotham et al., 2018):
L(x,λ) = J(x`) + `−1∑ k=0 λ>k (xk+1 − φk(xk)) + IX (x0). (14)
where x = (x0, · · · , x`) and λ = (λ0, · · · , λ`−1). This results in the dual problem g? ← maximize g(λ), where the dual function is
g(λ)= inf x`≤x`≤x̄` {J(x`) + λ>`−1x`}+ `−1∑ k=1 inf xk≤xk≤x̄k {λ>k−1xk−λ>k φk(xk)}+ inf x0∈X0 {−λ>0 φ0(x0)}.
By weak duality, g? ≤ J?. The inner minimization problems in g(λ) to compute the dual function for a given λ can often be solved efficiently or even in closed-form (Dvijotham et al., 2018). The resulting dual problem is unconstrained but non-differentiable; hence it is solved using dual subgradient method (Dvijotham et al., 2018). However, subgradient methods are known to be very slow with convergence rate O(1/ √ N) where N is the number of updates (Nesterov, 2003), making it inefficient to find exact solutions to the convex relaxation. On the other hand, this method can be stopped at any time to obtain a valid lower bound.

heading: 3.2 Lagrangian method via a non-separable splitting
text: Another related approach is the Lagrangian decomposition method from Bunel et al. (2020a). To decouple the constraints for the convex relaxation of (3), this approach can be viewed as introducing one set of intermediate variables yk as copies of xk to obtain
J?relaxed ←minimize J(y`) (15) subject to (yk, xk+1) ∈ Sφk k = 0, · · · , `− 1
yk = xk k = 0, · · · , ` x0 ∈ X
This splitting is in the spirit of the splitting introduced in Bunel et al. (2020a;b),2 and differs from our splitting which uses two sets of variables in (5). By relaxing the consensus constraints yk = xk, the Lagrangian is
L(x,y,µ) = J(y`) + ∑̀ k=0 µ>k (yk − xk) + `−1∑ k=0 ISφk (yk, xk+1) + IX (x0). (16)
Again the Lagrangian is separable and its minimization results in the following dual function
g(µ) = inf x`≤y`≤x̄` {J(y`)+µ>` y`}+ `−1∑ k=0 inf (yk,xk+1)∈Sφk {µ>k yk − µ>k+1xk+1}+ inf x0∈X {−µ>0 x0}.
Since the dual function is not differentiable, it must be maximized by a subgradient method, which again has an O(1/ √ N) rate. To induce differentiability in the dual function and improve speed, Bunel et al. (2020a) uses Augmented Lagrangian (AL). Since only one set of variables was introduced in (15), the AL is no longer separable across the primal variables. Therefore, for each update of the dual variable, the AL must be minimized iteratively. To this end, Bunel et al. (2020a) use the Frank-Wolfe Algorithm in a block-coordinate as an iterative subroutine. However, this slows down overall convergence and suffers from compounding errors when the sub-problems are not fully solved. When stopping early, the primal minimization must be solved to convergence in order to compute the dual function and produce a valid bound.
In contrast to the approach described above, in this paper we used a different variable splitting scheme in (5) that allows us to fully separate layers in a neural network. This subtle difference has a significant impact: we can efficiently minimize the corresponding AL in closed form, without resorting to any iterative subroutine. Specifically, we use the ADMM algorithm, which is known to converge at an O(1/N) rate (He & Yuan, 2012). In summary, our method enjoys an order of magnitude faster theoretical convergence, is more robust to numerical errors, and has minimal requirements for parameter tuning. To stop early, we can use a similar strategy as the Frank-Wolfe approach from Bunel et al. (2020a) and run the primal iteration to convergence with fixed dual variables to get a valid bound.

heading: 4 Experiments
text: The strengths of our method are (a) its ability to exactly solve LP relaxations and (b) do so at scales that were previously impossible. To evaluate this, we first demonstrate how solving the LP to optimality leads to tighter certified robustness guarantees in image classification and reinforcement learning tasks (Section 4.1). We then stress test our method in both speed and scalability against a commercial LP solver and in the large network setting (Section 4.2).
Setup In all the experiments, we focus on the setting of verification-agnostic networks, similar to Dathathri et al. (2020). However, we focus on significantly larger networks that
2If we define φk(xk) = Wk+1σ(xk) + bk+1, where Wk+1, bk+1 are the parameters of the affine layer and σ is a layer of activation functions, this splitting coincides with the one proposed in Bunel et al. (2020a;b)
have previously only been feasibly bounded with fast linear-based bounds of the LP relaxation (Wong et al., 2018; Xu et al., 2020). All the networks have been trained adversarially with the PGD attack (Madry et al., 2017). Full details about the network architecture, network training, and the parameters for our splitting method can be found in Appendix F.

heading: 4.1 Improved Bounds From Exact LP Solutions with ADMM
text: We first demonstrate how solving the LP exactly with our method results in tighter bounds than prior work. We consider two main settings: certifying the robustness of classifiers for CIFAR10 and deep Q-networks (DQNs) in Atari games. We defer additional analogous results in the smaller MNIST setting to Appendix F.3.
CIFAR10 Much work studying verification in CIFAR10 has focused primarily on a small CNN with only 6k hidden units—smaller than the classic LeNet from MNIST (LeCun et al., 2015). In this regime, tighter approaches such as SDP and branch-and-bound are applicable.
Instead, we consider an order of magnitude larger model, which cannot be feasibly solved by these alternatives. Up until this point, the only solutions for large networks were linear-based bounds to the LP relaxation (Wong et al., 2018; Xu et al., 2020). However, these bounds are known to be quite loose. How much better can we do if we solve the LP exactly?
In Table 1, we find that solving the LP exactly leads to significant gains in certified robustness for large networks, with up to 4% additional certified robustness over the best-performing alternative. All the methods in Table 1 are given the same time budget. Indeed, the better theoretical convergence guarantees of ADMM translate to better results in practice: when given a similar budget, the Lagrangian baselines have worse convergence and cannot verify as many examples. A complete description of these baselines is in Appendix F.3, with implementation details for the Lagrangian methods in Appendix G.
State-robust RL We demonstrate our approach on a non-classification benchmark from reinforcement learning: verifying the robustness of deep Q-networks (DQNs) to adversarial state perturbations (Zhang et al., 2020). Specifically, we verify whether a learned DQN policy outputs stable actions in the discrete space when given perturbed states. Similar to the large network considered in the CIFAR10 setting, this benchmark has only been demonstrably verified with fast but loose linear-based bounds (Xu et al., 2020).
Similar to the CIFAR10 setting, we observe consistent improvement in certified robustness of the DQN when solving the LP exactly with ADMM across multiple RL settings. We summarize the results using our method and the linear-based bounds on LP relaxations (Wong et al., 2018; Xu et al., 2020) in Table 2. Further details regarding the dataset, the DQN architectures, and this task can be found in Appendix F.4.

heading: 4.2 Speed and Scalability
text: In this section, we stress test several aspects of DeepSplit for solving LP relaxations. We first compare solving speeds of our ADMM-based LP solver to a commercial-grade LP solver. We then push the limits of architecture size and solve the LP relaxation for a standard ResNet18.
Speed comparison We first compare the solving speeds of our method with state-of-theart solvers for convex relaxations: a commercial-grade LP solver, Gurobi. Since Gurobi cannot handle large networks, we benchmark the approaches on a fully connected network that Gurobi can handle (see Appendix H for further details on experimental setup, as well as additional speed comparisons to alternative approaches).
We find that our method provides a nearly 7x speedup over Gurobi. On average, it takes ADMM 38 seconds per example to calculate these intermediate bounds using a single GeForce GTX 1080, in comparison to 258 seconds per example for Gurobi on an Intel Core i7-6700K.
Scalability Finally, to test the scalability and generality of our approach, we consider solving the LP relaxation for a ResNet18, which up to this point has simply not been possible. The only applicable method here is LiRPA (Xu et al., 2020)—a highly scalable implementation of the linear-based bounds that works for arbitrary networks but can be quite loose in practice. We defer specific experimental details to Appendix I.
We find that exact LP solutions with our ADMM solver can also produce substantial improvements in the bound at ResNet18 scales, as shown in Figure 2. For a substantial number of examples, we find that ADMM can find significantly tighter bounds (especially for lower bounds). A tabular presentation of the results is in Table 5 of Appendix I.

heading: 5 Conclusion
text: In this paper, we proposed DeepSplit, a scalable and modular operator splitting technique for solving convex relaxation-based verification problems for neural networks. The method can exactly solve large-scale LP relaxations with GPU acceleration with favorable convergence rates. Our approach leads to significantly tighter bounds across a range of classification and reinforcement learning benchmarks, and can scale to a standard ResNet18. We leave as future work a further investigation of variations of ADMM that can improve convergence rates in deep learning-sized problem instances, as well as extensions beyond the LP setting.

heading: A Convergence Analysis of DeepSplit
text: By defining x1 = x, x2 = (y, z) (the primal variables) and ν = (λ,µ) (the scaled dual variables), we can write the split convex relaxation in (6) as
minimize f1(x1) + f2(x2) (17) subject to A1x1 +A2x2 = 0
with the corresponding Augmented Lagrangian
L(x1,x2,ν) = f1(x) + f2(x2) + ρ
2 (‖A1x1 +A2x2 + ν‖22 − ‖ν‖22) (18)
where f1 : Rn → R∪{+∞} and f2 : R2n−2n` → R∪{+∞} are extended real-valued functions defined as f1(x1) := J(x`)+IX (x0), f2(x2) := ∑`−1 k=0 ISφk (yk, zk). Moreover, A1x1+A2x2 = 0 represents the set of equality constraints yk = xk and xk+1 = zk for k = 0, · · · , `. The dual function is given by
g(ν) = inf x1,x2
L(x1,x2,ν). (19)
By Danskin’s theorem (Bertsekas, 1997), the sub-differential of the dual function is given by ∂g(ν) = {A1x̄1 +A2x̄2 : (x̄1, x̄2) ∈ arg minL(x1,x2,ν)} (20)
Here x̄1, x̄2 are minimizers of the Lagrangian (not necessarily unique), which satisfy the optimality conditions
0 ∈ ∂f1(x̄1) + ρA>1 (A1x̄1 +A2x̄2 + ν) (21) 0 ∈ ∂f2(x̄2) + ρA>2 (A1x̄1 +A2x̄2 + ν)
First, we show that the sub-differential is a singleton, i.e., g is continuously differentiable. Suppose (x̄1, x̄2) ∈ arg minL(ξ,ν) and (w̄1, w̄2) ∈ arg minL(ξ,ν) are two distinct minimizers of the Lagrangian, hence satisfying
0 ∈ ∂f1(w̄1) + ρA>1 (A1w̄1 +A2w̄2 + ν) (22) 0 ∈ ∂f2(w̄2) + ρA>2 (A1w̄1 +A2w̄2 + ν)
By monotonicity of the sub-differentials, we can write (∂f1(x̄1)− ∂f1(w̄1))>(x̄1 − w̄1) ≥ 0 (23) (∂f2(x̄2)− ∂f2(w̄2))>(x̄2 − w̄2) ≥ 0
By substituting (21) and (22) in (23), we obtain
− (ρA>1 (A1x̄1 +A2x̄2 + ν)− ρA>1 (A1w̄1 +A2w̄2 + ν))>(x̄1 − w̄1) ≥ 0 (24) − (ρA>2 (A1x̄1 +A2x̄2 + ν)− ρA>2 (A1w̄1 +A2w̄2 + ν))>(x̄2 − w̄2) ≥ 0
By adding the preceding inequalities, we obtain −ρ‖A1x̄1 +A2x̄2 − (A1w̄1 +A2w̄2)‖22 ≥ 0 (25)
This implies that A1x̄1 + A2x̄2 = A1w̄1 + A2w̄2 and hence, the sub-differential ∂g is a singleton.
Convergence. When X is a closed nonempty convex set, IX (x0) is a convex closed proper (CCP) function. Assuming that J is also CCP, then we can conclude that f1 is CCP. Furthermore, since the sets Sφk are nonempty convex sets, we can conclude that f2 is CCP. Under these assumptions, the augmented Lagrangian has a minimizer (not necessarily unique) for each value of the dual variables. Finally, under the assumption that the Augmented Lagrangian has a saddle point (which produces a solution to (17)), the ADMM algorithm we have primal convergence ‖rp‖ → 0 (see (31)), dual residual convergence ‖rd‖ → 0, as well as objective convergence J(x`)→ J? (Boyd et al., 2011). We remark that the convergence guarantees of ADMM holds even if f1 and f2 assume the value +∞. This is the case for indicator functions resulting in projections in the first two updates of ADMM.

heading: B ADMM for General Computational Graphs
text: To handle general computational graphs, we follow the approach used in (Wong et al., 2018). Specifically, consider a generalized `-"layer" neural network given by the equations
xk+1 = k∑ i=1 φik(xi) ≡ Φk(x1:k) k = 0, · · · , `− 1 (26)
where x1:k denotes the concatenated vector of variables (x1, · · · , xk). Note that this view of a network reduces to the feedforward equations in (2) when φik(xi) = 0 for all i < k. Furthermore, this view subsumes arbitrary skip connections, such as residual connections. Then, the graph form of the verification problem in (4) can be rewritten for the generalized neural network as
J? ←minimize J(x`) (27) subject to (x1:k, xk+1) ∈ Gφk(x1:k) k = 0, · · · , `− 1
x0 ∈ X . Similar to before, let Sφk be a sound convex approximation of Gφk . Then, our convex relaxation introduces additional intermediate variables y1:k and zk to get the following convex relaxation (analogous to (6)):
p?relaxed ←minimize J(x`) (28) subject to y1:k = x1:k k = 0, · · · , `− 1
(y1:k, zk) ∈ Sφk k = 0, · · · , `− 1 xk+1 = zk k = 0, · · · , `− 1 x0 ∈ X .
Thus, producing an ADMM layer for a network with arbitrary skip connections reduces to deriving the corresponding projection operator onto Sφk . Since most residual connections simply add an affine transformation of a previous layer, these projections are straightforward and have a closed form (an example of a typical residual ADMM block is in Appendix C).

heading: C ADMM layers

heading: C.1 ReLU projection
text: Consider the ReLU function y = ϕ(x) = max(0, x) on the interval x ∈ [x, x̄]. The projection of a point (x(0), y(0)) onto the convex hull of G(ϕ) has a closed-form solution. We consider three cases.
1. Active neuron (x ≥ 0). In this case, the convex hull is given by Sϕ = {(x, y) | y = x, x ≤ x ≤ x̄}.
and the projection of (x(0), y(0)) onto Sϕ is given by x′ = y′ = min(max(x (0)+y(0)
2 , x), x̄).
2. Inactive neuron (x̄ ≤ 0). In this case, the convex hull is Sϕ = {(x, y) | y = 0, x ≤ x ≤ x̄}.
and the projection of (x(0), y(0)) onto Sϕ is given by x′ = min(max(x(0), x), x̄), y′ = 0. 3. Unknown neuron (x ≤ 0 ≤ x̄). The convex hull is given by
Sϕ = {(x, y) | y ≥ 0, y ≥ x, y ≤ y + s(x− x)}
where s = ȳ−yx̄−x , y = max(0, x), ȳ = max(0, x̄). We first project (x (0), y(0)) onto each facet of the triangle and then select the point with the minimal distance.
x(1) = min(max( x(0) + y(0)
2 , 0), x̄) y(1) = x(1)
x(2) = min(max( x(0) + sy(0) + s(sx− y)
s2 + 1 , x), x̄) y(2) = s(x(0) − x) + s2y(0) + y s2 + 1
x(3) = min(max(0, x(0)), x) y(3) = 0
The projected point is (x′, y′) = (x(i ?), y(i ?)), where i? = arg min1≤i≤3 √ (x(0) − x(i))2 + (y(0) − y(i))2.

heading: C.2 Convolutions
text: Although a convolution is technically a linear operator, it is impractical to directly form the inverse matrix from the projection step of the DeepSplit algorithm. Instead, we represent a typical convolutional layer fconv with stride, padding and bias as
fconv = fbias ◦ fds ◦ fcrop ◦ fcirc ◦ fpad = fpost ◦ fcirc ◦ fpad, (29)
where fpad is a padding step, fcirc is a circular convolution step, fcrop is a cropping step, fds is a downsampling step to handle stride greater than one, and fbias is a step that adds the bias. Then, the last three steps are combined into one post-processing step fpost to reduce the number of concensus constraints in the DeepSplit algorithm. The projection steps for all of these operators are presented next in Sections C.3 through C.6, and an efficient implementation of the projection step for the affine layer fcirc is in Appendix E.

heading: C.3 Padding
text: The padding layer fpad takes an image as input and adds padding to it. Denote the input image by yk and the padded image by zk. We can decompose the output zk into two vectors, z0k which is a copy of the input yk, and z 1 k which represents the padded zeros on the edges of image. Equivalently, the padding layer zk = φk(yk) can be written in an affine form
zk = [ z0k z1k ] = [ I 0 ] yk = Wkyk,
for which the projection operator reduces to the affine case.

heading: C.4 Cropping
text: The cropping layer fcrop crops the output of the circular convolution fcirc to the original size of the input image before padding. Denote yk the input image and zk the output image of the cropping layer. By decomposing the input image yk into the uncropped pixels y0k and the cropped pixels y1k, the cropping layer zk = φk(yk) has an affine formulation
zk = y 0 k = [I 0] [ y0k y1k ] = Wkyk
whose projection operator is given in Section 2.4.

heading: C.5 Down-sampling and bias
text: If the typical convolutional layer fconv has stride greater than one, a down-sampling layer is added in the DeepSplit algorithm, which essentially has the same affine form as the
cropping layer with different values of y0k and y 1 k. Therefore, the projection operator for the down-sampling layer reduces to the affine case as well.
The bias layer in the DeepSplit algorithm handles the case when the convolutional layer fconv has a bias bk and is implemented by zk = φk(yk) = yk + bk. This is an affine expression and its projection operator is given in Section 2.4.

heading: C.6 Convolutional post-processing layer
text: We combine the cropping, down-sampling and bias layers into one post-processing layer, i.e., fpost = fbias ◦ fds ◦ fcrop, as shown in (29). This reduces the total number of concensus constraints in the DeepSplit algorithm. Since all the three layers are in fact affine, the post-processing layer is also affine and its projection operator can be obtained correspondingly.

heading: C.7 Maxpooling
text: Consider a Maxpooling layer φ(x) = maxk xk, k = 1, · · · , n over x ≤ x ≤ x̄. (Bunel et al., 2020b, Appendix B.2) shows that the Maxpooling layer can be equivalently decomposed as a combination of linear and ReLU layers. This is motivated by the following fact that max(x1, x2, x3, x4) = max(max(x1, x2),max(x3, x4)) and max(x1, x2) = max(x1 − x2, 0) + max(x2 − x2, 0) + x2. Therefore, we can represent the Maxpooling layer by a composition of pairwise maximum functions which can themselves be decomposed into a combination of linear and ReLU layers. The projection operations onto the convex hull of the ReLU and linear layers are closed-form as shown in Appendix C.1 and Section 2.4, respectively.

heading: C.8 Residual connection
text: We consider a typical residual connection φk(yk, yi) = yk +Wiyi + bi for i < k, where yi is an arbitrary layer before yk and Wi, bi is any affine transformation (i.e. a 1 by 1 convolution for upsampling or downsampling layer frequently used in residual connections with differently sized feature maps). Then, we can write this as
φk(yki) = Wkiyki + bi, (30)
where yki = [ yk yi ] and Wki = [I Wi]. Then, the projection operator reduces to the affine case but with weights Wki, bi for input yki.
C.9 Batch normalization and average pooling
The batch normalization and average pooling layers that appear in the architecture of the residual networks, e.g., ResNet18 applied in this paper, are essentially affine mapping layers whose weights and biases can be extracted accordingly. Therefore, the projection onto the graph of these two layers has a closed-form solution as shown in Eq. (11).
D Implementation Details

heading: D.1 Convergence and Stopping Criterion
text: The DeepSplit algorithm converges to the optimal solution of the convex problem (6) under mild conditions. Specifically, when J is closed, proper and convex, and when the sets Sk (convex outer approximations of the graph of the layers) along with X are closed and convex, we can resort to the convergence results of ADMM (Boyd et al., 2011).
For the LP relaxation (6) of a feed-forward neural network, the primal and dual residuals are defined as
rp = `−1∑ k=0 {‖y+k − x + k ‖ 2 + ‖x+k+1 − z + k ‖ 2} (31)
rd = ρ `−1∑ k=1 ‖(y+k − yk) + (z + k−1 − zk−1)‖ 2 2 + ρ ( ‖y+0 − y0‖22 + ‖z + `−1 − z`−1‖ 2 2 ) .
These are the residuals of the optimality conditions for (6) and converge to zero as the algorithm proceeds. A reasonable termination criterion is that the primal and dual residuals must be small, i.e. rp ≤ p and rd ≤ d, where p > 0 and d > 0 are tolerance levels (Boyd et al., 2011, Chapter 3). These tolerances can be chosen using an absolute and relative criterion, such as
p= √ p abs+ rel max{(‖x0‖22+2 `−1∑ i=1 ‖xi‖22+‖x`‖22) 1 2 + ( `−1∑ i=0 (‖yi‖22 + ‖zi‖22)1/2}
d= √ n abs+ rel(‖λ0‖22+ `−1∑ i=1 ‖λi+µi−1‖22 + ‖µ`−1‖22) 1 2 ,
where p = n0 + 2 ∑`−1 i=1 ni + n`, n = ∑` i=0 ni, abs > 0 and rel > 0 are absolute and relative tolerances. Here n is the dimension of x, the vector of primal variables that are updated in the first step of the algorithm, and p is the total number of consensus constraints.
For neural networks with general computational graphs such as residual networks, the primal and dual residuals as well as the stopping criteria have different representations from what’s shown above. But with our proposed splitting method, these representations are easy to derive from (Boyd et al., 2011, Chapter 3).

heading: D.2 Parameter Selection
text: A proper selection of the augmentation constant ρ has a dramatic effect on the convergence of the method. Large values of ρ enforces consensus more quickly, yielding smaller primal residuals but larger dual ones. Conversely, smaller values of ρ leads to larger primal and smaller dual residuals. A commonly used heuristic to make this trade-off is residual balancing He et al. (2000), in which the penalty parameter varies adaptively based on the following rule:
ρ+ =  τρ if rp > µrd ρ/τ if rd > µrp ρ otherwise,
where µ, τ > 1 are parameters. In our experiments, we found this rule to be effective in speeding up the practical convergence.

heading: D.3 Convolutional Layers
text: The projection step for affine layers from (11) requires multiplication by the matrix (Ink + W>k Wk)
−1 for that layer, where nk is the number of neurons in the layer. When handling networks with convolutional layers on image data, nk can easily exceed tens of thousands, so the resulting matrix and its inversion can exceed reasonable memory and computational constraints.
To make the update step practical and in line with typical computational costs of deep convolutional networks, we replace the typical deep learning convolution with an equivalent circular convolution. Specifically, let fconv be a typical strided convolution with padding. We can rewrite fconv as three sequential updates using a circular convolution as fconv = fpost ◦ fcirc ◦ fpad, where fpad is a padding step, fcirc is a circular convolution, and fpost is a
post-processing step that performs cropping, downsampling, and adds any bias term from fconv.
We can now treat these three updates as separate, individual layers in our ADMM algorithm. The key observation is that we can use the convolution theorem to implement the ADMM update for the circular convolution fcirc efficiently. Specifically, for an input of size nk × nk, the projection step from (11) for a circular convolution can be calculated in O(n2k log nk) using the fast Fourier transform. We provide the specific details of this procedure in Appendix E.

heading: E FFT implementation for circular convolutions
text: In order to efficiently implement convolutional layers, recall that we decompose a convolution into the following three steps:
fconv = fpost ◦ fcirc ◦ fpad (32)
We now discuss in detail how to efficiently perform the (y, z) update for multi-channel, circular convolutions fcirc using Fourier transforms. We begin with the single-channel setting, and then extend our procedure to the multi-channel setting. See Appendix C for details regarding the ADMM projection step for fpad and fpost.

heading: E.1 Single-channel circular convolutions
text: Let U represent the discrete Fourier transform (DFT) as a linear operator, and let W be the weight matrix for the circular convolution fcirc(x) = W ∗ x. Then, using matrix notation, the convolution theorem states that
fcirc(x) = W ∗ x = U∗(UW · Ux) = U∗DUx (33)
where D = diag(UW ) is a diagonal matrix containing the Fourier transform of W and U∗ is the conjugate transpose of U . Then, we can represent the inverse operator from (11) as
(I + f>circfcirc) −1 = U∗(I +DD)−1U (34)
Since (I +DD) is a diagonal matrix, its inverse can be computed quickly by simply inverting the diagonal elements, and requires storage space no larger than the original kernel matrix. Thus, multiplication by the inverse matrix for a circular convolution reduces to two DFTs and an element-wise product. For an input of size n× n, this step has an overall complexity of O(n2 log n) when using fast Fourier transforms.

heading: E.2 Multi-channel circular convolutions
text: We now extend the operation for single-channel circular convolutions to multi-channel, which is typically used in convolutional layers found in deep vision classifiers. Specifically, for a circular convolution with n input channels and m output channels, we have
fcirc(x)j = n∑ i=1 Wij ∗ xi (35)
where fcirc(x)j is the jth output channel output of the circular convolution, Wij is the kernel of the ith input channel for the jth output channel, and xi is the ith channel of the input x. The convolutional theorem again tells us that
fcirc(x)j = n∑ i=1 U∗DijUxi (36)
where Dij = diag(UWij). This can be re-written more compactly using matrices as
fcirc(x) = Ū ∗D̄Ū x̄ (37)
where
• Ū =  U · · · 0... . . . ... 0 · · · U  is a block diagonal matrix with n copies of U along the diagonal
• Ū∗ =  U ∗ · · · 0 ... . . . ...
0 · · · U∗  is a block diagonal matrix with m copies of U along the diagonal
• D̄ =  D11 · · · Dn1... . . . ... D1m · · · Dnm  is a block matrix with diagonal blocks where the ijth block is Dij
• x̄ =  x1... xn  is a vertical stacking of all the input channels. Then, we can represent the inverse operator from (11) as
(I + f>circfcirc) −1 = Ū∗(I + D̄D̄)−1Ū (38)
where I + D̄D̄ is a block matrix, where each block is a diagonal matrix. The inverse can then be calculated by the inverting sub-matrices formed from sub-indexing the diagonal components. Specifically, let D̄j::p be a slice of D̄ containing elements spaced m elements apart in both column and row directions, starting with the jth item. For example, D̄0::p is the matrix obtained by taking the top-left most element along the diagonal of every block. Then, for j = 1 . . .m, we have
(I + D̄D̄)−1j::p = ( (I + D̄D̄)j::p )−1 (39) Thus, calculating this matrix amounts to inverting a batch of p matrices of size m×m. For typical convolutional networks, m is typically well below 1, 000, and so this can be calculated quickly. Further note that this only needs to be calculated once as a pre-computation step, and can be reused across different inputs and objectives for the network.
Memory and runtime requirements In practice, we do not store the fully-expanded block diagonal matrices; instead, we omit the zero entries and directly store the the diagonal entries themselves. Consequently, for an input of size p, the diagonal matrices require storage of size O(mnp), and the inverse matrix requires storage of size O(m2p). Since the discrete Fourier transform can be done in O(p log p) time with fast Fourier transforms, the overall runtime of the precomputation step to form the matrix inverse is the cost of the initial DFT and the batch matrix inverse, or O(nmp log p+m3p). Finally, the runtime of the projection step is O((n+m)p log p+n2mp), which is the respective costs of the DFT transformations Ū and Ū∗, as well as the multiplication by D̄. Since the number of channels in a deep network are typically much smaller than the size of the input to a convolution (i.e. n < p and m < p), the costs of doing the cyclic convolution with Fourier transforms are in line with typical deep learning architectures.

heading: F Experimental details
text: F.1 Architectures, datasets, and training specifics
MNIST We consider two ReLU based architectures. The first one is a fully connected network with layer sizes 784−600−400−200−100−10 which we denote MNIST-A. It is more than triple the size of that considered by (Dathathri et al., 2020) with one additional layer. It is, however, still small enough such that Gurobi is able to solve the LP relaxation, and allows us to compare our running time against Gurobi. The second one is a convolutional network
which uses the small convolutional architecture from (Wong et al., 2018) and consists of two convolutional layers of size 16− 32 with kernel sizes 4− 4, strides 2− 2, and padding 1− 1. We denote this network by MNIST-B. This architecture is comparable to the convolutional architecture considered by (Dathathri et al., 2020), and allows us to compare our running time against the gradient-based SDP solver.
We train both models with an `∞ PGD adversary at radius = 0.1, using 7 steps of size α = 0.02, with batch size 100 for 100 epochs. We use the Adam optimizer with a cyclic learning rate (maximum learning rate of 0.005), and both models achieve a clean accuracy of 99%.
CIFAR10 For CIFAR10, we use the large convolutional architectures from (Wong et al., 2018), which consists of four convolutional layers with 32 − 32 − 64 − 64 channels, with strides 1− 2− 1− 2, kernel sizes 3− 4− 3− 4, and padding 1− 1− 1− 1. This is followed by three linear layers of size 512− 512− 10. This is significantly larger than the CIFAR10 architectures considered by (Dathathri et al., 2020), and has sufficient capacity to reach 43% adversarial accuracy against an `∞ PGD adversary at = 8/255.
The model is trained against a PGD adversary with 7 steps of size α = 2/255 at a maximum radius of = 8.8/255, with batch size 128 for 200 epochs. We used the SGD optimizer with a cyclic learning rate (maximum learning rate of 0.23), momentum 0.9, and weight decay 0.0005. The model achieves a clean test accuracy of 71.8%.
State-robust RL We use the pretrained, adversarially trained, DQNs released by (Zhang et al., 2020) from
https://github.com/chenhongge/SA_DQN
which has not included any license information up to the date of submission of this paper. These models were trained to be robust at = 1/255 with a PGD adversary for the Atari games Pong, Roadrunner, Freeway, and BankHeist. Each input to the DQN is of size 1× 84× 84, which is more than double the size of CIFAR10. The DQN architectures are convolutional networks, with three convolutional layers with 32 − 64 − 64 channels, with kernel sizes 8− 4− 3, strides 4− 2− 1, and no padding. This is followed by two linear layers of size 512−K, where K is the number of discrete actions available in each game.

heading: F.2 Timing experiments
text: Runtime summary: We certify the robustness of classifiers for MNIST and CIFAR10 and DQNs in Atari games to `∞ perturbations over a range of radii. The specific verification setups can be found in Appendix F.3 and F.4. We report the approximate runtime for all experiments in Table 3.
Effects of ADMM parameters: In this subsection, we demonstrate the effects of the algorithm parameters on the convergence of ADMM through two networks: MNIST-A (fully-connected) and MNIST-B (convolutional). Specifically, we focus on the choice of ρ and the application of residual balancing.
We conduct our experiment on the 1938-th example which is randomly chosen from the MNIST test data set. For this example, both MNIST-A and MNIST-B predicts its class (number 4) correctly. We add an `∞ perturbation of radius = 0.02 to the input image and verify if the network outputs are robust with respect to class number 3. The maximum number of iterations is restricted to 3000. The objective values, primal and dual residuals of ADMM for MNIST-A and MNIST-B under different augmentation parameters ρ are plotted in Figure 4 and 5, respectively. The residual balancing in this experiment is applied with initial ρ = 10.0, τ = 2, and µ = 10 as described in Section D.2.
The update rule of ADMM suggests that a large value of ρ tends to produce small primal residuals since it puts a large penalty on the violation of the primal feasibility. However, the dual residuals for such ρ tends to be larger. Conversely, a small ρ tends to reduce the dual residuals at the cost of larger primal residuals. This phenomenon is illustrated empirically with the fully-connected network MNIST-A in Figure 4 and the convolutional network MNIST-B in Figure 5.
Since ADMM terminates when both the primal and dual residuals are small enough, in practice we prefer to choose the augmentation parameter ρ not too large or too small in order to balance the reduction in the primal and dual residuals. An effective way to choose a good ρ is residual balancing, which tries to keep the primal and dual residuals close during the ADMM updates by adjusting ρ online. In both Figure 4 and Figure 5, ADMM with residual balancing is initialized with ρ = 10.0 and shows significant improvement in convergence rate compared with the case of constant ρ = 10.0. As observed in our other experiments, with residual balancing, ADMM becomes insensitive to the initialization of ρ and usually achieves a good convergence rate.
In all of the test accuracy certification results reported in this paper, we initialize ρ = 1.0 and apply residual balancing when running ADMM. In different experiments, the stopping criterion parameters abs, rel of ADMM are chosen by trial-and-error to achieve a balance between the accuracy and the running time of the algorithm. Although ADMM has convergence guarantees for its objective values, primal and dual residuals, it may take too many iterations for ADMM to achieve a solution of high accuracy (Boyd et al., 2011).
F.3 Image classification
In this section, we report the details of certified test accuracy of our proposed method (ADMM) and the scalable linear-based bounds on LP relaxations (Linear) (Wong et al., 2018; Xu et al., 2020) in the image classification tasks of MNIST and CIFAR10.
MNIST: For both the fully-connected network MNIST-A and the convolutional network MNIST-B, we apply LP-ADMM with abs = 10−3, rel = 10−3 as the stopping criterion, and residual balancing with the initial ρ = 1.0. For each classifier network, we go through the 10000 examples in the MNIST test data set and a range of `∞ perturbation radii to count the number of certified examples. To make the counting more efficient, we (i) search over in a descending order since examples that are robust for a larger are also robust for a smaller , and (ii) only apply ADMM on examples that cannot be verified by Linear (Wong et al., 2018; Xu et al., 2020) since Linear (Wong et al., 2018; Xu et al., 2020) gives a more relaxed bound than the LP-relaxation. For each batch of test examples, ADMM solves ten optimization problems of the form (5) with the linear objective function given with respect to each prediction class.
For a range of ∈ {0.01, 0.02, · · · , 0.10}, the numbers of verified examples of ADMM and Linear (Wong et al., 2018; Xu et al., 2020) and their differences are shown in Figure 6 for MNIST-A and MNIST-B, respectively. Exact certified accuracy is reported in Table 4 for a range of .
CIFAR10: On the CIFAR10 data set, we compare ADMM and Linear (Wong et al., 2018; Xu et al., 2020) following the same process as on the MNIST data set. The stopping criterion for LP-ADMM is set as abs = 10−4, rel = 10−4 and the range of is set as {0.5/255, 1.0/255, · · · , 4.0/255}. The numbers of certified examples of ADMM and Linear (Wong et al., 2018; Xu et al., 2020) and their differences are shown in Figure 7.
Further details for baselines in Table 1 In Table 1, we report the certified accuracy of solving the LP exactly with ADMM in comparison to a range of baselines. We compare with methods that have previously demonstrated the ability to bound networks of this size: fast

heading: F.4 Reinforcement learning
text: We compare the tightness of ADMM and Linear (Wong et al., 2018; Xu et al., 2020) by verifying the robustness of DQNs on three Atari game benchmarks: BankHeist, Roadrunner, and Pong 3. The DQNs applied in these experiments are introduced in Appendix F.1.
For each benchmark, we collect 10, 000 frames (each with dimension 1× 84× 84) across 100 episodes using the natural policy with 100 frames sampled randomly from each episode as the test data set. Note that the input images to the DQNs are already pre-processed such that the pixel values are normalized to [0, 1] with a single channel.
In each game, for each frame in the sampled data set, we verify if the DQN does not change its actions when an `∞ perturbation of various radii is added to the frame which is the state observation of the agent. The number of actions in BankHeist, Roadrunner and Pong are 6, 6, 4, respectively. Essentially, we reduce the verification of the DQN with finite discrete action space to the same setting as verifying classifiers. Therefore, we apply the same verification process as descried in Appendix F.3 with ∈ {0.1/255, 0.2/255, · · · , 1.0/255} on
3(Zhang et al., 2020) consider one additional RL setting (Freeway). However, the released PGD-trained DQN is completely unverifiable for nearly all epsilons that we considered.
all three tasks. The numbers of certified examples of ADMM and Linear (Wong et al., 2018; Xu et al., 2020) and their differences on each task are plotted in Figure 8.

heading: G Additional details on comparison to Lagrangian-based methods
text: In Section 4.1, we compared the certified test accuracy of ADMM with a suite of Lagrangianbased methods (Bunel et al., 2020a) on CIFAR10. These Lagrangian-based baselines were implemented using the codes at
https://github.com/oval-group/decomposition-plnn-bounds
under the MIT license. As mentioned in Section 3.2, these baselines require solving an inner optimization problem through iterative algorithms. In the experiments of Table 1, the number of iterations of different iterative algorithms are bounded separately such that each Lagrangian-based method listed in Section 4.1 has runtime of 9 seconds to finish verifying one example.
Our ADMM solver averages 9 seconds runtime per example in this verification task, which is the same as the average runtime of the Lagrangian-based methods, with the stopping criterion of abs = 10−4, rel = 10−3 and ρ = 1.0.
Architecture size. We note that the architecture we verify in this paper is 10x larger than considered by (Bunel et al., 2020b). However, the released framework could not handle this model size due to their implementation of intermediate bounds. In order to handle this model size, we externally calculated the intermediate bounds (using (Wong et al., 2018)) and loaded these into the framework for the Lagrangian-based baselines.

heading: H Speed comparison
text: Comparison with LP solver: To demonstrate the effectiveness of GPU-acceleration in the DeepSplit algorithm, we compare the runtime of DeepSplit and Gurobi in solving LP relaxations that bound the output range of MNIST-A network (the fully connected network defined in Appendix F.1) with `∞ perturbations in the input. Specifically, for a given example in the MNIST test data set, we apply DeepSplit/Gurobi layer-by-layer to find the tightest pre-activation bounds under the LP-relaxation.
Recall that the MNIST-A network is a fully-connected netowrk with architecture 784− 600− 400− 200− 100− 10. With the Gurobi solver, we need to solve 2× 600 LPs sequentially to obtain the lower and upper bounds for the first activation layer, 2× 400 LPs for the second activation layer, and so forth. With DeepSplit, the pre-activation bounds can be computed in batch and allows GPU-acceleration.
In our experiment, we fix the radius of the `∞ perturbation at the input image as = 0.02. For the Gurobi solver, we randomly choose 10 samples from the test data set and compute the pre-activation bounds layer-by-layer. The LP relaxation is formulated in CVXPY and solved by Gurobi v9.1 on an Intel Core i7-6700K CPU, which has 4 cores and 8 threads. For each example, the total solver time of Gurobi is recorded with the average solver time being 275.9 seconds. For the DeepSplit method, we compute the pre-activation bounds layer-by-layer on 19 randomly chosen examples. The algorithm applies residual balancing with the initial ρ = 1.0 and the stopping criterion is given by abs = 10−4, rel = 10−3. The total running time of DeepSplit is 717.9 seconds, with 37.8 seconds per example on average. With the GPU-acceleration, out method achieves 7x speedup in verifying NN properties compared with the commercial-grade Gurobi solver.
Comparison with SDP solver: Our method is significantly faster than solving the convex relaxation through other approaches. For example, it took our method approximately 20 hours to verify all MNIST test set images for the convolutional architecture at 10 different epsilon values (See Table 3), taking on average 2 hours per epsilon to verify all 10,000 examples. See F.2 for the running time summary of each experiment. In contrast, the SDP relaxation (Raghunathan et al., 2018) takes 3 hours to verify 500 examples at one epsilon value for a similarly sized CNN. This highlights the difference in speed and scalability of solving the LP relaxation over the SDP one, albeit at the cost of looser verification guarantees.
Comparison with branch-and-bound solver: Branch-and-bound methods have exponential runtime as opposed to the polynomial runtime of LP-solvers. In practice, this leads to significantly longer solve times that has largely limited these methods to the small network regime. For example, the fastest branch-and-bound method from Bunel et al. (2020a) takes on average 6 minutes per example to verify a small convolutional network, whereas our solver takes on average 9 seconds per example (40x faster) to verify an order-of-magnitude larger network.

heading: I Scalability
text: To highlight the scalability of our approach, we consider a ResNet18 network trained on CIFAR10 whose max pooling layer is replaced by a down-sampling convolutional layer for comparison with LiRPA (Xu et al., 2020) 4 (codes available at
https://github.com/KaidiXu/auto_LiRPA
under the BSD 3-Clause "New" or "Revised" license), which is capable of computing provable linear bounds for the outputs of general neural networks and is the only method available so far that can handle ResNet18. The ResNet18 is adversarially trained using the fast adversarial training code from (Wong et al., 2020).
In our experiments, for the first 100 test examples in CIFAR10, we use LiRPA to compute the preactivation bounds for each ReLU layer in ResNet18 and then apply ADMM to compute the lower and upper bounds of ResNet18 outputs (there are 10 outputs corresponding to the 10 classes of the dataset). The ADMM is run with stopping criterion abs = 10−5, rel = 10−4 and augmentation parameter ρ = 1.0. With = 1/255, the lower/upper bounds of ADMM from these 100 examples are arranged in an ascending order in Figure 9 and are compared with those obtained by LiRPA. We observe that the bounds from ADMM are consistently tighter than those from LiRPA. The average running time of ADMM across the 100 examples is 2312 seconds. Overall bound statistics are shown in Table 5.
4The max pooling layer has not been considered in the implementation of LiRPA by the submission of this paper.

Review summaries:
This paper proposes an operator splitting method which solves a convex relaxation of a learge-scale nonconvex optimization problem for analyzing the worst-case performance of deep neural networks against input perturbations. The proposed method is modular, scalable and parallelizable. Experimentally, the authors demonstrate that the proposed method has tighter bounds on the worst-case performance of large CNNs in image classification and RL settings. 

The paper proposes a novel method for the neural network verification problem. The authors use variable splitting and Lagrangian relaxation to relax the verification problem, which is then solved with ADMM. The method is GPU friendly (similar to previous work like Dvijotham et al. 2018 and Bunel et al. 2020a), and claims to solve the relaxation to optimality.

This paper proposes an efficient solver for computing the convex relaxation of Neural Networks, which is an important component of verification methods. The method is based on ADMM, with the authors proposing a novel decomposition of the optimization problem, which allows them to get closed form solutions of all the intermediate steps of the algorithm, The proposed method also improves on the guaranteed convergence rate and is shown to scale to larger networks than the ones considered in other works.  Evaluation is performed on robustness verification on Cifar10 and on Reinforcement learning models.

The authors present DeepSplit, a novel solver for a popular neural network convex relaxation. Relying on ADMM and on a careful problem re-formulation, the authors achieve both a O(1/T) converge rate and closed-form solutions for the inner problems, differently from previously presented solvers for the same relaxation. Computational results are presented for incomplete neural network verification, showing that DeeSplit scales to a ResNet18 and achieves better bounds in the same time compared to relevant baselines.