Title: CURVATURE-BASED ROBUSTNESS CERTIFICATES AGAINST ADVERSARIAL EXAMPLES

Abstract: A robustness certificate against adversarial examples is the minimum distance of a given input to the decision boundary of the classifier (or its lower bound). For any perturbation of the input with a magnitude smaller than the certificate value, the classification output will provably remain unchanged. Computing exact robustness certificates for deep classifiers is difficult in general since it requires solving a nonconvex optimization. In this paper, we provide computationally-efficient robustness certificates for deep classifiers with differentiable activation functions in two steps. First, we show that if the eigenvalues of the Hessian of the network (curvatures of the network) are bounded, we can compute a robustness certificate in the l 2 norm efficiently using convex optimization. Second, we derive a computationallyefficient differentiable upper bound on the curvature of a deep network. We also use the curvature bound as a regularization term during the training of the network to boost its certified robustness against adversarial examples. Putting these results together leads to our proposed Curvature-based Robustness Certificate (CRC) and Curvature-based Robust Training (CRT). Our numerical results show that CRC outperforms CROWN's certificate when trained with our regularizer while CRT leads to higher certified accuracy compared to standard adversarial training.

1. INTRODUCTION:

Modern neural networks achieve high accuracy on tasks such as image classification and speech recognition, but are known to be brittle to small, adversarially chosen perturbations of their inputs (Szegedy et al., 2014) . A classifier which correctly classifies an image x, can be fooled by an adversary to misclassify an adversarial example x + δ, such that x + δ is indistinguishable from x to a human. Adversarial examples can also fool systems when they are printed out on a paper and photographed with a smart phone (Kurakin et al., 2016a) . Even in a black box threat model, where the adversary has no access to the model parameters, attackers could target autonomous vehicles by using stickers or paint to create an adversarial stop sign that the vehicle would interpret as a yield or another sign (Papernot et al., 2016) . This trend is worrisome and suggests that these vulnerabilities need to be appropriately addressed before neural networks can be deployed in security critical applications.

In the last couple of years, several empirical defenses have been proposed for training classifiers to be robust against adversarial perturbations (Madry et al., 2018; Samangouei et al., 2018; Zhang et al., 2019; Papernot et al., 2016; Kurakin et al., 2016b; Miyato et al., 2017; Zheng et al., 2016) Although these defenses robustify classifiers to particular types of attacks, they can be still vulnerable against stronger attacks (Athalye et al., 2018; Carlini & Wagner, 2017; Uesato et al., 2018; Athalye & Carlini, 2018) . For example, (Athalye et al., 2018) showed most of the empirical defenses proposed in ICLR 2018 can be broken by developing tailored attacks for each of them.

To end the cycle between defenses and attacks, a line of work on certified defenses has gained attention where the goal is to train classifiers whose predictions are provably robust within some given region (Huang et al., 2016; Katz et al., 2017; Ehlers, 2017; Carlini et al., 2017; Cheng et al., 2017; Lomuscio & Maganti, 2017; Dutta et al., 2018; Fischetti & Jo, 2018; Bunel et al., 2017; Wang et al., 2018a; Wong & Kolter, 2017; Wang et al., 2018b; Wong et al., 2018; Raghunathan et al., 2018b; a; Dvijotham et al., 2018a; b; Croce et al., 2018; Singh et al., 2018; Gowal et al., 2018; Gehr et al., 2018; Mirman et al., 2018; Zhang et al., 2018b; Weng et al., 2018) . These methods, however, do not scale to large and practical networks used in solving modern machine learning problems. Another line of defense work focuses on randomized smoothing where the prediction is robust within some region around the input with a user-chosen probability (Liu et al., 2017; Cao & Gong, 2017; Lécuyer et al., 2018; Li et al., 2018; Cohen et al., 2019; Salman et al., 2019) . Although these methods can scale to large networks, certifying robustness with probability close to 1 often requires generating a large number of noisy samples around the input which leads to high test-time computational complexity.

If the classifier f (.) was linear, the distance of an input point x to its decision boundary (i.e. the robustness certificate) can be computed efficiently using a convex optimization. For example, the l 2 robustness certificate in that case would be equal to f (x) ∇ x f (x) . However, modern classifiers based on neural networks are not linear and and can have non-zero curvatures in different parts of the input domain. The deviation of the classifier from the linear model makes the robustness certification problem to be a non-convex optimization which is often difficult to solve exactly. However, if we could compute global bounds on the maximum curvature values of the classification network, one may be able to compute computationally-efficient lower bounds on the robustness certificate even for non-linear deep classifiers. This is the key intuition of our results in this paper.

In this work, we derive a global bound on the Lipschitz constant of the gradient of deep neural networks with differentiable activation functions (such as sigmoid, tanh, softplus, etc.) . This provides an upper bound on the magnitude of the eigenvalues of the Hessian or the curvature values of the classification network. Using this global curvature bound and for the l 2 metric, we tackle both the certification and attack problems. In the certification problem and for a given pre-trained classifier, we provide a computationally-efficient lower bound on the distance of a point to the classification decision boundary. In the related attack problem, for a given input and a region around it, our goal is to find a perturbed input (an adversarial example) that maximizes the loss inside the given region. The outcome of the attack problem is then used in the adversarial training procedure (Madry et al., 2018) to further robustify the network. Furthermore, our global curvature bound is differentiable and we show that adding it to the loss function as a regularizer boosts certified robustness measures.

We note that other recent works (e.g. Moosavi Dezfooli et al. (2019) ; Qin et al. (2019) ) empirically show that using an estimate of curvature at inputs as a regularizer leads to empirical robustness on par with the adversarial training. In this work, however, we use a provable global upper bound on the curvature (and not an estimate) as a regularizer and show that it results in high certified robustness. Moreover, previous works have tried to certify robustness by bounding the Lipschitz constant of the neural network (Szegedy et al., 2014; Peck et al., 2017; Zhang et al., 2018c; Anil et al., 2018; Hein & Andriushchenko, 2017) . Our approach, however, is based on bounding the Lipschitz constant of the gradient of deep neural networks. We discuss existing works in more details in Appendix A.

Below, we state the key theoretical results of this paper informally while detailed statements of these results are presented in Section 4.

Theorem (informal) 1. Let z (L) i denotes the i th logit of an L layer fully-connected neural network with differentiable activation functions. Then, the curvature of the neural network function is globally bounded as follows:

mI ≼ ∇ 2 x z (L) i ≼ M I, ∀x ∈ R D

where m and M can be computed efficiently using parameters of the network.

This result along with the min-max theorem leads to the following curvature robustness certificate:

Theorem (informal) 2. Consider a network whose curvature values are bounded. For a given input x (0) with the true label y and the attack target t (t ≠ y), let p * cert denote the exact robustness certificate, i.e. the distance of x (0) to the decision boundary. We can efficiently compute d * cert such that d * cert ≥ p * cert . Moreover, if the solution x (cert) for d * cert satisfies z

(L) y = z (L) t , then d * cert = p * cert .

We have similar results for the attack problem. For simplicity, we summarize definitions of p * cert , d * cert , p * attack , d * attack in Table 1 . In summary, in this paper, we make the following contributions:

• We provide global bounds on the eigenvalues of the Hessian of a deep neural network with differentiable activation functions (Theorem 3 and Theorem 4). In addition to the adversarial robustness problem, these bounds may be of an independent interest for readers.

Certificate problem (-) = cert:

Attack problem (-) = attack primal problem, p * (-)

min f (x)=0 1 2 x -x (0) 2 min x-x (0) ≤ρ f (x) dual function, d (-) (η) min x 1 2 x -x (0) 2 + ηf (x) min x f (x) + η 2( x -x (0) 2 -ρ 2 )

When is dual solvable? 

-1 M ≤ η ≤ -1 m -m ≤ η dual problem, d * (-) max -1 M ≤η≤-1 m d cert (η) max -m≤η d attack (η) When primal = dual? f (x (cert) ) = 0 x (attack) -x (0) = ρ

(L) ∶ R D → R C

where D is the dimension of the input. For an input x, we use z (I) (x) ∈ R N I and a (I) (x) ∈ R N I to denote the input (before applying the activation function) and output (after applying the activation function) of neurons in the I th hidden layer of the network, respectively. To simplify notation and when no confusion arises, we make the dependency of z (I) and a (I) to x implicit. We define a (0) (x) = x and N 0 = D.

With a fully connected architecture, each z (I) and a (I) is computed using a transformation matrix W (I) ∈ R N I ×N I-1 , the bias vector b (I) ∈ R N I and an activation function σ(.) as follows:

z (I) (x) = W (I) a (I-1) (x) + b (I) , a (I) (x) = σ z (I) (x) .

We use (z

(L) i -z (L) j )(x) as a shorthand for z (L) i (x) -z (L) j (x).

We use [p] to denote the set {1, . . . , p} and [p, q], p ≤ q to denote the set {p, p + 1, . . . , q}. We use small letters i, j, k etc to denote the index over a vector or rows of a matrix and capital letters I, J to denote the index over layers of network. The element in the i th position of a vector v is given by v i , the vector in the i th row of a matrix A is A i while the element in the i th row and j th column of A is A i,j . We use v and A to denote the 2-norm and the operator 2-norm of the vector v and the matrix A, respectively. We use v and A to denote the vector and matrix constructed by taking the elementwise absolute values. We use λ max (A) and λ min (A) to denote the largest and smallest eigenvalues of a symmetric matrix A. We use diag(v) to denote the diagonal matrix constructed by placing each element of v along the diagonal. We use ⊙ to denote the Hadamard Product, I to denote the identity matrix. We use ≼ and ≽ to denote Linear Matrix Inequalities (LMIs) such that given two symmetric matrices A and B where A ≽ B means A -B Positive Semi-Definite (PSD).

3. USING DUALITY TO SOLVE THE ATTACK AND CERTIFICATE PROBLEMS:

Consider an input x (0) with true label y and the attack target t. In the certificate problem, our goal is to find a lower bound of the minimum l 2 distance between x (0) and the decision boundary, z

(L) y = z (L)

t . The problem for solving the exact distance (primal) can be written as:

p * cert = min z (L) y (x)=z (L) t (x) 1 2 x -x (0) 2 = min x max η 1 2 x -x (0) 2 + η z (L) y -z (L) t (x) . (1)

However, solving the above problem can be hard in general. Using the minimax theorem (primal ≥ dual), we can write the dual of the above problem as follows:

EQUATION

From the theory of duality, we know that d cert (η) for each value of η gives a lower bound on the exact certification value (the primal solution) p * cert . However, since z

(L) y -z (L) t

is non-convex, solving d cert (η) for every η can be difficult. In the next section, we will prove that the curvature of the function z

EQUATION

In this case, we have the following theorem:

Theorem 1. d cert (η) is a convex optimization problem for -1 M ≤ η ≤ -1 m. Moreover, If x (cert) is the solution to d * cert such that z (L) y (x (cert) ) = z (L) t (x (cert) ), then p * cert = d * cert .

Below, we briefly outline the proof while the full proof is presented in Appendix D.1. The Hessian of the objective function of the dual d cert (η), i.e the function inside the min x is given by:

∇ 2 x 1 2 x -x (0) 2 + η z (L) y -z (L) t (x) = I + η∇ 2 x z (L) y -z (L) t

From equation (3), we know that the eigenvalues of I + η∇ 2 x (z

(L) y -z (L) t ) are bounded between (1 + ηm, 1 + ηM ) if η ≥ 0, and in (1 + ηM, 1 + ηm) if η ≤ 0.

In both cases, we can see that for -1 M ≤ η ≤ -1 m, all eigenvalues will be non-negative, making the objective function convex. When x (cert) satisfies z 

(L) y = z (L) t , d * cert = 1 2 x (cert) -x (0)

such that x (attack) -x (0) = ρ, p * attack = d * attack .

The proof is presented in Appendix D.2. Both Theorems 1, 2 hold for any non-convex function with continuous gradients. They can also be of interest in problems such as optimization of neural nets.

Using Theorems 1 and 2, we have the following definitions for certification and attack optimizations:

Definition 1. (Curvature-based Certificate Optimization) Given an input x (0) with true label y, the false target t, we define (η (cert) , x (cert) ) as the solution of the following max-min optimization:

max -1 M ≤η≤-1 m min x 1 2 x -x (0) 2 + η z (L) y -z (L) t (x)

We refer to x (cert)x (0) as the Curvature-based Robustness Certificate (CRC).

Definition 2. (Curvature-based Attack Optimization) Given input x (0) with label y, false target t, and the l 2 ball radius ρ, we define (η (attack) , x (attack) ) as the solution of the following optimization:

max η≥-m min x η 2 x -x (0) 2 -ρ 2 + z (L) y -z (L) t (x) .

When x (attack) is used for training in an adversarial training framework, we call the method the Curvature-based Robust Training (CRT).

Since both curvature-based certificate and attack optimizations are convex optimization problems, any convex optimization solver can be used to solve them. In our implementation, we use majorizationminimization to solve the dual function for a given η and bisection method to maximize over η. More details are given in Appendix C.4 and C.5.

4. CURVATURE BOUNDS FOR DEEP NETWORKS:

In this section, we provide a computationally efficient approach to compute the curvature bounds for neural networks with differentiable activation functions. To the best of our knowledge, there is no prior work on finding provable bounds on the curvature values of deep neural networks. Our results rely on a closed form expression for the Hessian of the i th logit as a sum of matrix products (Section 4.1). After establishing this result, we first derive curvature bounds for a two-layer network in Section 4.2 and then extend the bounds to deeper networks in Section 4.3.

4.1. CLOSED FORM EXPRESSION FOR THE HESSIAN:

Using the chain rule of second derivatives, we can derive ∇ 2 x z (L) i as a sum of matrix products: Lemma 1. Given an L layer neural network, the Hessian of the i th hidden unit with respect to the input x, i.e ∇ 2 x z (L) i is given by the following formula:

∇ 2 x z (L) i = L-1 I=1 B (I) T diag F (L,I) i ⊙ σ ′′ z (I) B (I)

where B (I) is the Jacobian of z (I) with respect to x (dimensions N I × D), and F (L,I) is the Jacobian of z (L) with respect to a (I) (dimensions

N L × N I ).

The proof is presented in Appendix D.3. Using the chain rule of gradient, we can compute B (I) , F (L,I) matrices in Lemma 1 recursively as follows:

B (1) = W (1) B (I) = W (I) diag σ ′ z (I-1) B (I-1) I ∈ [2, L -1] (4) F (L,L-1) = W (L) F (L,I) = W (L) diag σ ′ z (L-1) F (L-1,I) I ∈ [L -2] (5)

This leads to a fast back-propagation like method that can be used to compute the Hessian. Note that Lemma 1 only assumes a matrix multiplication operation from a (I-1) to z (I) . Since a convolution operation can also be expressed as a matrix multiplication, we can directly extend this lemma to deep convolutional networks. Furthermore, Lemma 1 can also be of independent interest in other related problems such as higher-order interpretation methods for deep learning (e.g. Singla et al. (2019) ).

4.2. CURVATURE BOUNDS FOR TWO LAYER NETWORKS:

For a two-layer network and using Lemma 1, ∇ 2 x z

(2) yz

(2) t is given by:

∇ 2 x z (2) y -z (2) t = W (1) T diag W (2) y -W (2) t ⊙ σ ′′ z (1) W (1)

Note that only σ ′′ (z (1) ) depends on x. We can maximize and minimize each element in the diag term, (W

EQUATION

i ) independently subject to the constraint that σ ′′ (.) is bounded. Using this procedure, we construct matrices P and N that satisfy properties given in the following theorem: Theorem 3. Given a two layer network whose activation function has bounded second derivative:

h L ≤ σ ′′ (x) ≤ h U ∀x ∈ R (a)

We have the following linear matrix inequalities (LMIs):

N ≼ ∇ 2 x z (2) y -z (2) t ≼ P ∀x ∈ R D (b) If h U ≥ 0 and h L ≤ 0, P is a PSD matrix, N is a NSD matrix.

(c) This gives the following global bounds on the eigenvalues of the Hessian:

mI ≼ ∇ 2 x z (2) y -z (2) t ≼ M I,

where M = λ max (P), m = λ min (N)

P and N are independent of x and defined in equations ( 55) and (56) in Appendix D.4.

The proof is presented in Appendix D.4. Because power iteration finds the eigenvalue with largest magnitude, we can use it to find m and M only when P is PSD and N is NSD. We solve for h U and h L for sigmoid, tanh, softplus activation functions in Appendix E and show that this is in fact the case for them. Note that this result does not hold for ReLU networks since the ReLU function is not differentiable. However, in Appendix F , we devise a method to compute the certificate for a two layer ReLU network by finding a quadratic lower bound for z

(2)

EQUATION

t .

Using Lemma 1, we know that ∇ 2 x z (L) i is a sum product of matrices B (I) and F (L,I) i

. Thus, if we can find upper bounds for B (I) and F (L,I) i , we can get upper bounds for ∇ 2 x z (L) i

. Using this intuition (details are presented in Appendix D.5), we have the following result: Theorem 4. Given an L layer neural network whose activation function satifies:

σ ′ (x) ≤ g, σ ′′ (x) ≤ h ∀x ∈ R, the absolute value of eigenvalues of ∇ 2 x z (L) i

is globally bounded by the following quantity:

∇ 2 x z (L) i ≤ h L-1 I=1 r (I) 2 max j S (L,I) i,j , ∀x ∈ R D

where r (I) and S (L,I) are independent of x and defined recursively as:

EQUATION

The above expressions allows for an efficient computation of the curvature bounds for deep networks.

We consider simplification of this result for sigmoid, tanh, softplus activations in Appendix E.

Note that bounds for z

(L) y -z (L) t can be computed by replacing W (L) i with W (L) y -W (L) t

in Theorem 4. The resulting bound is independent of x, and only depends on network weights W (I) , the true label y, and the target t. We denote it with K(W, y, t). To simplify notation, when no confusion arises we denote it with K. In our experiments, for two layer networks, we use M , m from Theorem 3 (since it provides tighter curvature bounds). For deeper networks (L ≥ 3), we use M = K, m = -K.

5. ADVERSARIAL TRAINING WITH CURVATURE REGULARIZATION:

Using Theorem 2 (b), we know that if we solve the curvature-based attack optimization and obtain ρ = x (attack)x (0) , x (attack) is provably the closest adversarial example to x (0) . However, when Figure 1 : K ub and K lb are upper and lower curvature bounds of the network with Sigmoid activations (averaged over (y, t) pairs). When γ = 0 (no curvature regularization), networks adversarially trained with CRT or PGD both have high curvatures. However, CRT even with a small γ leads to a significant decrease in curvature bounds (note the log-scale of y-axis). Similar results hold for networks with Tanh activations (Appendix Figure 2 ) we performed adversarial training (with ρ = 0.5), we found that the curvature bound is loose and almost none of training inputs lead to zero primal-dual gap with ρ = x (attack)x (0) . To fix this issue, we use a regularizer that penalizes the curvature bound, K. Using equations ( 6) and ( 7), we can compute K using absolute value, matrix multiplications, and operator norm W (I) , I ∈ [L] . Since the gradient of operator norm does not exist in standard libraries, we created a new layer where the gradient of W (I) , i.e ∇ W (I) W (I) is given by:

∇ W (I) W (I) = u (I) v (I) T u (I) , v (I) satisfy W (I) v (I) = W (I) u (I)

This approach to compute the gradient of the largest singular value of a matrix has also been used in previous ICLR work (Miyato et al., 2018) . Implementation details are in Appendix C.1. Thus, the per-sample loss for training with curvature regularization is:

cross entropy z (L) (x (0) ), y + γK(W, y, t)

where y is the true label of the input x (0) , t is the target label and γ is the regularizer for penalizing large curvatures. Similar to the adversarial training, in CRT, we use x (attack) instead of x (0) .

6. EXPERIMENTS:

The empirical robust accuracy means the fraction of test samples that were correctly classified after running an l 2 bounded PGD attack (Madry et al., 2018) Comparison with existing certificates: In Table 3 , we compare CRC with CROWN-general (Zhang et al., 2018a) . For 2-layer networks, CRC outperforms CROWN significantly. For deeper networks, CRC works better only when the network is trained with curvature regularization. However, even with small γ = 0.005, we see a significant increase in CRC but a very small drop in the test accuracy (without any adversarial training). We can see that with γ = 0.01, non-trivial certified accuracies of 83.53%, 88.33%, 89.61% can be achieved on 2, 3, 4 layer sigmoid networks, respectively, without any adversarial training. Adversarial training using CRT further boosts certified accuracy to 95.59%, 94.99% and 93.41%, respectively. We observe similar results with Tanh networks in Appendix Section G.2. We show some results on CIFAR-10 dataset in Table 6 . We again observe improvements in the robustness certificate and certified robust accuracy using CRC and CRT.

In Figure 1 , we plot the effect of γ on the curvature upper bound K ub and a lower bound K lb of a 4-layer network with Sigmoid activations. K lb is computed by taking the maximum of the largest eigenvalue of the Hessian across all test images with label y and the second largest logit t, then averaging across different (y, t). Similarly, K ub is the mean of K over all pairs (y, t) (details in Appendix G.5). We observe that without any curvature regularization (when γ = 0), both standard adversarial training with PGD as well as the CRT lead to networks with high curvatures. However, CRT with even a small γ leads to a significant decrease in curvature bounds. Similar trends can be observed for networks with Tanh activations (Appendix Figure 2 ). Curvature bounds are higher for the Tanh networks compared to the Sigmoid ones due to having larger g and h parameters for Tanh in Theorem 4. Moreover, we report curvature bounds for networks with different depth in Appendix Table 12 . We observe that increasing depth increases curvature bounds.

Comparison with existing adversarial training methods: We compare CRT with adversarial training methods namely PGD (Madry et al., 2018) and TRADES (Zhang et al., 2019) in Table 4 .

We observe that none of the other methods give higher certified accuracy or robustness certificates than our proposed methods. We observe similar results with Tanh networks (Appendix Table 9 ). Moreover, in Appendix Table 10 , we observe that CRT outperforms Randomized Smoothing (Cohen et al., 2019) for 2 and 3 layer networks. Since TRADES and Randomized Smoothing were designed for untargeted attacks while CRT is for targeted attacks, to have a fair comparison, we modify the multi-class version of the cross entropy loss with its binary version (details in Appendix Section G.3). However, we emphasize that there is a fundamental difference between our certificate (and CROWN) with the smoothing based method since our certificate is deterministic while the smoothing-based certificate is probabilistic (with high probability).

We also compare against robust training technique proposed in Wong et al. (2018) in Table 5 with l 2 radius ρ = 1.58. We note that the activation functions used in this comparison are different because the method proposed in Wong et al. (2018) uses ReLU while our proposed CRT can only work with fully differentiable activation functions such as softplus, sigmoid, tanh. For CRT, we use a network with softplus activation function with curvature regularization. We observe that for two and three layer networks, CRT gives higher certified robust accuracy compared to that of Wong et al. (2018) while the difference between certified accuracy of the two approaches for a 4 layer network is small (around 0.22%). We emphasize that to get good certified accuracy using our method, we need to add curvature regularization during the training. Otherwise, the curvature bounds and therefore the robustness certificates are loose. To use local curvature bounds, we need to ensure that the gradient descent trajectory does not escape the local region where curvature values are bounded. To do this, we reduce the step size by a factor of two until the trajectory lies inside the l 2 ball of radius ρ around the input point.

It is straightforward to see that local bounds on curvature values can be tighter than global ones which can lead to better robustness certifications. In Table 7 , we show significant improvements for the CRC certificate and some improvements on the certified accuracy for two-layer sigmoid and tanh networks on the MNIST dataset. In this case, to compute local curvature bounds, we need to have a bound on h L and h U (Theorem 3) within the ball of radius ρ around the input sample. These bounds can be significantly tighter than the global bounds. For example, if a sigmoid neuron's output is bounded by -0.4 and 0.6 (those bounds can be obtained efficiently using CROWN), the second derivative of sigmoid is bounded between -0.048 and +0. For deeper networks, however, computing local curvature bounds is more challenging. We leave exploring this direction for the future work.

7.2. EXTENSION TO CONVOLUTIONAL NEURAL NETWORKS:

The formula derived in Lemma 1 is valid even for convolutional neural networks. However, to use our methods, we need to bound the singular values of the Jacobian of the convolution operation. In order to do this, one can use spectral bounds for convolution layers derived in Sedghi et al. (2018) . Furthermore, other heuristic techniques such as the one proposed in Miyato et al. (2018) can be used for the curvature regularization as well. We present some prelimiary results using these techniques for 2 layer convolutional networks on MNIST and CIFAR-10 in Appendix Section G.1.

8. CONCLUSION:

In this paper, we develop computationally-efficient convex relaxations for robustness certification and adversarial attack problems given the classifier has a bounded curvature. We also show that this convex relaxation is tight under some general conditions. To be able to use proposed certification and attack convex optimizations, we derive global curvature bounds for deep networks with differentiable activation functions. This result is a consequence of a closed-form expression that we derived for the Hessian of a deep network. Our empirical results indicate that our proposed curvature-based robustness certificate outperforms the CROWN certificate by an order of magnitude while being faster to compute as well. Furthermore, adversarial training using our attack method coupled with curvature regularization results in a significantly higher certified robust accuracy than the existing adversarial training methods. Scaling up our proposed curvature-based robustness certification and training methods as well as further tightening the derived curvature bounds are among interesting directions for the future work.

Empirical defenses:

The best known empirical defense is adversarial training (Kurakin et al., 2016b; Madry et al., 2018; Zhang et al., 2019) . In this method, a neural network is trained to minimize the worst-case loss over a region around the input. Although such defenses seem to work on existing attacks, there is no guarantee that a more powerful attack would not break them. In fact, most such defenses proposed in the literature were later broken by stronger attacks (Athalye et al., 2018; Carlini & Wagner, 2017; Uesato et al., 2018; Athalye & Carlini, 2018) . To end this arms race between defenses and attacks, a number of works have tried to focus on certified defenses that have formal robustness guarantees.

Certified defenses A classifier is said to be certifiably robust if one can easily obtain a guarantee that a classifier's prediction remains constant within some region around the input. Such defenses typically rely on certification methods which are either exact or conservative. Exact methods report whether or not there exists a adversarial perturbation inside some l p norm ball. In contrast, conservative methods either certify that no adversarial perturbation exists or decline to make a certification; they may decline even when no such perturbation exists. Exact methods are usually based on Satisfiability Modulo Theories (Huang et al., 2016; Katz et al., 2017; Ehlers, 2017; Carlini et al., 2017) and Mixed Integer linear programming (Cheng et al., 2017; Lomuscio & Maganti, 2017; Dutta et al., 2018; Fischetti & Jo, 2018; Bunel et al., 2017) . Unfortunately, they are computationally inefficient and difficult to scale up to even moderately sized neural networks. In contrast, conservative methods are more scalable and efficient which makes them useful for building certified defenses (Wang et al., 2018a; Wong & Kolter, 2017; Wang et al., 2018b; Wong et al., 2018; Raghunathan et al., 2018b; a; Dvijotham et al., 2018a; b; Croce et al., 2018; Singh et al., 2018; Gowal et al., 2018; Gehr et al., 2018; Mirman et al., 2018; Zhang et al., 2018b; Weng et al., 2018) . However, even these methods have not been shown to scale to practical networks that are large and expressive enough to perform well on ImageNet, for example. To scale to such large networks, randomized smoothing has been proposed as a probabilistically certified defense.

Randomized smoothing Randomized smoothing was previously proposed by several works (Liu et al., 2017; Cao & Gong, 2017) 

B THE ATTACK PROBLEM:

For a given input x (0) with true label y and attack target t, consider the attack problem. We are given that the eigenvalues of the Hessian ∇ 2 x z

(L) y -z (L) t

are bounded below i.e:

mI ≼ ∇ 2 x z (L) y -z (L) t ∀x ∈ R D Here m < 0 (since z (L) y -z (L) t is not convex in general).

The goal here is to find an adversarial example inside a l 2 ball of radius ρ such that (z

(L) y -z (L) t )(x) is minimized.

That is, we want to solve the following optimization:

p * attack = min x-x (0) ≤ρ z (L) y -z (L) t (x) = min x max η≥0 z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2 (8)

This optimization can be hard in general. Using the max-min inequality (primal ≥ dual), we have:

p * attack ≥ max η≥0 d attack (η), d attack (η) = min x z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2 (9)

We know that for every η ≥ 0, d attack (η) gives a lower bound to the primal solution p * attack . But solving d attack (η) for any η ≥ 0 can be hard unless the objective is convex. We prove that if the eigenvalues of the Hessian are bounded below i.e:

mI ≼ ∇ 2 x z (L) y -z (L) t ∀x ∈ R D

In general m < 0, since (z

(L) y -z (L)

t ) is non-convex. d attack (η) is a convex optimization problem for -m ≤ η. Equivalently the objective function, i.e the function inside the min x :

z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2 is a convex function in x for -m ≤ η

The Hessian of the above function is given by:

∇ 2 x z (L) y -z (L) t + ηI

Since we know that eigenvalues of ∇ 2 x (z

(L) y -z (L)

t ) ≽ mI, we know that eigenvalues of the above Hessian are ≥ η + m. For η ≥ -m, the eigenvalues are positive implying that the objective function is convex.

Since d attack (η) gives a lower bound to p * attack for every η ≥ 0, we get the following result:

EQUATION

Note that if x (attack) is the solution to d * attack such that: Our objective is to compute derivative of the largest singular value, i.e W (I) with respect to W (I) . I) are the singular vectors such that:

x (attack) -x (0) = ρ,

EQUATION

W (I) v (I) = W (I) u (I)

v (I) , W (I) 2 can be computed by running power iteration on W (I) T W (I) . u (I) can be computed using the identity:

u (I) = W (I) v (I) γ (I)

We use 25 iterations of the power method to compute the above quantities.

C.2 UPDATE EQUATION FOR THE CERTIFICATE PROBLEM:

Our goal is to minimize xx (0) such that z (L) y z (L) t (x) = 0. We know that the Hessian satisfies the following LMIs:

EQUATION

K is given by Theorem 4 for neural network of any depth (L ≥ 2). For 2 layer networks, M and m are given by Theorem 3. But for deeper networks (L ≥ 3), M = K, m = -K. In either case, K ≥ max( m , M ). Thus, we also have:

-KI ≼ ∇ 2 x z (L) y -z (L) t ≼ KI (12)

We will solve the dual (d * cert ) of the attack problem (p * cert ). The primal problem (p * cert ) is given by:

p * cert = min z (L) y (x)=z (L) t (x) 1 2 x -x (0) 2 = min x max η 1 2 x -x (0) 2 + η z (L) y -z (L) t (x)

Using inequality (11) and Theorem 1 part (a), we know that the dual of the above problem is convex when

-1 M ≤ η ≤ -1 m.

The corresponding dual problem (d * cert ) is given by:

d * cert = max -1 M ≤η≤-1 m d cert (η), d cert (η) = min x 1 2 x -x (0) 2 + η z (L) y -z (L) t (x)

For a given η, we have the following optimization:

d cert (η) = min x 1 2 x -x (0) 2 + η z (L) y -z (L) t (x)

We will use majorization-minimization to solve this optimization.

At a point x (k) , we aim to solve for the point x (k+1) that decreases the objective function. Using the Taylor's theorem at point x (k) , we have:

z (L) y -z (L) t (x) = z (L) y -z (L) t x (k) + g (k) T x -x (k) + 1 2 x -x (k) T H (ξ) x -x (k)

where g (k) is the gradient of (z k) and H (ξ) is the Hessian at a point ξ on the line connecting x and x (k) .

(L) y -z (L) t ) at x (

Multiplying both sides by η, we get the following equation:

EQUATION

Using inequality (12), we know that -KI ≼ H

EQUATION

Using equation ( 13) and inequality ( 14):

η z (L) y -z (L) t (x) ≤ η z (L) y -z (L) t (x (k) ) + η g (k) T (x -x (k) ) + ηK 2 x -x (k) 2

Adding 1 2 xx (0) 2 to both sides, we get the following inequality:

1 2 x -x (0) 2 + η z (L) y -z (L) t (x) ≤ 1 2 x -x (0) 2 + η z (L) y -z (L) t (x (k) ) + η g (k) T (x -x (k) ) + ηK 2 x -x (k) 2

LHS is the objective function of d cert (η) and RHS is an upper bound. In majorization-minimization, we minimize an upper bound on the objective function. Thus we set the gradient of RHS with respect to x to zero and solve for x:

∇ x 1 2 x -x (0) 2 + η z (L) y -z (L) t (x (k) ) + η g (k) T (x -x (k) ) + ηK 2 x -x (k) 2 = 0 x -x (0) + ηg (k) + ηK x -x (k) = 0 (1 + ηK )x -x (0) + ηg (k) -ηK x (k) = 0 x = -(1 + ηK ) -1 ηg (k) -ηK x (k) -x (0)

This gives the following iterative equation:

x (k+1) = -(1 + ηK ) -1 ηg (k) -ηK x (k) -x (0) (15) C.3 UPDATE EQUATION FOR THE ATTACK PROBLEM

Our goal is to minimize z

(L) y -z (L) t

within an l 2 ball of radius of ρ. We know that the Hessian satisfies the following LMIs:

EQUATION

K is given by Theorem 4 for neural network of any depth (L ≥ 2). For 2 layer networks, M and m are given by Theorem 3. But for deeper networks (L ≥ 3), M = K, m = -K. In either case, K ≥ max( m , M ). Thus, we also have:

EQUATION

We solve the dual (d * attack ) of the attack problem (p * attack ) for the given radius ρ. The primal problem (p * attack ) is given by:

p * attack = min x-x (0) ≤ρ z (L) y -z (L) t = min x max η≥0 z (L) y -z (L) t + η 2 x -x (0) 2 -ρ 2

Using inequality ( 16) and Theorem 2 part (a), we know that the dual of the above problem is convex when -m ≤ η.

The corresponding dual problem (d * cert ) is given by:

d * attack = max η≥-m d attack (η), d attack (η) = min x z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2

For a given η, we have the following optimization:

d attack (η) = min x z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2

We will use majorization-minimization to solve this optimization.

At a point x (k) , we have to solve for the point x (k+1) that decreases the objective function. Using the Taylor's theorem at point x (k) , we have:

EQUATION

where g (k) is the gradient of (z k) and H (ξ) is the Hessian at a point ξ on the line connecting x and x (k) .

(L) y -z (L) t ) at x (

Using inequality 17, we know that -KI ≼ H

EQUATION

Using equation ( 18) and inequality (19):

z (L) y -z (L) t (x) ≤ z (L) y -z (L) t (x (k) ) + g (k) T x -x (k) + K 2 x -x (k) 2

Adding η 2( xx (0) 2ρ 2 ) to both sides, we get the following inequality:

z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2 ≤ z (L) y -z (L) t (x (k) ) + g (k) T x -x (k) + K 2 x -x (k) 2 + η 2 x -x (0) 2 -ρ 2

LHS is the objective function of d attack (η) and RHS is an upper bound. In majorization-minimization, we minimize an upper bound on the objective function. Thus we set the gradient of RHS with respect to x to zero and solve for x:

∇ x z (L) y -z (L) t (x (k) ) + g (k) T x -x (k) + K 2 x -x (k) 2 + η 2 x -x (0) 2 -ρ 2 = 0 g (k) + K x -x (k) + η x -x (0) = 0 (K + η)x + g (k) -Kx (k) -ηx (0) = 0 x = -(K + η) -1 g (k) -Kx (k) -ηx (0)

This gives the following iterative equation:

EQUATION

C.4 ALGORITHM TO COMPUTE THE CERTIFICATE:

We start with the following initial values of x, η, η min , η max :

η min = -1 M, η max = -1 m η = 1 2 (η min + η max ), x = x (0)

To solve the dual for a given value of η, we run 20 iterations of the following update (derived in Appendix C.2):

Algorithm 1 Certificate optimization Require: input x (0) , label y, target t m, M, K ← compute bounds(z

(L) y -z (L) t ) η min ← -1 M η max ← -1 m η ← 1 2(η min + η max ) x ← x (0) for i in [1, . . . , 30] do for j in [1, . . . , 20] do g ← compute gradient(z (L) y -z (L) t , x) if ηg + (x -x (0) ) < 10 -5 then break end if x ← -(1 + ηK ) -1 ηg -ηK x -x (0) end for if (z (L) y -z (L) t )(x) > 0 then η min ← η else η max ← η end if η ← (η min + η max ) 2 end for return x

Algorithm 2 Attack optimization Require: input x (0) , label y, target t , radius ρ m, M, K ← compute bounds(z

(L) y -z (L) t ) η min ← -m η max ← 20(1 -m) η ← 1 2(η min + η max ) x ← x (0) for i in [1, . . . , 30] do for j in [1, . . . , 20] do g ← compute gradient(z (L) y -z (L) t , x) if g + η(x -x (0) ) < 10 -5 then break end if x ← -(K + η) -1 g -Kx -ηx (0) end for if x -x (0) < ρ then η max ← η else η min ← η end if η ← (η min + η max ) 2 end for return x D PROOFS D.1 PROOF OF THEOREM 1 (a) d cert (η) = min x 1 2 x -x (0) 2 + η z (L) y (x) -z (L) t (x) ∇ 2 x 1 2 x -x (0) 2 + η z (L) y (x) -z (L) t (x) = I + η∇ 2 x z (L) y -z (L) t

We are given that the Hessian ∇ 2 x (z

(L) y -z (L)

t ) satisfies the following LMIs:

mI ≼ ∇ 2 x z (L) y -z (L) t ≼ M I ∀x ∈ R n

The eigenvalues of I + η∇ 2 x (z

(L) y -z (L)

t ) are bounded between:

(1 + ηM, 1 + ηm), if η < 0 (1 + ηm, 1 + ηM ), if η > 0

We are given that η satisfies the following inequalities where m < 0, M > 0 since (z

(L) y - z (L)

t ) is neither convex, nor concave as a function of x:

-1

M ≤ η ≤ -1 m , m < 0, M > 0

We have the following inequalities: Let η (cert) , x (cert) be the solution of the above dual optimization (d * cert ) such that z (L) y (x

1 + ηM ≥ 0, 1 + ηm ≥ 0 Thus, I + η∇ 2 x (z (L) y -z (L) t ) is a PSD matrix for all x ∈ R D when -1 M ≤ η ≤ -1 m . Thus 1 2 x -x (0) 2 + η(z (L) y -z (L) t )(x)

(cert) ) = z (L) t (x (cert) ) (22) d *

cert is given by the following:

d * cert = 1 2 x (cert) -x (0) 2 + η (cert) z (L) y (x (cert) ) -z (L) t (x (cert) ) =0

Since we are given that z

(L) y (x (cert) ) = z (L) t (x (cert)

), we get the following equation for d * cert :

d * cert = 1 2 x (cert) -x (0) 2 (23)

Since p * cert is given by the following equation:

p * cert = min z (L) y (x)=z (L) t (x) 1 2 x -x (0) 2 (24)

Using equations ( 22) and ( 24), p * cert is the minimum value of 1 2 xx (0) 2 ∀x ∶ z

(L) y (x) = z (L) t (x): p * cert ≤ 1 2 x (cert) -x (0) 2 (25)

From equation ( 23), we know that d * cert = 1 2 x (cert)x (0) 2 . Thus, we get: 

EQUATION

d attack (η) = min x z (L) y -z (L) t (x) + η 2 x -x (0) 2 -ρ 2 ∇ 2 x z (L) y -z (L) t (x) + η 2 x -x (0) 2 = ∇ 2 x z (L) y -z (L) t + ηI

Since the Hessian ∇ 2 x (z

(L) y -z (L) t ) is bounded below: mI ≼ ∇ 2 x z (L) y -z (L) t ∀x ∈ R n

The eigenvalues of ∇ 2 x (z

(L) y -z (L)

t ) + ηI are bounded below: 

(m + η)I ≼ ∇ 2 x z (L) y -z (L) t + ηI Since η ≥ -m. η + m ≥ 0 Thus ∇ 2 x (z (L) y -z (L) t ) + ηI is a PSD matrix for all x ∈ R D when η ≥ -m. Thus (z (L) y -z (L) t )(x) + η 2( x -x (0) 2 -ρ 2 ) is a convex

Let η (attack) , x (attack) be the solution of the above dual optimization (d * attack ) such that

x (attack) -x (0) = ρ (28)

d * attack is given by the following:

d * attack = z (L) y -z (L) t (x (attack) ) + η (attack) 2 x (attack) -x (0) 2 -ρ 2 =0

Since we are given that x (attack)x (0) = ρ, we get the following equation for d * attack :

EQUATION

Since p * attack is given by the following equation:

EQUATION

Using equations ( 28) and (30), p * attack is the minimum value of (z

EQUATION

From equation ( 29), we know that d * attack = (z We have to prove that for an L layer neural network, the hessian of the i th hidden unit in the L th layer with respect to the input x, i.e ∇ 2 x z (L) i is given by the following formula:

(L) y -z (L) t )(x (attack)

EQUATION

where B (I) , I ∈ [L] is a matrix of size N I × D defined as follows:

EQUATION

and F (L,I) , I ∈ [L -1] is a matrix of size N L × N I defined as follows:

F (L,I) = ∇ a (I) z (L) 1 , ∇ a (I) z (L) 2 , . . . , ∇ a (I) z (L) N L T , I ∈ [L -1] (35) ∇ 2 x z (L) i

can be written in terms of the activations of the previous layer using the following formula:

∇ 2 x z (L) i = N I-1 j=1 W (L) i,j ∇ 2 x a (L-1) j (36)

Using the chain rule of the Hessian and a (I) = σ(z (I) ), we can write ∇ 2 x a

(L-1) j in terms of ∇ x z (L-1) j and ∇ 2 x z (L-1) j

as the following:

EQUATION

Replacing ∇ 2 x a (L-1) j using equation (37) into equation (36), we get:

EQUATION

For each I ∈ [2, L], i ∈ N I , we define the matrix A (I) i as the following:

∇ 2 x z (I) i = N I-1 j=1 W (I) i,j σ ′′ z (I-1) j ∇ x z (I-1) j ∇ x z (I-1) j T A (I) i + N I-1 j=1 W (I) i,j σ ′ z (I-1) j ∇ 2 x z (I-1) j (39) A (I) i = N I-1 j=1 W (I) i,j σ ′′ z (I-1) j ∇ x z (I-1) j ∇ x z (I-1) j T (40) Substituting A (L) i

using equation (40) into equation ( 38), we get:

EQUATION

)

We first simplify the expression for

A (L) i . Note that A (L) i

is a sum of symmetric rank one matrices

∇ x z (L-1) j ∇ x z (L-1) j T with the coefficient W (L) i,j σ ′′ z (L-1) j

for each j. We create a diagonal matrix for the coefficients and another matrix B (L-1) such that each j th row of B (L-1) is the vector

∇ x z (L-1) j

. This leads to the following equation:

EQUATION

B (I) where I ∈ [L] is a matrix of size N I × D defined as follows:

B (I) = ∇ x z (I) 1 , ∇ x z (I) 2 , . . . , ∇ x z (I) N I T , I ∈ [L]

Thus B (I) is the jacobian of z (I) with respect to the input x.

Using the chain rule of the gradient, we have the following properties of B (I) :

EQUATION

EQUATION

Similarly, F (I,J) where I ∈ [L], J ∈ [I -1] is a matrix of size N I × N J defined as follows:

F (I,J) = ∇ a (J) z (I) 1 , ∇ a (J) z (I) 2 , . . . , ∇ a (J) z (I) N I T , I ∈ [L], J ∈ [I -1]

Thus F (I,J) is the jacobian of z (I) with respect to the activations a (J) .

Using the chain rule of the gradient, we have the following properties for F (L,I) :

EQUATION

EQUATION

Recall that in our notation: For a matrix E, E i denotes the column vector constructed by taking the transpose of the i th row of the matrix E. Thus i th row of

W (L) is W (L) i T and F (L,I) is F (L,I) i T .

Equating the i th rows in equation ( 46), we get:

F (L,I) i T = W (L) i T diag σ ′ z (L-1) F (L-1,I)

Taking the transpose of both the sides and expressing the RHS as a summation, we get:

F (L,I) i = W (L) i T diag σ ′ z (L-1) F (L-1,I) T = N L-1 j=1 W (L) i,j σ ′ z (L-1) j F (L-1,I) j (47)

Substituting W (L) using equation ( 45) into equation ( 42), we get:

A (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) (48) Substituting A (L) i

using equation ( 48) into ( 41), we get:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + N L-1 j=1 W (L) i,j σ ′ z (L-1) j ∇ 2 x z (L-1) j (49)

Thus, equation ( 49) allows us to write the hessian of i th unit at layer L, i.e ∇ 2 x z (L) i in terms of the hessian of j th unit at layer L -1, i.e ∇ 2 x z (L-1) j . We will prove the following using induction:

EQUATION

Note that for L = 2, ∇ 2 x z (L-1) j = 0, ∀j ∈ N 1 . Thus using (49) we have:

∇ 2 x z (2) i = B (1) T diag F (2,1) i ⊙ σ ′′ z (1) B (1)

Hence the induction hypothesis ( 50) is true for L = 2. Now we will assume ( 50) is true for L -1. Thus we have:

∇ 2 x z (L-1) j = L-2 I=1 B (I) T diag F (L-1,I) j ⊙ σ ′′ z (I) B (I) ∀j ∈ N L-1 (51)

We will prove the same for L.

Using equation ( 49), we have:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + N L-1 j=1 W (L) i,j σ ′ z (L-1) j ∇ 2 x z (L-1) j

In the next set of steps, we will be working with the second term of the above equation, i.e

∑ N L-1 j=1 W (L) i,j σ ′ (z (L-1) j )(∇ 2 x z (L-1) j

):

Substituting ∇ 2 x z

(L-1) j using equation ( 51) we get:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + N L-1 j=1 W (L) i,j σ ′ z (L-1) j L-2 I=1 B (I) diag F (L-1,I) j ⊙ σ ′′ z (I) B (I) T

Combining the two summations in the second term, we get:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + N L-1 j=1 L-2 I=1 W (L) i,j σ ′ z (L-1) j B (I) T diag F (L-1,I) j ⊙ σ ′′ z (I) B (I)

Exchanging the summation over I and summation over j:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + L-2 I=1 N L-1 j=1 W (L) i,j σ ′ z (L-1) j B (I) T diag F (L-1,I) j ⊙ σ ′′ z (I) B (I)

Since B (I) is independent of j, we take it out of the summation over j:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + L-2 I=1 B (I) T N L-1 j=1 W (L) i,j σ ′ z (L-1) j diag F (L-1,I) j ⊙ σ ′′ z (I) B (I)

Using the property,

α (diag(u)) + β (diag(v)) = diag (αu + βv) ∀α, β ∈ R, u, v ∈ R n ;

we can move the summation inside the diagonal:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + L-2 I=1 B (I) T diag N L-1 j=1 W (L) i,j σ ′ z (L-1) j F (L-1,I) j ⊙ σ ′′ z (I) B (I)

Since σ ′′ z (I) is independent of j, we can take it out of the summation over j:

∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + L-2 I=1 B (I) T diag N L-1 j=1 W (L) i,j σ ′ z (L-1) j F (L-1,I) j ⊙ σ ′′ z (I) B (I)

Using equation ( 47), we can replace

∑ N L-1 j=1 W (L) i,j σ ′ z (L-1) j F (L-1,I) j with F (L,I) i : ∇ 2 x z (L) i = B (L-1) T diag F (L,L-1) i ⊙ σ ′′ z (L-1) B (L-1) + L-2 I=1 B (I) T diag F (L,I) i ⊙ σ ′′ z (I) B (I) ∇ 2 x z (L) i = L-1 I=1 B (I) T diag F (L,I) i ⊙ σ ′′ z (I) B (I)

D.4 PROOF OF THEOREM 3:

Using Lemma 1, we have the following formula for ∇ 2 x z

(2)

EQUATION

)

We are also given that the activation function σ satisfies the following property:

EQUATION

(a) We have to prove the following linear matrix inequalities (LMIs):

EQUATION

where P and N are given as following:

EQUATION

)

EQUATION

EQUATION

y,i -W (2) t,i ≥ 0 h U , W (2) y,i -W (2) t,i ≤ 0 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ (57)

We first prove:

N ≼ ∇ 2 x z (2) y -z (2) t ∀x ∈ R D :

We substitute ∇ 2 x z

(2)

y -z (2) t

and N from equations ( 52) and ( 56) respectively in

∇ 2 x z (2) y -z (2) t -N: ∇ 2 x z (2) y -z (2) t -N = N1 i=1 W (2) y,i -W (2) t,i σ ′′ z (1) i -W (2) y,i -W (2) t,i n i W (1) i W (1) i T = N1 i=1 W (2) y,i -W (2) t,i σ ′′ z (1) i -n i W (1) i W (1) i T Thus ∇ 2 x z (2) y -z (2) t

-N is a weighted sum of symmetric rank one matrices i.e,

W:

(1) i

(1) i T and it is PSD if and only if coefficient of each rank one matrix i.e,

(2)

y,i -W (2) t,i σ ′′ z (1) i

n i is positive. Using equations ( 53) and ( 57), we have the following:

W (2) y,i -W (2) t,i ≥ 0 ⇒ n i = h L ⇒ σ ′′ z (1) i -n i ≥ 0 ∀i ∈ [N 1 ], ∀x ∈ R D W (2) y,i -W (2) t,i ≤ 0 ⇒ n i = h U ⇒ σ ′′ z (1) i -n i ≤ 0 ∀i ∈ [N 1 ], ∀x ∈ R D ⇒ W (2) y,i -W (2) t,i σ ′′ z (1) i -n i ≥ 0 ∀i ∈ [N 1 ], ∀x ∈ R D (58) Thus ∇ 2 x z (2) y -z (2) t

-N is a PSD matrix i.e:

EQUATION

n i always positive using eq. ( 58)

EQUATION

Now we prove that ∇ 2 x z

(2)

EQUATION

yz

(2) t and P from equations ( 52) and ( 56) respectively in P -∇ 2 x z

(2)

EQUATION

EQUATION

EQUATION

(2) t is a weighted sum of symmetric rank one matrices i.e, W

(1) i T and it is PSD if and only if coefficient of each rank one matrix i.e, W

EQUATION

(1) i is positive. Using equations ( 53) and (57), we have the following:

EQUATION

yz

(2) t is PSD matrix i.e:

P -∇ 2 x z (2) y -z (2) t = N1 i=1 W (2) y,i -W (2) t,i p i -σ ′′ z

(1) i always positive using eq. ( 60)

EQUATION

Thus by proving the LMIs ( 59) and (61), we prove (54). (b) We have to prove that if h U ≥ 0 and h L ≤ 0, P is a PSD matrix, N is a NSD matrix.

We are given h U ≥ 0, h L ≤ 0. Using equation (57), we have the following:

EQUATION

Thus P is a weighted sum of symmetric rank one matrices i.e, W

(1) i T and each

coefficient p i W (2) y,i -W (2) t,i

is positive.

EQUATION

y,i -W

(2) t,i always positive using eq. ( 62)

EQUATION

Using equation ( 57), we have the following:

EQUATION

N = N1 i=1 n i W (2) y,i -W (2) t,i

always positive using eq. ( 63)

EQUATION

Thus P is a PSD and N is a NSD matrix if h U ≥ 0 and h L ≤ 0.

(c) We have to prove the following global bounds on the eigenvalues of ∇ 2 x (z

EQUATION

Let v * , x * be vectors such that:

(v * ) T ∇ 2 x * z (2) y -z (2) t v * = max x max v =1 v T ∇ 2 x z (2) y -z (2) t v

Thus using inequality (64):

EQUATION

Let v * , x * be vectors such that:

(v * ) T ∇ 2 x * z (2) y -z (2) t v * = min x min v =1 v T ∇ 2 x z (2) y -z (2) t v

Thus using inequality (66):

EQUATION

Using the inequalities ( 65) and ( 67), we get:

mI ≼ ∇ 2 x z (2) y -z (2) t ≼ M I, where M = max v =1 v T Pv, m = min v =1 v T Nv D.5 PROOF OF THEOREM 4

We are given that the activation function σ is such that σ ′ , σ ′′ are bounded, i.e:

EQUATION

We have to prove the following:

∇ 2 x z (L) i ≤ h L-1 I=1 r (I) 2 max j S (I) i,j ∀x ∈ R D

where S (L,I) is a matrix of size N L × N I defined as follows:

EQUATION

and r (I) is a scalar defined as follows:

EQUATION

We will prove the same in 3 steps.

In step (a), we will prove:

EQUATION

In step (b), we will prove:

EQUATION

In step (c), we will use (a) and (b) to prove:

EQUATION

Note that B (I) and F (L,I) are defined using (34) and ( 35) respectively.

(a) We have to prove that for

L ≥ 2, I ∈ [L -1], i ∈ N L , j ∈ N I : F (L,I) i,j ≤ S (L,I) i,j ∀x ∈ R D

where S (L,I) is a matrix of size N I × N J defined as follows:

S (L,I) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ W (L) I = L -1 g W (L) S (L-1,J) I ∈ [L -2]

We first prove the case when

I = L -1. Using equation (45), F (L,L-1) i,j = W (L) i,j . Since S (L,L-1) i,j = W (L) i,j : F (L,L-1) i,j = S (L,L-1) i,j

Hence for L ≥ 2, I = L -1, we have equality in (71). Hence proved. Now, we will use proof by induction.

To prove the base case L = 2, note that I = L -1 = 1 is the only possible value for I. Thus, using the result for I = L -1, the theorem holds for L = 2. This proves the base case. Now we assume the induction hypothesis is true for depth = L -1, I ∈ [L -2]. and prove for depth = L, I ∈ [L -1]. Since for I = L -1, we have proven already, we prove for I ≤ L -2.

Using equation ( 47), we have the following formula for F (L,I) i :

F (L,I) i = N L-1 k=1 W (L) i,k σ ′ z (L-1) k F (L-1,I) k

Taking the j th element of the vectors on both sides:

EQUATION

By induction hypothesis, we know that:

EQUATION

Using the absolute value properties for equation ( 74), we have:

F (L,I) i,j = N L-1 k=1 W (L) i,k σ ′ z (L-1) k F (L-1,I) k,j F (L,I) i,j ≤ N L-1 k=1 W (L) i,k σ ′ z (L-1) k F (L-1,I) k,j F (L,I) i,j ≤ N L-1 k=1 W (L) i,k σ ′ z (L-1) k F (L-1,I) k,j

Using σ ′ (x) ≤ g ∀x ∈ R (inequality (68)) :

F (L,I) i,j ≤ g N L-1 k=1 W (L) i,k F (L-1,I) k,j

Using the induction hypothesis (inequality (75)):

F (L,I) i,j ≤ g N L-1 k=1 W (L) i,k S (L-1,I) k,j

Using equation ( 69) for definition of S (L,I) i,j

:

F (L,I) i,j ≤ S (L,I) i,j

Hence we prove (71) for all L ≥ 2 and I ≤ L -1 using induction.

(b) We have to prove that for 1 ≤ I ≤ M -1:

B (I) ≤ r (I) , ∀x ∈ R D

where r (I) is a scalar given as follows:

r (I) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ W (1) I = 1 g W (I) r (I-1) I ∈ [2, L -1]

Using equation ( 43), for I = 1 we have:

EQUATION

Using equation ( 44), for I > 1, we have:

B (I) = W (I) diag σ ′ z (I-1) B (I-1) B (I) ≤ W (I) diag σ ′ z (I-1) B (I-1) Since diag σ ′ z (I-1) = max j σ ′ z (I-1) j

, using equation (68):

EQUATION

Using inequalities ( 76) and ( 77), the proof follows using induction.

(c) We have to prove that:

∇ 2 x z (L) i ≤ h L-1 I=1 r (I) 2 max j S (I) i,j

Using Lemma 1, we have the following equation for ∇ 2 x z (L) i :

∇ 2 x z (L) i = L-1 I=1 B (I) T diag F (L,I) i ⊙ σ ′′ z (I) B (I)

Using the properties of norm we have:

∇ 2 x z (L) i = L-1 I=1 B (I) T diag F (L,I) i ⊙ σ ′′ z (I) B (I) ≤ L-1 I=1 diag F (L,I) i ⊙ σ ′′ z (I) B (I) 2 ≤ L-1 I=1 max j F (L,I) i,j σ ′′ z (I) j B (I) 2

In the last inequality, we use the property that norm of a diagonal matrix is the maximum absolute value of the diagonal element. Using the product property of absolute value, we get:

∇ 2 x z (L) i ≤ L-1 I=1 max j F (L,I) i,j σ ′′ z (I) j B (I) 2 Since F (L,I) i,j

and σ ′′ z (I) j are positive terms:

∇ 2 x z (L) i ≤ L-1 I=1 max j F (L,I) i,j max j σ ′′ z (I) j B (I) 2

Since σ ′′ is bounded by h:

∇ 2 x z (L) i ≤ h L-1 I=1 max j F (L,I) i,j B (I) 2

Using inequality (71):

∇ 2 x z (L) i ≤ h L-1 I=1 max j S (I) i,j B (I) 2

Using inequality (72):

∇ 2 x z (L) i ≤ h L-1 I=1 r (I) 2 max j S (I) i,j ∀x ∈ R D E COMPUTING g, h, h U AND h L FOR DIFFERENT ACTIVATION FUNCTIONS E.1 SOFTPLUS ACTIVATION

For softplus activation, we have the following. We use S(x) to denote sigmoid:

σ(x) = log(1 + exp(x)) σ ′ (x) = S(x) σ ′′ (x) = S(x)(1 -S(x))

To bound S(x)(1 -S(x)), let α denote S(x). We know that 0 ≤ α ≤ 1:

α(1 -α) = 1 4 - 1 2 -α 2 Thus, S(x)(1 -S(x)

) is maximum at S(x) = 1 2 and minimum at S(x) = 0 and S(x) = 1. The maximum value is 0.25 and minimum value is 0.

0 ≤ S(x)(1 -S(x)) ≤ 0.25 ⇒ 0 ≤ σ ′′ (x) ≤ 0.25

Thus, h U = 0.25, h L = 0 (for use in Theorem 3) and g = 1, h = 0.25 (for use in Theorem 4).

E.2 SIGMOID ACTIVATION:

For sigmoid activation, we have the following. We use S(x) to denote sigmoid:

σ(x) = S(x) = 1 1 + exp(-x) σ ′ (x) = S(x)(1 -S(x)) σ ′′ (x) = S(x)(1 -S(x))(1 -2S(x))

The second derivative of sigmoid (σ ′′ (x)) can be bounded using standard differentiation. Let α denote S(x). We know that 0 ≤ α ≤ 1:

h L ≤ σ ′′ (x) ≤ h U h L = min 0≤α≤1 α(1 -α)(1 -2α) h U = max 0≤α≤1 α(1 -α)(1 -2α)

To solve for both h L and h U , we first differentiate α(1α)(1 -2α) with respect to α:

∇ α (α(1 -α)(1 -2α)) = ∇ α 2α 3 -3α 2 + α = 6α 2 -6α + 1

Solving for 6α 2 -6α + 1 = 0, we get the solutions:

α = 3 + √ 3 6 , 3 - √ 3 6

Since both (3 + √ 3 6), (3 -√ 3 6) lie between 0 and 1, we check for the second derivatives:

∇ 2 α (α(1 -α)(1 -2α)) = ∇ α 6α 2 -6α + 1 = 12α -6 = 6(2α -1) At α = (3 + √ 3) 6, ∇ 2 α = 6(2α -1) = 2 √ 3 > 0. At α = (3 - √ 3) 6, ∇ 2 α = 6(2α -1) = -2 √ 3 < 0. Thus α = (3 + √ 3) 6 is a local minima, α = (3 - √ 3

) 6 is a local maxima. Substituting the two critical points into α(1α)(1 -2α), we get h U = 9.623 × 10 -2 , h L = -9.623 × 10 -2 . Thus, h U = 9.623 × 10 -2 , h L = -9.623 × 10 -2 (for use in Theorem 3) and g = 0.25, h = 0.09623 (for use in Theorem 4).

E.3 TANH ACTIVATION:

For tanh activation, we have the following:

σ(x) = tanh(x) = exp(x) -exp(-x) exp(x) + exp(-x) σ ′ (x) = (1 -tanh(x)) (1 + tanh(x)) σ ′′ (x) = -2 tanh(x) (1 -tanh(x)) (1 + tanh(x))

The second derivative of tanh , i.e (σ ′′ (x)) can be bounded using standard differentiation. Let α denote tanh(x). We know that -1 ≤ α ≤ 1:

h L ≤ σ ′′ (x) ≤ h U h L = min 0≤α≤1 -2α(1 -α)(1 + α) h U = max 0≤α≤1 -2α(1 -α)(1 + α)

To solve for both h L and h U , we first differentiate -2α(1α)(1 + α) with respect to α:

∇ α (-2α(1 -α)(1 + α)) = ∇ α 2α 3 -2α = 6α 2 -2

Solving for 6α 2 -2 = 0, we get the solutions:

EQUATION

Since both -1 √ 3, 1 √ 3 lie between -1 and 1, we check for the second derivatives:

∇ 2 α (-2α(1 -α)(1 + α)) = ∇ α 6α 2 -2 = 12α At α = -1 √ 3, ∇ 2 α = 12α = -4 √ 3 < 0. At α = 1 √ 3, ∇ 2 α = 12α = 4 √ 3 > 0. Thus α = 1 √ 3 is a local minima, α = -1 √

3 is a local maxima. Substituting the two critical points into -2α(1α)(1 + α), we get h U = 0.76981, h L = -0.76981. Thus, h U = 0.76981, h L = -0.76981 (for use in Theorem 3) and g = 1, h = 0.76981 (for use in Theorem 4).

F QUADRATIC BOUNDS FOR TWO-LAYER RELU NETWORKS For a 2 layer network with ReLU activation, such that the input x lies in the ball xx (0) ≤ ρ, we can compute the bounds over z (1) directly:

W (1) i x (0) + b (1) i -ρ W (1) i ≤ z (1) i ≤ W (1) i x (0) + b (1) i + ρ W (1) i

Thus we can get a lower bound and upper bound for each z

(1) i . We define d i and u i as the following:

d i = W (1) i x (0) + b (1) i -ρ W (1) i (78) u i = W (1) i x (0) + b (1) i + ρ W (1) i (79)

We can derive the following quadratic lower and upper bounds for each a

(1) i :

a (1) i ≤ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -d i (u i -d i ) 2 z (1) i 2 + u 2 i + d 2 i (u i -d i ) 2 z (1) i - u 2 i d i (u i -d i ) 2 d i ≤ u i u i (u i -d i ) 2 z (1) i 2 - 2u i d i (u i -d i ) 2 z (1) i + u i d 2 i (u i -d i ) 2 d i ≥ u i a (1) i ≥ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 2 d i ≤ u i z (1) i d i ≥ 2 u i 1 u i -d i z (1) i 2 - d i u i -d i z

(1) i otherwise

The above steps are exactly the same as the quadratic upper and lower bounds used in (Zhang et al., 2018a) .

Using the above two inequalities and the identity:

z (2) y -z (2) t = N1 i=1 W (2) y,i -W (2)

t,i a

(1) i

we can compute a quadratic lower bound for z

(2) yz

(2) t in terms of z

(1) i by taking the lower bound for a

(1) i when W

(2)

y,i -W (2) t,i > 0 and upper bound when W (2) y,i -W (2) t,i <= 0. Furthermore since z (1) i = W (1) i x + b (1)

i , we can express the resulting quadratic in terms of x. Thus, we get the following quadratic function :

z (2) y -z (2) t ≥ 1 2 x T Px + q + r

The coefficients P, q and r can be determined using the above procedure. Note that unlike in (Zhang et al., 2018a) , RHS can be a non-convex function. Thus, it becomes an optimization problem where the goal is to minimize the distance 1 2 xx (0) 2 subject to RHS (which is quadratic in x) being zero. That is both our objective and constraint are quadratic functions. In the optimization literature, this is called the S-procedure and is one of the few non-convex problems that can be solved efficiently (Boyd & Vandenberghe, 2004 ).

We start with two initial values called ρ low (initialized to 0) and ρ high (initialized to 5). We start with an initial value of ρ, initialized at 1 2 (ρ low + ρ high ) to compute d i (eq. ( 78)) and u i (eq. ( 79)). If the final distance after solving the S-procedure is less than ρ, we set ρ low = ρ. if the final distance is greater than ρ, we set ρ high = ρ. Set new ρ = 1 2 (ρ low + ρ high ). Repeat until convergence.

G ADDITIONAL EXPERIMENTS:

Empirical accuracy means the fraction of test samples that were correctly classified after running a PGD attack (Madry et al., 2018) with an l 2 bound on the adversarial perturbations. Certified accuracy means the fraction of test samples that were classified correctly initially and had the robustness certificate greater than a pre-specified attack radius ρ. For both empirical and certified accuracy, we use ρ = 0.5. Unless otherwise specified, we use the class with the second largest logit as the attack target for the given input (i.e. the class t). All experiments were run on the MNIST dataset while noting that our results are scalable for more complex datasets 

G.1 RESULTS ON CONVOLUTIONAL NEURAL NETWORKS:

We use a 2 layer (1 hidden layer network). We use 64 filters in the convolution layer and softplus activation function and stride of length 2. This is followed by reshaping and a fully connected layer. 

G.3 COMPARING RANDOMIZED SMOOTHING AND TRADES WITH CRT:

Randomized smoothing is designed to work in untargeted attack settings while CRT is for targeted attacks. Thus, to do a fair comparison of CRT with randomized smoothing, we make the following changes in randomized smoothing.

First, we use n 0 = 100 initial samples to select the label class (l) and false target class (t). The samples for estimation were n = 100, 000 and failure probability was α = 0.001. Then we use the binary version of randomized smoothing for estimation, i.e classify between y and t. To find the adversarial example for adversarial training, we use the cross entropy loss for 2 classes (y and t).

For TRADES, we select the class with second highest logit as the target class t and use the 2 class version of the cross entropy loss for finding the adversarial example. In Table 12 , we measure how the standard accuracy, empirical accuracy, certified accuracy, upper bound on the curvature K ub , lower bound on the curvature K lb , changes as we increase the regularization parameter γ and the network is trained with CRT.

In Table 13 , we measure how the standard accuracy, empirical accuracy, certified accuracy, CROWN and CRC changes as we increase the regularization parameter γ and the network is trained without any adversarial training.

G.5 COMPUTING K lb AND K ub First, note that K does not depend on the input, but on network weights W (I) , label y and target t.

Different images may still have different K because label y and target t may be different.

To compute K lb in the table, first for each pair y and t, we find the largest eigenvalue of the Hessian of all test images that have label y and second largest logit of class t. Then we take the max of the largest eigenvalue across all test images. This gives a rough estimate of the largest curvature in the vicinity of test images with label y and target t. We can directly take the mean across all such pairs to compute K lb . However, we find that some pairs y and t were infrequent (with barely 1,2 test images in them). Thus, for all such pairs we cannot get a good estimate of the largest curvature in vicinity.

We select all pairs y and t that have at least 100 images in them and compute K lb by taking the mean across all such pairs.

To compute K ub in the table, we compute K for all pairs y and t that have at least 100 images, i.e at least 100 images should have label y and target t. And then we compute the mean across all K that satisfy this condition. This was done to do a fair comparison with K lb .

Figure 2 shows a plot of the K ub and K lb with increasing γ for a tanh network. When γ = 0 (no curvature regularization), networks adversarially trained with CRT or PGD both have high curvatures. However, CRT with even a small γ leads to a significant decrease in curvature bounds. Also, curvature bounds are higher with Tanh than with Sigmoid. Results are similar to Sigmoid in Figure 1 .

Table 12 : In this table, we measure the effect of increasing γ, when the network is trained with CRT on standard, empirical, certified robust accuracy, K lb and K ub (defined in subsection G.5) for different depths (2, 3, 4 layer) and activations (sigmoid, tanh). We find that for all networks γ = 0.01 works best. We find that the lower bound, K lb increases (for γ = 0) for deeper networks suggesting that deep networks have higher curvature. Furthermore, for a given γ (say 0.005), we find that the gap between K ub and K lb increases as we increase the depth suggesting that K is not a tight bound for deeper networks. 

:



Figures and Tables:

is a convex function in x and d cert (η) is a convex optimization problem. (b) For every value of η, d cert (η) is a lower bound for p * cert . Thus d * cert = max -1 M ≤ η ≤-1 m d cert (η) is a lower bound for p * cert , i.e:

. The notation (L × [1024], activation) denotes a neural network with L layers with the specified activation function, (γ = c) denotes standard training with γ set to c, (CRT, c) denotes CRT training with γ = c. Certificates CROWN and CRC are computed over 150 correctly classified images.

Figure 2: K ub and K lb are upper and lower curvature bounds of the network with Tanh activations (averaged over (y, t) pairs).When γ = 0 (no curvature regularization), networks adversarially trained with CRT or PGD both have high curvatures. However, CRT with even a small γ leads to a significant decrease in curvature bounds. Also, curvature bounds are higher with Tanh than with Sigmoid. Results are similar to Sigmoid in Figure1.

A summary of various primal and dual concepts used in the paper. f denotes the function of the decision boundary, i.e. z function of the neural network is z

Table content: <table/>

2 , using the duality theorem and definition of p * cert , we get p * cert = d * cert . Next, we consider the attack problem. The goal here is to find an adversarial example inside an l 2 ball of radius ρ such that z

Table content: <table><tr><td>(L) y theorem for the attack problem (p  *  (L) t is minimized. Using similar arguments, we can get the following -z attack , d</td></tr></table>

Comparison between CROWN-general(Zhang et al., 2018a)  and CRC. Note that both CROWN and CRC are computed on CPU. However, running time numbers are not directly comparable because our CRC implementation uses a batch of images while the CROWN implementation uses a single image at a time.

Table content: <table><tr><td>Network</td><td>Training</td><td>Standard Accuracy</td><td>Certified Robust Accuracy</td><td colspan="2">Certificate (mean) CROWN CRC</td><td>Time per image (seconds) CROWN CRC</td></tr><tr><td>2×[1024], sigmoid</td><td colspan="2">standard γ = 0.005 97.96% 98.37% γ = 0.01 98.08%</td><td>54.17% 82.68% 83.53%</td><td>0.28395 0.36125 0.32548</td><td colspan="2">0.48500 0.1818 0.83367 0.1599 0.84719 0.1732</td><td>0.1911 0.2229 0.2186</td></tr><tr><td/><td colspan="2">CRT, 0.01 98.57%</td><td colspan="2">95.59% 0.43061</td><td colspan="2">1.54673 0.1823</td><td>0.1910</td></tr><tr><td>3×[1024], sigmoid</td><td colspan="2">standard γ = 0.005 97.98% 98.37% γ = 0.01 97.71%</td><td>0.00% 88.66% 88.33%</td><td>0.24644 0.38030 0.39799</td><td colspan="2">0.06874 1.6356 0.99044 1.6220 1.07842 1.6342</td><td>0.5012 0.5319 0.5295</td></tr><tr><td/><td colspan="2">CRT, 0.01 97.23%</td><td colspan="2">94.99% 0.39603</td><td colspan="2">1.24100 1.5625</td><td>0.5013</td></tr><tr><td>4×[1024], sigmoid</td><td colspan="2">standard γ = 0.005 97.74% 98.39% γ = 0.01 97.41%</td><td>0.00% 88.95% 89.61%</td><td>0.19501 0.36863 0.40620</td><td colspan="2">0.00454 4.7814 0.91840 5.1667 1.05323 4.6296</td><td>0.8107 0.8567 0.8328</td></tr><tr><td/><td colspan="2">CRT, 0.01 97.83%</td><td colspan="2">93.41% 0.40327</td><td colspan="2">1.06208 4.1830</td><td>0.8088</td></tr><tr><td colspan="4">7 EXTENSIONS AND DISCUSSION</td><td/><td/></tr></table>

Comparison between CRT, PGD(Madry et al., 2018) and TRADES(Zhang et al., 2019).

Table content: <table><tr><td>Network 2×[1024], softplus 3×[1024], softplus</td><td>Training CRT, 0.01 CRT, 0.01 CRT, 0.03 CRT, 0.05</td><td>Standard Accuracy 98.68% 98.26% 97.82% 97.43%</td><td>Certified Robust Accuracy 69.79% 14.21% 50.72% 57.78%</td><td>Network 2×[1024], relu 3×[1024], relu</td><td>Standard Accuracy 89.33% 89.12%</td><td>Certified Robust Accuracy 44.29% 44.21%</td></tr><tr><td>4×[1024], softplus</td><td>CRT, 0.05 CRT, 0.01 CRT, 0.03</td><td>96.33% 97.80% 97.09%</td><td>44.44% 6.25% 29.64%</td><td>4×[1024], relu</td><td>90.17%</td><td>44.66%</td></tr></table>



Table content: <table><tr><td>Network</td><td>Training</td><td>Standard Accuracy</td><td>Empirical Robust Accuracy</td><td>Certified Robust Accuracy</td><td colspan="2">Certificate (mean) CROWN CRC</td></tr><tr><td>2 × [1024], sigmoid</td><td colspan="2">standard 46.23% γ = 0.01 45.42%</td><td>37.82% 38.17%</td><td>14.10% 26.50%</td><td>0.37219 0.40540</td><td>0.38173 0.55010</td></tr><tr><td>3 × [1024], sigmoid</td><td colspan="2">standard 48.57% γ = 0.01 50.31%</td><td>34.80% 39.87%</td><td>0.00% 18.28%</td><td>0.19127 0.24778</td><td>0.01404 0.37895</td></tr><tr><td>4 × [1024], sigmoid</td><td colspan="2">standard 46.04% γ = 0.01 48.28%</td><td>34.38% 40.10%</td><td>0.00% 21.07%</td><td>0.19340 0.29654</td><td>0.00191 0.40005</td></tr></table>

Results for CIFAR-10 dataset (only curvature regularization, no CRT training) Corollary 1. Both Curvature-based Certificate and Attack Optimizations are convex optimization problems within an l 2 ball of radius ρ around the input sample where curvature values are bounded.

Table content: <table/>

048, much better than the worst case bound -0.09623 and +0.09623 used in current global curvature bound. Comparison between Certified Robust accuracy and CRC for 2 layer sigmoid and tanh networks using global and local curvature bounds

Table content: <table><tr><td>Network</td><td>Training</td><td>Certified Accuracy (Global)</td><td>Certified Accuracy (Local)</td><td>CRC Certificate (Global)</td><td>CRC Certificate (Local)</td></tr><tr><td>2×[1024], sigmoid</td><td colspan="2">CRT, 0.0 CRT, 0.01 95.59% 95.04% CRT, 0.02 95.21%</td><td>95.31% 95.59% 95.21%</td><td>1.0011 1.5705 1.6720</td><td>1.3901 1.7262 1.7397</td></tr><tr><td>2×[1024], tanh</td><td colspan="2">CRT, 0.0 CRT, 0.01 95.00% 92.69% CRT, 0.02 94.77%</td><td>94.39% 95.00% 94.77%</td><td>0.8028 1.4832 1.5848</td><td>1.0511 1.6068 1.6592</td></tr></table>

function in x and d attack (η) is a convex optimization problem. (b) For every value of η, d attack (η) is a lower bound for p * attack . Thus d * attack = max -m≤η d attack (η) is a lower bound for p *

Table content: <table/>

Comparison between certified robust accuracy for different values of the regularization parameter γ for a single hidden layer convolutional neural network with softplus activation function G.2 RESULTS FOR TANH NETWORKS

Table content: <table><tr><td>Network</td><td>Training</td><td>Standard Accuracy</td><td>Empirical Robust Accuracy</td><td>Certified Robust Accuracy</td><td colspan="2">Certificate (mean) CROWN CRC</td></tr><tr><td>2×[1024], tanh</td><td colspan="2">PGD TRADES 98.63% 98.76% CRT, 0.01 98.52%</td><td>95.79% 96.20% 95.90%</td><td colspan="2">84.11% 93.72% 95.00% 0.37691 0.30833 0.40601</td><td>0.61340 0.86287 1.47016</td></tr><tr><td>3×[1024], tanh</td><td colspan="2">PGD TRADES 98.16% 98.78% CRT, 0.01 98.15%</td><td>94.92% 94.78% 95.00%</td><td colspan="2">0.00% 0.00% 94.16% 0.28004 0.12706 0.15875</td><td>0.03036 0.02983 1.14995</td></tr><tr><td>4×[1024], tanh</td><td colspan="2">PGD TRADES 97.08% 98.53% CRT, 0.01 97.24%</td><td>94.53% 92.85% 93.05%</td><td colspan="2">0.00% 0.00% 91.37% 0.33649 0.07439 0.11889</td><td>0.00140 0.00068 0.93890</td></tr></table>

Comparison between CRT, PGD(Madry et al., 2018) and TRADES(Zhang et al., 2019) for Tanh networks. CRC outperforms CROWN significantly for 2 layer networks and when trained with our regularizer for deeper networks. CRT outperforms TRADES and PGD giving higher certified accuracy.

Table content: <table/>

Comparison between CRT and Randomized Smoothing(Cohen et al., 2019). s denotes the standard deviation for smoothing. We use ρ = 0.5. For CRT, we use γ = 0.01

Table content: <table><tr><td>Network</td><td>Randomized Smoothing s = 0.25 s = 0.50 s = 1.0</td><td>CRT</td></tr><tr><td colspan="3">2 × [1024], sigmoid 93.75% 93.09% 88.91% 95.61%</td></tr><tr><td>2 × [1024], tanh</td><td colspan="2">94.61% 93.08% 82.26% 95.00%</td></tr><tr><td colspan="3">3 × [1024], sigmoid 94.00% 93.03% 86.58% 94.99%</td></tr><tr><td>3 × [1024], tanh</td><td colspan="2">93.69% 91.68% 80.55% 94.16%</td></tr><tr><td colspan="3">4 × [1024], sigmoid 93.68% 92.45% 84.99% 93.41%</td></tr><tr><td>4 × [1024], tanh</td><td colspan="2">93.57% 92.19% 83.90% 91.37%</td></tr></table>

Comparison between CRC and CROWN-general (CROWN-Ada for relu) for different targets. For CRT training, we use γ = 0.01. We compare CRC with CROWN-general for different targets for 150 correctly classified images. Runner-up means class with second highest logit is considered as adversarial class. Random means any random class other than the label is considered adversarial. Least means class with smallest logit is adversarial. For 2-layer networks, CRC outperforms CROWN-general significantly even without adversarial training. For deeper networks (3 and 4 layers), CRC works better on networks that are trained with curvature regularization.

Table content: <table><tr><td>Network</td><td>Training</td><td>Target</td><td colspan="2">Certificate (mean) CROWN CRC</td><td colspan="2">Time per Image (s) CROWN CRC</td></tr><tr><td/><td/><td colspan="2">runner-up 0.50110</td><td colspan="2">0.59166 0.1359</td><td>2.3492</td></tr><tr><td>2 × [1024], relu</td><td>standard</td><td>random</td><td>0.68506</td><td colspan="2">0.83080 0.2213</td><td>3.5942</td></tr><tr><td/><td/><td>least</td><td>0.86386</td><td colspan="2">1.04883 0.1904</td><td>3.0292</td></tr><tr><td/><td/><td colspan="2">runner-up 0.28395</td><td colspan="2">0.48500 0.1818</td><td>0.1911</td></tr><tr><td/><td>standard</td><td>random</td><td>0.38501</td><td colspan="2">0.69087 0.1870</td><td>0.1912</td></tr><tr><td>2 × [1024], sigmoid</td><td/><td colspan="2">least runner-up 0.43061 0.47639</td><td colspan="2">0.85526 0.1857 1.54673 0.1823</td><td>0.1920 0.1910</td></tr><tr><td/><td>CRT, 0.01</td><td>random</td><td>0.52847</td><td colspan="2">1.99918 0.1853</td><td>0.1911</td></tr><tr><td/><td/><td>least</td><td>0.62319</td><td colspan="2">2.41047 0.1873</td><td>0.1911</td></tr><tr><td/><td/><td colspan="2">runner-up 0.23928</td><td colspan="2">0.40047 0.1672</td><td>0.1973</td></tr><tr><td/><td>standard</td><td>random</td><td>0.31281</td><td colspan="2">0.52025 0.1680</td><td>0.1986</td></tr><tr><td>2 × [1024], tanh</td><td/><td colspan="2">least runner-up 0.37691 0.38964</td><td colspan="2">0.63081 0.1726 1.47016 0.1633</td><td>0.1993 0.1963</td></tr><tr><td/><td>CRT, 0.01</td><td>random</td><td>0.45896</td><td colspan="2">1.87571 0.1657</td><td>0.1982</td></tr><tr><td/><td/><td>least</td><td>0.52800</td><td colspan="2">2.21704 0.1697</td><td>0.1981</td></tr><tr><td/><td/><td colspan="2">runner-up 0.24644</td><td colspan="2">0.06874 1.6356</td><td>0.5012</td></tr><tr><td/><td>standard</td><td>random</td><td>0.29496</td><td colspan="2">0.08275 1.5871</td><td>0.5090</td></tr><tr><td>3 × [1024], sigmoid</td><td/><td colspan="2">least runner-up 0.39603 0.33436</td><td colspan="2">0.09771 1.6415 1.24100 1.5625</td><td>0.5056 0.5013</td></tr><tr><td/><td>CRT, 0.01</td><td>random</td><td>0.46808</td><td colspan="2">1.54622 1.6142</td><td>0.4974</td></tr><tr><td/><td/><td>least</td><td>0.51906</td><td colspan="2">1.75916 1.6054</td><td>0.4967</td></tr><tr><td/><td/><td colspan="2">runner-up 0.08174</td><td colspan="2">0.01169 1.4818</td><td>0.4908</td></tr><tr><td/><td>standard</td><td>random</td><td>0.10012</td><td colspan="2">0.01432 1.5906</td><td>0.4963</td></tr><tr><td>3 × [1024], tanh</td><td/><td colspan="2">least runner-up 0.28004 0.12132</td><td colspan="2">0.01757 1.5888 1.14995 1.4832</td><td>0.5076 0.4926</td></tr><tr><td/><td>CRT, 0.01</td><td>random</td><td>0.32942</td><td colspan="2">1.41032 1.5637</td><td>0.4957</td></tr><tr><td/><td/><td>least</td><td>0.38023</td><td colspan="2">1.65692 1.5626</td><td>0.4930</td></tr><tr><td/><td/><td colspan="2">runner-up 0.19501</td><td colspan="2">0.00454 4.7814</td><td>0.8107</td></tr><tr><td/><td>standard</td><td>random</td><td>0.21417</td><td colspan="2">0.00542 4.6313</td><td>0.8377</td></tr><tr><td>4 × [1024], sigmoid</td><td/><td colspan="2">least runner-up 0.40327 0.22706</td><td colspan="2">0.00609 4.7973 1.06208 4.1830</td><td>0.8313 0.8088</td></tr><tr><td/><td>CRT, 0.01</td><td>random</td><td>0.47038</td><td colspan="2">1.29095 4.3922</td><td>0.7333</td></tr><tr><td/><td/><td>least</td><td>0.52249</td><td colspan="2">1.49521 4.4676</td><td>0.7879</td></tr><tr><td/><td/><td colspan="2">runner-up 0.03554</td><td colspan="2">0.00028 5.7016</td><td>0.8836</td></tr><tr><td/><td>standard</td><td>random</td><td>0.04247</td><td colspan="2">0.00036 5.8379</td><td>0.8602</td></tr><tr><td>4 × [1024], tanh</td><td/><td colspan="2">least runner-up 0.33649 0.04895</td><td colspan="2">0.00044 5.8298 0.93890 3.8815</td><td>0.9045 0.8182</td></tr><tr><td/><td>CRT, 0.01</td><td>random</td><td>0.41617</td><td colspan="2">1.18956 4.0013</td><td>0.8215</td></tr><tr><td/><td/><td>least</td><td>0.47778</td><td colspan="2">1.41429 4.3856</td><td>0.8311</td></tr></table>

In this table, we measure the impact of increasing curvature regularization (γ) on accuracy, empirical robust accuracy, certified robust accuracy, CROWN-general and CRC when the network is trained without any adversarial training. We find that adding a very small amount of curvature regularization has a minimal impact on the accuracy but significantly increases CRC. Increase in CROWN certificate is not of similar magnitude. Somewhat surprisingly, we observe that even without any adversarial training, we can get nontrivial certified accuracies of 84.73%, 88.66%, 89.61% on 2,3,4 layer sigmoid networks respectively.

Table content: <table><tr><td>Network</td><td>γ</td><td>Standard Accuracy</td><td>Empirical Robust Accuracy</td><td>Certified Robust Accuracy</td><td colspan="2">Curvature bound (mean) K lb K ub</td></tr><tr><td/><td>0.0</td><td>98.77%</td><td>96.17%</td><td>95.04%</td><td>7.2031</td><td>72.0835</td></tr><tr><td>2×[1024], sigmoid</td><td colspan="2">0.005 98.82% 0.01 98.57% 0.02 98.59%</td><td>96.33% 96.28% 95.97%</td><td>95.61% 95.59% 95.22%</td><td>3.8411 2.8196 2.2114</td><td>8.2656 5.4873 3.7228</td></tr><tr><td/><td>0.03</td><td>98.30%</td><td>95.73%</td><td>94.94%</td><td>1.8501</td><td>2.9219</td></tr><tr><td/><td>0.0</td><td>98.65%</td><td>95.48%</td><td>92.69%</td><td>12.8434</td><td>107.5689</td></tr><tr><td>2×[1024], tanh</td><td colspan="2">0.005 98.71% 0.01 98.52% 0.02 98.35%</td><td>95.88% 95.90% 95.71%</td><td>94.76% 95.00% 94.77%</td><td>4.8116 3.4269 2.3943</td><td>10.1860 6.3529 4.1513</td></tr><tr><td/><td>0.03</td><td>98.29%</td><td>95.39%</td><td>94.54%</td><td>1.9860</td><td>3.933</td></tr><tr><td/><td>0.</td><td>98.52%</td><td>90.26%</td><td>0.00%</td><td>19.2131</td><td>3294.9070</td></tr><tr><td/><td colspan="2">0.005 98.41%</td><td>95.81%</td><td>94.91%</td><td>2.6249</td><td>13.4985</td></tr><tr><td>3×[1024], sigmoid</td><td>0.01 0.02 0.03</td><td>98.23% 97.99% 97.86%</td><td>95.70% 95.33% 94.98%</td><td>94.99% 94.64% 94.15%</td><td>1.9902 1.4903 1.2396</td><td>8.6654 5.4380 4.1409</td></tr><tr><td/><td>0.04</td><td>97.73%</td><td>94.60%</td><td>93.88%</td><td>1.0886</td><td>3.3354</td></tr><tr><td/><td>0.05</td><td>97.60%</td><td>94.45%</td><td>93.65%</td><td>0.9677</td><td>2.7839</td></tr><tr><td/><td>0.</td><td>98.19%</td><td>86.38%</td><td>0.00%</td><td colspan="2">133.7992 17767.5918</td></tr><tr><td/><td colspan="2">0.005 98.13%</td><td>94.56%</td><td>93.01%</td><td>3.2461</td><td>17.5500</td></tr><tr><td>3×[1024], tanh</td><td>0.01 0.02 0.03</td><td>98.15% 97.84% 97.70%</td><td>95.00% 94.79% 94.19%</td><td>94.16% 94.05% 93.42%</td><td>2.2347 1.6556 1.3546</td><td>10.8635 6.7072 5.0533</td></tr><tr><td/><td>0.04</td><td>97.57%</td><td>94.04%</td><td>92.95%</td><td>1.1621</td><td>4.0071</td></tr><tr><td/><td>0.05</td><td>97.31%</td><td>93.66%</td><td>92.65%</td><td>1.0354</td><td>3.3439</td></tr><tr><td/><td>0.</td><td>98.22%</td><td>83.04%</td><td>0.00%</td><td>86.9974</td><td>343582.3125</td></tr><tr><td/><td colspan="2">0.005 98.18%</td><td>95.02%</td><td>93.20%</td><td>2.1760</td><td>15.3358</td></tr><tr><td/><td>0.01</td><td>97.83%</td><td>94.65%</td><td>93.41%</td><td>1.6823</td><td>10.2289</td></tr><tr><td>4×[1024], sigmoid</td><td>0.02 0.03 0.04</td><td>97.33% 97.07% 96.70%</td><td>94.02% 93.52% 92.78%</td><td>92.94% 92.65% 91.95%</td><td>1.2089 1.0144 0.8840</td><td>6.5573 4.9576 3.9967</td></tr><tr><td/><td>0.05</td><td>96.38%</td><td>92.29%</td><td>91.33%</td><td>0.7890</td><td>3.4183</td></tr><tr><td/><td>0.06</td><td>96.29%</td><td>92.17%</td><td>91.11%</td><td>0.7128</td><td>3.0050</td></tr><tr><td/><td>0.07</td><td>96.08%</td><td>91.83%</td><td>90.67%</td><td>0.6614</td><td>2.6905</td></tr><tr><td/><td>0.</td><td>97.45%</td><td>75.18%</td><td>0.00%</td><td colspan="2">913.6984 37148156</td></tr><tr><td/><td colspan="2">0.005 97.48%</td><td>93.29%</td><td>89.98%</td><td>2.8690</td><td>18.8079</td></tr><tr><td/><td>0.01</td><td>97.24%</td><td>93.05%</td><td>91.37%</td><td>1.9114</td><td>12.2148</td></tr><tr><td>4×[1024], tanh</td><td>0.02 0.03 0.04</td><td>96.82% 96.27% 95.62%</td><td>92.65% 91.43% 90.69%</td><td>91.35% 90.09% 89.41%</td><td>1.3882 1.1643 0.9620</td><td>7.1771 5.1671 3.9061</td></tr><tr><td/><td>0.05</td><td>95.77%</td><td>90.69%</td><td>89.40%</td><td>0.9160</td><td>3.2909</td></tr><tr><td/><td>0.06</td><td>95.52%</td><td>90.00%</td><td>88.38%</td><td>0.8234</td><td>2.8808</td></tr><tr><td/><td>0.07</td><td>95.24%</td><td>89.51%</td><td>87.91%</td><td>0.7540</td><td>2.5635</td></tr></table>

Peer review (iclr2020):
This paper gives a global curvature bound for general neural networks, and used this bound for certifying robustness of neural networks. The basic techniques include writing the neural network verification problem as a constrained optimization problem (similar to [1]), and construct a Lagrangian dual for it, and discuss the situation where the dual problem is convex and can be solved efficiently. The global curvature bound is the crucial step to determine the condition that the optimization problem is convex. The benefits of the proposed methods include efficiency (bounds m and M only need to be solved for each target label pair, and the optimization process to get the certificate is fast), and tight certified accuracy for shallow networks. Also, it can be used during training to improve the tightness of CRC bounds. There are also a few weakness: the certificate unfortunately only works for non-ReLU networks, and only for L2 norm; for networks beyond 2 layers, CRC cannot outperform CROWN unless it is specially trained using the CRC loss. Questions and concerns: 1. For the curvature regularization training, some comparisons against other certified defenses are needed - PGD and TRADES are not good examples since they are not certified defense methods, and 0% certified accuracy are obtained. I understand many existing works like Wong et al., 2018 only work on ReLU networks, but I believe interval bound propagation (IBP) based methods [2] should be easily applicable to any monotonic activation functions. At least, the authors should compare verified accuracy at the same epsilon settings (pho=1.58) as in Wong et al., 2018. In page 22 (last page of appendix in arxiv version) of [3], you can find their L2 robustness training results. How does the certified accuracy of CRC compare to these results? 2. The paper mentions the Attack problem and proposes an attack algorithm (Algorithm 2), however I am not able to find any experiments on the attack. If the authors want to claim the contribution of curvature based attacks, some empirical results should be given, and at least compare it to a 200-step PGD baseline and see which one finds smaller adversarial examples. 3. The claim (in introduction and conclusions) that CRC is much faster than CROWN probably won*t hold under a fair comparison. I believe the gradient decent based certification algorithm (Algorithm 1) was computed efficiently on GPU, where the NN is defined. For a fair comparison, it should be also compared to a GPU implementation of CROWN ([4] provides such an implementation). In my experience, CROWN is an efficient algorithm that can be even used iteratively during training (as a certified defense [4]), so for the small networks used in this paper, it should compute bounds almost instantly on GPUs. I think it is better to revoke this claim, and in experiments the authors should clearly state that CROWN was computed on CPU so time is not comparable. 4. Also the claim *CRC outperforms CROWN*s certificate significantly* should be made clearer that it only holds for CRC trained models. According to Table 3, CRC is significantly worse than CROWN if the model is not CRC trained. Additionally, it should be made clear that the proposed method currently only applies to L2 norm setting. 5. (Minor) There are several duplicated references, e.g., on page 9, Cohen et al., (randomized smoothing) was cited twice as two different papers; the same is with Madry et al., on page 10, and Zhang et al. (CROWN) on page 11. These causes some confusions, e.g., on page 2, when talking about the bounded Lipschitz constant, I believe the correct citation for (Zhang et al., 2018b) should be another paper [5] from the same first author which is on bounding Lipschitz constant. Improvements and extensions: I think the proposed method can be greatly improved by using *local curvature*, where the curvature bounds m and M are computed within a local region near an input point. This is sufficient as long as our optimization does not escape this safe region, and the safe region can be naturally defined as the perturbation radius to be certified. The local curvature can be obtained by giving a tighter bounds on the second derivative of activation function, rather than considering just the worst case. Using CROWN, we can obtain pre-activation upper and lower bounds. These bounds can be used to bound the second derivative. For example, if a tanh neuron*s input is bounded by -0.1 and 0.1 (those bounds can be obtained efficiently by CROWN), the second derivative of tanh is bounded between -0.197356 and +0.197356, much better than the worst case bound -0.76981 and +0.76981 used in current global curvature bound. A similar technique to bound Jacobian matrix was used in [5]. Including some results for local curvature certificates will greatly increase the contribution of this paper, and make it a complete work. I strongly encourage the authors to do so, and feel free to discuss with me on any questions. Despite some concerns, the main contribution of giving global curvature bounds of neural networks is valid. Thus I vote for accepting this paper, however the authors should make sure to address all the concerns. [1] Salman, Hadi, et al. *A convex relaxation barrier to tight robust verification of neural networks.* arXiv preprint arXiv:1902.08722 (2019). [2] Gowal, Sven, et al. *On the effectiveness of interval bound propagation for training verifiably robust models.* arXiv preprint arXiv:1810.12715 (2018). [3] Wong, Eric, et al. *Scaling provable adversarial defenses.* Advances in Neural Information Processing Systems. 2018. https://arxiv.org/pdf/1805.12514.pdf [4] Zhang, Huan, et al. *Towards Stable and Efficient Training of Verifiably Robust Neural Networks.* arXiv preprint arXiv:1906.06316 (2019). [5] Zhang, H., Zhang, P., & Hsieh, C. J. Recurjac: An efficient recursive algorithm for bounding jacobian matrix of neural networks and its applications. arXiv preprint arXiv:1810.11783 (2018).