SCALING-UP DIVERSE ORTHOGONAL CONVOLU-
Year: 2021

Abstract:
Enforcing orthogonality in convolutional neural networks is an antidote for gradient vanishing/exploding problems, sensitivity to perturbation, and bounding generalization errors. Many previous approaches for orthogonal convolutions enforce orthogonality on its flattened kernel, which, however, do not lead to the orthogonality of the operation. Some recent approaches consider orthogonality for standard convolutional layers and propose specific classes of their realizations. In this work, we propose a theoretical framework that establishes the equivalence between diverse orthogonal convolutional layers in the spatial domain and the paraunitary systems in the spectral domain. Since 1D paraunitary systems admit a complete factorization, we can parameterize any separable orthogonal convolution as a composition of spatial filters. As a result, our framework endows high expressive power to various convolutional layers while maintaining their exact orthogonality. Furthermore, our layers are memory and computationally efficient for deep networks compared to previous designs. Our versatile framework, for the first time, enables the study of architectural designs for deep orthogonal networks, such as choices of skip connection, initialization, stride, and dilation. Consequently, we scale up orthogonal networks to deep architectures, including ResNet and ShuffleNet, substantially outperforming their shallower counterparts. Finally, we show how to construct residual flows, a flow-based generative model that requires strict Lipschitzness, using our orthogonal networks.

heading: 1 INTRODUCTION
text: Convolutional neural networks (CNNs), whose deployment has witnessed extensive empirical success, still exhibit a range of limitations that are not thoroughly studied. Firstly, deep convolutional networks are in general difficult to learn, and their high performance heavily relies on techniques that are not fully understood, such as skip-connections (He et al., 2016a), batch normalization (Ioffe & Szegedy, 2015), specialized initialization (Glorot & Bengio, 2010). Secondly, they are notoriously sensitive to imperceptible perturbations, including adversarial attacks (Goodfellow et al., 2014) or geometric transformations (Azulay & Weiss, 2019). Finally, a precise characterization of their generalization property is still under active investigation (Neyshabur et al., 2018; Jia et al., 2019).
Orthogonal networks, which have a “flat” spectrum with all singular values of the layer being 1 (each layer’s output norm ‖y‖ equals its input norm ‖x‖, ∀x), alleviate all problems above. As shown in recent works, by enforcing orthogonality in neural nets, we obtain (1) easier optimization (Zhang et al., 2018a; Qi et al., 2020): since each orthogonal layer preserves the gradient norm during backpropagation, an orthogonal network is free from gradient vanishing/ exploding problems; (2) robustness against adversarial perturbation (Anil et al., 2019; Li et al., 2019b; Trockman & Kolter, 2021): since each orthogonal layer is 1-Lipschitz, an orthogonal network can not amplify any perturbation to the input to flip the output prediction; (3) better generalizability as proved in (Jia et al., 2019): a network’s generalization error is positively related to the standard deviation of each linear layer’s singular values, thus encouraging orthogonality in the network lowers its generalization error.
Despite the benefits, enforcing orthogonality in convolutional networks is challenging. To avoid the strict constraint, orthogonal initialization (dynamical isometry) (Pennington et al., 2017; Xiao et al., 2018) and orthogonal regularization (Wang et al., 2019; Qi et al., 2020) are adopted to address the gradient vanishing/exploding problems. However, these methods are not suitable for applications that
require strict Lipschitzness, such as adversarial robustness (Anil et al., 2019) and residual flows (Chen et al., 2019), as they do not enforce strict orthogonality (and Lipschitzness) during training.
Our goal is to enforce exact orthogonality in state-of-the-art convolutional networks without expensive computations. We identify three main challenges in doing so. Challenge I: Achieving exact orthogonality throughout the training process. Prior works such as orthogonal regularization (Wang et al., 2019) and reshaped kernel orthogonality (Jia et al., 2017; Cisse et al., 2017), while enjoying algorithmic simplicity, do not meet the requirement of exact orthogonality. Note that enforcing orthogonality during training is necessary as a post-training orthogonalization of the weights can ruin the performance. Challenge II: Avoiding expensive computations. An efficient algorithm is crucial for scalability to large networks. Existing work based on projected gradient descent (Sedghi et al., 2019), however, requires expensive projection after each update. For instance, the projection step in Sedghi et al. (2019) computes an SVD and flattens the spectrum to enforce orthogonality, which costs O(size(feature) · channels3) for a convolutional layer. Challenge III: Scaling-up to state-of-the-art deep convolutional networks. There are many variants to the standard convolutional layer, including dilated, strided, group convolutions, which are essential for state-of-the-art deep convolutional networks. However, none of the existing methods proposes mechanisms to orthogonalize these variants. The lack of techniques, as a result, limits the broad applications of orthogonal convolutional layers to state-of-the-art deeper convolutional networks.
We resolve challenges I, II, and III by proposing complete parameterizations for orthogonal 1Dconvolutions and separable orthogonal 2D-convolutions. First, using the convolution theorem (Oppenheim et al., 1996) (convolution in the spatial domain is equivalent to multiplication in the spectral domain), we reduce the problem of designing orthogonal convolutions to constructing unitary matrices for all frequencies, i.e., a paraunitary system (Vaidyanathan, 1993). Therefore, we obtain a parameterization for all separable orthogonal 2D-convolutions, guaranteeing exact orthogonality and high expressive power. Note that there is no previous approach that achieves exact orthogonality (up to machine precision). For the first time, our versatile framework enables the study in the designs of deep orthogonal networks, such as skip connection, initialization, stride, and dilation. Consequently, we scale orthogonal networks to deep architectures, including ResNet and ShuffleNet, substantially outperforming their shallower counterparts. Our proposed method serves as gaining insight into the trade-off between orthogonality and expressiveness. We observe that exact orthogonality is essential in these deeper architectures with more than 10 layers (detailed discussion of exact orthogonality in Appendix E.3). Finally, we show how to deploy our orthogonal networks in Residual Flow (Chen et al., 2019), a flow-based generative model that requires strict Lipschitzness (Appendix F).

heading: Summary of Contributions:
text: 1. We establish the equivalence between orthogonal convolutions in the spatial domain and paraunitary systems in the spectral domain. Consequently, we can interpret existing approaches for orthogonal convolutions as implicit designs of paraunitary systems. 2. Based on a complete factorization of 1D paraunitary systems, we propose the first exact and complete design for orthogonal 1D-convolutions as well as separable orthogonal 2D-convolutions, ensuring exact orthogonality and high expressive power. None of the existing works achieve a complete parameterization of all orthogonal 2D-convolutions. 3. We prove that orthogonality for various convolutional layers (strided, dilated, group) are also entirely characterized by paraunitary systems. Consequently, our design easily extends to these variants, ensuring both completeness and exactness of the orthogonal convolutions. 4. We study the design considerations for orthogonal networks (choices of skip connection, initialization, depth, width, and kernel size), and show that orthogonal networks can scale to deep architectures including ResNet, WideResNet, and ShuffleNet.

heading: 2 ORTHOGONAL CONVOLUTIONS VIA PARAUNITARY SYSTEMS
text: Designing orthogonal convolutional layer {ht,s : yt = ht,s ∗ xs}T,St=1,s=1 (s, t index input/output channels) in the spatial domain is challenging. When a convolution layer is written as a matrix-vector product, where the matrix is block-circulant and its (t, s)th takes a circulant form Cir (ht,s) as:
Cir (ht,s) =  ht,s[1] ht,s[N ] · · · ht,s[2] ht,s[2] ht,s[1] ht,s[N ] · · ·
... . . .
. . . ...
ht,s[N ] · · · ht,s[2] ht,s[1]
 ∈ RN×N . (2.1)
Therefore, the layer is orthogonal if the block-circulant matrix [Cir (ht,s)] T,S t=1,s=1 is orthogonal. However, it is not obvious how to enforce orthogonality in a block-circulant matrix.

heading: 2.1 ACHIEVING ORTHOGONAL CONVOLUTIONS BY PARAUNITARY SYSTEMS
text: We propose a novel design of orthogonal convolutions from a spectral domain perspective, motivated by the convolution theorem (Theorem 2.1). For simplicity, we group the entries at the same locations into a vector/matrix, e.g., we denote {xs[n]}Ss=1 as x[n] ∈ RS and {ht,s[n]} T,S t=1,s=1 as h[n] ∈ RT×S . Theorem 2.1 (Convolution theorem (Oppenheim et al., 1996)). For a standard convolution layer h: y[i] = ∑ n h[n]x[i − n], the convolution in the spatial domain is equivalent to a matrixvector multiplication in the spectral domain, i.e., Y (z) = H(z)X(z),∀z ∈ C. Here, X(z) =∑N−1 n=0 x[n]z −n, Y (z) = ∑N−1 n=0 y[n]z −n, H(z) = ∑L n=−L h[n]z
−n denote the z-transforms of input, output, kernel respectively, where N is the length of x,y and [−L,L] is the span of the filter h.
The convolution theorem states that a standard convolution layer is a matrix-vector multiplication in the spectral domain. As long as the transfer matrix H(z) is unitary at z = ejω for all frequencies ∀ω ∈ R (j is the imaginary unit), the convolutional layer h is orthogonal. Therefore, as a major novelty of our paper, we design orthogonal convolutions through designing unitary transfer matrix H(ejω) at all frequencies ω ∈ R, which is known as a paraunitary systems (Vaidyanathan, 1993; Strang & Nguyen, 1996). We prove in Theorem B.5 that a convolutional layer is orthogonal in the spatial domain if and only if it is paraunitary in the spectral domain.
Benefits through paraunitary systems. (1) The spectral representation simplifies the designs of orthogonal convolutions and avoids the analysis of block-circulant structure. (2) Since paraunitary systems are necessary and sufficient for orthogonal convolutions, it is impossible to find an orthogonal convolution whose transfer matrix is not paraunitary. (3) There exists a complete factorization of paraunitary systems: any paraunitaryH(z) can be realized through multiplications in the spectral domain, as will show in Equation (2.2a). (4) Since multiplications in spectral domain correspond to convolutions in spatial domain, any orthogonal convolution can be realized as cascaded convolutions of multiple sub-layers, each parameterized by an orthogonal matrix. (5) There are mature methods that represent orthogonal matrices via unconstrained parameters. Consequently, we can learn orthogonal convolutions using standard optimizers on a model parameterized via our design.
Interpretation of existing methods. Since paraunitary system is a necessary and sufficient condition for orthogonal convolution, all existing approaches, including singular value clipping and masking (SVCM) (Sedghi et al., 2019), block convolution orthogonal parameterization (BCOP) (Li et al., 2019b), Cayley Convolution (CayleyConv) (Trockman & Kolter, 2021), skew orthogonal convolution (SOC) design paraunitary systems implicitly. We discuss these interpretations in Appendix C.3.3, and show that our SC-Fac has the lowest computational complexity among all approaches.

heading: 2.2 REALIZING PARAUNITARY SYSTEMS VIA RE-PARAMETERIZATION
text: After reducing the problem of orthogonal convolutions to paraunitary systems, we are left with the question of how to realize paraunitary systems. We use a complete factorization form of paraunitary systems to realize any paraunitary systems. According to Theorem C.7, any paraunitary systemH(z) can be written as the form in Equation (2.2a) and anyH(z) in this form is a paraunitary system:
H(z) = V (z;U (−L)) · · ·V (z;U (−1))QV (z−1;U (1)) · · ·V (z−1;U (L)), (2.2a)
where V (z;U (`)) = (I −U (`)U (`) > ) +U (`)U (`) > z, ∀` ∈ {−L, · · · ,−1} ∪ {1, · · · , L}. (2.2b)
Here Q is an orthogonal matrix and each U (`) is a column-orthogonal matrix whose number of columns is sampled uniformly from {1, · · · , T}. As multiplications in the spectral domain are equivalent to convolutions in the spatial domain, this complete spectral factorization of paraunitary systems in Equation (2.2a) allows us to parameterize any orthogonal convolution in the spatial domain as cascaded convolutions of V (z;U (`))’s spatial counterparts and the orthogonal matrixQ.
Model design in the spatial domain. We obtain a complete design of orthogonal 1D-convolutions: using learnable (column)-orthogonal matrices ({U (`)}−1`=−L,Q, {U (`)}L`=1), we parameterize a size (L+ L+ 1) convolution as cascaded convolutions of the following filters in the spatial domain{[
I −U (`)U (`) > , U (`)U (`) >]}−1 `=−L , Q, {[ U (`)U (`) > , I −U (`)U (`) >]}L `=1 . (2.3)
Figure 1 provides a visualization of our proposed design of orthogonal convolution layers; each block denotes a convolution and the form of the filter is displayed in each of the block. In practice, we compose all (L+ L+ 1) filters into one for orthogonal convolution, which not only increases the computational parallelism but also avoids storing intermediate outputs between filters.
Using this 1D-convolution, we obtain a complete design for separable orthogonal 2D-convolutions, which can be represented as a convolution of two orthogonal 1D-convolutional layers. With separability, we obtain a design of orthogonal 2D-convolutions by composing two complete designs of orthogonal 1D-convolutions. As a result, we can parameterize separable orthogonal 2D-convolution with filter size (L1+L1+1)×(L2+L2+1) as a convolution of two orthogonal 1D-convolutions with learnable (column-)orthogonal matrices ({U (`)1 } −1 `=−L1 ,Q1, {U (`)1 } L1 `=1) and ({U (`) 2 } −1 `=−L2 ,Q2, {U (`)2 } L2 `=1). With a complete factorization of paraunitary systems (1D and separable 2D), we reduce the problem of designing orthogonal convolutions to the one for orthogonal matrices.
Parameterization for orthogonal matrices. In Appendix C.4, we perform a comparative study on different parameterizations of orthogonal matrices, including the Björck orthogonalization (Anil et al., 2019; Li et al., 2019b), the Cayley transform (Helfrich et al., 2018; Maduranga et al., 2019), and the exponential map (Lezcano-Casado & Martínez-Rubio, 2019). In our implementation, we adopt a modified exponential map due to its efficiency, exactness, and completeness. The exponential map is a surjective mapping from a skew-symmetry matrixA to a special orthogonal matrix U (i.e., det(U) = 1) with U = exp(A) = I +A+A2/2 + · · · , where the infinite series is computed up to machine-precision (Higham, 2009). To parameterize all orthogonal matrices, we introduce an additional orthogonal matrixQ in U = Q exp(A), whereQ is (randomly) generated at initialization and fixed during training (Lezcano Casado, 2019) .
Now we have an end-to-end pipeline for orthogonal convolutions as shown in Figure 2 (See Algorithm 1 in Appendix E for a pseudo implementation). Since our method relies on separability and complete factorization for 1D paraunitary systems, we call it Separable Complete Factorization (SC-Fac).

heading: 3 UNIFYING ORTHOGONAL CONVOLUTIONS VARIANTS
text: Various convolutional layers (strided, dilated, and group convolution) are widely used in neural networks. However, it is not apparent how to enforce their orthogonality, as the convolution theorem (Theorem 2.1) only holds for standard convolutions. Previous approaches only deal with standard convolutions (Sedghi et al., 2019; Li et al., 2019b; Trockman & Kolter, 2021), thus orthogonality for state-of-the-art architectures has never been studied before.
We address this limitation by modifying convolution theorem for each variant of convolution layer, which allows us to design these variants using paraunitary systems. Theorem 3.1 (Convolution and paraunitary theorems for various convolutions). Strided, dilated, and group convolutions can be unified in the spectral domain as Y (z) =H(z)X(z), where Y (z), H(z),X(z) are modified Z-transforms of y, h, x. We instantiateH(z) for strided convolutions in Proposition C.4, dilated convolution in Proposition C.5, and group convolution in Proposition C.6. Furthermore, a convolution is orthogonal if and only ifH(z) is paraunitary.
In Table 1, we formulate strided, dilated, and group convolutions in the spatial domain, interpreting them as up-sampled or down-sampled variants of a standard convolution. Now, we introduce the concept of up-sampling and down-sampling precisely below.
Given a sequence x, we introduce its up-sampled sequence x↑R with sampling rate R as x↑R[n] , x[n/R] for n ≡ 0 (mod R). On the other hand, its (r,R)-polyphase component xr|R indicates the r-th down-sampled sequence with sampling rate R, defined
as xr|R[n] , x[nR + r]. We illustrated an example of x↑R and xr|R in Figure 3 when sampling rate R = 2. The Z-transforms of x↑R, xr|R are denoted as X↑R(z), Xr|R(z) respectively. Their relations toX(z) are studied in Appendix C.1.
Now we are ready to interpret convolution variants. (1) Strided convolution is used to adjust the feature resolution: a strided convolution (↓R-strided) decreases the resolution by down-sampling after a standard convolution, while a transposed strided convolutional layer (↑R-strided) increases the resolution by up-sampling before a standard convolution. (2) Dilated convolution increases the receptive field of a convolution without extra parameters: an R-dilated convolution up-samples its filters before convolution with the input. (3) Group convolution reduces the parameters and computations, thus widely used by efficient architectures: a G-group convolution divides the input/output channels into G groups and restricts the connections within each group. In Appendix C.2, we prove that a convolution is orthogonal if and only if its modified Z-transformH(z) is paraunitary.

heading: 4 LEARNING DEEP ORTHOGONAL NETWORKS WITH LIPSCHITZ BOUNDS
text: In this section, we switch our focus from layer design to network design. In particular, we aim to study how to scale-up deep orthogonal networks with Lipschitz bounds.
Lipschitz networks (Anil et al., 2019; Li et al., 2019b; Trockman & Kolter, 2021), whose Lipschitz bounds are imposed by their architectures, are proposed as competitive candidates to guarantee robustness in deep learning. A Lipschitz network consists of orthogonal layers and GroupSort activations (Anil et al., 2019) — both are 1-Lipschitz and gradient norm preserving. See Appendix D for more discussions on the properties of GroupSort and Lipschitz networks. Given a Lipschitz constant L, a Lipschitz network f can compute a certified radius for each input from its output margin. Formally, denote the output margin of an input x with label c as
Mf (x) , max(0, f(x)c −max i6=c f(x)i), (4.1)
i.e., the difference between the correct logit and the second largest logit. Then the output is robust to perturbation such that f(x+ ) = f(x) = c,∀ : ‖ ‖ <Mf (x)/ √ 2L.
Despite the benefit, existing architectures for Lipschitz networks remain simple and shallow, and a Lipschitz network is typically an interleaving cascade of orthogonal layers and GroupSort activations (Li et al., 2019b). More advanced architectures, such as ResNet and ShuffleNet, are still out of reach. While orthogonal layers supposedly substitute the role of batch normalization (Pennington et al., 2017; Xiao et al., 2018; Qi et al., 2020), other critical factors, including skip-connections (He et al., 2016a;b) and proper initialization (Glorot & Bengio, 2010) are lacking. In this section, we explore skip-connections and initialization methods toward addressing this problem.
Skip-connections. Two general types of skip-connections are widely used in deep networks, one based on addition and another on concatenation. The addition-based connection is proposed in ResNet (He et al., 2016a), and adopted in SE-Net (Hu et al., 2018) and EffcientNet (Tan & Le, 2019). The concatenation-based connection is proposed in flow-based generative models (Dinh et al., 2014; 2016; Kingma & Dhariwal, 2018), and adopted in DenseNet (Huang et al., 2017) and ShuffleNet (Zhang et al., 2018b; Ma et al., 2018). In what follows, we propose Lipschitz skip-connections with these two mechanisms, illustrated in Figure 6 (in Appendix D). Proposition 4.1 (Lipschitzness of residual blocks). Suppose f1, f2 are L-Lipschitz (for residual/shortcut branches) and α ∈ [0, 1] is a learnable scalar, then an additive residual block f : f(x) , αf1(x)+(1−α)f2(x) is L-Lipschitz. Alternatively, suppose g1, g2 are L-Lipschitz and P is a channel permutation, then a concatenative residual block g : g(x) , P [ g1(x1); g2(x2) ] is LLipschitz, where [·; ·] denotes channel concatenation, and x is split into x1 and x2, i.e., x = [x1,x2].
Initialization. Proper initialization is crucial in training deep networks (Glorot & Bengio, 2010; He et al., 2016a). Various methods are proposed to initialize orthogonal matrices, including the identical/permutation and torus initialization in the context of orthogonal RNNs (Henaff et al., 2016; Helfrich et al., 2018; Lezcano-Casado & Martínez-Rubio, 2019). However, initialization of orthogonal convolutions was not systematically studied, and all previous approaches inherit the initialization from the underlying parameterization (Li et al., 2019b; Trockman & Kolter, 2021). In Proposition D.1, we study the condition when a paraunitary system (represented in Equation (2.2a)) reduces to an orthogonal matrix. This reduction allows us to apply the initialization methods for orthogonal matrices (e.g., uniform, torus) to orthogonal convolutions.
In the experiments, we will evaluate the impact of different choices of skip-connections and initialization methods to the performance of deep Lipschitz networks.

heading: 5 RELATED WORK
text: Dynamical isometry (Pennington et al., 2017; 2018; Chen et al., 2018; Pennington et al., 2018) aims to address the gradient vanishing/exploding problems in deep vanilla networks with orthogonal initialization. These works focus on understanding the interplay between initialization and various nonlinear activations. However, these approaches do not guarantee Lipschitzness after training, thus are unsuitable for applications that require strict characterization of Lipschitz constants, such as adversarial robustness (Anil et al., 2019) and residual flows (Behrmann et al., 2019).
Learning orthogonality has three typical approaches: regularization, parameterization (representing the feasible set with unconstrained parameters), and projected gradient descent (PGD) / Riemannian gradient descent (RGD). While regularization is approximate, the latter two learn exact orthogonality. (1) For orthogonal matrices, various regularizations are proposed in Xie et al. (2017); Bansal et al. (2018). Alternatively, numerous parameterizations exist, including Householder reflections (Mhammedi et al., 2017), Given rotations (Dorobantu et al., 2016), Cayley transform (Helfrich et al., 2018), matrix exponential (Lezcano-Casado & Martínez-Rubio, 2019), and algorithmic unrolling (Anil et al., 2019; Huang et al., 2020). Lastly, Jia et al. (2017) propose PGD via singular value clipping, and Vorontsov et al. (2017); Li et al. (2019a) consider RGD. (2) For orthogonal convolutions, some existing works learn orthogonality for the flattened matrix (Jia et al., 2017; Cisse et al., 2017; Bansal et al., 2018) or each output channel (Liu et al., 2021). However, these methods do not lead to orthogonality (norm preserving) of the operation. Sedghi et al. (2019) propose to use PGD via singular value clipping and masking, which is expensive and can result in inexact orthogonality. Recent works adopt parameterizations, using block convolutions (Li et al., 2019b), Cayley transform (Trockman & Kolter, 2021), or convolution exponential (Singla & Feizi, 2021). Note that network deconvolution (Ye et al., 2020) aims to whiten the activations (i.e., orthogonalize their distribution), but the added whitening operations do not necessarily preserve orthogonality.
Paraunitary systems are extensively studied in filter banks and wavelets (Vaidyanathan, 1993; Strang & Nguyen, 1996; Lin & Vaidyanathan, 1996). Classic theory shows that 1D-paraunitary systems are fully characterized by a spectral factorization (see Chapter 14 of Vaidyanathan (1993)), but not all MD-paraunitary systems admit a factorized form (see Chapter 8 of Lin & Vaidyanathan (1996)). While the complete characterization of MD-paraunitary systems is known in theory (which incurs a difficult problem of solving a system of nonlinear equations) (Venkataraman & Levy, 1995; Zhou, 2005), most practical constructions use separable paraunitary systems (Lin & Vaidyanathan, 1996) and special classes of non-separable paraunitary systems (Hurley & Hurley, 2012). The equivalence between orthogonal convolutions and paraunitary systems thus opens the opportunities to apply these classic theories in designing orthogonal convolutions.

heading: 6 EXPERIMENTS
text: In the experiments, we achieve the following goals. (1) We demonstrate in Section 6.1 that our separable complete factorization (SC-Fac) achieves precise orthogonality (up to machine-precision), resulting in more accurate orthogonal designs than previous ones (Sedghi et al., 2019; Li et al., 2019b; Trockman & Kolter, 2021). (2) Despite the differences in preciseness, we show in Section 6.2 that different realizations of paraunitary systems only have a minor impact on the adversarial robustness of Lipschitz networks. (3) Due to the versatility of our convolutional layers and architectures, in Section 6.3, we explore the best strategy to scale Lipschitz networks to wider/deeper architectures. (4) In Appendix F, we further demonstrate in a successful application of orthogonal convolutions in residual flows (Chen et al., 2019). Training details are provided in Appendix E.1.

heading: 6.1 EXACT ORTHOGONALITY
text: We evaluate the orthogonality of our SC-Fac layer verse previous approaches, including CayleyConv (Trockman & Kolter, 2021), BCOP (Li et al., 2019b), SVCM (Sedghi et al., 2019), RKO (Cisse et al., 2017), OSSN (Miyato et al., 2018). Our experiments are based on a convolutional layer with 64 input channels and 16× 16 input size. We orthogonalize the layer using each approach, and evaluate it with Gaussian inputs. For our SC-Fac layer, We initialize all orthogonal matrices uniformly, while we use built-in initialization for others. We evaluate the difference between 1 and the ratio of the output norm to the input norm — a layer is exactly orthogonal if the number is close to 0.
(1) Standard convolution. We show in Table 2 (Left) that our SC-Fac is orders of magnitude more precise than all other approaches. The SC-Fac layer is in fact exactly orthogonal up to machine epsilon, which is 2−24 ≈ 5.96 × 10−8 for 32-bits floats. While RKO and OSSN are known not to be orthogonal, we surprisingly find that SVCM is far from orthogonal due to its masking step. (2) Convolutions variants. In Section 3, we show that various orthogonal convolutions can be constructed using paraunitary systems. We verify our theory in Table 2 (Right): our SC-Fac layers are exactly orthogonal (up to machine precision) for various convolutions.

heading: 6.2 ADVERSARIAL ROBUSTNESS
text: In this subsection, we evaluate the adversarial robustness of Lipschitz networks. Following the setup in (Trockman & Kolter, 2021), we adopt KW-Large, ResNet9, WideResNet10-10 as the backbone architectures, and evaluate their robust accuracy on CIFAR-10 with different designs of orthogonal
convolutions. We extensively perform a hyper-parameter search and choose the best hyper-parameters for each approach based on the robust accuracy. The details of the hyper-parameter search is in Appendix E. We run each model with 5 different seeds and report the best accuracy.
(1) Certified robustness. Following Li et al. (2019b), we use the raw images (without normalization) for network input to achieve the best certified accuracy. As shown in Table 3 (Top), different realizations of paraunitary systems, SC-Fac, CayleyConv and BCOP achieve comparable performance — CayleyConv performs < 1% better in clean accuracy, but the difference in robust accuracy are negligible. (2) Practical robustness. Trockman & Kolter (2021) shows that the certified accuracy is too conservative, and it is possible to increase the practical robustness (against PGD attacks) with a standard input normalization. Notice that the normalization increases the Lipschitz bound, thus lower the certified accuracy. Our experiments in Table 3 (Bottom) are based on ResNet9, WideResNet1010 (Trockman & Kolter, 2021) and a deeper WideResNet22. For the shallow architectures (ResNet9, WideResNet10-10), our SC-Fac, CayleyConv, and BCOP again achieve comparable performance — CayleyConv is slightly ahead in robust accuracy. For the deeper architecture, our SC-Fac has a clear advantage in both clean and robustness accuracy, and the clean accuracy to only 5% lower than a traditional ResNet 32 trained with batch normalization. Surprisingly, we find that RKO also performs well in robust accuracy while not exactly orthogonal. In summary, our experiments show that various paraunitary realizations provide different impacts on certified and practical robustness. While exact orthogonality provides tight Lipschitz bound, there is a trade-off between the exact orthogonality and the practical robustness (especially with the shallow architectures).

heading: 6.3 SCALING-UP DEEP ORTHOGONAL NETWORKS WITH LIPSCHITZ BOUNDS
text: All previously proposed Lipschitz networks (Li et al., 2019b; Trockman & Kolter, 2021) consider only shallow architectures (≤ 10 layers). In this part, we investigate various factors to scale Lipschitz networks to deeper architectures: skip-connection, depth-width, receptive field, and down-sampling.
(1) Skip-connections. Conventional wisdom suggests that skip-connections mainly address the gradient vanishing/exploding problem; thus, they are not needed for orthogonal networks. To understand their role, we perform an experiment that trains deep Lipschitz networks without skipconnection and with skip-connections based on addition/concatenation (see Section 4). As shown in Table 4 (left), the network with additive skip-connection substantially outperforms the other two, and the one without skip-connections performs the worst. Therefore, we empirically show that (additive) skip-connection is crucial in deep Lipschitz networks. (2) Depth and width. Exact orthogonality is criticized for harming the expressive power of neural networks. We show that the decrease of expressive power can be alleviated by increasing the network depth/width. In Table 3 (Bottom) and Table 7 (Appendix E), we observe that deeper/wider architectures increase both the clean and robust accuracy. (3) Initialization methods. We try different initialization methods, including identical, permutation, uniform, and torus (Henaff et al., 2016; Helfrich et al., 2018). We find that identical initialization works the best for deep Lipschitz networks (> 10 layers), while all methods perform similarly in shallow networks as shown in Table 6 (Appendix E). (4) Receptive field and down-sampling. Previous works (Li et al., 2019b; Trockman & Kolter, 2021) use larger kernel
size and no stride for Lipschitz networks. In Table 4 (Right), we perform a study on the effects of kernel/dilation size and down-sampling types for the orthogonal convolutions. We find that an average pooling as down-sampling consistently outperforms strided convolutions. Furthermore, a larger kernel size helps to boost the performance. (5) Run-time and memory comparison. We find that previously proposed orthogonal convolutions such as CayleyConv, BCOP, and RKO require more GPU memory and computation time than SC-Fac. Therefore, we could not to scale them due to memory constraints (for 22 and 32 layers using Tesla V100 32G). In order to scale up Lipschitz networks, economical implementation of orthogonal convolution is crucial. As shown in Figure 4, for deep and wide architectures, our SC-Fac is the most computationally and memory efficient method and the only method that scales to a width increase of 10 on WideResNet22. Missing numbers in Figure 4 and Table 7 (Appendix E) are due to the large memory requirement.
In summary, additive skip-connections are still essential for learning deep orthogonal networks. Due to the orthogonal constraints, it is helpful to increase the depth/width of the network. However, this significantly increases the memory requirement; thus, a cheap implementation (like SC-Fac) is desirable. Finally, we find that a larger kernel size and down-sampling based on average pooling is helpful, unlike standard practices in deep networks.

heading: 7 CONCLUSION
text: In this paper, we present a paraunitary framework for orthogonal convolutions. Specifically, we establish the equivalence between orthogonal convolutions in the spatial domain and paraunitary systems in the spectral domain. Therefore, any design for orthogonal convolutions is implicitly constructing paraunitary systems. We further show that the orthogonality for variants of convolution (strided, dilated, and group convolutions) is also fully characterized by paraunitary systems. In summary, paraunitary systems are all we need to ensure orthogonality for diverse types of convolutions.
Based on the complete factorization of 1D paraunitary systems, we develop the first exact and complete design of separable orthogonal 2D-convolutions. Our versatile design allows us to study the design principles for orthogonal convolutional networks. Consequently, we scale orthogonal networks to deeper architectures, substantially outperforming their shallower counterparts. In our experiments, we observe that exact orthogonality plays a crucial role in learning deep Lipschitz networks. In the future, we plan to investigate other use cases that exact orthogonality is essential.

heading: ETHICS STATEMENT
text: Our work lies in the foundational research on neural information processing. Specifically, we establish the equivalence between orthogonal convolutions in neural networks and paraunitary systems in signal processing. Furthermore, our presented orthogonal convolutional layers are plug-and-play modules that can replace various convolutional layers in neural networks. Consequently, our modules are applicable in Lipschitz networks for adversarial robustness, recurrent networks for learning long-term dependency, or flow-based networks for effortless reversibility.
The vulnerability of neural networks raises concerns about their deployment in security-sensitive scenarios, such as healthcare systems or self-driving cars. In our experiment, we demonstrate a successful application of orthogonal convolutions in learning robust networks. Furthermore, these networks achieve higher robust accuracy without additional techniques such as adversarial training or randomized smoothing. Therefore, our research contributes to the robust learning of neural networks and potentially leads to their broader deployment.
Furthermore, we reduce the inference cost of our layer to be the same as a traditional convolution layer, significantly lowering the expense for deployment compared with other methods for orthogonal networks. However, our layers are memory and computationally more expensive during training than traditional layers. The overhead to the already expensive training cost exacerbates the concerns on the efficacy of learning neural networks. Therefore, balancing between robustness and efficiency is an important research topic that requires more research in the future. We develop more efficient implementation than previous approaches, narrowing the gap between these two conflicting goals.

heading: REPRODUCIBILITY STATEMENT
text: We provide detailed proofs to our theoretical claims in Appendices B to D. We describe experimental setups and training details of our models in Appendix E. We include our code for experiments in the supplementary material.

heading: Appendices: Scaling-up Diverse Orthogonal Convolutional Networks by a Paraunitary Framework
text: Notations. We use non-bold letters for scalar (e.g., x) and bold ones for vectors or matrices (e.g., x). We denote sequences in the spatial domain using lower-case letters (e.g., x[n], x[n]) and their spectral representations using upper-case letters (e.g., X(z),X(z)). For a positive integer, say R ∈ Z+, we abbreviate the set {0, 1, · · · , R− 1} as [R], and whenever possible, we use its lower-case letter, say r ∈ [R], as the corresponding iterator.
Assumptions. For simplicity, we assume all sequences are L2 with range Z = {0,±1,±2, · · · } (a sequence x = {x[n], n ∈ Z} is L2 if ∑ n∈Z ‖x[n]‖2 < ∞). Such assumption is common in the literature, which avoids boundary conditions in signal analysis. To deal with periodic sequences (finite sequences with circular padding), one can either adopt the Dirac function in the spectral domain or use discrete Fourier transform to compute the spectral representations. In our implementation, we address the boundary condition case by case for each convolution type, with which we achieve exact orthogonality in the experiments (Section 6.1).

heading: A PSEUDO CODE FOR OUR SC-FAC ALGORITHM
text: We include the pseudo-code for separable complete factorization (Section 2) in Algorithm 1 and diverse orthogonal convolutions (Section 3) in Algorithm 2. The pseudo-code in Algorithm 1 consists of three parts: (1) First, we obtain orthogonal matrices from skew-symmetric matrices using matrix exponential. We use GeoTorch library (Lezcano Casado, 2019) for the function matrix_exp in our implementation; (2) Subsequently, we construct two 1D paraunitary systems using these orthogonal matrices; (3) Lastly, we compose two 1D paraunitary systems to obtain one 2D paraunitary systems The pseudo-code in Algorithm 2 consists of two parts: (1) First, we reshape each paraunitary system into an orthogonal convolution depending on the stride; and (5) second, we concatenate the orthogonal kernels for different groups and return the output.

heading: B ORTHOGONAL CONVOLUTIONS VIA PARAUNITARY SYSTEMS
text: In this section, we prove the convolution theorem and Parseval’s theorem for standard convolutional layers. Subsequently, we prove the paraunitary theorem which establishes the equivalence between orthogonal convolutional layers and paraunitary systems.

heading: B.1 SPECTRAL ANALYSIS OF STANDARD CONVOLUTION LAYERS
text: Standard convolutional layers are the default building blocks for convolutional neural networks. One such layer consists of a filter bank with T × S filters h = {hts[n], n ∈ Z}t∈[T ],s∈[S], where S, T are the number of input and output channels respectively. The layer maps an S-channel input x = {xs[i], i ∈ Z}s∈[S] to a T -channel output y = {yt[i], i ∈ Z}t∈[T ] according to
yt[i] = ∑ s∈[S] ∑ n∈Z hts[n]xs[i− n], (B.1)
where i indexes the output location to be computed, and n indexes the filter coefficients. Alternatively, we can rewrite Equation (B.1) in matrix-vector form as
y[i] = ∑ n∈Z h[n]x[i− n], (B.2)
where each h[n] ∈ RT×S is a matrix, and each x[i− n] ∈ RS or y[i] ∈ RT is a vector. Notice that in Equation (B.2), we group entries from all channels into a vector or matrix (e.g., from {x0[n]}s∈[S] to x[n]), different from a common notation that groups entries from all locations into a vector or matrix (e.g., from {xs[n], n ∈ Z} into xs). In the matrix-vector form, a standard convolutional layer computes a vector sequence y = {y[i] ∈ RT , i ∈ Z} with a convolution between a matrix sequence h = {h[n] ∈ RT×S , n ∈ Z} and a vector sequence x = {x[i] ∈ RS , i ∈ Z}.
Algorithm 1: Separable Complete Factorization (SC-Fac) Input: Number of channels C, kernel size K = 2L+ 1, and
Skew-symmetric matrices {A(`)d } withA (`) d ∈ RC×C ,∀` ∈ [−L,L], d ∈ {1, 2}.
Output: A paraunitary systemH ∈ RC×C×K×K . Initialization: Sample N (`)d from {1, · · · , C} uniformly ∀` ∈ [−L,L], d ∈ {1, 2} /* Iterate for vertical/horizontal dimensions */ for d = 1 to 2 do
/* 1) Compute orthogonal matrices from skew-symmetric matrices */ /* Iterate for filter locations */ for ` = −L to L do
if ` = 0 then Qd ← matrix_exp(A(0)d ) // use matrix_exp() in GeoTorch (Lezcano Casado, 2019) else
U (`) d ← select(matrix_exp(A (`) d ), cols = N (`) d ) // selects the first cols
columns of the matrix end if
end for /* 2) Compose 1D paraunitary systems from orthogonal matrices */ Hd ← Qd for ` = 1 to L do
Hd ← conv1d(Hd, [ U (`) d U (`) d > , I −U (`)d U (`) d > ] )
Hd ← conv1d( [ I −U (−`)d U (−`) d > , U (−`) d U (−`) d > ] ,Hd)
end for end for /* 3) Compose a 2D paraunitary systems from two 1D paraunitary */ H← Compose(H1,H2) // i.e., H:,:,i,j = (H2):,:,j(H1):,:,i where the 1D
paraunitary systems H1 and H2 are of size C × C ×K returnH
Algorithm 2: Construct Diverse Orthogonal Convolutions from Paraunitary Systems Input: Number of base channels C, kernel size K = R(2L+ 1), stride R, dilation D, number of groups G Output: An orthogonal kernelW ∈ RT×S×K×K Set K ′ ← K/R, number of input channels S ← GC/R2 and output channels T ← GC for g = 0 to G− 1 do
/* 1) Construct orthogonal convolutions from paraunitary systems */ Initialize skew-symmetric matrices {{A(`,g)d }L`=−L}2d=1 for the current g Hg ← Algorithm 1: SC-Fac(C,K ′, {{A(`,g)d }L`=−L}2d=1) Hg ← reshape(Hg, (C,C,K ′,K ′)→ (C/R2, C,K,K))
end for /* 2) Concatenate orthogonal convolutions from different groups */ W ← concatenate({Hg}G−1g=0 , dim = 0) returnW (where the filter for input channel s and output channel t isWt,s,:,: ∈ RK×K)
Let us first define the Z-transform and various types of Fourier transform in Definition B.1 before proving the convolution theorem (Theorem 2.1).
Definition B.1 (Z-transform and Fourier transforms). For a sequence (of scalars, vectors, or matrices) x = {x[n], n ∈ Z}, its Z-transformX(z) is defined as
X(z) = ∑ n∈Z x[n]z−n, (B.3)
where z ∈ C is a complex number such that the infinite sum is convergent. If z is restricted to the unit circle z = ejω (i.e., |z| = 1), the z-transformX(z) reduces to a discrete-time Fourier transform (DTFT) X(ejω). If ω is further restricted to a finite set ω ∈ {2πk/N, k ∈ [N ]}, the DTFT X(ejω) reduces to an N -points discrete Fourier transform (DFT)X(ej2πk/N ).
The celebrated convolution theorem states that the convolution in spatial domain (Equation (B.2)) is equivalent to a multiplication in the spectral domain, i.e., Y (z) =H(z)X(z),∀z ∈ C.
Proof of Theorem 2.1. The proof follows directly from the definitions of standard convolution (Equation (B.2)) and Z-transform (Equation (B.3)).
Y (z) = ∑ i∈Z y[i]z−i (B.4)
= ∑ i∈Z (∑ n∈Z h[n]x[i− n] ) z−i (B.5)
= ∑ n∈Z h[n]z−n (∑ i∈Z x[i− n]z−(i−n) )
(B.6)
= (∑ n∈Z h[n]z−n )(∑ k∈Z x[k]z−k ) (B.7)
=H(z)X(z), (B.8)
where Equations (B.4) and (B.8) use the definition of Z-transform, Equation (B.5) uses the definition of convolution, and Equation (B.7) makes a change of variable k = i− n.
Next, we introduce the concepts of inner product and Frobenius norm for sequences. We then prove Parseval’s theorem, which allows us to compute the sequence norm in the spectral domain.
Definition B.2 (Inner product and norm for sequences). Given two sequences x = {x[n], n ∈ Z} and y = {y[n], n ∈ Z} with x[n],y[n] having the same dimension for all n, the inner product of these two sequences is defined as
〈x,y〉 , ∑ n∈Z 〈x[n],y[n]〉 (B.9)
where 〈x[n],y[n]〉 denotes the Frobenius inner product between x[n] and y[n]. Subsequently, we can define the Frobenius norm of a sequence using inner product as
‖x‖ , √ 〈x,x〉 (B.10)
Theorem B.3 (Parsavel’s theorem). Given a sequence x = {x[n], n ∈ Z}, its sequence norm ‖x‖ can be computed byX(ejω) in the spectral domain as
‖x‖2 = ∑ n∈Z ‖x[n]‖2 = 1 2π ∫ π −π ‖X(ejω)‖2dω, (B.11)
where ‖X(ejω)‖2 =X(ejω)†X(ejω) is an inner product between two complex arrays.
Proof of Theorem B.3. The theorem follows from the definitions of convolution and discrete-time Fourier transform (DTFT).
1
2π ∫ π −π ∥∥X(ejω)∥∥2 dω = 1 2π ∫ π −π 〈 X(ejω),X(ejω) 〉 dω (B.12)
= 1
2π ∫ π −π 〈∑ n∈Z x[n]e−jωn, ∑ m∈Z x[m]e−jωm 〉 dω (B.13)
= ∑ n∈Z ∑ m∈Z 〈x[n],x[m]〉 1 2π ∫ π −π e−jω(m−n)dω (B.14)
= ∑ n∈Z ∑ m∈Z 〈x[n],x[m]〉1m=n (B.15)
= ∑ n∈Z 〈x[n],x[n]〉 = ∑ n∈Z ‖x[n]‖2 , (B.16)
where Equation (B.14) is due to the bi-linearity of inner products, and Equation (B.15) makes uses of the fact that ∫ π −π e −jωkdω = 0 for k 6= 0 and ∫ π −π e −jωkdω = ∫ π −π dω = 2π for k = 0.

heading: B.2 EQUIVALENCE BETWEEN ORTHOGONAL CONVOLUTIONS AND PARAUNITARY SYSTEMS
text: With the sequence norm introduced earlier, we formally define orthogonality for convolutional layers. Definition B.4 (Orthogonal convolutional layer). A convolution layer is orthogonal if the input norm ‖x‖ is equal to the output norm ‖y‖ for arbitrary input x, that is
‖y‖ , √∑ n∈Z ‖y[n]‖2 = √∑ n∈Z ‖x[n]‖2 , ‖x‖, (B.17)
where ‖x‖ (or ‖y‖) is defined as the squared root of ∑ n∈Z ‖x[n]‖2 (or ∑ n∈Z ‖y[n]‖2).
This definition of orthogonality not only applies to standard convolutions in Equation (B.2) but also variants of convolutions in Appendix C.3. In this section, however, we first establish the equivalence between orthogonality for standard convolutions and paraunitary systems. Theorem B.5 (Paraunitary theorem). A standard convolutional layer (in Equation (B.2)) is orthogonal (by Definition B.4) if and only if its transfer matrixH(z) is paraunitary, i.e.,
H(z) † H(z) = I, ∀|z| = 1⇐⇒H(ejω)†H(ejω) = I, ∀ω ∈ R. (B.18)

heading: In other words, the transfer matrixH(ejω) is unitary for all frequencies ω ∈ R.
text: Proof of Theorem B.5. We first prove that a convolutional layer is orthogonal if its transfer matrix H(z) is paraunitary (i.e.,H(ejω) is unitary for any frequency ω ∈ R).
‖y‖2 = 1 2π ∫ π −π ∥∥Y (ejω)∥∥2 dω (B.19) = 1
2π ∫ π −π ∥∥H(ejω)X(ejω)∥∥2 dω (B.20) = 1
2π ∫ π −π ∥∥X(ejω)∥∥2 dω (B.21) = ‖x‖2 (B.22)
where Equations (B.19) and (B.22) are due to Parseval’s theorem (Theorem B.3), Equations (B.19) and (B.20) follows from the convolution theorem (Theorem 2.1), and Equation (B.21) utilizes that H(ejω) is unitary for any frequency ω ∈ R (thus ‖H(ejω)X(ejω)‖ = ‖X(ejω)‖ for anyX(ejω)). The ‘only if’ part also holds in practice. We here prove by contradiction using periodic inputs (e.g., finite inputs with circular padding). Suppose there exists a frequency ω and H(ejω) is not unitary.
SinceH(ejω) is continuous (due to h ∈ L2 and dominated convergence theorem), there exist integers N , k, such that ω ≈ 2kπ/N andH(ej2πk/N ) is also not unitary. As a result, there exists a complex vector u such that v =H(ej2πk/N )u while ‖v‖ 6= ‖u‖. Therefore, we can construct two periodic sequences x = {x[n], n ∈ [N ]} and y = {y[n], n ∈ [N ]} such that
x[n] = uej2πnk/N =⇒ y[n] = vej2πnk/N . (B.23)
Now the input norm ‖x‖ = √∑
n∈[N ] ‖x[n]‖2 = √ N‖u‖ is not equal to the output norm ‖y‖ =√∑
n∈[N ] ‖y[n]‖2 = √ N‖v‖, i.e., the layer is not orthogonal, which leads to a contradiction.

heading: C A PARAUNITARY FRAMEWORK FOR ORTHOGONAL CONVOLUTIONS

heading: Spectral

heading: C.1 MULTI-RESOLUTION ANALYSIS
text: Multi-resolution operations are essential in various convolutional layers, in particular strided and dilated convolutions. In order to define and analyze these convolutions rigorously, we first review the concepts of up-sampling, down-sampling, and polyphase components.
(1) Up-sampling. Given a sequence (of scalars, vectors, matrices) x = {x[n], n ∈ Z}, its upsampled sequence x↑R = {x↑R[n], n ∈ Z} is defined as
x↑R[n] , { x[n/R] n ≡ 0 (mod R) 0 otherwise , (C.1)
where R ∈ Z+ is the up-sampling rate. Accordingly, we denote the Z-transform of x↑R as
X↑R(z) = ∑ n∈Z x↑R[n]z−n. (C.2)
The following proposition shows thatX↑R(z) is easily computed fromX(z).
Proposition C.1 (Z-transform of up-sampled sequence). Given a sequence x and its up-sampled sequence x↑R, their Z-transformsX(z) andX↑R(z) are related by
X↑R(z) =X(zR). (C.3)
Proof of Proposition C.1. The proof makes use of the definition of Z-transform (Equation (B.3)). X↑R(z) = ∑ n∈Z x↑R[n]z−n (C.4)
= ∑ m∈Z x↑R[mR]z−mR (C.5)
= ∑ m∈Z x[m](zR)−m =X(zR), (C.6)
where Equation (C.5) makes a change of variables m = n/R since x↑R[n] = 0,∀n 6= mR.
(2) Down-sampling and polyphase components. Different from the up-sampled sequence, there exist multiple down-sampled sequences, depending on the phase of down-sampling. These sequences are known as the polyphase components. Specifically, given a sequence (of scalars, vectors, or matrices) x = {x[n], n ∈ Z}, its rth polyphase component xr|R = {xr|R[n], n ∈ Z} is defined as
xr|R[n] , x[nR+ r], (C.7)
where R ∈ Z+ is the down-sampling rate. We further denote the Z-transform of xr|R as Xr|R(z) = ∑ n∈Z xr|R[n]z−n, (C.8)
Note that r ∈ Z is an arbitrary integer, which does not necessarily take values from [R]. In fact, we have x(r+kR)|R[n] = xr|R[n+ k] andXr+kR|R(z) = zkXr|R(z). In Proposition C.2, we establish the relation betweenH(z) and {Hr|R(z)}r∈[R], i.e., to representH(z) in terms of {Hr|R(z)}r∈[R]. Proposition C.2 (Polyphase decomposition). Given a sequence x and its polyphase components xr|R’s, the Z-transformX(z) can be represented by {Xr|R(z)}r∈[R] as
X(z) = ∑ r∈[R] Xr|R(zR)z−r. (C.9)
Proof of Proposition C.2. We start withX(z), and try to decompose it into its polyphase components X0|R(z), · · · ,XR−1|R(z).
X(z) = ∑ n∈Z x[n]z−n (C.10)
= ∑ r∈[R] ∑ m∈Z x[mR+ r]z−(mR+r) (C.11) = ∑ r∈[R] (∑ m∈Z x[mR+ r]z−mR ) z−r (C.12) = ∑ r∈[R] (∑ m∈Z xr|R[m](zR)−m ) z−r (C.13)
= ∑ r∈[R] Xr|R(zR)z−r, (C.14)
where Equation (C.11) makes a change of variables n = mR+r, and Equation (C.13) is the definition of polyphase components xr|R[m] = x[mR+ r].
For simplicity, we stack R consecutive polyphase components into polyphase matrices as
X [R](z) =  X 0|R(z)
... XR−1|R(z)  , X̃ [R](z) = [X−0|R(z); · · · ;X−(R−1)|R(z)] . (C.15) The following proposition extends the Parseval’s theorem in Theorem B.3 and shows that the sequence norm ‖x‖ can also be computed in terms of the polyphase matrixX [R](z) (or X̃ [R](z)).
Proposition C.3 (Parseval’s theorem for polyphase matrices). Given a sequence x, its sequence norm ‖x‖ can be computed byX [R](ejω) (or X̃ [R](ejω)) in the spectral domain as
‖x‖2 = 1 2π ∫ π −π ∥∥∥X [R](ejω)∥∥∥2 dω = 1 2π ∫ π −π ∥∥∥X̃ [R](ejω)∥∥∥2 dω. (C.16) Proof of Proposition C.3. The proof follows the standard Parseval’s theorem in Theorem B.3. We only prove the first part of the proposition (usingX [R](ejω)) as follows.
‖x‖2 = ∑ n∈Z ‖x[n]‖2 (C.17)
= ∑ r∈[R] ∑ m∈Z ‖x[mR+ r]‖2 (C.18)
= ∑ r∈[R] ∑ m∈Z ∥∥∥xr|R[m]∥∥∥2 (C.19) = 1
2π ∑ r∈[R] ∫ π −π ∥∥∥Xr|R(ejω)∥∥∥2 dω (C.20) = 1
2π ∫ π −π ∥∥∥X [R](ejω)∥∥∥2 dω, (C.21) where Equation (C.17) follows the definition of sequence norm, Equation (C.18) changes variables as n = mR + r, Equation (C.19) is the definition of polyphase components, and Equation (C.20) applies Parseval’s theorem to xr|R’s. The second part (using X̃ [R](ejω)) can be proved similarly.

heading: C.2 UNIFYING VARIOUS CONVOLUTIONAL LAYERS IN THE SPECTRAL DOMAIN
text: In Appendix B, the convolution theorem states that a standard convolutional layer is a matrix-vector product Y (z) =H(z)X(z) in the spectral domain, and the layer is orthogonal if and only ifH(z) is paraunitary (Theorem B.5). However, the canonical convolution theorem does not hold for variants of convolutions, thus enforcing a paraunitaryH(z) may not lead to orthogonal convolution. In this subsection, we address this limitation by showing that various convolutions can be uniformly written as Y (z) = H(z)X(z), where Y (z), H(z), X(z) are some spectral representations of y, h, x. Subsequently, we prove that any of these layers is orthogonal if and only if itsH(z) is paraunitary.
(1) Strided convolutional layers are widely used in neural networks to adjust the feature resolution: a strided convolution layer decreases the resolution by down-sampling after a standard convolution, while a transposed convolution increases the resolution by up-sampling before a standard convolution.
Formally, a strided convolutional layer with stride R (abbrev. as ↓R-strided convolution) computes its output following y[i] = ∑ n∈Z h[n]x[Ri− n]. (C.22) In contrast, a transposed strided convolutional layer with stride R (abbrev. as ↑R-strided convolution) computes its output according to
y[i] = ∑ n∈Z h[n]x↑R[i− n]. (C.23)
Proposition C.4 (Orthogonality of strided convolutional layers). For a ↓R-strided convolution, the spatial convolution in Equation (C.22) leads to the following spectral representation:
Y (z) = H̃ [R](z)X [R](z) (C.24)
And for an ↑R-strided convolution, the spatial convolution is represented in spectral domain as:
Y [R](z) =H [R](z)X(z) (C.25)
Furthermore, a ↓R-strided convolution is orthogonal if and only if H̃ [R](z) is paraunitary, and an ↑R-strided convolution is orthogonal if and only ifH [R][z] is paraunitary.
Proof of Proposition C.4. (1a) ↓R-strided convolutions. We first prove the spectral representation of ↓R-strided convolution in Equation (C.24).
Y (z) = ∑ i∈Z y[i]z−i = ∑ i∈Z (∑ n∈Z h[n]x[Ri− n] ) z−i (C.26)
= ∑ i∈Z ∑ r∈[R] ∑ m∈Z h[mR− r]x[(i−m)R+ r]  z−i (C.27) = ∑ r∈[R] (∑ m∈Z h[mR− r]z−m (∑ i x[(i−m)R+ r]z−(i−m) )) (C.28) = ∑ r∈[R] (∑ m∈Z h[mR− r]z−m )(∑ i′∈Z x[i′R+ r]z−i ′ ) (C.29)
= ∑ r∈[R] H−r|R(z)Xr|R(z), (C.30)
where Equation (C.26) follows from the definitions of the ↓R-strided convolution (Equation (C.22)) and the Z-transform (Equation (B.3)), Equation (C.27) makes a change of variables n = mR − r, Equation (C.29) further changes i′ = i−m, and Equation (C.30) is due to the definition of polyphase components (Equation (C.7)). Now We rewrite the last equation concisely as
Y (z) = [ H−0|R(z); · · · ;H−(R−1)(z) ]︸ ︷︷ ︸ H̃ [R](z)
 X 0|R(z)
... XR−1|R(z)  ︸ ︷︷ ︸ X [R](z) , (C.31)
which is the spectral representation of ↓R-strided convolutions in Equation (C.24). Now we prove the orthogonality condition for ↓R-strided convolutions.
‖y‖2 = 1 2π ∫ π −π ∥∥Y (ejω)∥∥2 dω (C.32) = 1
2π ∫ π −π ∥∥∥H̃ [R](ejω)X [R](ejω)∥∥∥2 dω (C.33) = 1
2π ∫ π −π ∥∥∥X [R](ejω)∥∥∥2 dω (C.34) = ‖x‖2 , (C.35)
where Equations (C.32) and (C.35) are due to Parseval’s theorems (Theorem B.3 and Proposition C.3), Equation (C.33) follows from the spectral representation of the ↓ R-strided convolution (Equation (C.24)), and Equation (C.34) utilizes that the transfer matrix is unitary at each frequency. The “only if” part can be proved by contradiction similar to Theorem B.5.
(1b) ↑R-strided convolutions. According to Proposition C.1, the Z-transform of x↑R is X(zR). Therefore, an application of the convolution theorem (Theorem 2.1) on Equation (C.23) leads us to
Y (z) =H(z)X(zR) (C.36) Expanding Y (z) andH(z) using polyphase decomposition (Proposition C.2), we have∑
r∈[R]
Y r|R(zR)z−r = ∑ r∈[R] Hr|R(zR)z−r X(zR) (C.37) ∑ r∈[R] Y r|R(zR)z−r = ∑ r∈[R] ( Hr|R(zR)X(zR) ) z−r (C.38)
Y r|R(zR) =Hr|R(zR)X(zR), ∀r ∈ [R] (C.39) Y r|R(z) =Hr|R(z)X(z), ∀r ∈ [R], (C.40)
where Equation (C.39) is due to the uniqueness of Z-transform, and Equation (C.40) changes the variables from zR to z. Again„we can rewrite the last equation in concisely as Y 0|R(z) Y 1|R(z)
... Y R−1|R(z)  ︸ ︷︷ ︸ Y [R](z) =  H0|R(z) H1|R(z) ... HR−1|R(z)  ︸ ︷︷ ︸ H [R](z) X(z), (C.41)
which is the spectral representation of ↑R-strided convolutions in Equation (C.25). Lastly, we prove the orthogonality condition for ↑R-strided convolutions.
‖y‖2 = 1 2π ∫ π −π ∥∥∥Y [R](ejω)∥∥∥2 dω (C.42) = 1
2π ∫ π −π ∥∥∥H [R](ejω)X(ejω)∥∥∥2 dω (C.43) = 1
2π ∫ π −π ∥∥X(ejω)∥∥2 dω (C.44) = ‖x‖2 , (C.45)
where Equations (C.42) and (C.45) are due to Parseval’s theorems (Theorem B.3 and Proposition C.3), Equation (C.43) follows from the spectral representation of the ↑ R-strided convolution (Equation (C.25)), and Equation (C.44) uses the fact that the transfer matrix is unitary for each frequency. The “only if” part can be proved by contradiction similar to Theorem B.5.
(2) Dilated convolutional layer is proposed to increase the receptive field of a convolutional layer without extra parameters and computation. The layer up-samples its filter bank before convolution with the input. R-dilated convolutional layer) computes its output with the following equation:
y[i] = ∑ n∈Z h↑R[n]x[i− n] (C.46)
Proposition C.5 (Orthogonality of dilated convolutional layer). For an R-dilated convolution, the spatial convolution in Equation (C.46) leads to a spectral representation as
Y (z) =H(zR)X(z), (C.47)
Furthermore, an R-dilated convolutional layer is orthogonal if and only ifH(zR) is paraunitary.
Proof of Proposition C.5. According to Proposition C.1, the Z-transform of h↑R isH(zR). Therefore, the “if” part follows directly from the convolution theorem. The “only if” part can be proved by constructing a counterexample similar to Theorem B.5.
Notice thatH(zR) is paraunitary if and only ifH(ejω) is unitary for all frequency ω ∈ R, which is the same asXHz being paraunitary. In other words, any filter bank that is orthogonal for a standard convolution is also orthogonal for a dilated convolution and vice versa.
(3) Group convolutional layer is proposed to reduce the parameters and computations and used in many efficient architectures, including MobileNet, ShuffleNet. The layer divides both input/output channels into multiple groups and restricts the connections within each group.
Formally, a group convolutional layer with G groups (abbrev. as G-group convolutions) is parameterized by G filter banks {hg}g∈[G], each consists of (T/G)× (S/G) filters. The layer maps an S channels input x to a T channels output y according to
y[i] = ∑ n∈Z blkdiag ( {hg[n]}g∈[G] ) x[i− n], (C.48)
where blkdiag({·}) computes a block diagonal matrix from a set of matrices.
Proposition C.6 (Orthogonality of group convolutional layer). For a G-group convolution, the spatial convolution in Equation (C.48) leads a spectral representation as , their z-transforms satisfy
Y (z) = blkdiag ( {Hg(z)}g∈[G] ) X(z), (C.49)
Furthermore, a G-group convolutional layer is orthogonal if and only if the block diagonal matrix is paraunitary, i.e., each hg(z) is paraunitary.
Proof of Proposition C.6. Due to the convolution theorem, it suffices to prove that the Z-transform of a sequence of block diagonal matrices is also block diagonal in the spectral domain.
∑ n∈Z
h 0[n]
. . . hG−1[n]  ︸ ︷︷ ︸ blkdiag ( {hg[n]}g∈[G] ) z−n =  ∑ n∈Z h 0[n]z−n . . . ∑ n∈Z h G−1[n]z−n  (C.50)
= H 0(z)
. . . HG−1(z)  ︸ ︷︷ ︸ blkdiag ( {Hg(z)}g∈[G] ) . (C.51)
As a result, we can write the orthogonality condition asH 0(z)
. . . hG−1(z)
 † H 0(z) . . .
HG−1(z)  = H 0(z)†H0(z)
. . . HG−1(z)†HG−1(z)  = I, ∀|z| = 1. (C.52) The equation impliesHg(z)†Hg(z) = I,∀|z| = 1,∀g ∈ [G], i.e., eachHg(z) is paraunitary.

heading: C.3 REALIZATIONS OF PARAUNITARY SYSTEMS
text: In this subsection, we first prove that all finite-length 1D-paraunitary systems can be represented in a factorized form. Next, we show how we can construct MD-paraunitary systems using 1D systems. Lastly, we study the relationship of existing approaches to paraunitary systems.

heading: C.3.1 COMPLETE FACTORIZATION OF 1D-PARAUNITARY SYSTEMS
text: The classic theorem for spectral factorization of paraunitary systems is traditionally developed for causal systems (Vaidyanathan, 1993; Kautsky & Turcajová, 1994). Given a causal paraunitary system of length L (i.e., polynomial in z−1), there always exists a factorization such that
H(z) = QV (z−1;U (1)) · · ·V (z−1;U (L−1)), (C.53) whereQ is an orthogonal matrix, U (`) is a column-orthogonal matrix, and V (z;U) is defined as
V (z;U) = (I −UU>) +UU>z. (C.54) In Theorem C.7, we extends this theorem from causal systems to finite-length (but non-causal) ones. Theorem C.7 (Complete factorization for 1D-paraunitary systems). Suppose that a paraunitary system H(z) is finite-length, i.e., it can be written as ∑ n h[n]z
−n for some sequence {h[n], n ∈ [−L,L]}, then it can be factorized in the following form:
H(z) = V (z;U (−L)) · · ·V (z;U (−1))QV (z−1;U (1)) · · ·V (z−1;U (L)), (C.55) where Q is an orthogonal matrix, U (`) is a column-orthogonal matrix, and V (z;U) is defined in Equation (C.54). Consequently, the paraunitary system H(z) is parameterized by L + L + 1 (column-)orthogonal matricesQ and U (`)’s.
Proof for Theorem C.7. Given a non-causal paraunitary systemH(z), we can always find a causal counterpart Ĥ(z) such that H(z) = zLĤ(z) (This can be done by shifting the causal system backward by L steps, which is equivalent to multiplying zL in the spectral domain). Since the causal system Ĥ(z) admits a factorization in Equation (C.55), we can write the non-causal systemH(z) as
H(z) = zLQV (z−1; Û (1)) · · ·V (z−1; Û (L+L)). (C.56)
Therefore, it suffices to show that for an orthogonal matrixQ and any column-orthogonal matrix Û , we can always find another column-orthogonal matrix U such that
zQV (z−1; Û) = V (z;U)Q. (C.57)
If the equation above is true, we can set U (`) = Û (`−1−L) for ` < 0 and U (`) = Û (`−L) for ` > 0, which will convert Equation (C.56) into Equation (C.55).
Now we start to prove Equation (C.57). Note that any column-orthogonal Û has a complement Ū such that [Û , Ū ] is orthogonal and I = ÛÛ> + ŪŪ>. We then rewrite Equation (C.57) as
zQV (z−1; Û) = zQ(I − ÛÛ> + ÛÛ>z−1) (C.58) = Q(I − ŪŪ> + ŪŪ>z) (C.59) = (I −QŪŪ>Q> +QŪŪ>Q>z)Q (C.60) = (I −UU> +UU>z)Q (C.61) = V (z;U)Q, (C.62)
where in Equation (C.61) we set U = QŪ . This completes the proof.

heading: C.3.2 MULTI-DIMENSIONAL (MD) PARAUNITARY SYSTEMS
text: If the data are multi-dimensional (MD), we will need MD-convolutional layers in neural networks. Analogously, we can prove the equivalence between orthogonal MD-convolutions in the spatial domain and MD-paraunitary systems in the spectral domain, i.e.,
H(z)†H(z) = I, z = (z1, · · · , zD), |zd| = 1,∀d ∈ [D], (C.63)
where D is the data dimension. In this work, we adopt a parameterization based on separable systems. Definition C.8 (Separable MD-paraunitary system). A MD-paraunitary systemH(z) is separable if there exists D 1D-paraunitary systemsH1(z1), · · · ,HD(zD) such that
H(z) =H(z1, · · · , zD) ,H1(z1) · · ·HD(zD). (C.64)
Therefore, we can construct an MD-paraunitary system with D number of 1D-paraunitary systems, each of which is represented in Equation (C.55). Notice that not all MD-paraunitary systems are separable, thus the parameterization in Equation (C.64) is not complete (see Section 5 for a discussion). However, we can guarantee that our parameterization realizes all separable MD-paraunitary systems — each separable paraunitary system admits a factorization in Equation (C.64), where each 1D-system admits a factorization in Equation (C.55).
C.3.3 INTERPRETATIONS OF PREVIOUS APPROACHES
In Theorem B.5, we have shown that a paraunitary transfer matrix is both necessary and sufficient for a convolution to be orthogonal. Therefore, we can interpret all approaches for orthogonal convolutions as implicit constructions of paraunitary systems, including singular value clipping and masking (SVCM) (Sedghi et al., 2019), block convolution orthogonal parameterization (BCOP) (Li et al., 2019b), Cayley convolution (CayleyConv) (Trockman & Kolter, 2021), skew orthogonal convolution (SOC) (Singla & Feizi, 2021). Furthermore, we prove how orthogonal regularization (Wang et al., 2019; Qi et al., 2020) encourages the transfer matrix to be unitary for all frequencies.
(1) Singular value clipping and masking (SVCM) (Sedghi et al., 2019) clips all singular values ofH(ejω) to ones for each frequency ω after gradient update. Since the clipping step can arbitrarily
enlarge the filter length, SVCM subsequently masks out the coefficients outside the filter length. However, the masking step breaks the orthogonality, as we have seen in the experiments (Section 6.1).
(2) Block convolution orthogonal parameterization (BCOP) (Li et al., 2019b) tries to generalize the spectral factorization of 1D-paraunitary systems in Equation (C.55) to 2D-paraunitary systems.
H(z1, z2) = z K−1 2 1 z K−1 2 2 QV (z1, z2;U (1) 1 ,U (1) 2 ) · · ·V (z1, z2;U (K−1) 1 ,U (K−1) 2 ), (C.65)
where V (z1, z2;U (`) 1 ,U (`) 2 ) = V (z1;U (`) 1 )V (z2;U (`) 2 ). In other words, this approach makes each V -block, instead of the whole paraunitary system, separable. This factorization in BCOP is incomplete for 2D-paraunitary system — unlike 1D-paraunitary systems, not every 2D-paraunitary system admits a factorized form (Lin & Vaidyanathan, 1996).
(2) Calyey convolution (CayleyConv) (Trockman & Kolter, 2021) aims to generalize the Cayley transform for orthogonal matrices in Equation (C.88) to 2D-paraunitary systemsH(z1, z2):
H(z1, z2) = (I −A(z1, z2)) (I +A(z1, z2))−1 , (C.66)
where A(z1, z2) is a skew-Hermitian matrix for |z1| = 1, |z2| = 1 (i.e., A(ejω1 , ejω2)† = A(ejω1 , ejω2) for any ω1, ω2). Since Cayley transform cannot parameterize a matrix with singular value −1 for any frequency, the CayleyConv is not a complete parameterization.
(4) Skew orthogonal convolution (SOC) (Singla & Feizi, 2021) aims to generalize the matrix exponential for orthogonal matrices (Equation (C.90)) to convolution exponential for 2D-paraunitary systemsH(z1, z2):
H(z1, z2) = exp (A(z1, z2)) , ∞∑ k=0 A(z1, z2) k! = I +A(z1, z2) + A(z1, z2) 2 2 + · · · , (C.67)
where A(z1, z2) is skew-Hermitian matrix for any matrix for |z1| = 1, |z2| = 1 (in other words, A(ejω1 , ejω2)† = A(ejω1 , ejω2) for any ω1, ω2). It is not resolved whether all 2D-paraunitary systems can be represented in terms of convolution exponential.
(5) Orthogonal regularization (Ortho-Reg) (Wang et al., 2019; Qi et al., 2020) is developed to encourage orthogonality in convolutional layers. We show that such orthogonal regularization is equivalent to a unitary regularization of the paraunitary system with uniform weights on all frequencies.
∑ i∈Z ∥∥∥∥∥∑ n∈Z h[n]h[n−Ri]> − δ[i] ∥∥∥∥∥ 2 = 1 2π ∫ π −π ∥∥∥H [R](ejω)†H [R](ejω)− I∥∥∥2 dω, (C.68) where δ[0] = I is an identity matrix, δ[n] = 0 is a zero matrix for n 6= 0. We prove the equivalence more generally in Proposition C.9. However, this approach cannot enforce exact orthogonality and in practice requires hyperparameter search for a proper regularizer coefficient. Proposition C.9 (Parseval’s theorem for ridge regularization). Given a sequence of matrices h = {h[n], n ∈ Z}, the following four expressions are equivalent:
1
2π ∫ π −π ∥∥∥H [R](ejω)†H [R](ejω)− I∥∥∥2 dω, (C.69a) ∑ i∈Z ∥∥∥∥∥∑ n∈Z h[n]>h[n−Ri]− δ[i] ∥∥∥∥∥ 2 , (C.69b)
1
2π ∫ π −π ∥∥∥H [R](ejω)H [R](ejω)† − I∥∥∥2 dω, (C.69c) ∑ i∈Z ∥∥∥∥∥∑ n∈Z h[n]h[n−Ri]> − δ[i] ∥∥∥∥∥ 2 , (C.69d)
where ‖ · ‖ denotes the Frobenius norm of a matrix.
Proof of Proposition C.9. We first prove the equivalence between Equations (C.69a) and (C.69c).∥∥∥H [R](ejω)†H [R](ejω)− I∥∥∥2 = tr (( H [R](ejω) † H [R](ejω)− I )† ( H [R](ejω) † H [R](ejω)− I )) (C.70)
= tr ( H [R](ejω) † H [R](ejω)H [R](ejω) † H [R](ejω) ) − 2tr ( H [R](ejω) † H [R](ejω) ) + I (C.71)
= tr ( H [R](ejω)H [R](ejω) † H [R](ejω)H [R](ejω) †) − 2tr ( H [R](ejω)H [R](ejω) †) + I (C.72)
= tr (( H [R](ejω)H [R](ejω) † − I )† ( H [R](ejω)H [R](ejω) † − I )) (C.73) = ∥∥∥H [R](ejω)H [R](ejω)† − I∥∥∥2 , (C.74)
where Equations (C.70) and (C.74) make use of ‖A‖2 = tr(A†A), Equations (C.71) and (C.73) are due to the linearity of tr(·), and Equation (C.72) utilizes tr(AB) = tr(BA). Next, we prove the equivalence between Equations (C.69a) and (C.69b).
H [R](ejω) † H [R](ejω)− I = ∑ r∈[R] Hr|R(ejω) † Hr|R(ejω)− I (C.75)
= ∑ r∈[R] ∑ i∈Z (∑ m∈Z hr|R[m]>hr|R[m− i]− δ[i] ) e−jωi (C.76) = ∑ r∈[R] ∑ i∈Z (∑ m∈Z h[Rm+ r]>h[R(m− i) + r]− δ[i] ) e−jωi (C.77) = ∑ i∈Z ∑ r∈[R] ∑ m∈Z h[Rm+ r]>h[Rm+ r −Ri]− δ[i]
 e−jωi (C.78) = ∑ i∈Z (∑ n∈Z h[n]>h[n−Ri]− δ[i] ) e−jωi, (C.79)
where Equation (C.75) follows from the definition of polyphase matrix in Equation (C.15). Equation (C.76) uses a number of properties of Fourier transform: a Hermitian in the spectral domain is a transposed reflection in the spatial domain, a frequency-wise multiplication in the spectral domain is a convolution in the spatial domain, and an identical mapping in the spectral domain is an impulse sequence in the spatial domain. Equation (C.77) follows from the definition of polyphase components in Equation (C.7), and Equation (C.79) makes a change of variables n = Rm+ r. In summary, we show that the LHS (denotedD(ejω)) is a Fourier transform of the RHS (denoted as d[i]):
H [R](ejω) † H [R](ejω)− I︸ ︷︷ ︸
D(ejω)
= ∑ i∈Z (∑ n∈Z h[n]>h[n−Ri]− δ[n] ) ︸ ︷︷ ︸
d[i]
e−jωi. (C.80)
Applying Parseval’s theorem (Theorem B.3) to the sequence d = {d[i], i ∈ Z}, we have
1
2π ∫ π −π ∥∥∥H [R](ejω)†H [R](ejω)− I∥∥∥2 dω =∑ i∈Z ∥∥∥∥∥∑ n∈Z h[n]>h[n−Ri]− δ[i] ∥∥∥∥∥ 2 , (C.81)
which proves the equivalence between Equations (C.69a) and (C.69b). With almost identical arguments, we can prove the equivalence between Equations (C.69c) and (C.69d), that is
1
2π ∫ π −π ∥∥∥H [R](ejω)H [R](ejω)† − I∥∥∥2 dω =∑ i∈Z ∥∥∥∥∥∑ n∈Z h[n]h[n−Ri]> − δ[i] ∥∥∥∥∥ 2 , (C.82)
which completes the proof.
Table 5: Computational complexities of different approaches for orthogonal convolutions.
Approach Computational Complexity Normal O ( L2N2C2 ) SVCM O ( L2N2C2 +N2C3
) Ortho-Reg O ( L2N2C2 + L4C2
) CayleyConv O ( N2 log(N)C2 +N2C3
) SOC O ( KL2N2C2
) BCOP O ( L2N2C2 +KL2C3
) SC-Fac O ( L2N2C2 + L2C3
) In Table 5, we compare the computational complexities (forward pass) of orthogonal convolutions against normal convolution. For simplicity, we assume the feature maps have size N × N , the convolution filters have size L× L, and the maximum of input/output channels is C. We use K to denote the number of iterations for Björck’s algorithm in BCOP, or Taylor’s series order in SOC.
• For SC-Fac, BCOP, SVCM, or Ortho-Reg, the first term O(L2N2C2) is the base cost of a normal convolution, and the second term is the overhead for reconstruction, projection, or regularization. Note that the overhead in SC-Fac is comparable to Ortho-Reg, which is lower than SVCM or BCOP. If C < N2, the overhead in SC-Fac is negligible compared to the base cost. • For CayleyConv, the first term O(N2 log(N)C2) is the cost for fast Fourier transform (FFT), and the second term O(N2C3) is the cost for matrix inversion for all frequencies. The cost of SOC is exactly K times as a normal convolution. The computational complexities of CayleyConv and SOC are significantly higher than the one of normal convolution.
• For those approaches whose filters can be explicitly obtained (SC-Fac, BCOP, SVCM, and OrthoReg), the inference time of an orthogonal convolution is no different from a normal convolution, i.e., O(L2N2C2). On the other hand, for those approaches whose filters are implicitly defined (CayleyConv and SOC), the inference time is the same as the forward pass in training.

heading: C.4 CONSTRAINED OPTIMIZATION OVER ORTHOGONAL MATRICES
text: In Theorem C.7, we have shown how orthogonal matrices characterize paraunitary systems. Our remaining goal, therefore, is to parameterize orthogonal matrices using unconstrained parameters. In this part, we review popular parameterization methods: Householder reflections (Mhammedi et al., 2017; Mathiasen et al., 2020), Givens rotations (Dorobantu et al., 2016; Jing et al., 2017), Björck orthogonalization (Anil et al., 2019), Cayley transform (Helfrich et al., 2018; Maduranga et al., 2019), and exponential map (Lezcano-Casado & Martínez-Rubio, 2019; Lezcano Casado, 2019).
(1) Householder reflections (Mhammedi et al., 2017; Mathiasen et al., 2020). The parameterization represents an orthogonal matrix using a product of Householder matrices. Given N ∈ N and n ∈ {1, · · · , N}, we defineH(n) as a mapping from a vector v ∈ Rn (a scalar v ∈ R for n = 1) to a block-diagonal matrixH(n)(v) ∈ RN×N :
H(n)(v) =
[ IN−n
In − 2 vv >
‖v‖2
] , H(1)(v) = [ IN−1
v
] , (C.83)
where In ∈ Rn×n denotes an identity matrix. For n ≥ 2, H(n)(v) is the Householder matrix that represents the reflection with respect to the hyperplane orthogonal to concat(0N−n,v) ∈ RN and passing through the origin. For n = 1, the scalar v takes values {−1,+1}, which makes H(1)(v) either an identity matrix (for v = 1) or a Householder matrix (for v = −1). This method parameterizes an orthogonal matrix as a product ofH(n)’s:
U =H(n)(v(n)) · · ·H(2)(v(2))H(1)(v(1)), (C.84)
where each v(n) ∈ Rn is a learnable vector for n ≥ 2, and v(1) is a fixed constant that is generated at initialization and fixed afterward. Observing that Equation (C.84) is serial (unfriendly to parallel computing), Mathiasen et al. (2020) proposes to increase its parallelism using WY transform.
(2) Givens rotations (Dorobantu et al., 2016; Jing et al., 2017). The parameterization represents a special orthogonal matrix using a product of Given rotations matrices. Given N ∈ N and i, j ∈ {1, · · · , N} with i 6= j, we define G(i,j) as a mapping from an angle θ ∈ R to the corresponding rotation matrixG(i,j) ∈ RN×N :
G(i,j)(θ) =  1 · · · 0 · · · 0 · · · 0 ... . . . ... ... ... 0 · · · cos θ · · · − sin θ · · · 0 ... ... . . . ... ... 0 · · · sin θ · · · cos θ · · · 0 ... ... ... . . . ...
0 · · · 0 · · · 0 · · · 1
 . (C.85)
This method parameterizes a special orthogonal matrix U as a product ofG(i,j)’s:
U = G(0,1)(θ(0,1))G(0,2)(θ(0,2))G(1,2)(θ(1,2)) · · ·G(N−2,N−1)(θ(N−2,N−1)), (C.86) where θ(i,j)’s are N(N − 1) unconstrained parameters. Since the determinant of each rotation matrix is +1, their product U also has +1 determinant, thus is a special orthogonal matrix. Again, Equation (C.86) is highly serial, thus both Dorobantu et al. (2016) and Jing et al. (2017) propose to group N/2 Given rotations into a packed rotation to increase the parallelism.
(2) Björck orthogonalization (Anil et al., 2019). The algorithm was first introduced in Björck & Bowie (1971) to compute the closest orthogonal matrix of a given matrix. Given an initial matrix U0, this algorithm iteratively approaches its closest orthogonal matrix as:
Uk+1 = Uk
( I + 1
2 Pk + · · ·+ (−1)p ( − 12 p ) P pk ) ,∀k ∈ [K], (C.87)
where K is the iterative steps, Pk = I −U>k Uk, and p controls the trade-off between efficiency and accuracy at each step. When the algorithm is used for parameterization, it maps an unconstrained matrix U0 to an approximately orthogonal matrix UK in K steps. Although Björck parameterization is complete (since any orthogonal matrixQ can be represented by U0 = Q), it is inexact due to the iterative approximation.
(3) Cayley transform (Helfrich et al., 2018; Maduranga et al., 2019). The transform provides a bijective parameterization of orthogonal matrices without −1 eigenvalue with skew-symmetric matrices (i.e.,A> = −A) U = (I −A)(I +A)−1, (C.88) where the skew-symmetric matrixA is represented by its upper-triangle entries. Since orthogonal matrices with −1 eigenvalue are out of consideration, the parameterization is incomplete. Helfrich et al. (2018); Maduranga et al. (2019) overcome this difficulty by a scaled Cayley transform: U =D(I −A)(I +A)−1, (C.89) whereD is a diagonal matrix with±1 non-zero entries, which is (randomly) generated at initialization and fixed during training.
(4) Exponential map (Lezcano-Casado & Martínez-Rubio, 2019; Lezcano Casado, 2019). The mapping provides a surjective parameterization of all special orthogonal matrices (with +1 determinant) with using skew-symmetry matrices (i.e.,A> = −A).
U = exp(A) , ∞∑ k=0 Ak k! = I +A+ 1 2 A2 + · · · , (C.90)
where the infinite sum can be computed exactly up to machine-precision (Higham, 2009). LezcanoCasado & Martínez-Rubio (2019) derives an efficient backpropagation algorithm for the mapping, making it an exact and efficient parameterization in neural networks. To support a complete parameterization for all orthogonal matrices, Lezcano Casado (2019) extends the mapping as:
U = Q exp(A), (C.91) whereQ is an orthogonal matrix, which is generated at initialization and fixed during training.
In principle, we can use any of these approaches to parameterize the orthogonal matrices. In this work, we choose exponential map due to its exactness, efficiency, and completeness.

heading: D LEARNING DEEP ORTHOGONAL NETWORKS WITH LIPSCHITZ BOUNDS
text: In this section, we first discuss the properties of GroupSort and Lipschitz networks. Subsequently, we prove Proposition 4.1, which exhibits two approaches to construct Lipschitz residual blocks. Lastly, we prove in Proposition D.1 when a paraunitary system (represented by a complete factorization as in Theorem C.7) reduces to an orthogonal matrix. The reduction allows us to apply the initialization methods for orthogonal matrices to paraunitary systems.
GroupSort and orthogonality. The GroupSort activation separates inputs into groups and sorts each group into ascending order (Anil et al., 2019). It guarantees two properties: (1) The activation is norm preserving in the forward pass — a sorting function does not change the norm of any input vector. (2) The activation is gradient norm preserving in the backward pass — a sorting function acts as a permutation, which is orthogonal. GroupSort with the group size of two is a spacial case that we use in the paper followed by (Chernodub & Nowicki, 2016; Anil et al., 2019).
A Lipschitz network with orthogonal layers and GroupSort activations is locally orthogonal. Given any input x to the network, there exists a neighborhood N(x, r) with radius r such that the sorting order in each activation does not change. In this neighborhoodN(x, r), each GroupSort layer operates as a constant permutation (thus orthogonal); consequently, the whole network is locally orthogonal.
Lipschitz residual blocks. In Proposition 4.1, we prove the Lipschitzness of two types of residual blocks, one based on additive skip-connection and another based on concatenative one (See Figure 6).
Proof for Proposition 4.1. We first prove the Lipschitzness for the additive residual block f : f(x) , αf1(x) + (1− α)f2(x). Let x, x′ be two inputs to f and f(x), f(x′) be their outputs, we have
‖f(x′)− f(x)‖ = ∥∥(αf1(x′) + (1− α)f2(x′))− (αf1(x) + (1− α)f2(x))∥∥ (D.1)
= ∥∥α (f1(x′)− f1(x))+ (1− α) (f2(x′)− f2(x))∥∥ (D.2)
≤ α ∥∥f1(x′)− f1(x)∥∥+ (1− α)∥∥f2(x′)− f2(x)∥∥ (D.3)
≤ αL ‖x′ − x‖+ (1− α)L ‖x′ − x‖ (D.4) = L‖x′ − x‖, (D.5)
where Equation (D.3) makes uses of the triangle inequality, and Equation (D.4) is due to the LLipschitzness of both f1, f2. Therefore, we have shown that ‖f(x′)− f(x)‖ ≤ L‖x′ − x‖. Similarly, we prove the Lipschitzness for the concatenative residual block g : g(x) , P [ g1(x1); g2(x2) ] . Let x, x′ be two inputs to g and g(x), g(x′) be their outputs, we have
‖g(x′)− g(x)‖2 = ∥∥∥P ([g1(x1′); g2(x2′)]− [g1(x1); g2(x2)])∥∥∥ (D.6)
= ∥∥[g1(x′); g2(x′)]− [g1(x); g2(x)]∥∥2 (D.7)
= ∥∥∥g1(x1′)− g1(x)∥∥∥2 + ∥∥∥g2(x2′)− g2(x2)∥∥∥2 (D.8)
≤ L2 ∥∥∥x1′ − x1∥∥∥2 + L2 ∥∥∥x2′ − x2∥∥∥2 (D.9)
= L2 ∥∥∥[x1′;x2′]− [x1;x2]∥∥∥2 (D.10)
= L2 ‖x′ − x‖2 , (D.11) where Equation (D.7) utilizes ‖Px‖ = ‖x‖,∀x, and Equation (D.9) is due to the L-Lipschitzness of g1, g2. The equations above implies that ‖g(x′)− g(x)‖ ≤ L‖x′ − x‖.
Reduction of a paraunitary system to an orthogonal matrix. In Proposition D.1, we prove a special case when a paraunitary system reduces to an orthogonal matrix. The reduction allows us to apply the initialization methods for orthogonal matrices to paraunitary systems. Proposition D.1 (Reduction of a paraunitary matrix to an orthogonal matrix). Suppose a paraunitary systemH(z) takes the complete factorization in Equation (C.55), and assume L = L with U (−`) = QU (`) for all `, then the paraunitary matrixH(z) reduces to an orthogonal matrixQ,
H(z) = V (z;U (−L)) · · ·V (z;U (−1))QV (z−1;U (1)) · · ·V (z−1;U (L)) = Q. (D.12)
Proof for Proposition D.1. In order to prove Equation (D.12), it suffice to show that
V (z;U (−`))QV (z−1;U (`)) = Q, (D.13)
and Equation (D.12) will reduce recursively to the orthogonal matrixQ. For simplicity, we rewrite U (−`) as L and U (`) as R, by which we have L = QR (or R = Q>L) and we aim to prove V (z;L)QV (z−1;R) = Q. By the definition of V (z; ·) in Equation (C.54), we expand it as
V (z;L)QV (z−1;R) = [( I −LL> ) +LL>z ] Q [( I −RR> ) +RR>z−1 ] =
LL>Q(I −RR>)︸ ︷︷ ︸ c[−1] z + (I −LL>)Q(I −RR>) +LL>QRR>︸ ︷︷ ︸ c[0] + (I −LL>)QRR>︸ ︷︷ ︸ c[1] z−1 (D.14)
Therefore, we will need to show that c[−1] = 0, c[1] = 0 and c[0] = Q. We first show that both c[−1] for z and c[1] for z−1 are zero matrices.
c[−1] = LL>Q ( I −RR> ) (D.15)
= LL>Q−LL>QRR> (D.16) = L(Q>L)> −L(Q>L)>RR> (D.17) = LR> −L(R>R)R> (D.18) = LR> −LR> = 0, (D.19)
c[1] = (I −LL>)QRR> (D.20) = QRR> −LL>QRR> (D.21) = (QR)R> −LL>(QR)R> (D.22) = LR> −L(L>L)R> (D.23) = LR> −LR> = 0. (D.24)
Lastly, we show that the constant coefficient c[0] is equal toQ.
c[0] = (I −LL>)Q(I −RR>) +LL>QRR> (D.25) = Q−LL>Q−QRR> + 2LL>QRR> (D.26) = Q−LR> −LR> + 2LR> = Q (D.27)
which completes the proof.

heading: E SUPPLEMENTARY MATERIALS FOR EXPERIMENTS

heading: E.1 EXPERIMENTAL SETUP
text: Network architectures. For fair comparisons, we follow the architectures by Trockman & Kolter (2021) for KW-Large, ResNet9, WideResNet10-10 (i.e., shallow networks). We set the group size for GroupSort activations as 2 in all experiments. For networks deeper than 10 layers, we implement their architectures modifying from the Pytorch official implementation of ResNet. It is crucial to replace the global pooling before fully-connected layers with an average pooling with a window size of 4. For the average pooling, we multiply the output with the window size to maintain its 1-Lipschitzness. Other architectures, including ShuffleNet and plain convolutional network (ConvNet), are further modified from the ResNet, where only the skip-connections are changed or removed. We use the widen factor to indicate the channel number: we set the number of channels at each layer as base channels multiplied by the widen factor. The base channels are 16, 32, 64 for three groups of residual blocks. More details of the ResNet architecture can be found in the official Pytorch implementation.1
Learning strategies. We use the CIFAR-10 dataset for all our experiments. We normalize all input images to [0, 1] followed by standard augmentation, including random cropping and horizontal flipping. We use the Adam optimizer with a maximum learning rate of 10−2 coupled with a piece-wise triangular learning rate scheduler. We initialize all our SC-Fac layers as permutation matrices: (1) we select the number of columns for each pair U (`),U (−`) uniformly from {1, · · · , T} at initialization (the number is fixed during training); (2) for ` > 0, we sample the entries in U (`) uniformly with respect to the Haar measure; (3) for ` < 0, we set U (−`) = QU (`) according to Proposition D.1.

heading: E.2 ADDITIONAL EMPIRICAL RESULTS ON HYPER-PARAMETERS SELECTION
text: Multi-class hinge loss. Following previous works on Lipschitz networks (Anil et al., 2019; Li et al., 2019b; Trockman & Kolter, 2021), we adopt the multi-class hinge loss in training. For each model, we perform a grid search on different margins 0 ∈ {1× 10−3, 2× 10−3, 5× 10−3, 1× 10−2, 2× 10−2, 5× 10−2, 0.1, 0.2, 0.5} and report the best performance in terms of robust accuracy. Notice that the margin 0 controls the trade-off between clean and robust accuracy, as shown in Figure 7.
Initialization methods. In Proposition D.1, we have shown how to initialize our orthogonal convolutional layers as orthogonal matrices. In Table 6, we perform a study on different initialization methods, including identical, permutation, uniform, and torus (Henaff et al., 2016; Helfrich et al., 2018). We find that permutation works the best for WideResNet22-10, while all methods are similar in shallower WideResNet10-10. Therefore, we use permutation initialization for all other experiments.
1 https://github.com/pytorch/vision/blob/master/torchvision/models/
Network depth and width. Exact orthogonality is criticized for harming the expressive power of neural networks, and we find that increasing network depth/width can partially compensate for such loss. In Table 7, we perform a study on the impact of network depth/width on the predictive performance. As shown, deeper/wider architectures consistently improve both the clean and robust accuracy for our implementation. However, the best robust accuracy is achieved by a 22-layer network since we can afford a wide architecture for 34-layer architecture.
Comparison against normal convolutional networks. In Table 8, We perform a comparison between our orthogonal networks and normal convolutional networks. Their architecture are identical except for the activation function (GroupSort for ours and ReLU for normal convolutional networks). Since batch normalization is common in normal convolutional networks but not in Lipschitz networks,
we provide both results for normal convolutional networks with or without batch normalization. In the table, we report the clean/robust accuracy, train time for epoch, and inference time for the test set.
Robustness against `∞ attacks using adversarial training. Since orthogonality only guarantees `2 Lipschitzness, Lipschitz networks with orthogonal layers are not naturally robust to `∞ perturbations. To further guard Lipschitz networks against `∞ attacks, we follow the approach of adversarial training in Wong et al. (2020). For training, we use a FGSM variant with step size 10/255; for evaluation, we use 50 PGD iterations with step size 2/255 and 10 random restarts. We report the experimental results in Table 9. We observe that different orthogonal convolution methods achieve similar `∞ robustness on WideResNet-22. Furthermore, the Lipschitz networks with WideResNet22 architecture is consistently better than ResNet-18, which is previously used in Singla & Feizi (2021).

heading: E.3 ON THE NECESSITY OF EXACT ORTHOGONALITY
text: Due to the benefits like generalizability and robustness, achieving exact orthogonality in convolutions is the primary goal of a current research line (Sedghi et al., 2019; Li et al., 2019b; Trockman & Kolter, 2021; Singla & Feizi, 2021). However, until our work, no previous approach achieves orthogonality up to machine precision. Therefore, our proposed method serves as an extreme case in gaining insight into the trade-off between orthogonality and expressiveness.
In Section 6, we have seen that exact orthogonality is not critical in shallow Lipschitz networks for robustness, and various orthogonal convolutions (with different precision) achieve comparable results. However, our method is more favorable in deeper networks (with more than 10 layers) — we show the results in Table 7. It indicates that exact orthogonality is crucial in learning deep Lipschitz networks. However, without our implementation of exact orthogonality (which does not exist before), it is unclear whether exact orthogonality up to machine precision is needed in Lipschitz networks.
Moreover, exact orthogonality is essential for other important and timely applications. For example, reversible networks/normalizing flows (Kingma & Dhariwal, 2018; Van Den Berg et al., 2018) require exact orthogonality to compute inverse transform and determinants accurately.

heading: F ORTHOGONAL CONVOLUTIONS FOR RESIDUAL FLOWS
text: In this section, we first review the class of flow-based generative models (Papamakarios et al., 2019; Kobyzev et al., 2020). We focus on invertible residual network (Behrmann et al., 2019), a flow-based model that relies on Lipschitz residual block, and its extended version Residual Flow (Chen et al., 2019). We then show how to construct improved Residual Flow using our orthogonal convolutions.
Flow-based models. Given an observable vector x ∈ RD and a latent vector z ∈ RD, we define a bijective mapping f : RD → RD from the latent vector z to an observation x = f(z). We further define the inverse of f as F = f−1, with which we represent the likelihood of x by the one of z as:
ln pX(x) = ln pZ(z) + ln|detJF (x)|, (F.1)
where pX is the data distribution, pZ is the base distribution (usually a normal distribution), and JF (x) is the Jacobian of F at x. In practice, the bijective mapping f is composed by a sequence of K bijective mapping such that f = fK ◦ · · · ◦ f1, where each fk is named as a flow. Since the inverse mapping F = F1 ◦ · · ·FK transforms the data distribution pX into a normal distribution pZ , flow-based models are also known as normalizing flows. Accordingly, we rewrite Equation (F.1) as:
ln pX(x) = ln pZ(z) + K∑ k=1 ln|detJFk(x)|, (F.2)
In a practical flow-based model, we require efficient computations of (a) each bijective mapping fk, (b) its inverse mapping Fk = f−1k , and (c) the corresponding log-determinant ln|detJF (·)|.
Invertible residual networks (i-ResNets). Behrmann et al. (2019) proposes a flow-based model based on residual network (ResNet). Note that a block in ResNet is defined as F (x) = x+ g(x), where g is a convolutional network. In (Behrmann et al., 2019), the authors prove that F is a bijective mapping if g is 1-Lipschitz, and its inverse mapping can be computed by fixed-point iterations:
xk+1 = y − g(xk), (F.3)
where y = g(x) is the output of F and the initialization of the iterative algorithm is x0 := y. From the Banach fixed-point theorem, we have
‖x− xk‖2 = Lip(g)k
1− Lip(g) ‖x1 − x0‖ , (F.4)
i.e., the convergence rate is exponential in the number of iterations and smaller Lipschitz constant will yield faster convergence. Furthermore, the log-determinant can be computed as:
ln pX(x) = ln pZ(z) + tr (ln (I + Jg(x))) (F.5)
= ln pZ(z) + tr ( ∞∑ k=1 (−1)k+1 k [Jg(x)] k ) , (F.6)
where the infinite sum is approximated by truncation and the trace is efficiently estimated using the Hutchinson trace estimator tr(A) = Ev∼N (0,I)[v>Av].
To constrain the Lipschitz constant, i-ResNet uses spectral normalization on each linear layer in the block. Moreover, to improve optimization stability, i-ResNet changes the activation function from ReLU to ELU, ensuring nonlinear activations have continous derivatives.
As summarized in the conclusion of Behrmann et al. (2019), there are two remaining problems in this model: (1) The estimator of the log-determinant is biased and inefficient; (2) Designing and learning networks with a Lipschitz constraint are challenging — one needs to constrain each linear layer in the block instead of being able to control the Lipschitz constant of a block.
Residual Flow. Chen et al. (2019) address problem {(1) by proposing an unbiased Russian roulette estimator for Equation (F.6):
tr (ln (I + Jg(x))) = En,v [ n∑ k=1 (−1)k k v [ Jg(x)k ] v P(N ≥ k) ] , (F.7)
where n ∼ p(N) and v ∼ N (0, I). Residual Flow further changes the activation function from ELU to LipSwish. The LipSwich activation avoids derivative saturation, which occurs when the second derivative is zero in a large region. However, problem (2) remains unresolved.
Residual flows with orthogonal convolutions. We propose to address problem (2) by replacing the spectral normalized layers by our orthogonal convolutional layers (SC-Fac). Note that orthogonal convolutions directly control the Lipschitz constant of a ResNet block. We keep all other components unchanged — in particular, we use LipSwish activation instead of GroupSort, as GroupSort suffers from derivative saturation. We experiment our model on MNIST dataset. As shown in Table 10, our model substantially improve the performance over the original Residual Flow. We display some images generated by our model in Figure 8.

Review summaries:
This paper presents a method for enforcing strict orthogonality of convolutional layers, by means of a factorization in the spectral domain. It shows that the technique can be extended to dealing with practical convolutions with strides, dilations and groups. The experiments demonstrated the superior performance in terms of adversarial robustness.

This paper proposed a theoretical framework for orthogonal convolutional layers based on the equivalency of orthogonal convolution in the spatial domain and the paraunitary systems in the spectral domain. The proposed method parametrizes the orthogonal convolution layer as compositions of multiple convolutions in the spatial domain, resulting in exact orthogonality. The layers are also more memory and computationally efficient than most previous methods. 

This paper suggests a framework for designing orthogonal convolutions using paraunitary systems. The proposed framework studies orthogonalization of different convolution operations. The proposed methods were examined on several datasets in comparison with state-of-the-art orthogonalization methods.