Year: 2021

Abstract:
Training convolutional neural networks (CNNs) with a strict Lipschitz constraint under the l2 norm is useful for provable adversarial robustness, interpretable gradients and stable training. While 1-Lipschitz CNNs can be designed by enforcing a 1-Lipschitz constraint on each layer, training such networks requires each layer to have an orthogonal Jacobian matrix (for all inputs) to prevent the gradients from vanishing during backpropagation. A layer with this property is said to be Gradient Norm Preserving (GNP). In this work, we introduce a procedure to certify the robustness of 1-Lipschitz CNNs by relaxing the orthogonalization of the last linear layer of the network that significantly advances the state of the art for both standard and provable robust accuracies on CIFAR-100 (gains of 4.80% and 4.71%, respectively). We further boost their robustness by introducing (i) a novel Gradient Norm preserving activation function called the Householder activation function (that includes every GroupSort activation) and (ii) a certificate regularization. On CIFAR-10, we achieve significant improvements over prior works in provable robust accuracy (5.81%) with only a minor drop in standard accuracy (−0.29%).

heading: 1 INTRODUCTION
text: Given a neural network f : Rd → Rk, the Lipschitz constant1 Lip(f) enforces an upper bound on how much the output is allowed to change in proportion to a change in the input. Previous work has demonstrated that a small Lipschitz constant is useful for improved adversarial robustness (Szegedy et al., 2014; Cissé et al., 2017), generalization bounds (Bartlett et al., 2017; Long & Sedghi, 2020), interpretable gradients (Tsipras et al., 2018) and Wasserstein distance estimation (Villani, 2008). Lip(f) also upper bounds the increase in the norm of gradient during backpropagation and can thus prevent gradient explosion during training, enabling us to train very deep networks (Xiao et al., 2018). While heuristic methods to enforce Lipschitz constraints (Miyato et al., 2018; Gulrajani et al., 2017) have achieved much practical success, they do not provably enforce a bound on Lip(f) globally and it remains challenging to achieve similar results when Lip(f) is provably bounded.
Using the property: Lip(g ◦ h) ≤ Lip(g) Lip(h), the Lipschitz constant of the neural network can be bounded by the product of the Lipschitz constant of all layers. While this allows us to construct 1-Lipschitz neural networks by constraining each layer to be 1-Lipschitz, Anil et al. (2018) identified a key difficulty with this approach. Because a 1-Lipschitz layer can only reduce the norm of gradient during backpropagation, backprop through each layer reduces the gradient norm, resulting in small gradient values for layers closer to the input, making training slow and difficult. To address this problem, they introduce Gradient Norm Preserving (GNP) architectures where each layer preserves the gradient norm during backpropagation. This involves constraining the Jacobian of each linear layer to be an orthogonal matrix and using a GNP activation function called GroupSort. GroupSort activation function (Anil et al., 2018) first separates the vector of preactivations z ∈ Rm into groups of pre-specified sizes, sorts each group in the descending order and then concatenates these sorted groups. When the group size is 2, the resulting activation function is called MaxMin.
For 1-Lipschitz CNNs, the robustness certificate for a sample x from class l is computed as Mf (x)/ √ 2 whereMf (x) = fl(x)−maxi6=l fi(x). Naturally, larger values of fl(x) and smaller values of maxj 6=l fj(x) will lead to larger certificates. This requires the last weight matrix in the
1Unless specified, we assume the Lipschitz constant under the l2 norm in this work.
network, denoted by W ∈ Rk×m (k is the number of classes, m is the dimension of the penultimate layer, m > k), to enforce the following constraints throughout training:
∀j, ‖Wj,:‖2 = 1, i 6= j, Wj,: ⊥Wi,: where Wi,: denotes the ith row of W. Now suppose that for some input image with label l, we want to update W to increase the logit for the lth class. Since ‖Wl,:‖2 is constrained to be 1, the logit can only be increased by changing the direction of the vector Wl,:. Because the other rows {Wi,:, i 6= l} are constrained to be orthogonal to Wl,:, this further requires an update for all the rows of W. Thus during training, any update made to learn some class must necessitate the forgetting of information relevant for the other classes. This can be particularly problematic when the number of classes k and thus the number of orthogonality constraints per row (i.e., k − 1) is large (such as in CIFAR-100). To address this limitation, we propose to keep the last weight layer of the network unchanged. But then the resulting function is no longer 1-Lipschitz and the certificateMf (x)/ √ 2 is not valid. Thus, we introduce a new certification procedure that does not require the last weight layer of the network W to be orthogonal. Our certificate is then given by the following equation:
min i 6=l fl(x)− fi(x) ‖Wl,: −Wi,:‖2
However, a limitation of using the above certificate is that because the weight layers are completely unconstrained, larger norms of rows (i.e., ‖Wi,:‖) can result in larger values of ‖Wl,: −Wi,:‖2 and thus smaller certificate values. To address this limitation, we normalize all rows to be of unit norm before computing the logits. While this still requires all the rows to be of unit norm, their directions are now allowed to change freely thus preventing the need to update other rows and forgetting of learned information. We show that this provides significant improvements when the number of classes is large. We call this procedure Last Layer Normalization (abbreviated as LLN). On CIFAR-100, this significantly improves both the standard (> 3%) and provable robust accuracy (> 4% at ρ = 36/255) across multiple 1-Lipschitz CNN architectures (Table 1). Here, ρ is the l2 attack radius.
Another limitation of existing 1-Lipschitz CNNs (Li et al., 2019b; Trockman & Kolter, 2021; Singla & Feizi, 2021) is that their robustness guarantees do not scale properly with the l2 radius ρ. For example, the provable robust accuracy of (Singla & Feizi, 2021) drops ∼ 30% at ρ = 108/255 compared to 36/255 on CIFAR-10 (Table 2). To address this limitation, we introduce a certificate regularization denoted by CR (Section 5) that when used along with the Householder activation results in significantly improved provable robust accuracy at larger radius values with minimal loss in standard accuracy. On the CIFAR-10 dataset, we achieve significant improvements in the provable robust accuracy for large ρ = 108/255 (min gain of +4.96%) across different architectures with minimal loss in the standard accuracy (max drop of −0.56%). Results are in Table 2. Additionally, we characterize the MaxMin activation function as a special case of the more general Householder (HH) activations. Recall that given z ∈ Rm, the HH transformation is a linear function reflecting z about the hyperplane vTx = 0 (‖v‖2 = 1), given by (I− 2vvT )z where I− 2vvT is orthogonal because ‖v‖2 = 1. The Householder activation function σv is defined below:
σv(z) =
{ z, vT z > 0,
(I− 2vvT )z, vT z ≤ 0. (1)
First, note that since z = (I − 2vvT )z along vT z = 0, σv is continuous. Moreover, the Jacobian ∇z σv is either I or I− 2vvT (both orthogonal) implying σv is GNP. Since these properties hold ∀ v : ‖v‖2 = 1, v can be learned during the training. In fact, we prove that any GNP piecewise linear function that changes from Q1z to Q2z (Q1,Q2 are square orthogonal matrices) along vT z = 0 must satisfy Q2 = Q1 ( I− 2vvT ) to be continuous (Theorem 1). Thus, this characterization proves that every GroupSort activation is a special case of the more general Householder activation function (example in Figure 1, discussion in Section 6).
In summary, in this paper, we make the following contributions:
• We introduce a certification procedure without orthogonalizing the last linear layer called Last Layer Normalization. This procedure significantly enhances the standard and provable robust accuracy when the number of classes is large. Using the LipConvnet-15 network on CIFAR-100, our modification achieves a gain of +4.71% in provable robust accuracy (at ρ = 36/255) with a gain of +4.80% in standard accuracy (Table 1).
• We introduce a Certificate Regularizer that significantly advances the provable robust accuracy with a small reduction in standard accuracy. Using LipConvnet-15 network on CIFAR-10, we achieve +5.81% improvement in provable robust accuracy (at ρ = 108/255) with only a −0.29% drop in standard accuracy over the existing methods (Table 2).
• We introduce a class of piecewise linear GNP activation functions called Householder or HH activations. We show that the MaxMin activation is a special case of the HH activation for certain settings. We prove that Householder transformations are necessary for any GNP piecewise linear function to be continuous (Theorem 1).

heading: 2 RELATED WORK
text: Provably Lipschitz convolutional neural networks: The class of fully connected neural networks (FCNs) which are Gradient Norm Preserving (GNP) and 1-Lipschitz were first introduced by Anil et al. (2018). They orthogonalize weight matrices and use GroupSort as the activation function to design each layer to be GNP. While there have been numerous works on enforcing Lipschitz constraints on convolution layers (Cissé et al., 2017; Tsuzuku et al., 2018; Qian & Wegman, 2019; Gouk et al., 2020; Sedghi et al., 2019), they either enforce loose Lipschitz bounds or are not scalable to large networks. To ensure that the Lipschitz constraint on convolutional layers is tight, multiple recent works try to construct convolution layers with an orthogonal Jacobian matrix (Li et al., 2019b; Trockman & Kolter, 2021; Singla & Feizi, 2021). These approaches avoid the aforementioned issues and allow the training of large, provably 1-Lipschitz CNNs while achieving impressive results.
Provable defenses against adversarial examples: A provably robust classifier is one for which we can guarantee that the classifier’s prediction remains constant within some region around the input. Most of the existing methods for provable robustness either bound the Lipschitz constant of the neural network or the individual layers (Weng et al., 2018; Zhang et al., 2019; 2018; Wong et al., 2018; Wong & Kolter, 2018; Raghunathan et al., 2018; Croce et al., 2019; Singh et al., 2018; Singla & Feizi, 2020; Zhang et al., 2021). However, these methods do not scale to large and practical networks on the ImageNet dataset (Deng et al., 2009). To scale to such large networks, randomized smoothing (Liu et al., 2018; Cao & Gong, 2017; Lécuyer et al., 2018; Li et al., 2019a; Cohen et al., 2019; Salman et al., 2019; Levine et al., 2019; Kumar et al., 2020a;b) has been proposed as a probabilistically certified defense. However, certifying robustness with high probability requires generating a large number of noisy samples leading to high inference-time computation. In contrast, the defense we propose is deterministic and hence not comparable to randomized smoothing. While Levine & Feizi (2021) provide deterministic robustness certificates using randomized smoothing, their certificates are in the l1 norm and not directly applicable for the l2 threat model studied in this work. We discuss the differences between l1 and l2 certificates in Appendix Section C.

heading: 3 PROBLEM SETUP AND NOTATION
text: For a vector v, vj denotes its jth element. For a matrix A, Aj,: and A:,k denote the jth row and kth column respectively. Both Aj,: and A:,k are assumed to be column vectors (thus Aj,: is the transpose of jth row of A). Aj,k denotes the element in jth row and kth column of A. A:j,:k denotes the matrix containing the first j rows and k columns of A. The same rules are directly extended to higher order tensors. I denotes the identity matrix, R to denote the field of real numbers. For θ ∈ R, J+(θ) and J−(θ) denote the orthogonal matrices with determinants +1 and −1 defined as follows:
J+(θ) = [ cos θ sin θ − sin θ cos θ ] J−(θ) = [ cos θ sin θ sin θ − cos θ ] (2)
We construct a 1-Lipschitz neural network, f : Rd → Rk (d is the input dimension, k is the number of classes) by composing 1-Lipschitz convolution layers and GNP activation functions. To certify robustness for some input x with prediction l, we first define the margin of prediction: Mf (x) = max (0, fl(x)−maxi 6=l fi(x)) where fi(x) is the logit for class i and l is the correct label. Using Theorem 7 in Li et al. (2019b), we can derive the robustness certificate (in the l2 norm)
asMf (x)/ √ 2. Thus, the l2 distance of x to the decision boundary is lower bounded byMf (x)/ √ 2:
min i 6=l min fi(x∗)=fl(x∗)
‖x∗ − x‖2 ≥ Mf (x)√
2 (3)
We often use the abbreviation fi − fj : RD → R to denote the function so that:
(fi − fj) (x) = fi(x)− fj(x), ∀ x ∈ RD
Our goal is to train the neural network f to achieve the maximum possible provably robust accuracy while also simultaneously improving (or maintaining) standard accuracy.

heading: 4 LAST LAYER NORMALIZATION
text: To ensure that the network is 1-Lipschitz so that the certificate in equation (3) is valid, existing 1-Lipschitz neural networks require the weight matrices of all the linear layers of the network to be orthogonal. For the weight matrix in the last layer of the network (that maps the penultimate layer neurons to the logits), W ∈ Rk×m (k is the number of classes, m is the dimension of the penultimate layer, m > k), this enforces the following constraints on each row Wi,: ∈ Rm:
∀j, ‖Wj,:‖2 = 1, ∀i 6= j, Wj,: ⊥Wi,: (4) Now suppose that for some input x with label l, we want to update W to increase the logit for the lth class. Since ‖Wl,:‖2 is constrained to be 1, the gradient update can only change the direction of the vector Wl,:. But now, because the other rows {Wi,:, i 6= l} are constrained to be orthogonal to Wl,:, this further requires an update for all the rows of W. This has the negative effect that during training, any update made to learn some class must necessitate the forgetting of information relevant for the other classes. This can be particularly problematic when the number of classes k is large (such as in CIFAR-100) and thus the number of orthogonality constraints per row i.e. k − 1 is large. To address this limitation, first observe that the neural network from the input layer to the penultimate layer (i.e excluding the last linear layer) is 1-Lipschitz. Let g : Rd → Rm be this function so that f(x) = Wg(x) + b. This equation suggests that even if W is not orthogonal, the Lipschitz constant of the function fl − fi, can be computed by multiplying the Lipschitz constant of g (which is 1) and that of (Wl,: −Wi,:) (which is ‖Wl,: −Wi,:‖2). The robustness certificate can then be computed asMf (x)/‖Wl,: −Wi,:‖2. This procedure leads to the following proposition: Proposition 1. Given 1-Lipschitz continuous function g : Rd → Rm and W ∈ Rk×m, b ∈ Rk (k is the number of classes), construct a new function f : Rd → Rk defined as: f(x) = Wg(x) + b. Let fl(x) > maxi 6=l fi(x). The robustness certificate (under the l2 norm) is given by:
min i 6=l min fl(x∗)=fi(x∗) ‖x∗ − x‖2 ≥ min i 6=l fl(x)− fi(x) ‖Wl,: −Wi,:‖2
Proof. We proceed by computing the Lipschitz constant of the function fl − fi. The gradient of the function: fl − fi at x can be computed using the chain rule:
∇x (fl − fi) = (Wl,: −Wi,:)T ∇xg Since g is given to be 1-Lipschitz, the Lipschitz constant of fl − fi can be computed using the above equation as follows:
‖∇x (fl − fi) ‖2 ≤ ‖ (Wl,: −Wi,:)T (∇xg) ‖2 ‖∇x (fl − fi) ‖2 ≤ ‖Wl,: −Wi,:‖2‖∇xg‖2 ≤ ‖Wl,: −Wi,:‖2
Thus, the distance of x to the decision boundary fl − fi = 0, is lower bounded by:
min fl(x∗)=fi(x∗) ‖x∗ − x‖2 ≥ fl(x)− fi(x) ‖Wl,: −Wi,:‖2
Thus, the distance to decision boundary across all classes i 6= l is lower bounded by:
min i 6=l min fl(x∗)=fi(x∗) ‖x∗ − x‖2 ≥ min i 6=l fl(x)− fi(x) ‖Wl,: −Wi,:‖2
However, in our experiments, we found that using this procedure directly i.e without any constraint on the weight matrix W often results in large norms of row vectors ‖Wi,:‖2, and thus large values of ‖Wl,: −Wi,:‖2 and smaller certificates (Theorem 1). To address this problem, we normalize all rows of the matrix to be of unit norm before computing the logits so that for the input x, the logit gi(x) can be computed as follows:
gi(x) = (Wi,:)
T f(x)
‖Wi,:‖2 + bi
The robustness certificate can then be computed as follows:
min i6=l gl(x)− gi(x) ‖W(n)l,: −W (n) i,: ‖2 , where ∀j, W(n)j,: = Wj,: ‖Wj,:‖2
, gj(x) = ( W (n) j,: )T f(x) + bj
While each row Wi,: is still constrained to be of unit norm, unlike with orthogonality constraints (equation (4)), their directions are allowed to change freely. This provides significant improvements when the number of classes is large. We call this procedure Last Layer Normalization (LLN).

heading: 5 CERTIFICATE REGULARIZATION
text: A limitation of using cross entropy loss for training 1-Lipschitz CNNs is that it is not explicitly designed to maximize the marginMf (x) and thus, the robustness certificate. That is, once the cross entropy loss becomes small, the gradients will no longer try to further increase the marginMf (x) even though the network may have the capacity to learn bigger margins. To address this limitation, we can simply subtract the certificate i.e −Mf (x)/ √
2 from the usual cross entropy loss function during training. Observe that we subtract the certificate because we want to maximize the certificate values while minimizing the cross entropy loss. However, in our experiments we found that this regularization term excessively penalizes the network for the misclassified examples and as a result, the certificate values for the correctly classified inputs are not large. Thus, we propose to use the following regularized loss function during training:
min Ω
E(x,l)∼D [ ` (fΩ(x), l)− γ relu ( Mf (x)√
2
)] (5)
In the above equation, fΩ denotes the 1-Lipschitz neural network parametrized by Ω, fΩ(x) denotes the logits for the input x, ` (fΩ(x), l) is the cross entropy loss for input x with label l and γ > 0 is the regularization coefficient for maximizing the certificate. We have the minus sign in front of the regularization term γ relu(Mf (x)/ √ 2) because we want to maximize the certificate while minimizing the cross entropy loss. For wrongly classified inputs,Mf (x)/ √ 2 < 0 =⇒ relu(Mf (x)/ √
2) = 0. This ensures that the optimization tries to increase the certificates only for the correctly classified inputs. We call the above mentioned procedure Certificate Regularization (abbreviated as CR).

heading: 6 HOUSEHOLDER ACTIVATION FUNCTIONS
text: Recall that given z ∈ Rm, the Householder (HH) transformation reflects z about the hyperplane vTx = 0 where ‖v‖2 = 1. The linear transformation is given by the equation (I− 2vvT )z where I− 2vvT is orthogonal because ‖v‖2 = 1. Now, consider the nonlinear function σv defined below: Definition 1. (Householder Activation of Order 1) The activation function σv : Rm → Rm, applied on z ∈ Rm, is called the m-dimensional Householder Activation of Order 1:
σv(z) =
{ z, vT z > 0,
(I− 2vvT )z, vT z ≤ 0. (6)
Since σv is linear when vT z > 0 or vT z < 0, it is also continuous in both cases. At the hyperplane separating the two cases (i.e., vT z = 0) we have: (I − 2vvT )z = z − 2(vT z)v = z (both linear functions are equal). Thus, σv is continuous ∀ z ∈ Rm. Moreover, the Jacobian is either I or I− 2vvT which are both square orthogonal matrices. Thus, σv is also GNP and 1-Lipschitz. Since these properties hold for all v satisfying ‖v‖2 = 1, v can be made a learnable parameter.
While the above arguments suggest that HH transformations are sufficient to ensure such functions are continuous, we also prove that they are necessary. That is, we prove that if a GNP piecewise linear function g : Rm → Rm transitions between different linear functions Q1z and Q2z (in an open set S ⊂ Rm) along a hyperplane vT z = 0 (where ‖v‖2 = 1), then g is continuous in S if and only if Q2 = Q1(I− 2vvT ). This theoretical result provides a general principle for designing piecewise linear GNP activation functions. The formal result is stated in the following Theorem: Theorem 1. Given an open set S ⊂ Rm, orthogonal square matrices Q1 6= Q2, and vector v ∈ Rm (‖v‖2 = 1) such that S ∩ {z : vT z = 0} 6= ∅, the function g defined as follows:
g(z) = { Q1z, z ∈ S,vT z > 0, Q2z, z ∈ S,vT z ≤ 0
(7)
is continuous in S if and only if Q2 = Q1(I− 2vvT ).
Proof of Theorem 1 is in Appendix A.1. Note that since the matrix I− 2vvT has determinant −1, the above theorem necessitates that det(Q1) = −det(Q2) i.e the determinant of the Jacobian must change sign whenever the Jacobian of a piecewise linear GNP activation function changes.
Recall that for the MaxMin activation function, MaxMin(z1, z2) = (z1, z2) if z1 > z2 and (z2, z1) otherwise. Thus, the Jacobian of MaxMin for z1 > z2 case is I = J+(0) while for z1 ≤ z2 is J−(π/2). Using Theorem 1, we can easily prove that MaxMin is a special case of the more general Householder activation functions where the Jacobian J−(π/2) is replaced with J−(θ) and the conditions z1 > z2 are replaced with z1 sin(θ/2) > z2 cos(θ/2) (similarly for ≤). The construction of Householder activations in 2 dimensions, denoted by σθ, is given in the following corollary: Corollary 1. The function σθ : R2 → R2 defined as
σθ (z1, z2) =  [ 1 0 0 1 ] [ z1 z2 ] if z1 sin (θ/2)− z2 cos (θ/2) > 0[ cos θ sin θ
sin θ − cos θ ] [ z1 z2 ] if z1 sin (θ/2)− z2 cos (θ/2) ≤ 0
(8)
is continuous and is called 2D Householder Activation of Order 1.
The two cases are demonstrated in Figure 1a and Figure 1b, respectively. Since σθ is continuous, GNP and 1-Lipschitz ∀ θ ∈ R, θ is a learnable parameter. For θ = π/2 in equation (8), σθ is equivalent to MaxMin. Thus, σθ is at least as expressive as MaxMin.
A major limitation of using σv (equation (6)) directly is that it has only 2 linear regions and is thus limited in its expressive power. In contrast, MaxMin first divides the preactivation z ∈ Rm (assuming m is divisible by 2) into m/2 groups of size 2 each. Since each group has 2 linear regions, we get 2m/2 linear regions from the m/2 groups. Thus, to increase the expressive power, we similarly divide
z into m/2 groups of size 2 each and apply the 2-dimensional Householder activation function of Order 1 (σθ) to each group resulting in 2m/2 linear regions (same as MaxMin).
To apply σθ to the output of a convolution layer z ∈ Rm×n×n (m is the number of channels and n×n is the spatial size), we first split z into 2 tensors along the channel dimension giving the tensors: z:m/2,:,: and zm/2:,:,:. Each of these tensors is of size m/2× n× n giving n2m/2 groups. We use the same θ for each pair of channels (irrespective of spatial location) resulting in m/2 learnable parameters. We initialize each θ = π/2 so that σθ is equivalent to MaxMin at initialization.

heading: 7 EXPERIMENTS
text: Our goal is to evaluate the effectiveness of the three changes proposed in this work: (a) Last Layer Normalization, (b) Certificate regularization and (c) Householder activation functions. We perform experiments under the setting of provably robust image classification on CIFAR-10 and CIFAR-100 datasets using the same 1-Lipschitz CNN architectures used by Singla & Feizi (2021) (LipConvnet-5, 10, 15, . . . , 40) due to their superior performance over the prior works. We compare with the three orthogonal convolution layers in the literature: SOC (Singla & Feizi, 2021), BCOP (Li et al., 2019b) and Cayley (Trockman & Kolter, 2021) using MaxMin as the activation function.
We use SOC with MaxMin as the primary baseline for comparison in the maintext due to their superior performance over the prior works (BCOP, Cayley). Results using BCOP and Cayley convolutions are given in Appendix Sections H and I for completeness. We use the same implementations for these convolution layers as given in their respective github repositories. We compare the provable robust accuracy using 3 different values of the l2 perturbation radius ρ = 36/255, 72/255, 108/255. In both Tables 1 and 2, for all networks, we use SOC as the convolution layer. Using SOC, we achieve the same bound on the approximation error of an orthogonal Jacobian as achieved in Singla & Feizi (2021) i.e. 2.42× 10−6 across all convolution layers in all networks. Thus, even for the network with 40 layers, this results in maximum Lipschitz constant of (1 + 2.42× 10−6)40 = 1.000097 ≤ 1.0001. Thus, the Lipschitz constant across all our networks is bounded by 1.0001. The symbol HH (in Tables 1, 2) is for the 2D Householder Activation of order 1 or σθ (defined in equation (8)).
All experiments were performed using 1 NVIDIA GeForce RTX 2080 Ti GPU. All networks were trained for 200 epochs with initial learning rate of 0.1, dropped by a factor of 0.1 after 100 and 150 epochs. For Certificate Regularization (or CR), we set the parameter γ = 0.5.

heading: 7.1 RESULTS ON CIFAR-100
text: In Table 1, for each architecture, the row "SOC + MaxMin" uses the MaxMin activation, "+ LLN" adds Last Layer Normalization (uses MaxMin), "+ HH" replaces MaxMin with σθ (also uses LLN) , "+ CR" also adds Certificate Regularization with γ = 0.1 (uses both σθ and LLN). The column, "Increase (Standard)" denotes the increase in standard accuracy relative to "SOC + MaxMin".
By adding LLN (the row "+ LLN"), we observe gains in standard (min gain of 1.10%) and provable robust accuracy (min gain of 1.71% at ρ = 36/255) across all the LipConvnet architectures (gains relative to "SOC + MaxMin"). These gains are smallest for the LipConvnet-40 network with the maximum depth. However, replacing MaxMin with the σθ activation further improves the standard (min gain of 3.65%) and provable robust accuracy (min gain of 4.46% at ρ = 36/255) across all networks (again relative to "SOC + MaxMin"). We observe that replacing MaxMin with σθ significantly improves the performance of the deeper LipConvnet-35, 40 networks.
Adding CR further improves the provable robust accuracy while only slightly reducing the standard accuracy. Because LLN significantly improves the standard accuracy, we compare the standard accuracy numbers between rows "+ CR" and "+ LLN" to evaluate the drop due to CR. We observe a small drop in standard accuracy (−0.04%,−0.11%) only for LipConvnet-5 and LipConvnet-15 networks. For the other networks, the standard accuracy actually increases.

heading: 7.2 RESULTS ON CIFAR-10
text: In Table 2, for each architecture, the row "SOC + MaxMin" uses the MaxMin activation, the row "+ HH" uses σθ activation (replacing MaxMin) and the row "+ CR" also adds Certificate Regularization
with γ = 0.1 (again using σθ as the activation). Due to the small number of classes in CIFAR-10, we do not observe significant gains using Last Layer Normalization or LLN (Appendix Table 10). Thus, we do not include LLN for any of the results in Table 2. The column, "Increase (108/255)" denotes the increase in provable robust accuracy with ρ = 108/255 relative to "SOC + MaxMin".
For LipConvnet-25, 30, 35, 40 architectures, we observe gains in both the standard and provable robust accuracy by replacing MaxMin with the HH activation (i.e σθ). Similar to what we observe for CIFAR-100 in Table 1, the gains in provable robust accuracy (ρ = 108/255) are significantly higher for deeper networks: LipConvnet-35 (3.65%) and LipConvnet-40 (4.35%) with decent gains in standard accuracy (1.71 and 1.61% respectively).
Again similar to CIFAR-100, adding CR further boosts the provable robust accuracy while slightly reducing the standard accuracy. Comparing "+ CR" with "SOC + MaxMin", we observe small drops in standard accuracy for LipConvnet-5, 10, . . . , 30 networks (max. drop of −0.56%), and gains for LipConvnet-35 (+0.52%) and LipConvnet-40 (+0.96%). For provable robust accuracy (ρ = 108/255), we observe very significant gains of > 4.96% for all networks and > 8% for the deeper LipConvnet-35, 40 networks.

heading: 8 CONCLUSION
text: In this work, we introduce a procedure to certify robustness of 1-Lipschitz convolutional neural networks without orthogonalizing of the last linear layer of the network. Additionally, we introduce a certificate regularization that significantly improves the provable robust accuracy for these models at higher l2 radii. Finally, we introduce a class of GNP activation functions called Householder (or HH) activations and prove that the Jacobian of any Gradient Norm Preserving (GNP) piecewise linear function is only allowed to change via Householder transformations for the function to be continuous which provides a general principle for designing piecewise linear GNP functions. These ideas lead to improved deterministic `2 robustness certificates on CIFAR-10 and CIFAR-100.

heading: REPRODUCIBILITY
text: The code for the experiments is provided in the supplementary materials.

heading: ETHICS STATEMENT
text: This paper introduces novel set of techniques for improving the provable robustness of 1-Lipschitz CNNs. We do not see any foreseeable negative consequences associated with our work.

heading: A PROOFS

heading: A.1 PROOF OF THEOREM 1
text: Proof. We first prove that if Q2 = (I− 2vvT )Q1, then the function g is continuous. First, observe that for vT z > 0, g(z) = Q1z which is continuous. Similarly, for vT z < 0, g(z) = Q2z which is again continuous. This proves that the function g is continuous when vT z > 0 or vT z < 0. Thus, to prove continuity ∀ z ∈ S, we must prove that:
Q1z = Q2z ∀ z : vT z = 0 (9) Since Q2 = Q1(I− 2vvT ), we have:
Q2 −Q1 = −2Q1vvT (Q2 −Q1) z = −2 ( Q1vv T ) z = −2Q1v ( vT z ) The above equation directly proves (9).
Now, we prove the other direction i.e if g is continuous in S then, Q2 = Q1(I− 2vvT ). Since g is continuous for all z : vT z = 0, we have:
Q2z = Q1z ∀ z : vT z = 0 (Q2 −Q1) z = 0 ∀ z : vT z = 0
Since z ∈ Rm, we know that the set of vectors: {z : vT z = 0} spans a m− 1 dimensional subspace. Thus, the null space of Q2 −Q1 is of size m− 1. This in turn implies that Q2 −Q1 is a rank one matrix given by the following equation:
Q2 −Q1 = uvT (10) where the vector u remains to be determined. Since Q1 and Q2 are orthogonal matrices, we have the following set of equations:
QT2 Q2 = ( Q1 + uv T )T ( Q1 + uv T )
(11)
Q2Q T 2 = ( Q1 + uv T ) ( Q1 + uv T )T
(12) We first simplify equation (11):
QT2 Q2 = ( QT1 + vu T ) ( Q1 + uv T )
I = I + v ( QT1 u )T + ( QT1 u ) vT + ( uTu ) vvT
0 = v ( QT1 u )T + ( QT1 u ) vT + ( uTu ) vvT
− ( uTu ) vvT = v ( uTQ1 ) + ( QT1 u ) vT
Right multiplying both sides by v and using ‖v‖2 = 1, we get: − ( uTu ) v = ( uTQ1v ) v + QT1 u
QT1 u = − ( uTu + uTQ1v ) v = λv
u = λQ1v, where λ = − ( uTu + uTQ1v ) (13)
Substituting u using equation (13) in equation (10), we get:
Q2 −Q1 = λQ1vvT Q2 = Q1 ( I + λvvT ) Since QT2 Q2 = I, we have:
QT2 Q2 = ( Q1 ( I + λvvT ))T Q1 ( I + λvvT ) QT2 Q2 = ( I + λvvT ) QT1 Q1 ( I + λvvT
) I = ( I + λvvT ) ( I + λvvT
) I = I + 2λvvT + λ2vvT
=⇒ λ = 0 or λ = −2
Since λ = 0 would imply Q1 = Q2 which is not allowed by the assumption of the Theorem that Q1 6= Q2. λ = −2 is the only possibility allowed. This proves the other direction i.e:
Q2 = Q1 ( I− 2vvT )

heading: A.2 PROOF OF COROLLARY 1
text: Proof. Subsitute Q1,v as follows in Theorem 1.
Q1 = [ 1 0 0 1 ] v = [ + sin (θ/2) − cos (θ/2)
] Q2 = I− 2vvT
Q2 = [ 1 0 0 1 ] − 2 [ sin (θ/2) − cos (θ/2) ] [sin (θ/2) − cos (θ/2)]
Q2 = [ 1 0 0 1 ] − 2 [ sin2 (θ/2) − sin (θ/2) cos (θ/2) − sin (θ/2) cos (θ/2) cos2 (θ/2) ] Q2 = [ 1− 2 sin2 (θ/2) 2 sin (θ/2) cos (θ/2)
2 sin (θ/2) cos (θ/2) 1 + 2 cos2 (θ/2) ] Q2 = [ cos(θ) sin(θ) sin(θ) − cos(θ)
] Theorem 1 then directly implies the corollary.

heading: A.3 PROOF OF THEOREM 2
text: Theorem 2. Given: 0 ≤ θ0 < θ1 · · · < θ2n = 2π + θ0 such that ∑n−1 i=0 (θ2i+1 − θ2i) = π and
αi = 2 ∑i j=0 θi−j(−1)j , The function σΘ : R2 → R2 is continuous, GNP and 1-Lipschitz where Θ = [θ0, θ1, . . . , θ2n−1] (also called 2D Householder Activation of order n):
σΘ(z1, z2) =
[ cosαi sinαi
(−1)i sinαi (−1)i+1 cosαi ] [ z1 z2 ] θi ≤ ϕ < θi+1 (14)
where ϕ ∈ [θ0, θ2n = 2π + θ0) and cos(ϕ) = z1/ √ z21 + z 2 2 , sin(ϕ) = z2/ √ z21 + z 2 2 .
Proof. We are given the following:
n−1∑ i=0 (θ2i+1 − θ2i) = π, αi = 2 i∑ j=0 θi−j(−1)j (15)
Note that by definition (equation (14)), the function is linear for θi < ϕ < θi+1 and hence continuous. Furthermore, since ϕ ∈ [θ0, θ2n), we proceed to prove continuity for the following two cases: Case 1: θi − < ϕ < θi + , > 0, i ≥ 1 Case 2: θ0 < ϕ < θ0 + , θ2n − < ϕ < θ2n, > 0

heading: Proof for Case 1:
text: Using equation (14), we know that σΘ realizes different linear functions for θi − < ϕ < θi and θi < ϕ < θi + . Thus, for σΘ to be continuous, we require that the two linear functions be the same at the boundary i.e ϕ = θi.
We first write the input (z1, z2) in terms of shifted polar coordinates i.e: (r cos(ϕ), r sin(ϕ)) where r = √ z21 + z 2 2 and cosϕ = z1/ √ z21 + z 2 2 , sinϕ = z2/ √ z21 + z 2 2 , ϕ ∈ [θ0, θ0 + 2π) Thus, the function for θi − < ϕ < θi is given by:
σΘ (r cosϕ, r sinϕ) =
[ cosαi−1 sinαi−1
(−1)i−1 sinαi−1 (−1)i cosαi−1 ] [ r cosϕ r sinϕ ] (16)
Similarly, the function for θi < ϕ < θi + is given by:
σΘ (r cosϕ, r sinϕ) =
[ cosαi sinαi
(−1)i sinαi (−1)i+1 cosαi
] [ r cosϕ r sinϕ ] (17)
The difference in the function values at the boundary i.e ϕ = θi, obtained by subtracting equations (17) and (16) is given as follows:[
cosαi sinαi (−1)i sinαi (−1)i+1 cosαi ] [ r cos θi r sin θi ] − [ cosαi−1 sinαi−1 (−1)i−1 sinαi−1 (−1)i cosαi−1 ] [ r cos θi r sin θi ] = r [ cosαi − cosαi−1 sinαi − sinαi−1
(−1)i sinαi − (−1)i−1 sinαi−1 (−1)i+1 cosαi − (−1)i cosαi−1 ] [ cos θi sin θi ] = r [ cosαi − cosαi−1 sinαi − sinαi−1
(−1)i(sinαi + sinαi−1) (−1)i+1(cosαi + cosαi−1) ] [ cos θi sin θi ] Using sum-to-product trigonometric identities, the above equals:
= r  2 sin(αi−1−αi2 ) sin(αi−1+αi2 ) 2 sin(αi−αi−12 ) cos(αi+αi−12 ) 2(−1)i sin ( αi+αi−1
2
) cos ( αi−αi−1
2
) 2(−1)i+1 cos ( αi−1−αi
2
) cos ( αi−1+αi
2
)[cos θi sin θi ]
= 2r  sin(αi−1−αi2 ) sin(αi−1+αi2 ) − sin(αi−1−αi2 ) cos(αi+αi−12 ) (−1)i sin ( αi+αi−1
2
) cos ( αi−αi−1
2
) (−1)i+1 cos ( αi−1−αi
2
) cos ( αi−1+αi
2
)[cos θi sin θi ]
= 2r  sin(αi−1−αi2 ) (−1)i cos ( αi−1−αi
2
)[sin(αi−1+αi2 ) − cos(αi−1+αi2 )] [cos θisin θi ]
Using equation (15), we directly have: αi = 2θi − αi−1. Thus, the above equation reduces to:
= 2r
[ sin (θi − αi)
(−1)i cos (θi − αi)
] [sin (θi) − cos (θi)] [ cos θi sin θi ] = [ 0. 0. ] Hence, the linear functions given by equations (16) and (17) are equal at ϕ = θi. This proves that the function is continuous for Case 1.

heading: Proof for Case 2:
text: Using equation (14), we know that σΘ realizes different linear functions for θ0 < ϕ < θ0 + and θ2n − < ϕ < θ2n. Thus, for σΘ to be continuous, we require that the two linear functions be the same at the boundary i.e ϕ = θ0. As before, we first write the input (z1, z2) in terms of shifted polar coordinates i.e: (r cos(ϕ), r sin(ϕ)). Thus, the function for θ0 < ϕ < θ0 + is given by:
σΘ (r cosϕ, r sinϕ) = [ cosα0 sinα0 sinα0 − cosα0 ] [ r cosϕ r sinϕ ] (18)
Similarly, the function for θ2n − < φ < θ2n is given by:
σΘ (r cosϕ, r sinϕ) =
[ cosα2n−1 sinα2n−1
(−1)i sinα2n−1 (−1)i+1 cosα2n−1
] [ r cosϕ r sinϕ ] (19)
Using equation (15), α2n−1 is given as follows:
α2n−1 = 2 2n−1∑ i=0 θ2n−1−i(−1)i = 2 n−1∑ i=0 (θ2i+1 − θ2i) = 2π
Thus, equation (19) reduces to:
σΘ (r cosϕ, r sinϕ) = [ 1 0 0 1 ] [ r cosϕ r sinϕ ] (20)
The difference in the function values at the boundary i.e ϕ = θ0, obtained by subtracting equations (20) and (18) is given as follows:[
1 0 0 1 ] [ r cos θ0 r sin θ0 ] − [ cosα0 sinα0 sinα0 − cosα0 ] [ r cos θ0 r sin θ0 ] = r [ 1− cosα0 − sinα0 − sinα0 1 + cosα0 ] [ cos θ0 sin θ0
]

heading: Using the trigonometric identities: 1 − cos(θ) = 2 sin2(θ/2), 1 + cos(θ) = 2 cos2(θ/2) and
text: sin(θ) = 2 sin(θ/2) cos(θ/2), we have:
= r
[ 2 sin2(α02 ) −2 sin( α0 2 ) cos( α0 2 )
−2 sin(α02 ) cos( α0 2 ) 2 cos 2(α02 )
] [ cos θ0 sin θ0 ] = 2r
[ sin(α02 ) − cos(α02 ) ] [ sin(α02 ) − cos( α0 2 ) ] [cos θ0 sin θ0 ] Using equation (15), we have: α0 = 2θ0. Thus, the above equation reduces to:
= 2r [ sin(θ0) − cos(θ0) ] [sin(θ0) − cos(θ0)] [ cos θ0 sin θ0 ] = [ 0. 0. ] Hence, the linear functions given by equations (18) and (20) are equal at ϕ = θ0. This proves that the function is continuous for Case 2.
B SELECTION OF γ USING CROSS VALIDATION
Using 5000 held out samples from CIFAR-10, we tested 6 different values of γ shown in Table 3 and selected γ = 0.5 because it resulted in less than 0.5% decrease in standard accuracy while 4.96% increase in provably robust accuracy. We used the LipConvnet-5 network with the 2D Householder activation function of order 2 i.e σθ.
C DIFFERENCES BETWEEN l1 AND l2 CERTIFICATE
By the equivalence of norms, we have the following result:
x ∈ Rd, 1√ d ‖x‖1 ≤ ‖x‖2 ≤
√ d‖x‖1 (21)
Let us assume that we have an l1 certificate for the input x so that the prediction of a neural network remains constant in a region of l1 radius ρ1 around the input x i.e:
‖x∗ − x‖1 ≤ ρ1 (22)
We want to compute the l2 certificate implied by the certificate given in equation (22). Let ρ2 be the l2 certificate so that:
‖x∗ − x‖2 ≤ ρ2 =⇒ ‖x∗ − x‖1 ≤ ρ1 Using equation (21), we have the following set of equations:
‖x∗ − x‖2 ≤ ρ2 =⇒ 1√ d ‖x∗ − x‖1 ≤ ρ2 =⇒ ‖x∗ − x‖1 ≤
√ dρ2
Using equation (22), we have the following: √ dρ2 ≤ ρ1 =⇒ ρ2 ≤
ρ1√ d
Hence, the l2 norm certificate induced by the l1 norm certificate can be significantly smaller for high dimensional inputs. Using the l1 certificate used in Levine & Feizi (2021) i.e ρ1 = 4 and d = 32 √ 3
for CIFAR-10 and CIFAR-100, we get an implied certificate of ρ2 = 4/(32 √
3) = 0.07225. In contrast, in this work we study l2 robustness for much higher certificate values i.e 36/255 = 0.14118.

heading: D RESULTS USING REVISED LIPSCHITZ CONSTANTS
text: In Tables 4 and 5, we show results where the Lipschitz constant was computed using product of Lipschitz constant of all layers. The Lipschitz constant of each layer was computed using direct power iteration on the linear layer (using 50 iterations) and not using the approximation bound provided in Singla & Feizi (2021).
E VERIFICATION THAT σθ(z1, z2) ALWAYS LIES ON ONE SIDE OF THE HYPERPLANE
Consider the case: z1 sin (θ/2)− z2 cos (θ/2) > 0 In this case σθ(z1, z2) = (z1, z2) and the result follows directly. Consider the other case: z1 sin (θ/2)− z2 cos (θ/2) ≤ 0[
a1 a2
] = [ cos θ sin θ sin θ − cos θ ] [ z1 z2 ] = [ z1 cos θ + z2 sin θ z1 sin θ − z2 cos θ ] a1 sin (θ/2)− a2 cos (θ/2) = (z1 cos θ + z2 sin θ) sin (θ/2)− (z1 sin θ − z2 cos θ) cos (θ/2) = z1 (cos θ sin (θ/2)− sin θ cos (θ/2)) + z2 (sin θ sin (θ/2) + cos θ cos (θ/2)) = −z1 sin (θ/2) + z2 cos (θ/2)
Since z1 sin (θ/2)− z2 cos (θ/2) ≤ 0, we have −z1 sin (θ/2) + z2 cos (θ/2) ≥ 0.

heading: F HIGHER ORDER HOUSEHOLDER ACTIVATION FUNCTIONS
text: We know that MaxMin(z1, z2) = (max(z1, z2), min(z1, z2)) where z1, z2 ∈ R. Because max(z1, z2) > min(z1, z2), applying MaxMin again gives the same result i.e MaxMin◦MaxMin = MaxMin. Now consider the function σθ (discussed in the maintext, given below for convenience):
σθ (z1, z2) =  [ 1 0 0 1 ] [ z1 z2 ] if z1 sin (θ/2)− z2 cos (θ/2) > 0[ cos θ sin θ
sin θ − cos θ ] [ z1 z2 ] if z1 sin (θ/2)− z2 cos (θ/2) ≤ 0
From Figure 1, we observe that if (u1, u2) = σθ(z1, z2), then (u1, u2) always lies on the right side of the hyperplane (pink colored region). In other words, u1 sin (θ/2) − u2 cos (θ/2) > 0 ∀ z1, z2 (proof in Appendix E). This further implies that σθ ◦ σθ = σθ. However from Figure 2a in the maintext, we observe that if we use a different angle φ where φ 6= θ + 2nπ for some n ∈ Z, then σφ(u1, u2) 6= (u1, u2) for all (u1, u2) in the pink colored region (u1 sin (θ/2) − u2 cos (θ/2) > 0). This motivates us to construct the function σ(n) : R2 → R2 defined as follows:
σ(n) = σθ ◦ σθ ◦ σθ . . . ◦ σθ︸ ︷︷ ︸ n times, θ’s can be different
(23)
Clearly, σ(n) has a larger number of linear regions than σθ and thus expected to have more expressive power. However, a drawback of using σ(n) is that it requires a sequential application of σθ which can
be expensive if the number of iterations n is large. To address this limitation, first observe that σθ realizes the same linear function for (z1, z2) and (cz1, cz2) when c > 0 i.e∇(z1,z2) σθ = ∇(cz1,cz2) σθ. Since σθ is piecewise linear, σθ(cz1, cz2) = cσθ(z1, z2). Thus, the input of the next function in the iteration is scaled by c as well and its linear function (or the Jacobian) remains unchanged. By induction, same holds for all the subsequent iterations. By chain rule, the Jacobian of composition of functions is equal to the product of Jacobian of each individual function. Since the Jacobian of each function is unchanged on scaling by c > 0, the Jacobian ∇(z1,z2) σ(n) also remains unchanged: ∇(z1,z2) σ(n) = ∇(cz1,cz2) σ(n). This suggests that it is possible to determine the Jacobian ∇(z1,z2) σ(n) for the input (z1, z2) by first converting to the polar coordinates ( √ z21 + z 2 2 , ϕ) and
then using the phase angle ϕ alone (where cos(ϕ) = z1/ √ z21 + z 2 2 , sin(ϕ) = z2/ √ z21 + z 2 2).
This motivates us to construct another GNP piecewise linear activation that only depends on the phase of the input but unlike σθ, it is allowed to have more than 2 linear regions without requiring a sequential application. This construction is given in the following theorem (example in Figure 2b):
Theorem 2. Given: 0 ≤ θ0 < θ1 · · · < θ2n = 2π + θ0 such that ∑n−1 i=0 (θ2i+1 − θ2i) = π and
αi = 2 ∑i j=0 θi−j(−1)j , The function σΘ : R2 → R2 is continuous, GNP and 1-Lipschitz where Θ = [θ0, θ1, . . . , θ2n−1] (also called Householder Activation of order n in 2 dimensions):
σΘ(z1, z2) =
[ cosαi sinαi
(−1)i sinαi (−1)i+1 cosαi ] [ z1 z2 ] θi ≤ ϕ < θi+1
where ϕ ∈ [θ0, θ2n = 2π + θ0) and cos(ϕ) = z1/ √ z21 + z 2 2 , sin(ϕ) = z2/ √ z21 + z 2 2 .
Using the definition of αi, α2n−1 can be computed as follows:
α2n−1 = 2 2n−1∑ j=0 θ2n−1−j(−1)j = 2 n−1∑ j=0 (θ2n−1−2j − θ2n−2−2j) = 2 n−1∑ j=0 (θ2j+1 − θ2j) = 2π
σΘ(z1, z2) = [ cosα2n−1 sinα2n−1 − sinα2n−1 cosα2n−1 ] [ z1 z2 ] = (z1, z2) θ2n−1 ≤ ϕ < θ2n(= θ0 + 2π)
By continuity, σΘ(z1, z2) = (z1, z2) for ϕ = θ0. Thus if we set θ0 = 0, σΘ is fixed to be identity when ϕ = 0 (or z2 = 0). However, a learnable θ0 offers the flexibility of choosing arbitrary intervals around the ϕ = θ0 to be the identity function (while of course allowing θ0 = 0). Since we can choose any interval of 2π for the phase angle, we choose ϕ ∈ [θ0, θ2n = θ0 + 2π) instead of the usual [0, 2π) to allow this possibility. We call ( √ z21 + z 2 2 , ϕ), the shifted polar coordinates.
Additionally, we make the following observations about Theorem 2. First, by construction σΘ has 2n linear regions. Second, since ∑n−1 i=0 (θ2i+1− θ2i) = π, the sum of angles subtended by linear regions with determinant −1 equals π. This in turn implies that sum of angles subtended by linear regions with determinant +1 must also equal π. Third, again using ∑n−1 j=0 (θ2j+1 − θ2j) = π, we know that only 2n− 1 of the 2n parameters in [θ0, θ1, θ2, . . . , θ2n−1] can be chosen independently implying that σΘ has 2n− 1 learnable parameters. In contrast, σ(n) has only n learnable parameters. Fourth, when n = 1, σΘ reduces to the original σθ activation.
Because every function of the form σ(n) (equation (23)) can have potentially 2n linear regions, while σΘ has only 2n linear regions, σΘ cannot express every function of the form σ(n). The primary benefit of using σΘ is that it can be easily applied by first determining the angle ϕ (using shifted polar coordinates), the region [θi, θi+1) to which ϕ belongs and the Jacobian for this region. This requires 1 multiplication with the Jacobian instead of n required for σ(n).

heading: G EXTENSION TO HIGHER DIMENSIONS
text: We introduced Householder activation function of Order 1 in m dimensions in maintext Definition 1. However, it suffers from the limitation that it has only 2 linear regions thus limiting its expressive power. The construction given in Appendix Section F allows more than 2 linear regions but is valid only for 2 dimensional inputs. This motivates us to construct Householder activations that depend on all the m components of input z ∈ Rm, (m ≥ 3) while allowing for more than 2 linear regions.
A straightforward way of constructing such an activation function is to apply an orthogonal matrix Q ∈ Rm×m, followed by dividing the output Qz into groups of size 2 each and then applying σθ to each group. However, since 1-Lipschitz neural networks involve multiplication with an orthogonal weight matrix followed by GNP activation anyway, this construction is trivial because it does not lead to additional gains in expressive power over using 2 dimensional σΘ activation functions.
Recall that the function σv is given by the following equation:
σv(z) =
{ z, vT z > 0,
(I− 2vvT )z, vT z ≤ 0. (24)
By a similar analysis as for the 2-dimensional case (Figure 2 in maintext), a repeated application of σv leads to increased number of linear regions and thus higher expressive power. This motivates us to construct the function σ(m,n) : Rm → Rm by applying the function σv (equation (24)) n times iteratively with different learnable parameter v at each iteration:
σ(m,n) = σv ◦ σv ◦ σv . . . ◦ σv︸ ︷︷ ︸ n times, v’s can be different
Since σv realizes the same linear function for both the inputs z and cz i.e ∇z σv = ∇cz σv when c > 0, σ(m,n) satisfies this property as well. This suggests that it is possible to determine the Jacobian of σ(m,n) for the given input z by projecting z onto a unit sphere: z/‖z‖2. Moreover, we want our constructed function to have at least 2n linear regions while requiring k iterations of σv where k is independent of n. This motivates the following open question: Open Problem. Can non-trivial order-n (n > 1) householder activation functions with 2n linear regions be constructed for m dimensional input (m > 2) using k iterations of σv where k is independent of n (but may depend on m)?

heading: H ADDITIONAL RESULTS ON CIFAR-10
text: The rows "BCOP", "Cayley" and "SOC (baseline)" all use the MaxMin activation function. "SOC + HH" replaces MaxMin with 2D Householder activation of order 1 (σθ), "+ CR" adds Certificate Regularization (CR) with γ = 0.1 (while using σθ as the activation function).
In Table 9, the row "SOC + HH(2)" uses Householder activation of order 2 (σΘ defined in Theorem 2), "+ CR" adds Certificate Regularization (CR) with γ = 0.1 (while using the HH activation of order 2 i.e σΘ as the activation function).
None of the results in Tables 7, 8 and 9 include Last Layer Normalization (LLN).

heading: Architecture Running times (seconds) MaxMin HH

heading: I ADDITIONAL RESULTS ON CIFAR-100
text: All results in Tables 13, 14 and 15 include Last Layer Normalization (LLN).
The rows "BCOP", "Cayley" and "SOC (baseline)" all use the MaxMin activation function. "SOC + HH" replaces MaxMin with 2D Householder activation of order 1 (σθ), "+ CR" adds Certificate Regularization (CR) with γ = 0.1 (while using σθ as the activation function).
In Table 15, the row "SOC + HH(2)" uses Householder activation of order 2 (σΘ defined in Theorem 2), "+ CR" adds Certificate Regularization (CR) with γ = 0.1 (while using the HH activation of order 2 i.e σΘ as the activation function).

Review summaries:
This paper addresses the problem of provable L2 robustness in image classification. Three techniques are proposed as additions on top of the SOC work (Singla & Feizi, 2021): 1) let the last linear layer be non-orthogonal and adjust provable L2 robustness calculation accordingly; 2) add a loss term that encourages larger margins in logits; 3) a new nonlinearity.  Comparisons against the original SOC work are reported on CIFAR-10 and CIFAR-100.

The authors propose three methods to improve deterministic l2 robustness certificates of 1-Lipschitz convolutional neural networks.   1. Last Layer Normalization (LLN):      (1) Relax the orthogonalization of the last weight layer of the network W.      (2) The direction of the rows in W can change freely; do not need to update other rows in W and forgot the learned information.     (3) Efficient than original robustness certificate when the number of classes K is large.  2. Certificate Regularizer (CR):      (1) Enable the gradients to continue increasing the margin M_f(x) while entropy loss becomes small by adding the regularization term.      (2) The relu function in the maximization term ensures that the optimization tries to increase the certificates only for the correctly classified inputs. 3. Householder Activation Functions (HH):      (1) Introduce a general linear GNP activation function named HH for the special case MaxMin used in Anil et al., 2018. and other works.      (2) Prove this HH transformation is necessary for any GNP piecewise linear function to be continuous.      Their experimental results demonstrate the effectiveness of these proposed methods compared to the baseline method SOC. 

This work is following the line of (Anil et al 2018), its goal is to propose criteria that allows to train models wirth orthogonal weights and to certify they are. The authors generalize (Anil et al 2018) by considering a more expressive non-linearity which leads to unitary jacobians, with a softer criterium during supervision which allows to obtain better accuracies. Overall, I feel the paper is well motivated, and the results are convincing.