My stuff
Efficient and Accurate Estimation of Lipschitz Constants for Deep Neural Networks

Tight estimation of the Lipschitz constant for deep neural networks (DNNs) is useful in many applications ranging from robustness certification of classifiers to stability analysis of closed-loop systems with reinforcement learning controllers. Existing methods in the literature for estimating the Lipschitz constant suffer from either lack of accuracy or poor scalability. In this paper, we present a convex optimization framework to compute guaranteed upper bounds on the Lipschitz constant of DNNs both accurately and efficiently. Our main idea is to interpret activation functions as gradients of convex potential functions. Hence, they satisfy certain properties that can be described by quadratic constraints. This particular description allows us to pose the Lipschitz constant estimation problem as a semidefinite program (SDP). The resulting SDP can be adapted to increase either the estimation accuracy (by capturing the interaction between activation functions of different layers) or scalability (by decomposition and parallel implementation). We illustrate the utility of our approach with a variety of experiments on randomly generated networks and on classifiers trained on the MNIST and Iris datasets. In particular, we experimentally demonstrate that our Lipschitz bounds are the most accurate compared to those in the literature. We also study the impact of adversarial training methods on the Lipschitz bounds of the resulting classifiers and show that our bounds can be used to efficiently provide robustness guarantees.

A function  f : Rn → Rmis globally Lipschitz continuous on  X ⊆ Rnif there exists a nonnegative constant L ≥0 such that


The smallest such L is called the Lipschitz constant of f. The Lipschitz constant is the maximum ratio between variations in the output space and variations in the input space of f and thus is a measure of sensitivity of the function with respect to input perturbations.

When a function f is characterized by a deep neural network (DNN), tight bounds on its Lipschitz constant can be extremely useful in a variety of applications. In classification tasks, for instance, L can be used as a certificate of robustness of a neural network classifier to adversarial attacks if it is estimated tightly [33]. In deep reinforcement learning, tight bounds on the Lipschitz constant of a DNN-based controller can be directly used to analyze the stability of the closed-loop system. Lipschitz regularity can also play a key role in derivation of generalization bounds [5]. In these applications and many others, it is essential to have tight bounds on the Lipschitz constant of DNNs. However, as DNNs have highly complex and nonlinear structures, estimating the Lipschitz constant both accurately and efficiently has remained a significant challenge.

Our contributions. In this paper we propose a novel convex programming framework to derive tight bounds on the global Lipschitz constant of deep feed-forward neural networks. Our framework yields significantly more accurate bounds compared to the state-of-the-art and lends itself to a distributed implementation, leading to efficient computation of the bounds for large-scale networks.

Our approach. We use the fact that nonlinear activation functions used in neural networks are gradients of convex functions; hence, as operators, they satisfy certain properties that can be abstracted as quadratic constraints on their input-output values. This particular abstraction allows us to pose the Lipschitz estimation problem as a semidefinite program (SDP), which we call LipSDP. A striking feature of LipSDP is its flexibility to span the trade-off between estimation accuracy and computational efficiency by adding or removing extra decision variables. In particular, for a neural network with  ℓlayers and a total of n hidden neurons, the number of decision variables can vary from  ℓ(least accurate but most scalable) to O(n) (most accurate but least scalable). As such, we derive several distinct yet related formulations of LipSDP that span this trade-off. To scale each variant of LipSDP to larger networks, we also propose a distributed implementation.

Our results. We illustrate our approach in a variety of experiments on both randomly generated networks as well as networks trained on the MNIST [22] and Iris [10] datasets. First, we show empirically that our Lipschitz bounds are the most accurate compared to all other existing methods of which we are aware. In particular, our experiments on neural networks trained for MNIST show that our bounds almost coincide with the true Lipschitz constant and outperform all comparable methods. For details, see Figure 2a. Furthermore, we investigate the effect of two robust training procedures [20, 23] on the Lipschitz constant for networks trained on the MNIST dataset. Our results suggest that robust training procedures significantly decrease the Lipschitz constant of the resulting classifiers. Moreover, we use the Lipschitz bound for two robust training procedures to derive non-vacuous lower bounds on the minimum adversarial perturbation necessary to change the classification of any instance from the test set. For details, see Figure 3.

Related work. The problem of estimating the Lipschitz constant for neural networks has been studied in several works. In [33], the authors estimate the global Lipschitz constant of DNNs by the product of Lipschitz constants of individual layers. This approach is scalable and general but yields trivial bounds. We are only aware of two other methods that give non-trivial upper bounds on the global Lipschitz constant of fully-connected neural networks and can scale to networks with more than two hidden layers. In [9], Combettes and Pesquet derive bounds on Lipschitz constants by treating the activation functions as nonexpansive averaged operators. The resulting algorithm scales well with the number of hidden units per layer, but very poorly (in fact exponential) with the number of layers. In [36], Virmaux and Scaman decompose the weight matrices of a neural network via singular value decomposition and approximately solve a convex maximization problem over the unit cube. Notably, estimating the Lipschitz constant using the method in [36] is intractable even for small networks; indeed, the authors of [36] use a greedy algorithm to compute a bound, which may underapproximate the Lipschitz constant. Bounding Lipschitz constants for the specific case of convolutional neural networks (CNNs) has also been addressed in [4,5,42].

Using Lipschitz bounds in the context of adversarial robustness and safety verification has also been addressed in several works [30, 37, 38]. In particular, in [38], the authors convert the robustness analysis problem into a local Lipschitz constant estimation problem, where they estimate this local constant by a set of independently and identically sampled local gradients. This algorithm is scalable but is not guaranteed to provide upper bounds. In a similar work, the authors of [37] exploit the piece-wise linear structure of ReLU functions to estimate the local Lipschitz constant of neural networks. In [12], the authors use quadratic constraints and semidefinite programming to analyze local (point-wise) robustness of neural networks. In contrast, our Lipschitz bounds can be used as a global certificate of robustness and are agnostic to the choice of the test data.

1.1 Motivating applications

We now enumerate several applications that highlight the importance of estimating the Lipschitz constant of DNNs accurately and efficiently.

Robustness certification of classifiers. In response to fragility of DNNs to adversarial attacks, there has been considerable effort in recent years to improve the robustness of neural networks against adversarial attacks and input perturbations [15,20,21,23,25,41]. In order to certify and/or improve the robustness of neural networks, one must be able to bound the possible outputs of the neural network over a region of input space. This can be done either locally around a specific input [6,11,12,14,19,20,28,29,32,34,39,40], or globally by bounding the sensitivity of the function to input perturbations, i.e., the Lipschitz constant [17,27,33,38]. Indeed, tight upper bounds on the Lipschitz constant can be used to derive non-vacuous lower bounds on the magnitudes of perturbations necessary to change the decision of neural networks. Finally, an efficient computation of these bounds can be useful in either assessing robustness after training [12,28,29] or promoting robustness during training [20,35]. In the experiments section, we explore this application in depth.

Stability analysis of closed-loop systems with learning controllers. A central problem in learning-based control is to provide stability or safety guarantees for a feedback control loop when a learning-enabled component, such as a deep neural network, is introduced in the loop [3,7,18]. The Lipschitz constant of a neural network controller bounds its gain. Therefore a tight estimate can be useful for certifying the stability of the closed-loop system.

Notation. We denote 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. We denote by  Sn, Sn+, and  Sn++the sets of n-by-n symmetric, positive semidefinite, and positive definite matrices, respectively. The p-norm (p ≥1) is denoted by  ∥ · ∥p : Rn → R+. The  ℓ2-norm of a matrix  W ∈ Rm×nis the largest singular value of W. We denote the i-th unit vector in  Rnby  ei. We write diag(a1, ..., an) for a diagonal matrix whose diagonal entries starting in the upper left corner are  a1, · · · , an.

2.1 Problem statement

Consider an  ℓ-layer feed-forward neural network  f(x): Rn0 → Rnℓ+1described by the following recursive equations:


Here  x ∈ Rn0is an input to the network and  W k ∈ Rnk+1×nkand  bk ∈ Rnk+1are the weight matrix and bias vector for the k-th layer. The function  φis the concatenation of activation functions at each layer, i.e., it is of the form  φ(x) = [ϕ(x1) · · · ϕ(xn)]⊤. In this paper, our goal is to find tight bounds on the Lipschitz constant of the map  x �→ f(x) in  ℓ2-norm. More precisely, we wish to find the smallest constant  L2 ≥0 such that  ∥f(x) − f(y)∥2 ≤ L2∥x − y∥2for all  x, y ∈ Rn0.

The main source of difficulty in solving this problem is the presence of the nonlinear activation functions. To combat this difficulty, our main idea is to abstract these activation functions by a set of constraints that they impose on their input and output values. Then any property (including Lipschitz continuity) that is satisfied by our abstraction will also be satisfied by the original network.

2.2 Description of activation functions by quadratic constraints

In this section, we introduce several definitions and lemmas that characterize our abstraction of nonlinear activation functions. These results are crucial to the formulation of an SDP that can bound the Lipschitz constants of networks in Section 2.3.

Definition 1 (Slope-restricted non-linearity)  A function ϕ: R → R is slope-restricted on [α, β] where0  ≤ α < β < ∞ if


The inequality in (3) simply states that the slope of the chord connecting any two points on the curve of the function  x �→ ϕ(x) is at least  αand at most  β(see Figure 1). By multiplying all sides of (3) by (y − x)2, we can write the slope restriction condition as  α(y − x)2 ≤ (ϕ(y) − ϕ(x))(y − x) ≤ β(y − x)2. By the left inequality, the operator  ϕ(x) is strongly monotone with parameter  α[31], or equivalently the anti-derivative function�ϕ(x)dxis strongly convex with parameter  α. By the right-hand side inequality,  ϕ(x) is one-sided Lipschitz with parameter  β. Altogether, the preceding inequalities state that the anti-derivative function �ϕ(x)dxis  α-strongly convex and  β-smooth.

Note that all activation functions used in deep learning satisfy the slope restriction condition in (3) for some 0  ≤ α < β < ∞. For instance, the ReLU, tanh, and sigmoid activation functions are all slope restricted with  α = 0 and β = 1. More details can be found in [12].

Definition 2 (Incremental Quadratic Constraint [1]) A function  φ: Rn → Rnsatisfies the incremental quadratic constraint defined by  Q ⊂ S2nif for any  Q ∈ Qand  x, y ∈ Rn,


In the above definition, Q is the set of all multiplier matrices that characterize  φ, and is a convex cone by definition. As an example, the softmax operator  φ(x) = (�ni=1exp(xi))−1[exp(x1) · · ·exp(xn)]⊤is the gradient of the convex function  ψ(x) = log(�ni=1exp(xi)). This function is smooth and strongly convex with paramters  α = 0 and β = 1 [8]. For this class of functions, it is known that the gradient functionφ(x) = ∇ψ(x) satisfies the quadratic inequality [24]


Therefore, the softmax operator satisfies the incremental quadratic constraint defined by  Q = {λM | λ ≥ 0}, where M the middle matrix in the above inequality.

To see the connection between incremental quadratic constraints and slope-restricted nonlinearities, note that (3) can be equivalently written as the single inequality


which, in view of Definition 2, is an incremental quadratic constraint for  ϕ. From this perspective, incremental quadratic constraints generalize the notion of slope-restricted nonlinearities to multi-variable vector-valued nonlinearities.

Lemma 1 Suppose  ϕ: R → Ris slope-restricted on [α, β]. Define the set


We will see in the next section that the matrix T that parameterizes the multiplier matrix in (9) appears as a decision variable in an SDP, in which the objective is to find an admissible T that yields the tightest bound on the Lipschitz constant.

2.3 LipSDP for single-layer neural network

To develop an optimization problem to estimate the Lipschitz constant of a fully-connected feed-forward neural network, the key insight is that the Lipschitz condition in (1) is in fact equivalent to an incremental quadratic constraint for the map  x �→ f(x) characterized by the neural network. By coupling this to the incremental quadratic constraints satisfied by the cascade combination of the activation functions [13], we can develop an SDP to minimize an upper bound on the Lipschitz constant of f. This result is formally stated in the following theorem.


Figure 1: An illustrative description of encoding activation functions by quadratic constraints.

Theorem 1 (Lipshitz certificates for single-layer neural networks) Consider a single-layer neural network described by  f(x) = W 1φ(W 0x + b0) +  b1. Suppose  φ(x): Rn → Rn = [ϕ(x1) · · · ϕ(xn)], where ϕis slope-restricted in the sector [α, β]. Define  Tnas in (8). Suppose there exists a  ρ >0 such that the matrix inequality


holds for some  T ∈ Tn. Then  ∥f(x) − f(y)∥2 ≤ √ρ∥x − y∥2for all  x, y ∈ Rn0.

Theorem 1 provides us with a sufficient condition for  L2 = √ρto be an upper bound on the Lipschitz constant of  f(x) = W 1φ(W 0x+b0)+b1. In particular, we can find the tightest bound by solving the following optimization problem:


where the decision variables are (ρ, T) ∈ R+ × Tn. Note that  M(ρ, T) is linear in  ρand T and the set  Tnis convex. Hence, (11) is an SDP, which can be solved numerically for its global minimum.

2.4 LipSDP for multi-layer neural networks

We now consider the multi-layer case. Assuming that all the activation functions are the same, we can write the neural network model in (2) compactly as


where  x = [x0⊤ x1⊤ · · · xℓ⊤]⊤is the concatenation of the input and the activation values, and the matrices b, A, B and C are given by [12]


The particular representation in (12) facilitates the extension of LipSDP to multiple layers, as stated in the following theorem.

Theorem 2 (Lipschitz certificates for multi-layer neural networks) Consider an  ℓ-layer fully connected neural network described by (2). Let  n = �ℓk=1 nkbe the total number of hidden neurons and suppose the activation functions are slope-restricted in the sector [α, β]. Define  Tnas in (8). Define A and B as in (13). Consider the matrix inequality


If (14) is satisfied for some (ρ, T) ∈ R+ × Tn, then  ||f(x) − f(y)||2 ≤ √ρ ||x − y||2, ∀x, y ∈ Rn0.

In a similar way to the single-layer case, we can find the best bound on the Lipschitz constant by solving the SDP in (11) with  M(ρ, T) defined as in (14).

Remark 1 We have only considered the  ℓ2norm in our exposition. By using the inequality  ∥x∥p ≤


or, equivalently,  ∥f(y)−f(x)∥p ≤ n( 12 − 1q )m( 1p − 12 )L2∥y−x∥q.Hence,  n( 12 − 1q )m( 1p − 12 )L2is a Lipschitz constant of f when  ℓqand  ℓpnorms are used in the input and output spaces, respectively. We can also extend our framework to accommodate quadratic norms  ∥x∥P =√x⊤Px, where  P ∈ Sn++.

2.5 Variants of LipSDP: reconciling accuracy and efficiency

In LipSDP, there are n + 1 decision variables  ρ ≥ 0, λii, 1 ≤ i ≤ n, where n is the total number of hidden neurons. Using all these decision variables would provide the tightest convex relaxation in our formulation. However, solving this SDP with all the decision variables included is impractical for large networks. Nevertheless, we can consider a hierarchy of relaxations of LipSDP by removing a subset of the decision variables. Below, we give a brief description of the efficiency and accuracy of each variant. Throughout, we let n be the total number of neurons and  ℓthe number of hidden layers.

1. LipSDP-Neuron has n + 1 decision variables. For this case, we have  T = diag(λ11, · · · , λnn).

2. LipSDP-Layer considers only one constraint per layer, resulting in  ℓ+ 1 decision variables. It is the most scalable and least accurate method. For this variant, we have  T = blkdiag(λ1In1, · · · , λℓInℓ).

Parallel implementation by splitting. The Lipschitz constant of the composition of two or more functions can be bounded by the product of the Lispchtiz constants of the individual functions. By splitting a neural network up into small sub-networks, one can first bound the Lipschitz constant of each sub-network and then multiply these constants together to obtain a Lipschitz constant for the entire network. Because sub-networks do not share weights, it is possible to compute the Lipschitz constants for each sub-network in parallel. This greatly improves the scalability of of each variant of LipSDP with respect to the total number of activation functions in the network.

Repeated nonlinearities. The quadratic constraint in Lemma 1 is valid for functions of the form  φ(x) =[ϕ1(x1) · · · ϕn(xn)]⊤as along as all  ϕi’s are slope-restricted on [α, β]. Therefore, this quadratic constraint does not exploit the fact that the same activation functions are used. In Appendix A.8 we derive more refined quadratic constraints that exploit this additional structure when one of the points in the definition of Lipschitz continuity is fixed. The resulting SDP would be tighter but more expensive to solve.


Table 1: Computation time in seconds for evaluating Lipschitz bounds of one-hidden-layer neural networks with a varying number of hidden units. A plot showing the Lipschitz constant for each network tested in this table has been provided in the Appendix.

Table 2: Computation time in seconds for computing Lipschitz bounds of  ℓ-hidden-layer neural networks with 100 activation functions per layer. For LipSDP-Neuron and LipSDP-Layer, we split each network up into 5-layer sub-networks.


In this section we describe several experiments that highlight the key aspects of this work. In particular, we show empirically that our bounds are much tighter than any comparable method, we study the impact of robust training on our Lipschitz bounds, and we analyze the scalability of our methods.

Experimental setup. For our experiments we used MATLAB, the CVX toolbox [16] and MOSEK [2] on a 9-core CPU with 16GB of RAM to solve the SDPs. All classifiers trained on MNIST used an 80-20 train-test split.

Training procedures. Several training procedures have recently been proposed to improve the robustness of neural network classifiers. Two prominent procedures are the LP-based method in [20] and projected gradient descent (PGD) based method in [23]. We refer to these training methods as LP-Train and PGD-Train, respectively. Both procedures take as input a parameter  ϵthat defines the  ℓ∞perturbation of the training data points.

Baselines. Throughout the experiments, we will often show comparisons to what we call the naive lower and upper bounds. As has been shown in several previous works [9, 36], trivial lower and upper bounds on the Lipschitz constant of a feed-forward neural network with  ℓhidden-layers are given by  L2, lower =����W ℓW ℓ−1 · · · W 0����2and  L2, upper = �ℓi=0����W i����2, which we refer to them as naive lower and upper bounds, respectively. We are aware of only two methods that bound the Lipschitz constant and can scale to fully-connected networks with more than two hidden layers; these methods are [9], which we will refer to as CPLip, and [36], which is called SeqLip.

We compare the Lipschitz bounds obtained by LipSDP-Neuron, LipSDP-Layer, CPLip, and SeqLip in Figure 2a. It is evident from this figure that the bounds from LipSDP-Neuron are tighter than CPLip and SeqLip. Our results show that the true Lipschitz constants of the networks shown above are very close to the naive lower bound.

To demonstrate the scalability of the LipSDP formulations, we split a 100-hidden layer neural network into sub-networks with six hidden layers each and computed the Lipschitz bounds using LipSDP-Neuron and LipSDP-Layer. The results are shown in Figure 2b. Furthermore, in Tables 1 and 2, we show the computation time for scaling the LipSDP methods in the number of hidden units per layer and in the number of layers. In particular, the largest network we tested in Table 2 had 50,000 hidden neurons; SDPLip-Neuron took approximately 12 minutes to find a Lipschitz bound, and SDPLip-Layer took approximately 4 minutes.

Impact of robust training. In Figure 3, we empirically demonstrate that the Lipschitz bound of a neural network is directly related to the robustness of the corresponding classifier. This figure shows that LP-train and PGD-Train networks achieve lower Lipschitz bounds than standard training procedures. Figure 3a indicates that robust training procedures yield lower Lipschitz constants than networks trained with standard training procedures such as the Adam optimizer. Figure 3b shows the utility of sharply estimating


Figure 2: Comparison of the accuracy LipSDP methods to other methods that compute the Lipschitz constant and scalability analysis of all three SeqLip methods.


Figure 3: Analysis of impact of robust training on the Lipschitz constant and the distance to misclassification for networks trained on MNIST

the Lipschitz constant; a lower value of  L2guarantees that a neural network is more locally robust to input perturbations; see Proposition 1 in the Appendix.

In the same vein, Figure 4 shows the impact of varying the robustness parameter  ϵused in LP-Train and PGD-Train on the test accuracy of networks trained for a fixed number of epochs and the corresponding Lipschitz constants. In essence, these results quantify how much robustness a fixed classifier can handle before accuracy plummets. Interestingly, the drops in accuracy as  ϵincreases coincide with corresponding drops in the Lipschitz constant for both LP-Train and PGD-Train.

Robustness for different activation functions. The framework proposed in this work allows us to examine the impact of using different activation functions on the Lipschitz constant of neural networks. We trained two sets of neural networks on the MNIST dataset. The first set used ReLU activation functions,


Figure 4: Trade-off between accuracy and Lipschitz constant for different values of the robustness parameter used for LP-Train and PGD-Train. All networks had one hidden layer with 50 hidden neurons.

Figure 5: Lipschitz constants for topologically identical three-hidden-layer networks with ReLU and leaky ReLU activation functions. All classifiers were trained until they reached 97% test accuracy.

while the second set used leaky ReLU activations. Figure 5 shows empirically that the networks with the leaky ReLU activation function have larger Lipschitz constants than networks of the same architecture with the ReLU activation function.

In this paper, we proposed a hierarchy of semidefinite programs to derive tight upper bounds on the Lipschitz constant of feed-forward fully-connected neural networks. Some comments are in order. First, our framework can be directly used to certify convolutional neural networks (CNNs) by unrolling them to a large feed-forward neural network. A future direction is to exploit the special structure of CNNs in the resulting SDP. Second, we only considered one application of Lipschitz bounds in depth (robustness certification). Having an accurate upper bound on the Lipschitz constant can be useful in domains beyond robustness analysis, such as stability analysis of feedback systems with control policies updated by deep reinforcement learning. Furthermore, Lipschitz bounds can be utilized during training as a heuristic to promote out-of-sample generalization [35]. We intend to pursue these applications for future work.

[1] Beh¸cet A¸cıkme¸se and Martin Corless. Observers for systems with nonlinearities satisfying incremental quadratic constraints. Automatica, 47(7):1339–1348, 2011.

[2] MOSEK ApS. The MOSEK optimization toolbox for MATLAB manual. Version 8.1., 2017.

[3] Anil Aswani, Humberto Gonzalez, S Shankar Sastry, and Claire Tomlin. Provably safe and robust learning-based model predictive control. Automatica, 49(5):1216–1226, 2013.

[4] Radu Balan, Maneesh Singh, and Dongmian Zou. Lipschitz properties for deep convolutional networks. arXiv preprint arXiv:1701.05217, 2017.

[5] Peter L Bartlett, Dylan J Foster, and Matus J Telgarsky. Spectrally-normalized margin bounds for neural networks. In Advances in Neural Information Processing Systems, pages 6240–6249, 2017.

[6] Osbert Bastani, Yani Ioannou, Leonidas Lampropoulos, Dimitrios Vytiniotis, Aditya Nori, and Anto- nio Criminisi. Measuring neural net robustness with constraints. In Advances in neural information processing systems, pages 2613–2621, 2016.

[7] Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause. Safe model-based rein- forcement learning with stability guarantees. In Advances in neural information processing systems, pages 908–918, 2017.

[8] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.

[9] Patrick L. Combettes and Jean-Christophe Pesquet. Lipschitz certificates for neural network structures driven by averaged activation operators. arXiv preprint arXiv:1903.01014, 2019.

[10] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017.

[11] Souradeep Dutta, Susmit Jha, Sriram Sankaranarayanan, and Ashish Tiwari. Output range analysis for deep feedforward neural networks. In NASA Formal Methods Symposium, pages 121–138. Springer, 2018.

[12] Mahyar Fazlyab, Manfred Morari, and George J Pappas. Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming. arXiv preprint arXiv:1903.01287, 2019.

[13] Mahyar Fazlyab, Alejandro Ribeiro, Manfred Morari, and Victor M Preciado. Analysis of optimiza- tion algorithms via integral quadratic constraints: Nonstrongly convex problems. SIAM Journal on Optimization, 28(3):2654–2689, 2018.

[14] Timon Gehr, Matthew Mirman, Dana Drachsler-Cohen, Petar Tsankov, Swarat Chaudhuri, and Martin Vechev. Ai2: Safety and robustness certification of neural networks with abstract interpretation. In 2018 IEEE Symposium on Security and Privacy (SP), pages 3–18. IEEE, 2018.

[15] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy. Explaining and harnessing adversarial examples (2014). arXiv preprint arXiv:1412.6572.

[16] Michael Grant, Stephen Boyd, and Yinyu Ye. Cvx: Matlab software for disciplined convex programming, 2008.

[17] Todd Huster, Cho-Yu Jason Chiang, and Ritu Chadha. Limitations of the lipschitz constant as a defense against adversarial examples. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 16–29. Springer, 2018.

[18] Ming Jin and Javad Lavaei. Stability-certified reinforcement learning: A control-theoretic perspective. arXiv preprint arXiv:1810.11505, 2018.

[19] Matt Jordan, Justin Lewis, and Alexandros G Dimakis. Provable certificates for adversarial examples: Fitting a ball in the union of polytopes. arXiv preprint arXiv:1903.08778, 2019.

[20] J Zico Kolter and Eric Wong. Provable defenses against adversarial examples via the convex outer adversarial polytope. arXiv preprint arXiv:1711.00851, 1(2):3, 2017.

[21] Alexey Kurakin, Ian Goodfellow, and Samy Bengio. Adversarial machine learning at scale. arXiv preprint arXiv:1611.01236, 2016.

[22] Yann LeCun. The mnist database of handwritten digits. http://yann. lecun. com/exdb/mnist/, 1998.

[23] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. Towards deep learning models resistant to adversarial attacks. arXiv preprint arXiv:1706.06083, 2017.

[24] Yurii Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science & Business Media, 2013.

[25] Nicolas Papernot, Patrick McDaniel, Xi Wu, Somesh Jha, and Ananthram Swami. Distillation as a defense to adversarial perturbations against deep neural networks. In 2016 IEEE Symposium on Security and Privacy (SP), pages 582–597. IEEE, 2016.

[26] Jonathan Peck, Joris Roels, Bart Goossens, and Yvan Saeys. Lower bounds on the robustness to adversarial perturbations. In Advances in Neural Information Processing Systems, pages 804–813, 2017.

[27] Haifeng Qian and Mark N Wegman. L2-nonexpansive neural networks. arXiv preprint arXiv:1802.07896, 2018.

[28] Aditi Raghunathan, Jacob Steinhardt, and Percy Liang. Certified defenses against adversarial examples. arXiv preprint arXiv:1801.09344, 2018.

[29] Aditi Raghunathan, Jacob Steinhardt, and Percy S Liang. Semidefinite relaxations for certifying robust- ness to adversarial examples. In Advances in Neural Information Processing Systems, pages 10900–10910, 2018.

[30] Wenjie Ruan, Xiaowei Huang, and Marta Kwiatkowska. Reachability analysis of deep neural networks with provable guarantees. arXiv preprint arXiv:1805.02242, 2018.

[31] Ernest K Ryu and Stephen Boyd. Primer on monotone operator methods. Appl. Comput. Math, 15(1):3–43, 2016.

[32] Gagandeep Singh, Timon Gehr, Matthew Mirman, Markus P¨uschel, and Martin Vechev. Fast and effective robustness certification. In Advances in Neural Information Processing Systems, pages 10802– 10813, 2018.

[33] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus. Intriguing properties of neural networks. arXiv preprint arXiv:1312.6199, 2013.

[34] Vincent Tjeng, Kai Xiao, and Russ Tedrake. Evaluating robustness of neural networks with mixed integer programming. arXiv preprint arXiv:1711.07356, 2017.

[35] Yusuke Tsuzuku, Issei Sato, and Masashi Sugiyama. Lipschitz-margin training: Scalable certification of perturbation invariance for deep neural networks. In Advances in Neural Information Processing Systems, pages 6541–6550, 2018.

[36] Aladin Virmaux and Kevin Scaman. Lipschitz regularity of deep neural networks: analysis and efficient estimation. In Advances in Neural Information Processing Systems, pages 3835–3844, 2018.

[37] Tsui-Wei Weng, Huan Zhang, Hongge Chen, Zhao Song, Cho-Jui Hsieh, Duane Boning, Inderjit S Dhillon, and Luca Daniel. Towards fast computation of certified robustness for relu networks. arXiv preprint arXiv:1804.09699, 2018.

[38] Tsui-Wei Weng, Huan Zhang, Pin-Yu Chen, Jinfeng Yi, Dong Su, Yupeng Gao, Cho-Jui Hsieh, and Luca Daniel. Evaluating the robustness of neural networks: An extreme value theory approach. arXiv preprint arXiv:1801.10578, 2018.

[39] Eric Wong, Frank Schmidt, Jan Hendrik Metzen, and J Zico Kolter. Scaling provable adversarial defenses. In Advances in Neural Information Processing Systems, pages 8400–8409, 2018.

[40] Huan Zhang, Tsui-Wei Weng, Pin-Yu Chen, Cho-Jui Hsieh, and Luca Daniel. Efficient neural network robustness certification with general activation functions. In Advances in Neural Information Processing Systems, pages 4939–4948, 2018.

[41] Stephan Zheng, Yang Song, Thomas Leung, and Ian Goodfellow. Improving the robustness of deep neural networks via stability training. In Proceedings of the ieee conference on computer vision and pattern recognition, pages 4480–4488, 2016.

[42] Dongmian Zou, Radu Balan, and Maneesh Singh. On lipschitz bounds of general convolutional neural networks. arXiv preprint arXiv:1808.01415, 2018.

A.1 Robustness certification of DNN-based classifiers

Consider a classifier described by a feed-forward neural network  f : Rn → Rk, where n is the number of input features and k is the number of classes. In this context, the function f takes as input an instance or measurement x and returns a k-dimensional vector of scores – one for each class. The classification rule is based on assigning x to the class with the highest score. That is, we define the classification C(x) :  Rn →{1, . . . , k} to be  C(x) = argmax1≤i≤k fi(x). Now suppose that  x⋆is an instance that is classified correctly by the neural network. To evaluate the local robustness of the neural network around  x⋆, we consider a bounded set  A(x⋆) = {x | ∥x − x⋆∥2 ≤ ϵ} ⊆ Rnthat represents the set of all possible  ℓ2-norm perturbations of  x⋆. Then the classifier is locally robust at  x⋆against  A(x⋆) if it assigns all the perturbed inputs to the same class as the unperturbed input, i.e., if


In the following proposition, we derive a sufficient condition to guarantee local robustness around  x⋆for the perturbation set  A(x⋆).

Proposition 1 Consider a neural-network classifier  f : Rn → Rkwith Lipschitz constant  L2in the  ℓ2-norm. Let  ϵ >0 be given and consider the inequality


Then (16)  implies C(x) := argmax1≤i≤kfi(x) = C(x⋆) for all x ∈ A(x⋆).

The inequality in (16) provides us with a simple and computationally efficient test for assessing the point-wise robustness of a neural network. According to (16), a more accurate estimation of the Lipschitz constant directly increases the maximum perturbation  ϵthat can be certified for each test example. This makes the framework suitable for model selection, wherein one wishes to select the model that is most robust to adversarial perturbations from a family of proposed classifiers [26].

A.2 Proof of Proposition 1

Let  i⋆ = argmax1≤i≤k fi(x⋆) be the class of  x⋆. Define the polytope in the output space of f:


which is the set of all outputs whose score is the highest for class  i⋆; and, of course,  f(x⋆) ∈ Pi⋆. The distance of  f(x⋆) to the boundary  ∂Piof the polytope is the minimum distance of  f(x⋆) to all edges of the polytope:


Note that the Lipschitz condition implies  dist(f(x⋆), f(x)) = ∥f(x) − f(x⋆)∥2 ≤ L2ϵfor all  ∥x − x⋆∥2 ≤ ϵ. The condition in (16) then implies


and for all 1  ≤ j ≤ k, j ̸= i⋆. Therefore, the output of the classification would not change for any ∥x − x⋆∥2 ≤ ϵ.


Figure 6: Illustration of local robustness certification using the Lipschitz bound.

A.3 Proof of Lemma 1

Let  T = diag(λ11, · · · , λnn) with  λii ≥0 for all 1  ≤ i ≤ n. By expanding the quadratic form on the right-hand side of (9), we can write


where the last inequality follows from (7).

A.4 Proof of Theorem 1

Define  x1 = φ(W 0x + b0) ∈ Rnand  y1 = φ(W 0y + b0) ∈ Rnfor two arbitrary inputs  x, y ∈ Rn0. Using Lemma 1, we can write the quadratic inequality


where  T ∈ Tnand  Tnis defined as in (8). The preceding inequality can be simplified to


By left and right multiplying  M(ρ, T) in (10) by�(x−y)⊤ (x1−y1)⊤�and�(x−y)⊤ (x1−y1)⊤�⊤, respectively, and rearranging terms, we obtain


By adding both sides of the preceding inequalities, we obtain


or, equivalently,



Figure 7: Bounds on the image of an  ℓ2 0.01-ball under a neural network trained on the Iris dataset. The network was trained using the Adam optimizer and reached a test accuracy of 99%.

Finally, note that by definition of  x1and  y1, we have  f(x) = W 1x1+  b1and  f(y) = W 1y1+  b1. Therefore, the preceding inequality implies


A.5 Proof of Theorem 2

For two arbitrary inputs  x0, y0 ∈ Rn0, define  x = [x0⊤ · · · xℓ⊤]⊤and y = [y0⊤ · · · yℓ⊤]⊤Using the compact notation in (12), we can write


Multiply both sides of the first matrix in (14) by (x − y)⊤and (x − y), respectively and use the preceding identities to obtain


where the last inequality follows from Lemma 1. On the other hand, by multiplying both sides of the second matrix in (14) by (x − y)⊤and (x − y), respectively, we can write


where we have used the fact that  f(x) = W ℓxℓ+  bℓand  f(y) = W ℓyℓ+  bℓ. By adding both sides of (20) and (21), we get


When the LMI in (14) holds, the left-hand side of (22) is non-positive, implying that the right-hand side is non-positive.

A.6 Bounding the output set of a neural network classifier

As we have shown, accurate estimation of the Lipschitz constant can provide tighter bounds on adversarial examples drawn from a perturbation set around a nominal point ¯x. Figure 7 shows these bounds for a neural network trained on the Iris dataset. Because this dataset has three classes, we project the output sets onto the coordinate axes. This figure clearly shows that our bound is quite close to the naive lower bound. Indeed, the naive upper bound and the constant computed by CPLip are considerably looser.


Figure 8: Lipschitz bounds for networks in Table 1.

A.7 Further analysis of scalability

Figure 8 shows the Lipschitz bounds found for the networks used in Table 1.

A.8 LipSDP and repeated nonlinearities

The incremental quadratic constraint in Lemma 1 is valid for functions of the form  φ(x) = [ϕ1(x1) · · · ϕn(xn)]⊤as along as all  ϕi’s are slope-restricted on [α, β]. Therefore, this quadratic constraint does not exploit the fact that the same activation functions are used. In this subsection, we derive more refined quadratic constraints that exploit this additional structure when one of the points in the definition of Lipschitz continuity is fixed. Explicity, for a given  y ∈ Rn0, our goal is to compute the smallest  ρ(depending on y) such that


This inequality can be equivalently written as


To certify the quadratic inequality above, we make use of the following definition from [12].

Definition 3 Let  φ: Rn → Rnand suppose  Q ⊂ S2n+1is the set of all symmetric and indefinite matrices Q such that the inequality


holds for all  x ∈ X ⊆ Rn. Then we say  φsatisfies the quadratic constraint defined by Q on X.

Note that Quadratic Constraints as defined above are more general than Incremental Quadratic Constraints, and would allow us to capture more information about  φ. For example, for any pair (xi, xj) i ̸= j, since ϕ(xi) and  ϕ(xj) are obtained from the same nonlinearity, we can write the incremental quadratic constraint


This inequality can be equivalently represented as


This quadratic constraints couples neuron i and  j, i ̸= j, resulting in  O(n2) quadratic constraints for an activation layer with n neurons. In the following lemma, we characterize  φ(x) = [ϕ(x1) · · · ϕ(xn)]⊤with all possible incremental quadratic constraints.

Lemma 2 (Lemma 2 of [12]) Suppose  ϕ: R → Ris slope-restricted on [α, β]. Define the set


where  eiis the i-th standard basis vector in  Rn. Then for any  S ∈ Snthe vector-valued function  φ : Rn → Rndefined by  φ(x) = [ϕ(x1) · · · ϕ(xn)]⊤satisfies


In the above Lemma, we are only capturing the quadratic constraints implied by slope restriction that are imposed between any pair of neurons. However, depending on what activation function is being considered, more quadratic constraints can be included. For a detailed account of Quadratic Constraints, see [12].

LipSDP with repeated nonlinearities. Consider the compact representation (12) of a multi-layer network f defined by the recursion (2). Suppose  φsatisfies the quadratic constraint defined by Q on  Rn. Utilizing Definition 3, we can write the following quadratic constraints for any  Q ∈ Q,


Now consider the following LMI,


where  R(ρ) is defined in (23). By left- and right- multiplying both sides by�x⊤ 1�and�x⊤ 1�⊤, respectively, we obtain


The first term is non-negative by (29), implying that the second term is non-positive. Therefore,


Therefore, the feasibility of the LMI in (30) implies the inequality in (23). We remark that the LMI in (30) has now  O(n2) decision variables. For future work, we will invesitigate the effectiveness of these additional decision variables on the computed Lipschitz bound.

Designed for Accessibility and to further Open Science