My stuff
Lipschitz constant estimation of Neural Networks via sparse polynomial optimization

We introduce LiPopt, a polynomial optimization framework for computing increasingly tighter upper bounds on the Lipschitz constant of neural networks. The underlying optimization problems boil down to either linear (LP) or semidefinite (SDP) programming. We show how to use the sparse connectivity of a network, to significantly reduce the complexity of computation. This is specially useful for convolutional as well as pruned neural networks. We conduct experiments on networks with random weights as well as networks trained on MNIST, showing that in the particular case of the  ℓ∞-Lipschitz constant, our approach yields superior estimates, compared to baselines available in the literature.

We consider a neural network  fddefined by the recursion:


for an integer d larger than 1, matrices  {Wi}di=1of appropriate dimensions and an activation function σ, understood to be applied element-wise. We refer to d as the depth, and we focus on the case where fdhas a single real value as output.

In this work, we address the problem of estimating the Lipschitz constant of the network  fd. A function f is Lipschitz continuous with respect to a norm  ∥·∥if there exists a constant L such that for all x, y we have  |f(x) − f(y)| ≤ L∥x − y∥. The minimum over all such values satisfying this condition is called the Lipschitz constant of f and is denoted by L(f).

The Lipschitz constant of a neural network is of major importance in many successful applications of deep learning. In the context of supervised learning, Bartlett et al. (2017) show how it directly correlates with the generalization ability of neural network classifiers, suggesting it as model complexity measure. It also provides a measure of robustness against adversarial perturbations (Szegedy et al., 2014) and can be used to improve such metric (Cisse et al., 2017). Moreover, an upper bound on  L(fd)provides a certificate of robust classification around data points (Weng et al., 2018).

Another example is the discriminator network of the Wasserstein GAN (Arjovsky et al., 2017), whose Lipschitz constant is constrained to be at most 1. To handle this constraint, researchers have proposed different methods like heuristic penalties (Gulrajani et al., 2017), upper bounds (Miyato et al., 2018), choice of activation function (Anil et al., 2019), among many others. This line of work has shown that accurate estimation of such constant is key to generating high quality images.

Lower bounds or heuristic estimates of  L(fd)can be used to provide a general sense of how robust a network is, but fail to provide true certificates of robustness to input perturbations. Such certificates require true upper bounds, and are paramount when deploying safety-critical deep reinforcement learning applications (Berkenkamp et al., 2017; Jin & Lavaei, 2018). The trivial upper bound given by the product of layer-wise Lipschitz constants is easy to compute but rather loose and overly pessimistic, providing poor insight into the true robustness of a network (Huster et al., 2018).

Indeed, there is a growing need for methods that provide tighter upper bounds on  L(fd), even at the expense of increased complexity. For example Raghunathan et al. (2018a); Jin & Lavaei (2018); Fa- zlyab et al. (2019) derive upper bounds based on semidefinite programming (SDP). While expensive to compute, these type of certificates are in practice surprisingly tight. Our work belongs in this vein of research, and aims to overcome some limitations in the current state-of-the-art.

Our Contributions.

▷We present LiPopt, a general approach for upper bounding the Lipschitz constant of a neural network based on a relaxation to a polynomial optimization problem (POP) (Lasserre, 2015). This approach requires only that the unit ball be described with polynomial inequalities, which covers the common  ℓ2- and  ℓ∞-norms.

▷Based on a theorem due to Weisser et al. (2018), we exploit the sparse connectivity of neural network architectures to derive a sequence of linear programs (LPs) of considerably smaller size than their vanilla counterparts. We provide an asymptotic analysis of the size of such programs, in terms of the number of neurons, depth and sparsity of the network.

▷Focusing on the  ℓ∞-norm, we experiment on networks with random weights and networks trained on MNIST (Lecun et al., 1998). We evaluate different configurations of depth, width and sparsity and we show that the proposed sequence of LPs can provide tighter upper bounds on  L(fd)compared to other baselines available in the literature.

Notation. We denote by  nithe number of columns of the matrix  Wiin the definition (1) of the network. This corresponds to the size of the i-th layer, where we identify the input as the first layer. We let  n = n1 + . . . + ndbe the total number of neurons in the network. For a vector x, Diag(x) denotes the square matrix with x in its diagonal and zeros everywhere else. For an array X, vec(X) is the flattened array. The support of a sequence  supp(α)is defined as the set of indices j such that αjis nonzero. For  x = [x1, . . . , xn]and a sequence of nonnegative integers  γ = [γ1, . . . , γn]we denote by  xγthe monomial  xγ11 xγ22 . . . xγnn. The set of nonnegative integers is denoted by N.

Remark. The definition of network (1) covers typical architectures composed of dense and convolutional layers. In general, our proposed approach can be readily extended with minor modifications to any directed acyclic computation graph e.g., residual network architectures (He et al., 2016).

In this section we derive an upper bound on  L(fd)given by the value of a POP, i.e. the minimum value of a polynomial subject to polynomial inequalities. Our starting point is the following theorem, which casts L(f) as an optimization problem:

Theorem 1. Let f be a differentiable and Lipschitz continuous function on an open, convex subset X of an euclidean space. Let  ∥·∥∗be the dual norm. The Lipschitz constant of f is given by


For completeness, we provide a proof in appendix A. In our setting, we assume that the activation function  σis Lipschitz continuous and differentiable. In this case, the assumptions of Theorem 1 are fulfilled because  fdis a composition of activations and linear transformations. The differentiability assumption rules out the common ReLU activation  σ(x) = max{0, x}, but allows many others such as the exponential linear unit (ELU) (Clevert et al., 2015) or the softplus.

Using the chain rule, the compositional structure of  fdyields the following formula for its gradient:


For every  i = 1, . . . , d − 1we introduce a variable  si = σ′(fi(x))corresponding to the derivative of  σat the i-th hidden layer of the network. For activation functions like ELU or softplus, their derivative is bounded between 0 and 1, which implies that  0 ≤ si ≤ 1. This bound together with the definition of the dual norm  ∥x∥∗ := sup∥t∥≤1 tT ximplies the following upper bound of  L(fd):


We will refer to the polynomial objective of this problem as the norm-gradient polynomial of the network  fd, a central object of study in this work.

For some frequently used  ℓp-norms, the constraint  ∥t∥p ≤ 1can be written with polynomial inequalities. In the rest of this work, we use exclusively the  ℓ∞-norm for which  ∥t∥∞ ≤ 1is equivalent to the polynomial inequalities  −1 ≤ ti ≤ 1, for  i = 1, . . . , n1. However, note that when  p ≥ 2is a positive even integer,  ∥t∥p ≤ 1is equivalent to a single polynomial inequality  ∥t∥pp ≤ 1, and our proposed approach can be adapted with minimal modifications.

In such cases, the optimization problem in the right-hand side of (4) is a POP. Optimization of polynomials is a NP-hard problem and we do not expect to have efficient algorithms for solving (4) in this general form. In the next sections we describe LiPopt: a systematic way of obtaining an upper bound on  L(fd)via tractable approximation methods of the POP (4).

Local estimation. In many practical escenarios, we have additional bounds on the input of the network. For example, in the case of image classification tasks, valid input is constrained in a hypercube. In the robustness certification task, we are interested in all possible input in a  ϵ-ball around some data point. In those cases, it is interesting to compute a local Lipschitz constant, that is, the Lipschitz constant of a function restricted to a subset.

We can achieve this by deriving tighter bounds  0 ≤ li ≤ si ≤ ui ≤ 1, as a consequence of the restricted input (see for example, Algorithm 1 in Wong & Kolter (2018)). By incorporating this knowledge in the optimization problem (4) we obtain bounds on local Lipschitz constants of  fd. We study this setting and provide numerical experiments in section 7.3.

Choice of norm. We highlight the importance of computing good upper bounds on  L(fd)with respect to the  ℓ∞-norm. It is one of the most commonly used norms to assess robustness in the adversarial examples literature. Moreover, it has been shown that, in practice,  ℓ∞-norm robust networks are also robust in other more plausible measures of perceptibility, like the Wasserstein distance (Wong et al., 2019). This motivates our focus on this choice.


For ease of exposition, we rewrite (4) as a POP constrained in  [0, 1]nusing the substitution  s0 := (t+1)/2. Denote by p the norm-gradient polynomial, and let  x = [s0, . . . , sd−1]be the concatenation of all variables. Polynomial optimization methods (Lasserre, 2015) start from the observation that a value  λis an upper bound for p over a set K if and only if the polynomial  λ − pis positive over K.

In LiPopt, we will employ a well-known classical result in algebraic geometry, the so-called Krivine’s positivity certificate1, but in theory we can use any positivity certificate like sum-of-squares (SOS). The following is a straightforward adaptation of Krivine’s certificate to our setting:

Theorem 2. (Adapted from Krivine (1964); Stengle (1974); Handelman (1988)) If the polynomial λ − pis strictly positive on  [0, 1]n, then there exist finitely many positive weights  cαβsuch that


By truncating the degree of Krivine’s positivity certificate (Theorem 2) and minimizing over all possible upper bounds  λwe obtain a hierarchy of LP problems (Lasserre, 2015, Section 9):


where  N2nkis the set of nonnegative integer sequences of length 2n adding up to at most k. This is indeed a sequence of LPs as the polynomial equality constraint can be implemented by equating coefficients in the canonical monomial basis. For this polynomial equality to be feasible, the degree of the certificate has to be at least that of the norm-gradient polynomial p, which is equal to the depth d. This implies that the first nontrivial bound (θk < ∞) corresponds to k = d.

The sequence  {θk}∞k=1is non-incresing and converges to the maximum of the upper bound (4). Note that for any level of the hierarchy, the solution of the LP (6) provides a valid upper bound on  L(fd).

An advantage of using Krivine’s positivity certificate over SOS is that one obtains an LP hierarchy (rather than SDP), for which commercial solvers can reliably handle a large instances. Other positivity certificates offering a similar advantage are the DSOS and SDSOS hierarchies (Ahmadi & Majumdar, 2019), which boil down to LP or second order cone programming (SOCP), respectively.

Drawback. The size of the LPs given by Krivine’s positivity certificate can become quite large. The dimension of the variable c is  |N2nk | = O(nk). For reference, if we consider the MNIST dataset and a one-hidden-layer network with 100 neurons we have��N2n2 �� ≈ 1.5 × 106while��N2n3 �� ≈ 9.3 × 108. To make this approach more scalable, in the next section we exploit the sparsity of the polynomial p to find LPs of drastically smaller size than (6), but with similar approximation properties.

Remark. In order to compute upper bounds for local Lipschitz constants, first obtain tighter bounds 0 ≤ li ≤ si ≤ uiand then perform the change of variables  �si = (si − li)/(ui − li)to rewrite the problem (4) as a POP constrained on  [0, 1]n.

Many neural network architectures, like those composed of convolutional layers, have a highly sparse connectivity between neurons. Moreover, it has been empirically observed that up to 90% of network weights can be pruned (set to zero) without harming accuracy (Frankle & Carbin, 2019). In such cases their norm-gradient polynomial has a special structure that allows polynomial positivity certificates of smaller size than the one given by Krivine’s positivity certificate (Theorem 2).

In this section, we describe an implementation of LiPopt (Algorithm 1) that exploits the sparsity of the network to decrease the complexity of the LPs (6) given by the Krivine’s positivity certificate. In this way, we obtain upper bounds on  L(fd)that require less computation and memory. Let us start with the definition of a valid sparsity pattern:

Definition 1. Let I = {1, . . . , n} and p be a polynomial with variable  x ∈ Rn. A valid sparsity pattern of p is a sequence  {Ii}mi=1of subsets of I, called cliques, such that  �mi=1 Ii = Iand:


When the polynomial objective p in a POP has a valid sparsity pattern, there is an extension of Theorem 2 due to Weisser et al. (2018), providing a smaller positivity certificate for  λ − pover [0, 1]n. We refer to it as the sparse Krivine’s certificate and we include it here for completeness:

Theorem 3 (Adapted from Weisser et al. (2018)). Let a polynomial p have a valid sparsity pattern {Ii}mi=1. Define  Nias the set of sequences  (α, β) ∈ N2nwhere the support of both  αand  βis contained in  Ii. If  λ − pis strictly positive over  K = [0, 1]n, there exist finitely many positive weights  cαβsuch that


where the polynomials  hαβare defined as in (5).

The sparse Krivine’s certificate can be used like the general version (Theorem 2) to derive a sequence of LPs approximating the upper bound on  L(fd)stated in (4). However, the number of different polynomials  hαβof degree at most k appearing in the sparse certificate can be drastically smaller, the amount of which determines how good the sparsity pattern is.

We introduce a graph that depends on the network  fd, from which we will extract a sparsity pattern for the norm-gradient polynomial of a network.


Figure 1: Sparsity pattern of Proposition 1 for a network of depth three.

Figure 2: Structure of one set in the sparsity pattern from Proposition 1 for a network with 2D convolutional layers with  3 × 3filters.

Definition 2. Let  fdbe a network with weights  {Wi}di=1. Define a directed graph  Gd = (V, E)as:


which we call the computational graph of the network  fd.

In the graph  Gdthe vertex  s(i,j)represents the j-th neuron in the i-th layer. There is a directed edge between two neurons in consecutive layers if they are joined by a nonzero weight in the network. The following result shows that for fully connected networks we can extract a valid sparsity pattern from this graph. We relegate the proof to appendix B.

Proposition 1. Let  fdbe a dense network (all weights are nonzero). The following sets, indexed by i = 1, . . . , nd, form a valid sparsity pattern for the norm-gradient polynomial of the network  fd:


We refer to this as the sparsity pattern induced by  Gd. An example is depicted in in Figure 1.

Remark. When the network is not dense, the the second condition (Definition 1) for the sparsity pattern (9) to be valid might not hold. In that case we lose the guarantee that the values of the corresponding LPs converge to the maximum of the POP (4). Nevertheless, it still provides a valid positivity certificate that we use to upper bound  L(fd). In Section 7 we show that in practice it provides upper bounds of good enough quality. If needed, a valid sparsity pattern can be obtained via a chordal completion of the correlative sparsity graph of the POP (Waki et al., 2006).

We now quantify how good this sparsity pattern is. Let s be the size of the largest clique in a sparsity pattern, and let  Ni,kbe the subset of  Ni(defined in Theorem 3) composed of sequences summing up to k. The number of different polynomials for the k-th LP in the hierarchy given by the sparse Krivine’s certificate can be bounded as follows:


We immediately see that the dependence on the number of cliques m is really mild (linear) but the size of the cliques as well as the degree of the hierarchy can really impact the size of the optimization problem. Nevertheless, this upper bound can be quite loose; polynomials  hαβthat depend only on variables in the intersection of two or more cliques are counted more than once.

The number of cliques given in the sparsity pattern induced by  Gdis equal to the size of the last layer  m = ndand the size of each clique depends on the particular implementation of the network. We now study different architectures that could arise in practice, and determine the amount of polynomials in their sparse Krivine’s certificate.

Fully connected networks. Even in the case of a network with all nonzero connections, the sparsity pattern induced by  Gddecreases the size of the LPs when compared to Krivine’s certificate. In this case the cliques have size  n1 +. . .+nd−1 +1but they all have the same common intersection equal to all neurons up to the second-to-last hidden layer. A straightforward counting argument shows that the total number of polynomials is  O(n(n1 +. . .+nd−1 +1)k−1), improving the upper bound (10). Unstructured sparsity. In the case of networks obtained by pruning (Hanson & Pratt, 1989) or generated randomly from a distribution over graphs (Xie et al., 2019), the sparsity pattern can be arbitrary. In this case the size of the resulting LPs varies at runtime. Under the layer-wise assumption that any neuron is connected to at most r neurons in the previous layer, the size of the cliques in (9) is bounded as  s = O(rd). This estimate has an exponential dependency on the depth but ignores that many neurons might share connections to the same inputs in the previous layer, thus being potentially loose. The bound (10) implies that the number of different polynomials is  O(ndrdk).

2D Convolutional networks. The sparsity in the weight matrices of convolutional layers has a certain local structure; neurons are connected to contiguous inputs in the previous layer. Adjacent neurons also have many input pixels in common (see Figure 2). Assuming a constant number of channels per layer, the size of the cliques in (9) is  O(d3). Intuitively, such number is proportional to the volume of the pyramid depicted in Figure 2 where each dimension depends linearly on d. Using (10) we get that there are  O(ndd3k)different polynomials in the sparse Krivine’s certificate. This is a drastic decrease in complexity when compared to the unstructured sparsity case.

The use of sparsity in polynomial optimization preceeds Theorem 3 (Weisser et al., 2018). First studied in the context of sum-of-squares by Kojima et al. (2005) and further refined in Waki et al. (2006); Lasserre (2006) (and references therein), it has found applications in safety verification (Yang et al., 2016; Zhang et al., 2018), sensor localization Wang et al. (2006), optimal power flow (Ghaddar et al., 2015) and many others. Our work fits precisely into this set of important applications.


Another way of upper bounding  L(fd)comes from a further relaxation of (4) to an SDP. We consider the following equivalent formulation where the variables  siare normalized to lie in the interval [−1, 1], and we rename  t = s0:


Any polynomial optimization problem like (11) can be cast as a (possibly non-convex) quadratically constrained quadratic program (QCQP) by introducing new variables and quadratic constraints. This is a well-known procedure described in Park & Boyd (2017, Section 2.1). When d = 2 problem (11) is already a QCQP (for the  ℓ∞and  ℓ2-norm cases) and no modification is necessary.

QCQP reformulation. We illustrate the case d = 3 where we have the variables  s1, s2corresponding to the first and second hidden layer and a variable  s0corresponding to the input. The norm-gradient polynomial in this case is cubic, and it can be rewritten as a quadratic polynomial by introducing new variables corresponding to the product of the first and second hidden layer variables.

More precisely the introduction of a variable  s1,2with quadratic constraint  s1,2 = vec(s1sT2 )allows us to write the objective (11) as a quadratic polynomial. The problem then becomes a QCQP with variable  y = [1, s0, s1, s2, s1,2]of dimension  1 + n + n1n2.

SDP relaxation. Any quadratic objective and constraints can then be relaxed to linear constraints on the positive semidefinite variable  yyT = X ≽ 0yielding the so-called Shor’s relaxation of (11) (Park & Boyd, 2017, Section 3.3). When d = 2 the resulting SDP corresponds precisely to the one studied in Raghunathan et al. (2018a). This resolves a common misconception (Raghunathan et al., 2018b) that this approach is only limited to networks with one hidden layer.

Note that in our setting we are only interested in the optimal value rather than the optimizers, so there is no need to extract a solution for (11) from that of the SDP relaxation.

Drawback. This approach includes a further relaxation step from (11), thus being fundamentally limited in how tightly it can upper bound the value of  L(fd). Moreover when compared to LP solvers, off-the-shelf semidefinite programming solvers are, in general, much more limited in the number of variables they can efficiently handle.

In the case d = 2 this relaxation provides a constant factor approximation to the original QCQP (Ye, 1999). Further approximation quality results for such hierarchical optimization approaches to NP-hard problems are out of the scope of this work.

Relation to sum-of-squares. The QCQP approach might appear fundamentaly different to the hierarchical optimization approaches to POPs, like the one described in Section 3. However, it is known that Shor’s SDP relaxation corresponds exactly to the first degree of the SOS hierarchical SDP solution to the QCQP relaxation (Lasserre, 2000). Thus, the approach in section 3 and the one in this section are, in essence, the same; they only differ in the choice of polynomial positivity certificate.

Estimation of  L(fd)with  ℓ2-norm is studied by Virmaux & Scaman (2018); Combettes & Pesquet (2019); Fazlyab et al. (2019); Jin & Lavaei (2018). The method SeqLip proposed in Virmaux & Scaman (2018) has the drawback of not providing true upper bounds. It is in fact a heuristic method for solving (4) but which provides no guarantees and thus can not be used for robustness certification. In contrast the LipSDP method of Fazlyab et al. (2019) provides true upper bounds on  L(fd)and in practice shows superior performance over both SeqLip and CPLip (Combettes & Pesquet, 2019).

Despite the accurate estimation of LipSDP, its formulation is limited to the  ℓ2-norm. The only estimate available for other  ℓp-norms comes from the equivalence of norms in euclidean spaces. For instance, we can obtain an upper bound for the  ℓ∞-norm after multiplying the  ℓ2Lipschitz constant upper bound by the square root of the input dimension. The resulting bound can be rather loose and our experiments in section 7 confirm the issue. In contrast, our proposed approach LiPopt can acommodate any norm whose unit ball can be described via polynomial inequalities.

Let us point to one key advantage of LiPopt, compared to LipSDP (Jin & Lavaei, 2018; Fazlyab et al., 2019). In the context of robustness certification we are given a sample  x♮and a ball of radius ϵaround it. Computing an upper bound on the local Lipschitz constant in this subset, rather than a global one, can provide a larger region of certified robustness. Taking into account the restricted domain we can refine the bounds in our POP (see remark in section 1). This potentially yields a tighter estimate of the local Lipschitz constant. On the other hand, it is not clear how to include such additional information in LipSDP, which only computes one global bound on the Lipschitz constant for the unconstrained network.

Raghunathan et al. (2018a) find an upper bound for  L(fd)with  ℓ∞metric starting from problem (4) but only in the context of one-hidden-layer networks (d = 2). To compute such bound they use its corresponding Shor’s relaxation and obtain as a byproduct a differentiable regularizer for training networks. They claim such approach is limited to the setting d = 2 but, as we remark in section 5, it is just a particular instance of the SDP relaxation method for QCQPs arising from a polynomial optimization problem. We find that this method fits into the LiPopt framework, using SOS certificates instead of Krivine’s. We expect that the SDP-based bounds described in 5 can also be used as regularizers promoting robustness.

Weng et al. (2018) provide an upper bound on the local Lipschitz constant for networks based on a sequence of ad-hoc bounding arguments, which are particular to the choice of ReLU activation function. In contrast, our approach applies in general to activations whose derivative is bounded.

We consider the following estimators of  L(fd)with respect to the  ℓ∞norm:



We compare the bounds obtained by the algorithms described above on networks with random weights and either one or two hidden layers. We define the sparsity level of a network as the maximum number of neurons any neuron in one layer is connected to in the next layer. For example, the network represented on Figure 1 has sparsity 2. The non-zero weights of network’s i-th layer are sampled uniformly in  [− 1√ni , 1√ni ]where  niis the number of neurons in layer i.

For different configurations of width and sparsity, we generate 10 random networks and average the obtained Lipschitz bounds. For better comparison, we plot the relative error. Since we do not know the true Lipschitz constant, we cannot compute the true relative error. Instead, we take as reference the lower bound given by LBS. Figures 3 and 5 show the relative error, i.e.,  (ˆL − LLBS)/LLBSwhere  LLBSis the lower bound computed by LBS and ˆLis the estimated upper bound. Figures 9 and 10 in Appendix C we show the values of the computed Lipschitz bounds for 1 and 2 hidden layers respectively.

When the chosen degree for LiPopt-k is the smallest as possible, i.e., equal to the depth of the network, we observe that the method is already competitive with the SDP method, especially in the case of 2 hidden layers. When we increment the degree by 1, LiPopt-k becomes uniformly better than SDP over all tested configurations. We remark that the upper bounds given by UBP are too large to be shown in the plots. Similarly, for the 1-hidden layer networks, the bounds from LipSDP are too large to be plotted.

Finally, we measured the computation time of the different methods on each tested network (Figures 4 and 6). We observe that the computation time for LiPopt-k heavily depends on the network sparsity, which reflects the fact that such structure is exploited in the algorithm. In contrast, the time required for SDP does not depend on the sparsity, but only on the size of the network. Therefore as the network size grows (with fixed sparsity level), LipOpt-k obtains a better upper bound and runs faster. Also, with our method, we see that it is possible to increase the computation power in order to compute tighter bounds when required, making it more flexible than SDP in terms of computation/accuracy tradeoff. LiPopt uses the Gurobi LP solver, while SDP uses Mosek. All methods run on a single machine with Core i7 2.8Ghz quad-core processor and 16Gb of RAM.


Figure 3: Lipschitz approximated relative error for 1-hidden layer networks


Figure 4: Computation times for 1-hidden layer networks (seconds


Figure 5: Lipschitz approximated relative error for 2-hidden layer networks


Figure 6: Computation times for 2-hidden layer networks (seconds


Similarly, we compare these methods on networks trained on MNIST. The architecture we use is a fully connected network with two hidden layers with 300 and 100 neurons respectively, and with one-hot output of size 10. Since the output is multi-dimensional, we restrict the network to a single output, and estimate the Lipschitz constant with respect to label 8.

Moreover, in order to improve the scalability of our method, we train the network using the pruning strategy described in Han et al. (2015)2. After training the full network using a standard technique, the weights of smallest magnitude are set to zero. Then, the network is trained for additional iterations, only updating the nonzero parameters. Doing so, we were able to remove 95% of the weights, while preserving the same test accuracy. We recorded the Lipschitz bounds for various methods in Table 7.2. We observe clear improvement of the Lipschitz bound obtained from LiPopt-k compared to SDP method, even when using k = 3. Also note that the input dimension is too large for the method Lip-SDP to provide competitive bound, so we do not provide the obtained bound for this method.



In the of section 7.1, we study the improvement on the upper bound obtained by LiPopt, when we incorporate tighter upper and lower bounds on the variables  siof the polynomial optimization problem (4). Such bounds arise from the limited range that the pre-activation values of the network can take, when the input is limited to an  ℓ∞-norm ball of radius  ϵcentered at an arbitrary point  x0. The algorithm that computes upper and lower bounds on the pre-activation values is fast (it has the same complexity as a forward pass) and is described, for example, in Wong & Kolter (2018). The variables  sicorrespond to the value of the derivative of the activation function. For activations like ELU or ReLU, their derivative is monotonically increasing, so we need only evaluate it at the upper and lower bounds of the pre-activation values to obtain corresponding bounds for the variables  si.

We plot the local upper bounds obtained by LiPopt-3 for increasing values of the radius  ϵ, the bound for the global constant (given by LiPopt-3) and the lower bound on the local Lipschitz constant obtained by sampling in the  ϵ-neighborhood (LBS). We sample 15 random networks and plot the average values obtained. We observe clear gap between both estimates, which shows that larger certified balls could be obtained with such method in the robustness certification applications.


Figure 7: Global vs local Lipschitz constant bounds for 1-hidden layer networks


Figure 8: Global vs local Lipschitz constant bounds for 2-hidden layer networks

In this work, we have introduced a general approach for computing an upper bound on the Lipschitz constant of neural networks. This approach is based on polynomial positivity certificates and generalizes some existing methods available in the literature. We have empirically demonstrated that it can tightly upper bound such constant. The resulting optimization problems are computationally expensive but the sparsity of the network can reduce this burden.

In order to further scale such methods to larger and deeper networks, we are interested in several possible directions: (i) divide-and-conquer approaches splitting the computation on sub-networks in the same spirit of Fazlyab et al. (2019), (ii) exploiting parallel optimization algorithms leveraging the structure of the polynomials, (iii) custom optimization algorithms with low-memory costs such as Frank-wolfe-type methods for SDP (Yurtsever et al., 2019) as well as stochastic handling of constraints (Fercoq et al., 2019) and (iv), exploting the symmetries in the polynomial that arise from weight sharing in typical network architectures to further reduce the size of the problems.

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 725594 - time-data) and from the Swiss National Science Foundation (SNSF) under grant number 200021 178865. FL is supported through a PhD fellowship of the Swiss Data Science Center, a joint venture between EPFL and ETH Zurich. VC acknowledges the 2019 Google Faculty Research Award.

A. Ahmadi and A. Majumdar. Dsos and sdsos optimization: More tractable alternatives to sum of squares and semidefinite optimization. SIAM Journal on Applied Algebra and Geometry, 3(2):193–230, 2019. doi: 10.1137/18M118935X. URL https://doi.org/10.1137/ 18M118935X.

Cem Anil, James Lucas, and Roger Grosse. Sorting out Lipschitz function approximation. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 291–301, Long Beach, California, USA, 09–15 Jun 2019. PMLR.

Martin Arjovsky, Soumith Chintala, and L´eon Bottou. Wasserstein generative adversarial networks. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 214–223, International Convention Centre, Sydney, Australia, 06–11 Aug 2017. PMLR.

Peter L Bartlett, Dylan J Foster, and Matus J Telgarsky. Spectrally-normalized margin bounds for neural networks. In Advances in Neural Information Processing Systems 30, pp. 6240–6249. Curran Associates, Inc., 2017.

Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause. Safe model-based reinforcement learning with stability guarantees. In Advances in Neural Information Processing Systems 30, pp. 908–918. Curran Associates, Inc., 2017.

Moustapha Cisse, Piotr Bojanowski, Edouard Grave, Yann Dauphin, and Nicolas Usunier. Parseval networks: Improving robustness to adversarial examples. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 854–863, International Convention Centre, Sydney, Australia, 06–11 Aug 2017. PMLR.

Djork-Arn´e Clevert, Thomas Unterthiner, and Sepp Hochreiter. Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs). arXiv e-prints, art. arXiv:1511.07289, Nov 2015.

Patrick L. Combettes and Jean-Christophe Pesquet. Lipschitz Certificates for Neural Network Struc- tures Driven by Averaged Activation Operators. arXiv e-prints, art. arXiv:1903.01014, Mar 2019.

Mahyar Fazlyab, Alexander Robey, Hamed Hassani, Manfred Morari, and George J. Pappas. Effi- cient and Accurate Estimation of Lipschitz Constants for Deep Neural Networks. arXiv e-prints, art. arXiv:1906.04893, Jun 2019.

Olivier Fercoq, Ahmet Alacaoglu, Ion Necoara, and Volkan Cevher. Almost surely constrained con- vex optimization. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 1910–1919, Long Beach, California, USA, 09–15 Jun 2019. PMLR.

Jonathan Frankle and Michael Carbin. The lottery ticket hypothesis: Finding sparse, trainable neural networks. In International Conference on Learning Representations, 2019.

Bissan Ghaddar, Jakub Marecek, and Martin Mevissen. Optimal power flow as a polynomial opti- mization problem. IEEE Transactions on Power Systems, 31(1):539–546, 2015.

Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville. Im- proved training of wasserstein gans. In Advances in Neural Information Processing Systems 30, pp. 5767–5777. Curran Associates, Inc., 2017.

Song Han, Huizi Mao, and William J Dally. Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding. arXiv preprint arXiv:1510.00149, 2015.

David Handelman. Representing polynomials by positive linear functions on compact convex polyhedra. Pacific J. Math., 132(1):35–62, 1988. URL https://projecteuclid.org: 443/euclid.pjm/1102689794.

Stephen Jose Hanson and Lorien Y. Pratt. Comparing biases for minimal network construction with back-propagation. In Advances in Neural Information Processing Systems 1, pp. 177–185. Morgan-Kaufmann, 1989.

Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Conference on Computer Vision and Pattern Recognition, pp. 770–778, 06 2016. doi: 10.1109/CVPR.2016.90.

Todd Huster, Cho-Yu Jason Chiang, and Ritu Chadha. Limitations of the Lipschitz constant as a defense against adversarial examples. arXiv e-prints, art. arXiv:1807.09705, Jul 2018.

Ming Jin and Javad Lavaei. Stability-certified reinforcement learning: A control-theoretic perspec- tive. arXiv e-prints, art. arXiv:1810.11505, Oct 2018.

Masakazu Kojima, Sunyoung Kim, and Hayato Waki. Sparsity in sums of squares of polynomials. Mathematical Programming, 103(1):45–62, May 2005. ISSN 1436-4646. doi: 10.1007/s10107-004-0554-3.

Jean-Louis Krivine. Anneaux pr´eordonn´es. Journal d’analyse math´ematique, 12:p. 307–326, 1964.

J. B. Lasserre. Convergent lmi relaxations for nonconvex quadratic programs. In Proceedings of the 39th IEEE Conference on Decision and Control (Cat. No.00CH37187), volume 5, pp. 5041–5046 vol.5, Dec 2000. doi: 10.1109/CDC.2001.914738.

Jean B Lasserre. Convergent sdp-relaxations in polynomial optimization with sparsity. SIAM Journal on Optimization, 17(3):822–843, 2006.

Jean Bernard Lasserre. An Introduction to Polynomial and Semi-Algebraic Optimization. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2015. doi: 10.1017/ CBO9781107447226.

Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied to document recog- nition. Proceedings of the IEEE, 86(11):2278–2324, Nov 1998. doi: 10.1109/5.726791.

Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral normalization for generative adversarial networks. In International Conference on Learning Representations, 2018.

Jaehyun Park and Stephen Boyd. General heuristics for nonconvex quadratically constrained quadratic programming. arXiv preprint arXiv:1703.07870, 2017.

Aditi Raghunathan, Jacob Steinhardt, and Percy Liang. Certified defenses against adversarial exam- ples. In International Conference on Learning Representations, 2018a.

Aditi Raghunathan, Jacob Steinhardt, and Percy S Liang. Semidefinite relaxations for certifying robustness to adversarial examples. In Advances in Neural Information Processing Systems, pp. 10877–10887, 2018b.

Gilbert Stengle. A nullstellensatz and a positivstellensatz in semialgebraic geometry. Mathematische Annalen, 207(2):87–97, Jun 1974. ISSN 1432-1807. doi: 10.1007/BF01362149.

Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian J. Goodfel- low, and Rob Fergus. Intriguing properties of neural networks. In 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.

Aladin Virmaux and Kevin Scaman. Lipschitz regularity of deep neural networks: analysis and efficient estimation. In Advances in Neural Information Processing Systems 31, pp. 3835–3844. Curran Associates, Inc., 2018.

H. Waki, S. Kim, M. Kojima, and M. Muramatsu. Sums of squares and semidefinite program relaxations for polynomial optimization problems with structured sparsity. SIAM Journal on Optimization, 17(1):218–242, 2006. doi: 10.1137/050623802.

Zizhuo Wang, Song Zheng, Stephen Boyd, and Yinyu Ye. Further relaxations of the sdp approach to sensor network localization. Tech. Rep., 2006.

Tillmann Weisser, Jean B Lasserre, and Kim-Chuan Toh. Sparse-bsos: a bounded degree sos hierar- chy for large scale polynomial optimization with sparsity. Mathematical Programming Computation, 10(1):1–32, 2018.

Lily Weng, Huan Zhang, Hongge Chen, Zhao Song, Cho-Jui Hsieh, Luca Daniel, Duane Boning, and Inderjit Dhillon. Towards fast computation of certified robustness for ReLU networks. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 5276–5285, Stockholmsmssan, Stockholm Sweden, 10–15 Jul 2018. PMLR.

Eric Wong and Zico Kolter. Provable defenses against adversarial examples via the convex outer adversarial polytope. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 5286–5295, Stockholmsmssan, Stockholm Sweden, 10–15 Jul 2018. PMLR.

Eric Wong, Frank Schmidt, and Zico Kolter. Wasserstein adversarial examples via projected Sinkhorn iterations. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 6808–6817, Long Beach, California, USA, 09–15 Jun 2019. PMLR.

Saining Xie, Alexander Kirillov, Ross Girshick, and Kaiming He. Exploring Randomly Wired Neural Networks for Image Recognition. International Conference on Computer Vision, Apr 2019.

Zhengfeng Yang, Chao Huang, Xin Chen, Wang Lin, and Zhiming Liu. A linear programming relaxation based approach for generating barrier certificates of hybrid systems. In John Fitzgerald, Constance Heitmeyer, Stefania Gnesi, and Anna Philippou (eds.), FM 2016: Formal Methods, pp. 721–738, Cham, 2016. Springer International Publishing. ISBN 978-3-319-48989-6.

Yinyu Ye. Approximating quadratic programming with bound and quadratic constraints. Mathematical Programming, 84(2):219–226, Feb 1999. ISSN 1436-4646. doi: 10.1007/s10107980012a.

Alp Yurtsever, Olivier Fercoq, and Volkan Cevher. A conditional-gradient-based augmented la- grangian framework. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 7272–7281, Long Beach, California, USA, 09–15 Jun 2019. PMLR.

Y. Zhang, Z. Yang, W. Lin, H. Zhu, X. Chen, and X. Li. Safety verification of nonlinear hybrid systems based on bilinear programming. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 37(11):2768–2778, Nov 2018. doi: 10.1109/TCAD.2018.2858383.


Theorem. Let f be a differentiable and Lipschitz continuous function on an open, convex subset X

of a euclidean space. Let  ∥ · ∥be the dual norm. The Lipschitz constant of f is given by


Proof. First we show that  L(f) ≤ supx∈X ∥∇f(x)∥∗.


were we have used the convexity of X.

Now we show the reverse inequality  L(f) ≥ supx∈X ∥∇f(x)∥∗. To this end, we show that for any

positive  ϵ, we have that  L(f) ≥ supx∈X ∥∇f(x)∥∗ − ϵ.

Let  z ∈ Xbe such that  ∥∇f(z)∥∗ ≥ supx∈X ∥∇f(x)∥∗ − ϵ. Because X is open, there exists a

sequence  {hk}∞k=1with the following properties:

1. hk, ∇f(z)⟩ = ∥hk∥∥∇f(z)∥∗

2. z  + hk ∈ X

3. limk→∞ hk= 0.

By definition of the gradient, there exists a function  δsuch that  limh→0 δ(h) = 0and the following



For our previously defined iterates  hkwe then have


Dividing both sides by  ∥hk∥and using the definition of L(f) we finally get



Proposition. Let  fdbe a dense network (all weights are nonzero). The following sets, indexed by

i = 1, . . . , nd, form a valid sparsity pattern for the norm-gradient polynomial of the network  fd:


Proof. First we show that  ∪mi=1Ii = I. This comes from the fact that any neuron in the network is

connected to at least one neuron in the last layer. Otherwise such neuron could be removed from the

network altogether.

Now we show the second property of a valid sparsity pattern. Note that the norm-gradient polyno-

mial is composed of monomials corresponding to the product of variables in a path from input to a

final neuron. This imples that if we let  pibe the sum of all the terms that involve the neuron  s(d−1,i)

we have that  p = �i pi, and  pionly depends on the variables in  Ii.

We now show the last property of the valid sparsity pattern. This is the only part where we use that

the network is dense. For any network architecture the first two conditions hold. We will use the fact

that the maximal cliques of a chordal graph form a valid sparsity pattern (see for example Lasserre


Because the network is dense, we see that the clique  Iiis composed of the neuron in the last layer

s(d−1,i)and all neurons in the previous layers. Now consider the extension of the computational

graph ˆGd = (V, ˆE)where


which consists of adding all the edges between the neurons that are not in the last layer. We show

that this graph is chordal. Let  (a1, . . . , ar, a1)be a cycle of length at least 4 (r ≥ 4). notice that

because neurons in the last layer are not connected between them in ˆG, no two consecutive neurons

in this cycle belong to the last layer. This implies that in the subsequence  (a1, a2, a3, a4, a5)at most

three belong to the last layer. A simple analysis of all cases implies that it contains at least two

nonconsecutive neurons not in the last layer. Neurons not in the last layer are always connected in

ˆG. This constitutes a chord. This shows that  ˆGdis a chordal graph. Its maximal cliques correspond

exactly to the sets in proposition.











Figure 9: Lipschitz bound comparison for 1-hidden layer networks








Figure 10: Lipschitz bound comparison for 2-hidden layer networks

Designed for Accessibility and to further Open Science