Year: 2021

Abstract:
The convolution with orthogonal input-output Jacobian matrix, i.e., orthogonal convolution, has attracted substantial attention due to their excellent properties. A convolution layer with an orthogonal Jacobian matrix is 1-Lipschitz in the 2-norm, making the output robust to the perturbation in input. Meanwhile, an orthogonal Jacobian matrix preserves the gradient norm in back-propagation, which is critical for stable training deep networks. Nevertheless, existing orthogonal convolutions are burdened by high computational costs for preserving orthogonality. In this work, we exploit the relation between the singular values of the convolution layer’s Jacobian and the structure of the convolution kernel. To achieve the orthogonality, we explicitly construct the convolution kernel for enforcing all singular values of the convolution layer’s Jacobian to be 1s. After training, the explicitly constructed orthogonal (ECO) convolution can be constructed only once, and their weights are stored. Then, in evaluation, we only need to load the stored weights of the trained ECO convolution, and the computational cost of ECO convolution is the same as the standard dilated convolution. It is significantly more efficient than the recent state-of-the-art approach, skew orthogonal convolution (SOC) in evaluation. Our extensive experiments on CIFAR-10 and CIFAR-100 datasets demonstrate that the proposed ECO convolution is faster than SOC in evaluation while leading to competitive standard and certified robust accuracies.

heading: 1 INTRODUCTION
text: A layer with an orthogonal input-output Jacobian matrix is 1-Lipschitz in the 2-norm, robust to the perturbation in input. Meanwhile, it preserves the gradient norm in back-propagating the gradient, which effectively overcomes the gradient explosion and attenuation issues in training deep neural networks. In the past years, many studies have shown that exploiting orthogonality of the Jacobian matrix of layers in neural networks can achieve provable robustness to adversarial attacks (Li et al., 2019), stabler and faster training (Arjovsky et al., 2016; Xiao et al., 2018), and improved generalization (Cogswell et al., 2016; Bansal et al., 2018; Sedghi et al., 2019).
In a fully-connected layer y = Wx where W ∈ Rcout×cin is the weight matrix, the layer’s inputoutput Jacobian matrix J = ∂y∂x is just W. Thus, preserving the orthogonality of J can be accomplished through imposing the orthogonal constraint on W, which has been extensively studied in previous works (Mhammedi et al., 2017; Cissé et al., 2017; Anil et al., 2019). In contrast, in a convolution layer, the Jacobian matrix is no longer the weight matrix (convolution kernel). Instead, it is the circulant matrix composed of convolution kernel (Karami et al., 2019; Sedghi et al., 2019). Thus, generally, simply constructing an orthogonal convolution kernel cannot achieve an orthogonal convolution. Achieving an orthogonal Jacobian matrix in the convolution layer is more challenging than that in a fully-connected layer. It is plausibly straightforward to expand the convolution kernel to the doubly block circulant Jacobian matrix and directly impose the orthogonal constraint on the Jacobian matrix. Nevertheless, it is extremely difficult to construct a Jacobian matrix that is both doubly block-circulant and orthogonal.
Block Convolutional Orthogonal parameterization (BCOP) (Li et al., 2019) is one of the pioneering works for constructing the orthogonal convolution neural networks. It adopts the construction algorithm (Xiao et al., 2018) which decomposes a 2-D convolution into a stack of 1-D convolutions and a channel-wise orthogonal transformation. Trockman & Kolter (2021) maps the convolution kernel and the feature tensor into the frequency domain using Fast Fourier Transform (FFT), and achieves
the orthogonality of the Jacobian matrix through Cayley transform on the weight matrix in the frequency domain. They devise the convolution kernel of the same size as the input feature map, which takes more parameters than standard convolution layers with a local reception field. Meanwhile, FFT maps the real-value feature map and the weight matrix into matrices of complex values, increasing the computational cost to 4 times its counterpart in the real-value domain. Meanwhile, the Cayley transform requires computing the matrix inverse, which is not friendly for GPU computation. Recently, Skew Orthogonal Convolutions (SOC) (Singla & Feizi, 2021b) devises skew-symmetric filter and exploits matrix exponential (Hoogeboom et al., 2020) to attain the orthogonality of the Jacobian matrix. But SOC is slow in evaluation since it needs to apply a convolution filter multiple times sequentially on the feature map to obtain the Taylor expansion. Observing the efficiency limitations of the existing methods, in this work, we propose an explicitly constructed orthogonal (ECO) convolution, which is fast in both training and evaluation.
It is well known that the Jacobian matrix of a layer is orthogonal if and only if each of its singular values is 1. Thus we convert the problem of ensuring the orthogonality of the Jacobian matrix into making every singular value as 1. Based on the relation between the singular values of the Jacobian matrix for the convolution layer and the weight matrix of the convolution layer discovered by Sedghi et al. (2019), we construct the convolution kernel so that it ensures every singular value of the Jacobian matrix to be 1. Compared with the recent state-of-the-art method, SOC (Singla & Feizi, 2021b) implicitly approximating orthogonal convolution by multiple times convolution operations, ours explicitly builds the orthogonal convolution. Thus, we can directly deploy the constructed orthogonal convolution in evaluation, taking the same computational cost as the standard convolution. It is more efficient than SOC with multiple times convolution operations in evaluation. Experiments on CIFAR10 and CIFAR100 show that, in evaluation, taking less time, ours achieves competitive standard and robust accuracy compared with SOC.

heading: 2 RELATED WORK
text: Weight orthogonalization. In a fully-connected layer, the input-output Jacobian matrix is the transform itself. Thus many efforts are devoted to orthogonalizing the weights. Some early works exploit an orthogonal weight initialization to speed up the training (Saxe et al., 2014; Pennington et al., 2017) and build extremely deep neural networks (Xiao et al., 2018). Recently, more efforts are devoted to exploiting the orthogonality not only in initialization but also in the whole training process. Some approaches propose “soft” constraints on the weight matrix. For example, Bansal et al. (2018); Xiong et al. (2016) introduce mutual coherence and spectral restricted isometry as orthogonal regularization on the weight matrix. Parseval Networks (Cissé et al., 2017) adapts a regularizer to encourage the orthogonality of the weight matrices. However, these methods can not guarantee the exact orthogonality of the weight matrix. Some approaches orthogonalize features during the forward pass. For example, Huang et al. (2018b) extends Batch Normalization (Ioffe & Szegedy, 2015) with ZCA; (Huang et al., 2018a) solves a suboptimization problem to decorrelate deep features. These methods can not avoid computationally expensive operations like SVD, thus slow in practice. Some approaches enforce orthogonality by using the Riemannian optimization on the Stiefel manifold. For example, Casado & Martı́nez-Rubio (2019) uses Pade approximation as an alternative to matrix exponential mapping to update the weight matrix on the Stifel manifold. Li et al. (2020) uses an iterative Cayley transform to enforce weight matrix orthonormal. Anil et al. (2019) also attains orthogonal weights by orthogonal parameterization based on Björck orthonormalization (BjoRick & Bowie, 1971). Lastly, some approaches incorporate orthogonal constraint into their network architecture. For example, Mhammedi et al. (2017) builds orthgonal layers in RNNs with householder reflections.
Orthogonal convolution. In a convolution layer, the Jacobian matrix is no longer the weight matrix. Thus, the above-mentioned weight orthogonalization methods cannot be trivially applied in convolution layers. Wang et al. (2020); Qi et al. (2020) encourage the orthogonality of the Jacobian matrix of convolution layers through a regularizer, but they cannot achieve the strict orthogonality and cannot attain the provable robustness. Block Convolutional Orthogonal parameterization (BCOP) (Li et al., 2019) is a pioneering work for enforcing the orthogonality of the Jacobian matrix of a convolution layer. It conducts parameterization of orthogonal convolutions by adapting the construction algorithm proposed by Xiao et al. (2018). But BCOP is slow in training. Trockman & Kolter (2021) applies the Cayley transform to a skew-symmetric convolution weight in the Fourier domain so that the convolution recovered from the Fourier domain has an orthogonal Jacobian matrix. Though it achieves higher efficiency in some applications than BCOP, it is still significantly more costly than
standard convolutional layers. Skew Orthogonal Convolution (Singla & Feizi, 2021b) achieves orthogonal convolution through Taylor expansion of the matrix exponential. It is faster than BCOP in training but slower in evaluation. In contrast, our method is efficient in both training and evaluation. Recently, Su et al. (2021) achieves orthogonal convolutions via paraunitary systems in an elegant manner. Nevertheless, it also suffers from inefficiency issue.

heading: 3 PRELIMINARY
text: Notation. For a matrix M, M[i, j] denotes the element in the i-th row and the j-th column, M[i, :] denotes the i-th row and M[:, j] denotes the j-th column. We use the similar notation for indexing tensors. For n ∈ N, [n] denotes the set of n numbers from 0 to n − 1, that is, {0, 1, · · · , n − 1}. We denote the function to obtain the singular values by σ(·). Let vec(·) denotes the function which unfolds a matrix or higher-order tensor into a vector. Let ωk = exp(2πi/k) where i = √ −1. Let Fk ∈ Rk×k be matrix for the discrete Fourier transform (DFT) for k-element vectors and each entry of Fk is defined as Fk[p, q] = ω pq k .

heading: 3.1 RELATION BETWEEN THE JACOBIAN MATRIX AND THE WEIGHT MATRIX
text: We denote the input tensor of a convolution layer by X ∈ Rc×n×n, where n is the spatial size, and c is the number of input channels. We denote the convolution operation on X by convŴ(X ), where Ŵ ∈ Rk×k×c×c size where k is the size of the reception field and k is normally smaller than size of the input image, n. In this definition, we set the number of output channels identical to that of input channels. Meanwhile, we set the convolution stride as 1 and adapt the cyclic padding to make the output of the same size as the input.
We denote the output tensor by Y ∈ Rc×n×n and each element of Y is obtained by ∀l, s, t ∈ [c]× [n]× [n], Y [l, s, t] = ∑ r∈[c] ∑ p∈[k] ∑ q∈[k] X [l + r, s+ p, t+ q]Ŵ [p, q, l, r]. (1)
Since the convolution operation is a linear transformation, there exists a linear transformation M ∈ Rn2c×n2c which satisfies:
Y = convŴ(X )⇔ vec(Y) = Mvec(X ). (2)
To simplify illustration, Ŵ ∈ Rk×k×c×c is normally expanded into W ∈ Rn×n×c×c through zeropadding. We term W as the expanded convolution kernel. In W , only k2c2 elements are non-zero. Based on the above definition, we review the following theorem in Sedghi et al. (2019): Theorem 1 (see Sedghi et al. (2019) Section 2.2). For any expanded convolution kernel W ∈ Rn×n×c×c, let M is the matrix encoding the linear transformation of a convolution with W . For each p, q ∈ [n]× [n], let P(p,q) be the c× c matrix computed by
∀s, t ∈ [c]× [c], P(p,q)[s, t] = (F>nW [:, :, s, t]Fn)[p, q]. (3) Then
σ(M) = ⋃
p∈[n],q∈[n]
σ(P(p,q)). (4)
The above theorem gives the relation between the singular values of the transformation matrix M and the expanded convolution W . Theorem 2. For any W ∈ Rn×n×c×c, the Jacobian matrix of the convolution layer, J, is orthogonal if and only if ∀p, q ∈ [n]× [n], each singular value of the matrix P(p,q) is 1.
Proof. A real matrix is orthogonal if and only if each singular value is 1. Meanwhile, the Jacobian matrix J = ∂vec(Y)∂vec(X ) = M. Based on Theorem 1, each singular value of M is 1 if and only if ∀p, q ∈ [n]× [n], each singular value of the matrix P(p,q) is 1, thus we complete the proof.

heading: 3.2 MATRIX ORTHOGONALIZATION
text: Carley transform. Orthogonal parameterization has been extensively studied in the previous works. Many works (Absil et al., 2008; Helfrich et al., 2018; Maduranga et al., 2019) utilize Cayley transform on skew-symmetric matrices (A>= -A) to obtain the orthogonal weight matrices. Cayley
transform is a bijection between skew-symmetric matrices and orthogonal matrices without −1 eigenvalues. Specifically, it maps a skew-symmetric matrix A to an orthogonal matrix Q by
Q = (I−A)(I+A)−1. (5)
Li et al. (2020) proposes an iterative approximation of the Cayley transform for orthogonallyconstrained optimizers and achieves higher speed for training CNNs and RNNs.
Matrix exponential. Casado & Martı́nez-Rubio (2019) exploits the matrix exponential of the skewsymmetric matrix for orthogonal parameterization. To be specific, if A is a skew-symmetric matrix, exp(A) is orthogonal. Skew Orthogonal Convolutions (SOC) (Singla & Feizi, 2021b) approximates exp(A) through Taylor series exp(A) ≈ ST (A) = ∑T i=0 Ai
i! , and derives a bound on the approximation error:
‖exp(A)− ST (A)‖2 ≤ ‖A‖T2 T ! , (6)
where ‖ · ‖2 denotes the spectral norm, which is the maximum singular value of a matrix. When ‖A‖2 is small, the approximation can achieve high precision, thus SOC further adopts the spectral normalization on A to make ‖A‖2 small. Orthogonality of convolution. In a fully-connected layer, its Jacobian matrix is just the transpose of the weight matrix. Thus we can easily achieve the orthogonality of the Jacobian matrix through orthogonal parameterization on the weight matrix through the well-studied orthogonal parameterization approaches such as Cayley transform or matrix exponential of the skew-symmetric matrix. Nevertheless, in a convolution layer, the Jacobian matrix is no longer the transpose of the weight matrix. Thus we can not trivially apply the orthogonal parameterization approaches on the weight matrix to achieve the orthogonality of the Jacobian matrix.

heading: 3.3 LIPSCHITZNESS AND PROVABLY ROBUSTNESS
text: Lipschitz constant. A function f : Rm → Rn is L-Lipschitz under the `2 norm iff ‖f(x) − f(y)‖2 ≤ L‖x− y‖2, ∀(x,y) ∈ Rm × Rm, where L is a positive constant. We define the smallest value of L satisfying L-Lipschitz for f as the Lipschitz constant of f , denoted by Lip(f). Given two Lipschitz continuous functions f and g, the Lipschitz constant of the a composition of them, Lip(f ◦ g), is upper-bounded by the product of the Lipschitz constants of f and g (Weaver, 2018):
Lip(f ◦ g) ≤ Lip(f)Lip(g). (7)
Thus, given a network n(·) consisting of L layers, if the Lipschitz constant of each layer is 1, the Lipschitz constant of the whole network, Lip(n) is less than 1.
Provably robustness. Considering a neural network n(·) for c-class classification. Given an input x with the ground-truth label y ∈ [c], the network generates a c-dimensional prediction vector n(x) = [s1, · · · , sc] where the i-th item in n(x) denotes the confidence score for the class i. The margin of prediction for x by n(·) is defined as
Mn(x) = max(0, sy −maxi 6=ysi). (8)
As proved by Tsuzuku et al. (2018) and Li et al. (2019), for a network n(·) with a Lipschitz constant of L, x is provably robustly classified by n(·) under perturbation with a `2-norm of Mn(x)√2L , that is,
argmax i
n(x+ δ) = y, ∀‖δ‖2 < Mn(x)√
2L . (9)

heading: 4 METHOD
text: Based on Theorem 2, to achieve an orthogonal Jacobian matrix J for a convolution layer, we need to ensure that, ∀p, q ∈ [n] × [n], each singular value of the matrix P(p,q) is 1. When P(p,q) is a real matrix, this condition is equivalent to the orthogonality of P(p,q). To simplify the notation, we assemble the matrices P(p,q) (p, q ∈ [n]× [n]) into a tensor P ∈ Rn×n×c×c such that P [p, q, :, :] = P(p,q). Based on the definition, we rewrite the transformation from W to P(p,q) in Eq. (3) into
∀s, t ∈ [c]× [c],P [:, :, s, t] = F>nW [:, :, s, t]Fn. (10)
Since Fn is the full-rank square matrix, the mapping from W to P is invertible. Thus, to achieve an orthogonal Jacobian matrix, we can directly construct orthogonal P [p, q, :, :] (p, q ∈ [n]× [n]) using Carley transform or matrix exponential of skew-symmetric matrices (introduced in Section 3.2) to make each singular value of them identical to 1. Then the convolution kernel weight W can be determined uniquely from the P through the inverse transformation in Eq. (10). Since Fn is the Discrete Fourier Transform (DFT) matrix, F>nW [:, :, s, t]Fn in Eq. (10) is equivalent to conducting a 2-dimensional DFT on W [:, :, s, t]. Thus, we can recover W [:, :, s, t] through the inverse 2-dimensional DFT (IDFT2D) on the constructed orthogonal matrix P [:, :, s, t]. That is,
W [:, :, s, t] = IDFT2D(P [:, :, s, t]), ∀s, t ∈ [c]× [c]. (11)
Based on Theorem 2, the Jacobian matrix of the convolution layer with the kernel W is orthogonal. Nevertheless, the above approach has many drawbacks. First of all, for a generally constructed P , the recovered expanded convolution kernel W computed based on the inverse 2-dimensional DFT might contain n2c2 non-zeros elements. It significantly increases the number of parameters in the convolution layer from k2c2 to n2c2. Meanwhile, the convolution with kernel W has a global reception field, which is also slower than the original convolution kernel with a local reception field. Besides, the above approach needs to compute n2 times matrix orthogonalization (Carley transform or matrix exponential) on P [p, q, :, :] (p, q ∈ [n] × [n]). Considering the spatial size of the input feature tensor, n and the input/out channel number c are not a small value, n2 times orthogonalization on c× c matrices might be computationally expensive. Thus, to achieve a high efficient orthogonal convolution, we devise a 2-dimensional periodic P . We set the period as k and assume n is divisible by k. Specifically, P satisfies
P [p, q, :, :] = P [p+ i ∗ k, q + j ∗ k, :, :], ∀[p, q] ∈ [k]× [k], [i, j] ∈ [n/k]× [n/k]. (12)
We define P0 = P [0:k, 0:k, :, :] ∈ Rk×k×c×c. That is, P can be obtained by repeating P0 by n/k × n/k times in the first two dimensions. It has many benefits by the 2-dimensional periodic configuration. First of all, based on the property of DFT, the convolution kernel W [:, :, s, t] obtained from an inverse 2D DFT on the 2-dimensional periodic P [:, :, s, t] on contains only k2 non-zero elements. That is, the number of non-zero elements in the recovered W is only k2c2, which is the same as the number of elements in original convolution kernel Ŵ with a local reception field of k × k size. Besides, the non-zero elements in W can be efficiently obtained. To be specific, we define W0 ∈ Rk×k×c×c obtained from W0[:, :, s, t] = IDFT2D(P0[:, :, s, t]). Based on the property of DFT (Nussbaumer, 1981), W and W0 have the following relation:
W [i, j, :, :] = { W0[ i∗kn , j∗k n , :, :], if i% n k = 0 and j% n k = 0,
0c×c, otherwise. (13)
That is, we can obtain the non-zero elements of W ∈ Rn×n through computing W0 ∈ Rk×k, which just takes only k2 times 2D FFT on c × c matrices and is efficient. Moreover, since P is periodic along the first two dimensions, to achieve an orthogonal Jacobian matrix based on Theorem 2, we only need to explicitly construct k2 orthogonal matrices P0[p, q, :, :] ((p, q) ∈ [k] × [k]) to ensure that n2 matrices P [p, q, :, :] ((p, q) ∈ [n]× [n]) are orthogonal, which is also efficient. Moreover, as the recovered W [i, j, :, :] is non-zero only if i%nk = 0 and j% n k = 0. This property corresponds to the dilated convolution (Yu & Koltun, 2016). Specifically, using the relation between W and W0 in Eq. (13), the following equation holds:
ConvW(X ) = ConvW0(X ,dilation = n
k ). (14)
That is, we can achieve the orthogonal convolution through a dilated convolution with the kernel W0 and the dilation n/k. The dilated convolution has been well supported by mainstream deep learning frameworks like Pytorch and Tensorflow. Based on above analysis, the process of achieving the orthogonal convolution is decomposed into three steps:
1. Construct the orthogonal matrices P0[p, q, :, :], ∀(p, q) ∈ [k]× [k]. 2. Recover W0[:, :, s, t] from P0[:, :, s, t] using inverse 2D DFT, ∀(s, t) ∈ [c]× [c]. 3. Conduct the dilated convolution using the recovered kernel W0 and the dilation n/k.
Algorithm 1: The proposed efficient orthogonal convolution. In training, the operations in lines 1-14 are involved in orthogonal convolution. In contrast, in evaluation, only the operation in line 14 is conducted, the operations in lines 1-13 have been conducted before evaluation and the obtained W has been stored in the model for further evaluation. Input: Input tensor X ∈ Rn×n×c where c is the number of input channels, the spatial size of the convolution, k, the parameter tensor P̂0 ∈ RL×c×c (if k is odd, L is k 2+1 2 else k2+4 2 ) Output: Output tensor Y ∈ Rn×n×c from conducting orthogonal convolution on X . 1 for all i ∈ [L] do 2 P̂0[i, :, :]← P̂0[i, :, :]− P̂0[i, :, :]> 3 P̂0[i, :, :]← SpectralNorm(P̂0[i, :, :]) 4 P̂0[i, :, :]← ∑T i=0 P̂0[i,:,:]i i! 5 end 6 I = ConstructMap(k) % construct the mapping matrix I ∈ Rk×k 7 P0 = zeros(k, k, c, c) ; W0 = zeros(k, k, c, c) 8 for all (i, j) ∈ [k]× [k] do 9 P0[p, q, :, :] = P̂0[I[p, q], :, :]
10 end 11 for all (s, t) ∈ [c]× [c] do 12 W0[:, :, s, t] = IDFT2D(P0[:, :, s, t]) 13 end 14 Y = ConvW0(X ,dilation = n/k) 15 return Y
Below we introduce the details in these three steps individually.
Orthogonal parameterization for P0[p, q, :, :]. This step can be implemented by any matrix orthogonalization approach such as Carley Transform on skew-symmetric matrix or matrix exponential of the skew-symmetric matrix. By default, we adopt the Taylor expansion for the exponential of the skew-symmetric matrix proposed in SOC (Singla & Feizi, 2021b) in this step considering its simplicity and efficiency. To be specific, we start from a parameter tensor P0 ∈ Rk×k×c×c. ∀ (p, q) ∈ [k]× [k], we compute
P0[p, q, :, :]← P0[p, q, :, :]−P0[p, q, :, :]>,
P0[p, q, :, :]← T∑
i=0
P0[p, q, :, :]i
i! . (15)
Recall from Eq. (6) that the approximation error between exp(A) and the Taylor expansion Sk(A), ‖exp(A) − Sk(A)‖2, is bounded by ‖A‖ k 2
k! . That is, when ‖A‖2 < 1, it can achieve a good approximation. Thus, following SOC (Singla & Feizi, 2021b), we conduct spectral normalization on P0[p, q, :, :] so that ‖exp(P0[p, q, :, :])− ∑T i=0 P0[p,q,:,:]i i! ‖2 ≤ ‖P0[p,q,:,:]‖k2 k! = 1 k! .
Recover W0 from P0. ∀(s, t) ∈ [c]× [c], inverse 2D DFT is on P0[:, :, s, t] to obtain W0[:, :, s, t]: W0[:, :, s, t] = (F>k )−1P0[:, :, s, t]F−1k = IDFT2D(P0[:, :, s, t]). (16) As W0 is the convolution kernel for the further convolution, to make the convolution compatible with mainstream real-value convolutional neural networks, we should keep W0 real. Based on the property of DFT, for inverse 2D discrete Fourier transform, the output Y ∈ Rk×k if purely real if and only if the input X is conjugate-symmetry, that is, X[i, j] = X∗[(k−i)%k, (k−j)%k] for each (i, j) ∈ [k] × [k]. Thus, to make W0 real, we need set P0[i, j, :, :] := P∗0[(k−i)%k, (k−j)%k, : , :] for each (i, j) ∈ [k] × [k]. Since P0 is real-value, we just need to set P0[i, j, :, :] and P0[(k−i)%k, (k−j)%k, :, :] identical. Due to the conjugate-symmetry settings, among these k2 matrices {P0[(i, j, :, :]}k,ki=1,j=1, there are only L different matrices. L is (k2 + 1)/2 when k is odd and (k2+4)/2 when k is even. In the implementation, we encode all the information of P0 through a L× c× c tensor P̂0. Then we build P0 from P̂0 by a mapping matrix I ∈ Rk×k:
P0[i, j, :, :] = P̂0[I[i, j], :, :],∀(i, j) ∈ [k]× [k], (17)
where I is a mapping matrix, which maps a pair (i, j) ∈ [k]× [k] to a scalar l ∈ [1, L] and it satisfies I[i, j] = I[(k−i)%k, (k−j)%k],∀(i, j) ∈ [k]× [k]. (18)
The details of constructing the mapping matrix I is provided in Algorithm 2 of Appendix A.
Dilation convolution. In this step, we conduct the standard dilated convolution using the recovered kernel W0 obtained in Eq. (16) with a dilation n/k on the input X and obtain the output Y . We summarize the proposed model in Algorithm 1. It is worth noting that, in each training iteration, given a batch of B samples, the lines 1-13 are only conducted once, and the obtainedW0 is shared across B samples. When the training finishes, we store the computed W0 in line 12 in the model. In evaluation, we directly load the stored W0 and only need to compute the operation in line 14, which is just vanilla dilated convolution and thus is very efficient.
In the current configuration of Algorithm 1, we set P0[p, q, :, :] ∈ Rc×c to be real matrix. In fact, P0[p, q, :, :] can also be complex matrix. To make the output of the proposed orthogonal convolution compatible with existing mainstream real-value neural network, the only constraint is to make W0[:, :, i, j] = IDFT2DP0[:, :, i, j] a real matrix. That is, the only constraint is to set P0[:, :, i, j] to be conjugate symmetric and thus P0[:, :, i, j] is unnecessarily real. When considering P0[p, q, :, :] as a complex matrix, it makes the constructed orthogonal convolution kernel W0 more complete. Meanwhile, it brings more computational cost in training. In Appendix B, we compare with using a complex P0[p, q, :, :].

heading: 5 IMPLEMENTATION
text: Spectral normalization. For a small approximation error from Taylor expansion, we conduct spectral normalization on P0[p, q, :, :] to make ‖P0[p, q, :, :]‖2 ≤ 1. Following Miyato et al. (2018); Singla & Feizi (2021a;b), we use the power method to estimate the approximated spectral norm and then divide P0[p, q, :, :] by the estimated spectral norm. We use the mix-precision training supported by the APEX package for faster training and less GPU memory cost.
Padding. Following SOC (Singla & Feizi, 2021b), we use cyclic padding to substitute zeropadding since zero-padded orthogonal convolutions must be size 1 (see proof in Appendix D.2 of SOC (Singla & Feizi, 2021b)). To be specific, for a dilated convolution of k× k size with a dilation d, we add d(k − 1)/2 paddings on each side in a cyclic manner. Input and output channels. The proposed orthogonal convolution is devised under the condition that the number of channels in input is equal to that in output, that is, cin = cout = c. When cin 6= cout, following SOC (Singla & Feizi, 2021b), we additionally conduct the following operations:
1. if cin > cout, we apply our orthogonal convolution with a kernel size of k × k × cin × cin, generating the output with cin channels. We only select the first cout channels in the output.
2. if cin < cout, we use zero pad the input with cout − cin channels. Then we apply our orthogonal convolution with a kernel size of k × k × cout × cout.
Strided convolution. Given an input X ∈ Rn×n×cin , a convolution with stride s generates the output Y ∈ Rns×ns×cout . The proposed orthogonal convolution is designed for 1-stride convolution.
For a more general s-stride convolution, we decompose it into a down-sampling module and an orthogonal convolution, following BCOP (Li et al., 2019) and SOC (Singla & Feizi, 2021b). Specifically, the down-sampling module adapts the invertible down-sampling (Maduranga et al., 2019), which simply reshapes the input feature map X from Rn×n×cin to Rns×ns×cins2 . The invertible down-sampling decreases the spatial size of the feature map from k×k to k/s×k/s and meanwhile strictly preserves the gradient norm. Then we conduct the proposed 1-stride orthogonal convolution on the reshaped X . 1-Lipschitz architecture. Following SOC (Singla & Feizi, 2021b), we build a provably 1-Lipschitz architecture based on the proposed orthogonal convolution. For a fair comparison with SOC, our structure is the same as SOC, and we only replace the orthogonal convolution in SOC with the proposed orthogonal convolution. The details of the architecture are given in Table 1. It consists of L convolution layers in 5 blocks, and each block contains L/5 layers. Each block consists of two groups of convolutions layers. The first group repeats L/5−1 convolutions layers with stride 1, and the second group contains a single convolution layer with stride 2, and the second group contains a single convolution layer with stride 2. Following SOC, a MaxMin activation function (Anil et al., 2019) is added after each convolution layer.

heading: 6 EXPERIMENTS
text: Settings. The training is conducted on a single NVIDIA V100 GPU with 32G memory. Following Singla & Feizi (2021b), the whole training process takes 200 epochs for all models. The initial learning rate is 0.1 when the number of convolution layers is larger than 25 and 0.05 otherwise, and it decreases by a factor of 0.1 at the 50th and the 150th epoch. We set the weight decay as 0 when training BCOP (Li et al., 2019) and 5e-4 when training SOC (Singla & Feizi, 2021b) and our ECO. By default, we set the batch size as 256 on all models in training.
Evaluation metrics. We evaluate the performance proposed orthogonal convolution by two metrics: standard recognition accuracy and the robust recognition accuracy. The standard recognition
accuracy is the same as that used in classic image classification. In contrast, the robust recognition accuracy is the accuracy under an l2 perturbation. As shown in Section 3.3, L-Lipschitz network has the robustness certificate when the `2-norm perturbation is less than Mn(x)/ √ 2L. Thus, for a 1-Lipschitz we built based on the proposed ECO convolution has the robustness certificate as Mn(x)/ √ 2 ≥ . Following Li et al. (2019) and Singla & Feizi (2021b), when testing the robust accuracy, we set the default norm of perturbation, , as 36/255. In this case, the robust accuracy is the ratio of samples satisfying Mn(x)/ √ 2 ≥ 36/255.
Standard and robust recognition accuracy. Table 2 compares the recognition accuracy of the proposed ECO convolution with BCOP (Li et al., 2019) and SOC (Singla & Feizi, 2021b) on CIFAR10 and CIFAR100 datasets. As we can observe from the table that our ECO convolution achieves competitive performance compared with BCOP and SOC in both standard and robust recognition accuracy. Specifically, on CIFAR-10 dataset, when the number of convolution layers, L, is larger than 25, ours outperforms both SOC and BCOP in standard and robust recognition accuracy. Meanwhile, on CIFAR-100 dataset, our ECO convolution achieves higher robust accuracy than SOC in all models from LipConvnet-5 to LipConvnet-40.
Training and evaluation time. As shown in the last column of Table 2, SOC is slightly more efficient than our ECO convolution in training. It is because the dilation convolution in our ECO adds more padding, making the input feature map larger than that in SOC. Meanwhile, as shown in Table 3, our ECO convolution is significantly faster than SOC in evaluation. It is due to that we have already constructed the ECO convolution before evaluation for once. Thus, in evaluation, the computational cost of our ECO convolution is the same as standard convolution. In contrast, SOC needs to take T times convolution operations per convolution layer to approximate the orthogonal convolution, which takes T times the computational cost as standard convolution. By default, SOC sets T = 5 in training and T = 10 in evaluation. As shown in Table 3, in LipConvnet-40, ours ECO convolution achieves an around 6× speed-up ratio over SOC. Robust accuracy with respect to distortion. As mentioned, the 1-Lipschitz network has the robustness certificate as Mn(x)/ √ 2 ≤ and the testing accuracy is the ratio of samples satisfying
Mn(x)/ √ 2 ≤ . We evaluate the robust accuracy of LipConvnet-25 when the distortion norm adapts different values. We compare with SOC and the robust accuracies of SOC under different are copied from Singla et al. (2021) using the same MaxMin activation function. As shown in Table 4, our ECO convolution consistently outperforms ECO under different and the advantage is more significant for a larger .

heading: 7 CONCLUSION
text: This work investigates building the convolution layer with an orthogonal Jacobian matrix, which is 1-Lipschitz in the 2-norm and preserves the gradient. By exploiting the relation between the convolution kernel and the Jacobian matrix, we explicitly construct the convolution kernel to achieve the orthogonal convolution. The explicitly constructed orthogonal convolution (ECO) takes the exact computational cost as a standard dilated convolution in evaluation, which is more efficient than the existing state-of-the-art method, SOC, approximating the orthogonal convolution through multiple times convolution operations. Our experiments on CIFAR-10 and CIFAR-100 show that our ECO achieves competitive standard and robust accuracy as SOC using much less evaluation time.

heading: A CONSTRUCTING THE MAPPING MATRIX I
text: In Algorithm 2, we give the details of constructing the mapping matrix I.

heading: B COMPARISONS WITH THE COMPLEX SETTINGS
text: As we mentioned in Section 4, to obtain a real-value convolution kernel W0, the entries in P0 are unnecessarily real. We compare the performance when P0 is constrained to be real as Algorithm 1 and that when P0 is a complex tensor. As shown in Table 5, on CIFAR-10 benchmark, using complex P0, it consistently outperforms its real counterpart on both standard accuracy and robust accuracy. The accuracy gap between the complex P0 and its real counterpart is large for shallow networks and turns smaller as the network goes deeper. It is worth noting that the computational cost of multiplication between complex matrices is four times as that between real matrices of the same sizes. But the complex P0 is only used for constructing the real orthogonal kernel W0. After kernel construction, the computation cost of convolution is the same for the layer with complex P0 and that with real P0. Thus, as shown in Table 5, the training time per epoch when using complex P0 is around 2 to 3 times as that using real P0. In evaluation, since we have constructed the real W0 from a real or complex P0 before evaluation, the computational cost is the same no matter P0 is complex and real.

heading: C COMPARISONS WITH THE STANDARD CONVOLUTION
text: Here, we compare the proposed orthogonal convolution with the standard convolution without orthogonal constraint. To make a fair comparison, we only replace the orthogonal convolution layers of LipConvnet-5 by the standard convolution and keep the other layers unfixed. Since the network built upon the standard convolution cannot preserve the 1-Lipschitz property, we cannot directly obtain the robust accuracy based on Mn(x) as the original LipConvnet-5 network. Thus, we utilize Projected Gradient Descent (PGD) attack (Madry et al., 2017) to obtain the robust accuracy upper bound under the distortion with the `2 norm less than 36255 . The results are shown in Table 6. As shown in the table, on both CIFAR10 and CIFAR10 datasets, the network based on the standard con-
volution achieves higher standard accuracy than the network based on our ECO convolution. In the meanwhile, ours achieves higher robust accuracy on both datasets under the same scale of distortion.

heading: D STANDARD TRAINING ON RESNET
text: We conduct experiments based on the ResNet-18, which is one of mainstream backbones in computer vision. To be specific, we replace the standard convolution in ResNet-18 with the orthogonal convolutions based on BCOP, SOC and ECO. As shown in Table 7, the orthogonal convolution can not achieve the high recognition accuracy as its standard counterpart. But the advantage of orthogonal convolution layers is their provable robustness in LipConvNet, making them still attractive in scenarios when the recognition robustness is important.

heading: E TRUNCATED TAYLOR EXPANSION
text: Note that, the exact matrix exponential (Casado & Martı́nez-Rubio, 2019) exp(A) relies on singular value decomposition (SVD), which is not friendly for GPU computing platform. Thus, we adopt the truncated Taylor Expansion (Singla & Feizi, 2021b), ST (A) = ∑T i=0 Ai
i! , to approximate the exact matrix exponential exp(A). The truncation error R = exp(A)− ST (A) satisfies
‖R‖2 = ‖exp(A)− ST (A)‖2 ≤ ‖A‖T2 T ! . (19)
Below we first compare the time cost in GPU platform using the exact matrix exponential exp(A) and that based on the truncated Taylor Expansion ST (A). To be specific, we implement exp(A) simply based on the torch.matrix exp. We measure the time cost on GPU using the function torch.cuda.Event and we set T = 10 in truncated Taylor Expansion ST (A). As shown in Table 8, the latency based on the truncated Taylor Expansion is significantly less than that based on the exact matrix exponential.
Then we show the scale of the average factual truncated error R with respect to T . It is based on the trained LipConvnet-5 network built upon the proposed ECO convolution. By default, we set T = 5 in training and T = 10 in testing. As shown in Table 9, the scale of the truncation error is relatively small when T is 5 or 10. Meanwhile, we also show the influence of the truncation error on the orthogonality of the convolution. Following Trockman & Kolter (2021), we measure the orthogonality of a convolution layer through ‖vec(conv(X ))‖2‖vec(X )‖2 −1, and a strict orthogonal convolution layer should satisfy
‖vec(conv(X ))‖2 ‖vec(X )‖2 − 1 = 0. (20)
We estimate the deviation of the proposed ECO convolution from the strict orthogonal convolution through this measurement. As shown in Table 9, when T is 5 or 10, the deviation of our ECO convolution from the strict orthogonal convolution is relative marginal.

Review summaries:
Summarizing the paper -------------------------------      The paper proposes a method to construct a  convolutional layer with the orthogonal Jacobian matrix. Such a layer is 1-Lipschits: this property is an important one for building robust neural networks. The authors conduct experiments on classification tasks using CIFAR10/CIFAR100 datasets. They report standard and robust accuracies, training, and evaluation time. They claim the proposed method outperforms existing state-of-the-art approaches.      

This paper presents a method for enforcing strict orthogonality for convolutional layers. The method is based on considering the spectral domain, where the orthogonality of the 4D conv kernel (in spatial domain) is characterized as the orthogonality of 2D matrices, for which orthogonality can be enforced by existing techniques such as Carley transform. It is shown in experiments that this method is significantly faster than the recent skew orthogonal convolution (at ICML'21) method.

This paper studies how to construct orthogonal convolutional networks in an efficient way. To this end, this paper builds the connection between the DFT-transformed kernel with the common dilated convolution. During training, the forward pass can be done by a sequence of inverse DFT and dilated convolution. During testing, all the convolution kernels only need to be transformed once so that evaluation time is significantly reduced.