{'fraction_non_alphanumeric': 0.04491683025105899, 'fraction_numerical': 0.01242380411199504, 'mean_word_length': 4.194753790419965, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Orthogonal Deep Neural Networks

In this paper, we introduce the algorithms of Orthogonal Deep Neural Networks (OrthDNNs) to connect with recent interest of spectrally regularized deep learning methods. OrthDNNs are theoretically motivated by generalization analysis of modern DNNs, with the aim to find solution properties of network weights that guarantee better generalization. To this end, we first prove that DNNs are of local isometry on data distributions of practical interest; by using a new covering of the sample space and introducing the local isometry property of DNNs into generalization analysis, we establish a new generalization error bound that is both scale- and range-sensitive to singular value spectrum of each of networks' weight matrices. We prove that the optimal bound w.r.t. the degree of isometry is attained when each weight matrix has a spectrum of equal singular values, among which orthogonal weight matrix or a non-square one with orthonormal rows or columns is the most straightforward choice, suggesting the algorithms of OrthDNNs. We present both algorithms of strict and approximate OrthDNNs, and for the later ones we propose a simple yet effective algorithm called Singular Value Bounding (SVB), which performs as well as strict OrthDNNs, but at a much lower computational cost. We also propose Bounded Batch Normalization (BBN) to make compatible use of batch normalization with OrthDNNs. We conduct extensive comparative studies by using modern architectures on benchmark image classification. Experiments show the efficacy of OrthDNNs.


INTRODUCTION
Deep learning or deep neural networks (DNNs) have been achieving great success on many machine learning tasks, with image classification [45] as one of the prominent examples. Key design that supports success of deep learning can date at least back to Neocognitron [15] and Convolutional neural networks (CNNs) [36], which employ hierarchial, compositional design to facilitate learning target functions that approximately capture statistical properties of natural signals. Modern DNNs are usually over-parameterized and have very high model capacities, yet practically meaningful solutions can be obtained via simple back-propagation training of stochastic gradient descent (SGD) [37], where regularization methods such as early stopping, weight decay, and data augmentation are commonly used to alleviate the issue of overfitting.
Over the years, new technical innovations have been introduced to improve DNNs in terms of architectural design [20], [25], optimization [12], [17], [32], and also regularization [23], [27], which altogether make efficient and effective training of extremely over-parameterized models possible. While many of these innovations are empirically proposed, some of them are justified by subsequent theoretical studies that explain their practical effectiveness. For example, dropout training [23] is explained as an approximate regularization of adaptive weight decay in [3], [52]. Theoretically characterizing global optimality conditions of DNNs are also presented in [29], [58].
The above optimization and regularization methods aim to explain and address the generic difficulties of training * indicates equal contribution.
DNNs, and to improve efficient use of network parameters; they do not have designs on properties of solutions to which network training should converge. In contrast, there exist other deep learning methods that have favored solution properties of network parameters, and expect such properties to guarantee good generalization at inference time. In this work, we specially focus on DNN methods that impose explicit regularization on weight matrices of network layers [48], [53]. For example, Sokolic et al. [48] propose by theoretical analysis a soft regularizer that penalizes spectral norms of weight matrices. More recently, methods that regularize the whole spectrum of singular values and its range for each of networks' weight matrices are also proposed [2], [4], [11], [28], [54], [55]. They achieve clearly improved performance over those without imposing such a regularization. However, many of these methods are empirically motivated, with no theoretical justification on its effect on generalization. We aim to study this theoretical issue in this work.
Motivated by geometric intuitions from isometric mappings [26], we introduce a term of local isometry into the framework of generalization analysis via algorithmic robustness [57]. We use an intuitive and also formal definition of instance-wise variation space to characterize data distributions of practical interest, and prove that DNNs are of local isometry on such data distributions. More specifically, we prove that for a DNN trained on such a data distribution, a covering based on a linear partition (induced by the DNN) of the input space can be found such that DNN is locally linear in each covering ball, where we give bound on the diameters of covering balls in terms of spectral norms of the DNN's weight matrices. Based on a further proof that for a mapping induced by a linear DNN, degree of isometry is fully controlled by singular value spectrum of each of its weight matrices, we establish our generalization error (GE) bound for (nonlinear) DNNs, and show that it is both scaleand range-sensitive to singular value spectrum of each of their weight matrices. An illustration of our proofs is given in fig. 1.

arXiv:1905.05929v1 [cs.LG] 15 May 2019
Derivation of our bound is based on a new covering of the sample space, as illustrated in fig. 2, which enables explicit characterization of GEs caused by both the distance expansion and distance contraction of locally isometric mappings.
To attain an optimal GE bound w.r.t. the degree of isometry, we prove that the optimum is achieved when each weight matrix of a DNN has a spectrum of equal singular values, among which orthogonal weight matrix or a nonsquare one with orthonormal rows or columns is the most straightforward choice, suggesting the algorithms of Orthogonal Deep Neural Networks (OrthDNNs). Training to obtain a strict OrthDNN amounts to optimizing the weight matrices over their respective Stiefel manifolds, which, however, is very costly for large-sized DNNs. To achieve efficient learning, we propose a simple yet effective algorithm of approximate OrthDNNs called Singular Value Bounding (SVB). SVB periodically bounds, in the SGD based training iterations, all singular values of each weight matrix in a narrow band around the value of 1, thus achieving near orthogonality (row-or column-wise orthonormality) of weight matrices. In this work, we also discuss alternative schemes of soft regularization [11], [55] to achieve approximate OrthDNNs, and compare with our proposed SVB. Batch Normalization (BN) [27] is commonly used in modern DNNs, yet it has a potential risk of ill-conditioned layer transform, making it incompatible with OrthDNNs. We propose Degenerate Batch Normalization (DBN) and Bounded Batch Normalization (BBN) to remove such a potential risk, and to enable its use with strict and approximate OrthDNNs respectively.
To investigate the efficacy of OrthDNNs, we conduct extensive experiments of benchmark image classification [34], [45] on modern architectures [21], [25], [47], [56], [59]. These experiments show that OrthDNNs consistently improve generalization by providing regularization to training of these architectures. Interestingly, approximate Or-thDNNs perform as well as strict ones, but at a much lower computational cost. For approximate OrthDNNs, we also compare hard regularization via our proposed SVB and BBN with the alternative of soft regularization; our method consistently outperforms the alternative on modern architectures. In some of these studies, we investigate behaviors of our method under learning regimes from small to large sizes of training samples; results confirm the empirical strength of our method, especially for learning problems of smaller sample sizes. We also investigate robustness of our method against corruptions that are commonly encountered in natural images; our results demonstrate better robustness against such corruptions, and the robustness stands gracefully with increase of corruption severity levels.

Generalization analysis of DNNs
Classical theories of DNNs show that they are universal approximators [5], [24]. However, recent results from Zhang et al. [60] show an apparent puzzle that over-parameterized DNNs are able to shatter randomly labeled training data, suggesting worst-case generalization since test performance can only be at a chance level, while at the same time they perform well on practical learning tasks (e.g., ImageNet classification); the puzzle suggests that traditional analysis of data-independent generalization does not readily apply. They further conjecture [61] that over-parameterized DNNs, when trained via SGD, tend to find local solutions that fall in, with high probability, flat regions in the highdimensional solution space, which is even obvious when learning tasks are on natural signals; flat-region solutions may imply robustness in the parameter space of DNNs, which may further implies robustness in the input data space. Similar argument of flat-region solutions is also presented in [31], although Dinh et al. [13] argue that these flat minima can be equivalently converted as sharp minima without affecting network prediction. Generalization of DNNs is also explained by stochastic optimization. In [19], the notion of uniform stability [9] is extended to characterize the randomness of SGD, and a generalization bound in expectation is established for learning with SGD. The distribution-free stability bound of [19] is improved in [35] via the notion of on-average stability, revealing data-dependent behavior of SGD. To understand practical generalization of DNNs, Kawaguchi et al. [30] argue that independent of the hypothesis set and algorithms used, the learned model itself, possibly selected via a validation set, is the most important factor that accounts for good generalization; a generalization bound w.r.t. validation error is also presented in [30].
To further characterize generalization of DNNs with their weight matrices, Sokolic et al. [48] study DNNs as robust large-margin classifiers via the algorithmic robustness framework [57]. They introduce a notion of average Jacobian, and use spectral norm of the Jacobian matrix to locally bound the distance expansion from the input to the output space of a DNN; spectral norm of the Jacobian is further relaxed as the product of spectral norms of the network's weight matrices, which is used to establish the robustness based generalization bound. Bartlett et al. [6] use a scalesensitive measure of complexity to establish a generalization bound. They derive a margin-normalized spectral complexity, i.e., the product of spectral norms of weight matrices divided by the margin, via covering number approximation of Rademacher complexity; they further show empirically that such a bound is task-dependent, suggesting that SGD training learns parameters of a DNN whose complexity scales with the difficulty of the learning task.
While both of our bound and that of [48] are developed under the framework of algorithmic robustness [57], our bound is controlled by the whole spectrum of singular values, rather than spectral norm (i.e., the largest singular value) of each of the network's weight matrices, by introducing a term of local isometry into the framework. This also means that in contrast to [6], our bound is both scale-and range-sensitive to singular values of weight matrices. The fact that our bound is scale-sensitive in the sense of [6] implies that for difficult learning tasks, e.g., randomly labeled CIFAR10 [60], spectral norms of weight matrices would go extremely large, causing that diameters of covering balls go extremely small and correspondingly the second term of our bound (cf. theorem 3.2) that characterizes distribution mismatch between training and test samples dominates, and that the bound becomes vacuous. In contrast, for learning tasks of practical interest, e.g., standard CIFAR10 [34], the spectra of singular values of weight matrices are potentially in a benign range, and the bound is of practical use to inspire design of improved learning algorithms.

Compositional computations and isometries of DNNs
Montúfar et al. [41] characterize complexity of functions computable by DNNs and establish a lower bound on the maximal number of linear regions into which a DNN (with ReLU activation) can partition the input space, where the bound is derived by compositional replication of layer-wise space partitioning and grows exponentially with depth of the DNN. Our derivation of the analytic form of regionwise linear mapping (cf. lemma 3.2) borrows ideas from [41]. Similar compositional derivations for the number of computational paths from the network input to a hidden unit are also presented in [29], [30].
A property of dynamic isometry is studied in [46] to understand learning dynamics of deep linear networks. Pennington et al. [44] extend such studies to DNNs by employing powerful tools from free probability theory; they show that with orthogonal weight initialization, sigmoid activation function can preserve the max singular value to be 1 as layers go deeper, and isometry of DNNs can preserve for a large amount of time during training.
Geometric intuition of isometric mappings has been introduced to improve robustness of deep feature transformation [26], where DNNs are studied as a form of transformation functions. However, their development of robustness bound only uses explicitly the distance expansion constraint of isometric mappings; moreover, their studies are in the context of metric learning and for DNNs, they stay on a general function form, with no indications on how layerwise weight matrices affect generalization.

Regularization on weight matrices
Wang et al. [53] propose Extended Data Jacobian Matrix (EDJM) as a network analyzing tool, and study how the spectrum of EDJM affects performance of different networks of varying depths, architectures, and training methods. Based on these observations, they propose a spectral soft regularizer that encourages major singular values of EDJM to be closer to the largest one (practically implemented on weight matrix of each layer). As discussed above, a related notion of average Jacobian is used in [48] to motivate a soft regularizer that penalizes spectral norms of weight matrices.
There exist other recent methods [4], [11], [28], [55] that improve empirical performance of DNNs by regularizing the whole spectrum of singular values for each of networks' weight matrices. This is implemented in [11], [55] as soft regularizers that encourage the product between each weight matrix and its transpose to be close to an identity one. Different from [11], [55], we propose a hard regularization method termed Singular Value Bounding (SVB), which periodically bounds in the training process all singular values of each weight matrix in a narrow band around the value of 1, so that orthonormality of rows or columns of weight matrices can be approximately achieved.

Contributions
There exists a growing recent interest on using spectral regularization to improve training of DNNs. These methods impose explicit regularization on weight matrices of network layers by penalizing either their spectral norms [48], [53] or the whole spectrums of their singular values [4], [11], [28], [55]. The SVB algorithm proposed in our preliminary work [28] is among the later approach. Most of these methods are empirically motivated with no theoretical guarantees. In the present paper, we focus on theoretical analysis of these methods from the perspective of generalization analysis, and prove a novel GE bound for data distributions of practical interest. We also intensively compare empirical performance of these methods, and present their empirical strengths under various learning scenarios. We summarize our technical contributions as follows.
• We present in this paper a new generalization error bound for DNNs. We first prove that DNNs are of local isometry on data distributions of practical interest, where the degree of isometry is fully controlled by singular value spectrum of each of their weight matrices. By using a new covering of the sample space and introducing the local isometry property of DNNs into an algorithmic robustness framework, we establish our GE bound and show that it is both scaleand range-sensitive to singular value spectrum of each of networks' weight matrices. • We prove that the optimal bound w.r.t. the degree of isometry is attained when each weight matrix of a DNN has a spectrum of equal singular values, among which orthogonal weight matrix or a non-square one with orthonormal rows or columns is the most straightforward choice, suggesting the algorithms of Orthogonal Deep Neural Networks (Or-thDNNs). In this paper, we also present the algorithmic details of OrthDNNs.
• To address the heavy computation of strict Or-thDNNs, we propose a novel algorithm called Singular Value Bounding (SVB), which achieves approximate OrthDNNs via a simple scheme of hard regularization. We discuss alternative schemes of soft regularization, and compare with our proposed SVB. Batch normalization has a potential risk of ill-conditioned layer transform, making it incompatible with Or-thDNNs. We propose Degenerate Batch Normalization (DBN) and Bounded Batch Normalization (BBN) to remove such a potential risk, and to enable its use with strict and approximate OrthDNNs.

PROBLEM STATEMENT
We start by describing the formalism of classification problems that jointly learn a representation and a classifier, e.g., via Deep Neural Networks (DNNs).

The classification-representation-learning problem and its generalization error
Assume a sample space Z = X × Y, where X is the instance space and Y is the label space. We restrict ourselves to classification problems in this paper, and have x ∈ X as vectors in R n and y ∈ Y as a positive integer less than |Y| ∈ N. We use S m = {s i = (x i , y i )} m i=1 to denote the training set of size m whose examples are drawn independent and identically distributed (i.i.d.) according to an unknown distribution P .
We also denote S . Given a loss function L, the goal of learning is to identify a function f Sm : X → Y in a hypothesis space (a class F of functions) that minimizes the expected risk where z = (x, y) ∈ Z is sampled i.i.d. according to P . Since P is unknown, the observable quantity serving as a proxy to the expected risk R(f ) is the empirical risk One of the primary goals in statistical learning theory is to characterize the discrepancy between R(f Sm ) and R m (f Sm ), which is termed as generalization error -it is sometimes termed as generalization gap in the literature In this paper, we are interested in using DNNs to solve classification problems. It amounts to learning a map T , which extracts feature characteristic to a classification task, and minimizing R m simultaneously. We denote classification with this approach as a Classification-Representation-Learning (CRL) problem. We single out the map T because most of the theoretical analysis in this paper resolves around it. Rewriting the two risks by incorporating a map T (we write T when it is instantiated by a DNN), we have

Generalization analysis for robust algorithms with isometric mapping
The upper bounds of GE are generally established by leveraging on certain measures related to the capacity of hypothesis space F, such as Rademacher complexity and VCdimension [40]. These complexity measures capture global properties of F; however, GE bounds based on them ignore the specifically used learning algorithms. To establish a finer bound, one may resort to algorithm-dependent analysis [39], [57]. Our analysis of GE bound in this work is based on the algorithmic robustness framework [57] that has the advantage of conveying information of local geometry. We begin with the definition of robustness used in [57].
Definition 1 ((K, (·))-robustness). An algorithm is (K, (·))robust for K ∈ N and (·) : Z m → R, if Z can be partitioned into K disjoint sets, denoted by C = {C k } K k=1 , such that the following holds for all s i = (x i , y i ) ∈ S m , z = (x, y) ∈ Z, C k ∈ C: The gist of the definition is to constrain the variation of loss values on test examples w.r.t. those of training ones through local property of the algorithmically learned function. Intuitively, if s ∈ S m and z ∈ Z are "close" (e.g., in the same partition C k ), their loss should also be close, due to the intrinsic constraint imposed by f . For any algorithm that is robust, [57] proves Theorem 2.1 ( [57]). If a learning algorithm is (K, (·))-robust and L is bounded, a.k.a. L(f (x), y) ≤ M ∀z ∈ Z, for any η > 0, with probability at least 1 − η we have To control the first term, a natural approach is to constrain the variation of the loss function. Covering number [33] provides a way to bound the variation of the loss function, and more importantly, it conceptually realizes the actual number K of disjoint partitions.
In [26], they propose δ-isometry as a desirable property in CRL problem to help control the variation, where δisometry is a geometric property of mapping functions.
Definition 3 (δ-isometry). Given a map T that maps a metric space (P, ρ P ) to another metric space (Q, ρ Q ), it is called δisometry if the following inequality holds When T in eq. (1) and eq. (2) is of δ-isometry, by using ρ Q (T x, T x ) ≤ ρ P (x, x ) + δ a realization of algorithmic robustness (or GE bound in the form of Theorem 2.1) similar to [26] can be established for DNNs as follows.

Theorem 2.2.
Given an algorithm in a CRL problem, if the Lipschtiz constant of L w.r.t. T x is bounded by A, T is of δisometry, and X is compact with a covering number N γ/2 (X , ρ), then it is (|Y|N γ/2 (X , ρ), A(γ + δ))-robust.
Remark. The result in [26] is (|Y|N γ/2 (X , ρ), 2A(γ + δ))robust; the factor of 2 in the second term is dropped here due to the fact that in CRL problems, we are not doing metric learning as in [26], which involves two pairs of examples, and we only compare one pair of examples.
Remark. Denote ρ Q (T x, T x ) ≤ ρ P (x, x ) + δ as the expansion property of the δ-isometry, and ρ Q (T x, T x ) ≥ ρ P (x, x )−δ as its contraction property. We note that the above theorem is established by only exploiting the expansion property. After proving that DNNs achieve locally isometric mappings in section 3.1, we will show that a better generalization can be derived by considering the both properties.

Notations of deep neural networks
We study the map T as a neural network. We present the definition of Multi-Layer Perceptron (MLP) here, which captures all ingredients for theoretical analysis and enables us to convey the analysis without unnecessary complications, though we note that the analysis extends to Convolutional Neural Networks (CNNs) almost equally.
A MLP is a map that takes an input x ∈ R n from the space X , and builds its output by recursively applying a linear map W l followed by a pointwise non-linearity g where l ∈ {1, . . . , L} indexes the layer, x l ∈ R n l , x 0 = x, W l ∈ R n l ×n l−1 , and g denotes the activation function, which throughout the paper is the Rectifier Linear Unit (ReLU) [18]. Optionally g may include max pooling operator [7] after applying ReLU. We also denote the intermediate feature space g(W l x l−1 ) as X l and X 0 = X . Each X l is a metric space, and throughout the paper the metric is taken as the 2 norm || · || 2 , shortened as || · ||. We compactly write the map of a MLP as We denote the spectrum of singular values of a matrix W by σ(W ), and σ max and σ min are the maximum and minimum (nonzero) singular values of W respectively. We denote the rank of a matrix W by r(W ), and the null space of W by N (W ). We write the complement of N (W ) as X − N (W ).

NETWORKS
In this section, we develop GE bounds for CRL problems instantiated by DNNs. We identify two quantities that help control a bound, i.e., δ-isometry of T and the diameter γ of covering balls of X . We show that both of the two quantities can be controlled by constraining the spectrum of singular values of the weight matrix associated with each network layer, i.e., spectrums of singular values of To proceed, we consider in this paper variations of instances in X that are of practical interest -more specifically, those that output nonzero vectors after passing through a DNN. We first prove in lemma 3.1 that in such a variation subspace, mapping induced by a linear neural network is of δ-isometry, where δ is specified by the maximum and minimum singular values of weight matrices of all the network layers. For a nonlinear neural network, where we assume ReLU activation and optionally with max pooling, we consider the fact that it divides the input space X into a set of regions and within each region, it induces a linear mapping. In lemma 3.2, we specify the explicit form of region-wise mapping T q , associated with any linear region q, with submatrices of {W q i } i=1,...,L , based on which we prove in lemma 3.3 that a covering set for X can be found with a diameter γ of covering balls that is upper bounded by a quantity inversely proportional to the product of maximum singular values of weight matrices of some network layers. With the δ and γ specified in lemma 3.1 and lemma 3.3, we further prove in lemma 3.4 that in the instance-wise variation subspaces considered in this paper, a nonlinear neural network T is of local δ-isometry within each covering ball. The proofs are illustrated in fig. 1.
To develop a GE bound, we propose a covering scheme that includes instances of different labels into the same balls, thus reducing the size of covering set when compared with that in theorem 2.1. The covering scheme is illustrated in fig. 2. We correspondingly characterize both the errors caused by distance contraction between instances of different labels and those by distance expansion between instances of the same labels. Based on such characterization, we come with our main result of theorem 3.2.
Given the bound in theorem 3.2, we prove in lemma 3.5 that the optimal bound w.r.t. δ is obtained when all singular Black dots represent training instances in the sample space, and regions coded with different colors represent regions in the instance space that are hierarchically specified by layers of the DNN. The illustration depicts proofs of lemmas 3.1, 3.2, 3.3, and 3.4. Lemma 3.1 proves that the mapping induced by a linear DNN is of δ-isometry, where δ specified the expansion and contraction properties of the mapping and is determined by singular value spectrums of weight matrices of all the layers. Lemma 3.2 proves that a nonlinear DNN partitions the instance space into increasingly refined regions. As illustrated in the figure, the space is firstly partitioned into coarser regions (i.e., the center triangle and other three regions color coded as red, blue, and green), an additional layer further partitions some of the coarser regions into sets of smaller regions (e.g., those inside the region of center triangle), and the process goes recursively. Suppose that a region q is created by layer l of the DNN, and in region q, the nonlinear mapping defined by the matrix W l and activation function reduces to a linear mapping of W q l = diag(τ l (q))W l , where τ l (q) is a binary vector roughly indicating active neurons, and diag(·) diagonalizes τ l (q). Suppose q is created by layer l + 1 at the bottom part of the center triangle, and in region q , the nonlinear mappings of layer l and layer l + 1 are reduced to a linear mapping of W q l+1 = diag(τ l+1 (q ))W l+1 W q l , where symbols have similar meanings as described above. The phenomenon enables to find a covering for the sample space, such that in each covering ball that contains training instances, e.g., x, the DNN T defines a transformation that can be characterized as a linear mapping, e.g., T |x x. This is proved in lemma 3.3. Radius of the covering ball is illustrated as γ/2 in the figurea radius is acceptable as long as it is less than the smallest distance from any of the training instances to their respective region boundaries. Lastly, in lemma 3.4, we prove that by Cauchy interlacing law, a nonlinear DNN is of local δ-isometry within each covering ball specified above.
values of the weight matrix of each network layer are of equal ones, which inspires a straightforward choice of enforcing all singular values to have the value of 1, and thus the algorithms of OrthDNNs.

δ-isometry in deep neural networks
We begin with a few definitions necessary for the subsequent analysis.
Definition 4 (Variation subspace of an instance). Given an instance x ∈ X , suppose we are interested in the variation of a set X ⊆ X w.r.t. x. We call the linear vector space the variation subspace w.r.t. the instance x of X , shorten as variation subspace of instance x. Illustration of our used covering scheme that includes instances of different labels into same balls. Colored squares, circles, and stars represent instances of different labels, and unfilled (large) circles represent covering balls. Top: existing works [48], [57] usually separate instances of different labels into different covering balls, either by assuming that distances between instances of different labels in the sample space are infinite, or by using covering balls that are small enough not to contain instances of different labels; such a scheme can characterize the generalization errors caused by improper expansion of intraclass variations, but it cannot characterize the errors caused by improper contraction of inter-class differences. Bottom: we use a covering scheme that includes instances into covering balls regardless of their labels; it enables characterization of both types of the aforementioned errors, and leads to a possibly tighter generalization error bound in some cases.
The definition is to formalize variations of interest of particular instances, thus enabling us to discuss what variations a DNN is able to constrain. Correspondingly, we have the following definition of isometry.
Definition 5 (δ-isometry w.r.t. variation subspace of an instance). Given a map T that maps a metric space (P, ρ P ) to another metric space (Q, ρ Q ), it is called δ-isometry w.r.t. the variation subspace P x of instance x, if the following holds The following lemma specifies for a linear DNN the δisometry w.r.t. data variations of practical interest.

Lemma 3.1. Given a linear neural network T and an in-
The above lemma shows that as long as x varies within the complement of the null space of a linear DNN T , we can constrain variations induced by the mapping by the specified δ. Outside the space, it is unlikely of practical interest since T discards all the information about the variations.
To further proceed for nonlinear DNNs, we introduce some terminologies from hyperplane arrangement [42], and give definitions to describe the objects of interest exactly. 1. Note that the space X − N (T ) does not depend on x, and this is a trivial case where the variation subspaces of all instances are the same. We will see nontrivial cases later when dealing with nonlinear neural networks. Definition 6 ((Finite) Hyperplane arrangement). A finite hyperplane arrangement A is a finite set of affine hyperplanes in some vector space X ≡ K n , where K is a field and is taken as R in this paper.

Definition 7 (Region).
Denote H ∈ A an element of the arrangement, a region of the arrangement is a connected component of the complement R n − H∈A H. The set of all regions is denoted as R(A), shortened as R when no confusion exists.
We set up a labeling scheme for r ∈ R. Choosing a linear order in A, and α i is the normal of hyperplane H i . Let J = {1, 0}, and π i : J n → J be the projection onto the i-th coordinate. Define a map τ : X → J n by With the scheme, for any A, we would have an index set T (A, τ ), shortened as T , such that for any r ∈ R, it corresponds to a unique element in J n , denoted as τ (r). We will use τ (x) --labeling on elements, and τ (r) -labeling on regions, interchangely.

Definition 8 (Neuron).
A neuron a lk of a neural network T is a functional defined by where l ∈ {1, . . . , L} and k ∈ {1, . . . , n l }. All the neurons at layer l define a map, denoted as The lemma that follows is mostly an analysis of the domain of a DNN T . We begin with the following definition.
Definition 9 (Support of Neuron/DNN). Given a neuron a lk , the support of a lk is the set of instances in X that satisfy Similarly, the support of a neural network T is the set of instances in X that satisfy

Lemma 3.2.
A nonlinear neural network T divides X into a set of regions Q, and within each region q ∈ Q, T is linear w.r.t. variations of instances as long as they vary within q. We denote the linear mapping at q as T q and have where when layer l does not contain max pooling, τ l (q) is defined as and when layer l does contain max pooling, τ l (q) is defined as In the above definition, τ relu is defined as before, and where K is the set of indices of neurons being pooled over.
See the proof in appendix B.
Remark. The function τ l is intuitively a selection function that sets some rows of W l to zero, and selects a submatrix from it.
Remark. For q ∈ Q and q ∈ supp(T ), by the definition of linearity, T is still linear over q, i.e., the special case of T x = 0, ∀x ∈ q. In this case, x ∈ N (T q ).
In the proof of appendix B for lemma 3.2, we prove that within each region q ∈ Q, each neuron of a DNN is a linear functional. We summarize the result in the following corollary.
Corollary 3.1. A nonlinear neural network T divides X into a set of regions Q, and within each region q ∈ Q, the k th neuron a lk of the layer l is linear w.r.t. variations of x within q, and we have With the above lemma, we define the behavior of a DNN at a local area around x ∈ X or a local area around a set B ⊂ X as below.
Definition 10 (Linear neural network and neuron induced at x ∈ X from a nonlinear neural network). For any given x ∈ X with x ∈ q ∈ Q, we call the linear neural network T q the linear neural network induced by a nonlinear neural network T at x -denoting it as T |x , the linear neuron a q lk the linear neuron induced by the nonlinear neuron a lk at x -denoting it as a lk|x , and the submatrix of weight matrix of each layer l the submatrix induced by nonlinearity -denoting it as W l|x .
Definition 11 (Linear neural network and neuron induced at subset B ⊂ X from a nonlinear neural network). For any given B ⊂ X with B ⊂ q ⊂ Q, we call the linear neural network T q the linear neural network induced by a nonlinear neural network T at B -denoting it as T |B , the linear neuron a q lk the linear neuron induced by the nonlinear neuron a lk at B -denoting it as a lk|B , and the submatrix of weight matrix of each layer l the submatrix induced by nonlinearity -denoting it as W l|B .

Lemma 3.3.
For any nonlinear neural network T of L layers, a covering set for X can be found with a diameter γ = is a value depending on the training data and network weights, so is l(S m , T ) with 1 ≤ l(S m , T ) ≤ L (the dependence is specified in the proof), and σ i max is the maximum singular value of weight matrix W i of the i th layer.
See the proof in appendix C. We come with the local isometry property of DNNs after one more definition.
Definition 12 (γ-cover δ-isometry w.r.t. variation subspace of instance). Given an map T that maps a metric space (P, ρ P ) to another metric space (Q, ρ Q ), and an instance x ∈ P, it is called γ-cover δ-isometry w.r.t. variation space P x of x, if a γcover exists such that the following inequality holds where B denotes a ball given by the γ-cover.
m , where δ and γ are respectively specified in lemma 3.1 and lemma 3.3.
See the proof in appendix D.

Main results of generalization bound
With the local δ-isometry property of DNNs established, we derive in this section our main results of generalization bound.
Theorem 3.1. Given a CRL problem, the algorithm to learn is a nonlinear neural network of L layers, denoted as T . Suppose the following assumptions hold: Then, for any η > 0, with probability at least 1 − η we have where o(S m , T ) and 1 ≤ l(S m , T ) ≤ L are values depending on the training set S m and learned network T .
Proof. By lemma 3.4, we have that T is of γ-cover δisometry w.r.t. variation space of each training instance. The expansion property of δ-isometry gives Note that the proof differs from theorem 2.2 subtly, for that the loss difference needs to stay in the variation space of each covering ball. Since the tricky part is also present in the proof of theorem 3.2 (the condition x −x j ∈ P xj in eq. (7) ), to avoid tautology, we do not write the full proof here. The proof is finished by applying the robustness conclusion into theorem 2.1.
Regarding the 4 th assumption in theorem 3.1, it intuitively states that the local variation of interest w.r.t. each i.e., the variation vanishes after passing through the network. In practice, we are not interested in such a trivial case of vanishing local variations. Instead, it is the variation in the complement X − N (T |B ) that we want to constrain.
We have assumed that X is a regular k-dimensional manifold, whose covering number is ( C X γ/2 ) k , where C X is a constant that captures the "intrinsic" properties of X , and γ is the diameter of the covering ball. Such an assumption is general enough to accommodate at least visual data such as natural images and has been widely used [51].
In theorem 3.1, |Y|( C X γ/2 ) k corresponds the covering number of the joint space X × Y. We now show that with proper use of the contraction property of isometric mapping of DNNs, the constant |Y| in the second term of the upper bound can be removed while a modified first term has the potential to be small as well, indicating a better bound. Our idea is to directly exploit the covering number of the instance space X and measure the differences of the loss function L(f (T x), y) w.r.t. both arguments. Such a measure deals with instances of different labels but are close enough in X , thus characterizing both the errors that are caused by erroneously contracting the distance between instances from different classes and erroneously expanding the distance between instances from the same classes, instead of those of the same classes alone. However, it will cause the issue of infinite Lipschitz constant of loss function L(f (x), y). To see this, consider a binary classification problem whose loss function L is where (x, y) is an example, f (x) is a function that maps x to probability, and 1 is an indicator function. We provide an example case to illustrate the influence of metrics on Z.
Case. Let the metric on Z be ρ((x, y), (x , y )) = ||(x−x , y− y )|| = ||x − x || + |y − y |. Suppose that we have a pair of examples (x, y = 1) and (x, y = 0) that only differ in labels, we have When there exists an x such that f (x) → 1, we have |L(f (x), y) − L(f (x), y )| → +∞. This happens because as y changes values, due to its discreteness, it could induce a jump discontinuity on L(f (x), y), even though L(f (x), y) is Lipschitz continuous w.r.t. x. To avoid this, [57] and [26] employ a large covering number to ensure that examples in the same ball have the same label. We note that derivation of generalization bounds for robust algorithms concerns with the loss difference |L(f (x), y) − L(f (x ), y )| between example pairs. To address the aforementioned issue, we consider two separate cases for the loss difference: the cases of y = y and y = y . For the case y = y , we exploit the bounded Lipschitz constant of L w.r.t. T x. For the case y = y , we introduce the following pairwise error function to characterize the loss difference.
Definition 13 (Pairwise error function). Given a CRL problem, of which L is bounded for any compact set in Z, a.k.a. for z in any compact subset of Z, L(f (x), y) ≤ M , X is a regular k-dimensional manifold with a γ-cover, and T is of γcover δ-isometry, a pairwise error function (PE) of the tuple (L, f, T , Z, γ) is defined as It characterizes the largest loss difference for examples in Z that may arise due to the contraction and expansion properties of δisometry mapping. Note that PE(δ) is a monotonously increasing function of δ -a larger δ means more feasible examples in X and possibly larger distance contraction/expansion, leading to a possibly larger value of PE(δ). Theorem 3.2. Given a CRL problem, the algorithm to learn is a nonlinear neural network of L layers, denoted as T . Suppose the following assumptions hold: 1) ||x|| ≤ b ∀ x ∈ X , i.e., instances are norm bounded; 2) the loss function L is bounded, a.k.a. ∀z ∈ Z, L(f (T x), y) ≤ M , and the Lipschitz constant of L w.r.t T x is bounded by A; 3) X is a regular k-dimensional manifold with a covering number Then, for any η > 0, with probability at least 1 − η we have where , and γ is the same as that of theorem 3.1.
Proof. Similar to the proof of theorem 2.1, we partition the space Z via the assumed γ-cover. Since X is a k-dimensional manifold, its covering number is upper bounded by C k X /(γ/2) k . Let K be the overall number of covering set, which is upper bounded by C k X /(γ/2) k . Denote C i the i th covering ball and let N i be the set of indices of training examples that fall into C i . Note that (|N i |) i=1,...,K is an IDD multimonial random variable with parameters m and (|µ(C i )|) i=1,...,K . Then Remember that z = (x, y). We consider the two cases of y = y j and y = y j to bound eq. (7). When y = y j , by the assumption that T is of γ-cover δ-isometry w.r.t. P x of x ∈ S (x) m and the Lipschitz constant of L w.r.t. T x is A, suppose the maximum is achieved at x k and x k ∈ C p , we have where the second inequality holds since the γ-cover δ- When y = y j , given any training x k , we have where the inequality holds since by γ-cover δ-isometry of T , Thus eq. (7) is less than or equal to max{A(γ+δ), PE(δ)}. By Breteganolle-Huber-Carol inequality, eq. (10) is less than or equal to M log(2)2 k+1 C k X γ k m The proof is finished.
We now specify a case where the obtained bound in theorem 3.2 is tighter than that in theorem 3.1; for example, in the ball that covers (x k , y k ), few examples with y = y k are misclassified. In this case, , and the covering number is shrunken by a factor of |Y|. Our result only incrementally improves the generalization bound. However, it clearly shows that the contraction property of isometric mapping plays an important role in bounding the generalization error.

Suggestion of new algorithms
Many quantities exist in the GE bound established in theorem 3.2. Except γ and δ, all others are independent of the neural network. Although both γ and δ are controlled by singular values of weight matrices, we note that γ, which specifies the size of covering balls for a covering of X , is more of a trade-off parameter that balances between the first and second term of the GE bound, than of a variable used to control GE, as long as its values satisfy the condition of To control the bound via δ, we note that the minimum value of the bound w.r.t. δ is achieved when δ = 0, which implies L i=1 σ i max = 1 and L i=1 σ i min = 1. In the following lemma, we show that the condition is achieved only when σ i max = σ i min , ∀i = 1, . . . , L. Lemma 3.5. In theorem 3.2, δ = 0 is achieved only when Proof. It is straightforward to see that In the following, we show the two conditions hold only when σ i max = σ i min , ∀i = 1, . . . , L.
We show in lemma 3.5 that the optimal GE bound of theorem 3.2 w.r.t. δ is achieved only when all singular values of each of weight matrices of a DNN are equal. Among various solutions, the most straightforward one is that all singular values are equal to 1; in other words, each weight matrix has orthonormal rows or columns. This inspires a new set of algorithms that we generally term as Orthogonal Deep Neural Networks (OrthDNNs).

ALGORITHMS OF ORTHOGONAL DEEP NEURAL NETWORKS
In this section, we first present the algorithm of strict Or-thDNNs by enforcing strict orthogonality of weight matrices during network training. It amounts to optimizing weight matrices on their respective Stiefel manifolds, which however, is computationally prohibitive for large-sized networks. To achieve efficient OrthDNNs, we propose a novel algorithm called Singular Value Bounding (SVB), which achieves approximate OrthDNNs via a simple scheme of hard regularization. We discuss alternative schemes of soft regularization for approximate OrthDNNs, and compare with our proposed SVB. Batch Normalization [27] is commonly used to accelerate training of modern DNNs, yet it has a potential risk of ill-conditioned layer transform, causing its incompatibility with OrthDNNs. In fact, direct use of BN in OrthDNNs makes it ineffective to enforce strict orthogonality of weight matrices. We propose Degenerate Batch Normalization (DBN) to enable its use with strict OrthDNNs. We also propose Bounded Batch Normalization (BBN) to remove the potential risk of ill-conditioned layer transform. We finally explain how OrthDNNs are used for convolutional kernels.
Denote parameters of a DNN collectively as Θ = {W l , b l } L l=1 , where {b l } L l=1 are bias terms. We discuss algorithms of strict or approximate OrthDNNs in the following context. Given a training set {x i , y i } m i=1 , we write the training objective as L ({x i , y i } m i=1 ; Θ). Training is based on SGD (or its variants [49]), which updates Θ via a simple rule of Θ t+1 ← Θ t − η ∂L ∂Θ t , where η is the learning rate, and the gradient ∂L ∂Θ t is usually computed from a mini-batch of training examples. Network training proceeds by sampling for each iteration t a mini-batch from {x i , y i } m i=1 , until a specified number of iterations or the training loss plateaus.

The case of strict orthogonality
Enforcing orthogonality of weight matrices during network training amounts to solving the following constrained optimization problem where O stands for the set of matrices whose row or column vectors are orthonormal. For W l of any l th layer, problem (13) in fact constrains its solution set as a Riemannian manifold called Stiefel manifold, which is defined as M l = {W l ∈ R n l ×n l−1 |W l W l = I} assuming n l ≥ n l−1 , and is an embedded submanifold of the space R n l ×n l−1 , where I is an identity matrix. In literature, optimization of a differentiable cost function on such a matrix manifold and its convergence analysis have been intensively studied [1], [8]. For completeness, we briefly present the solving algorithm of (13) as follows. Denote T W l M l as the tangent space to M l at the current W l ∈ M l . First-order methods such as SGD first find a tangent vector Ω W l ∈ T W l M l that describes the steepest descent direction for the cost, and update W l as W l − ηΩ W l with the step size η that satisfies conditions of convergence, and then perform a retraction R W l (−ηΩ W l ) that defines a mapping from the tangent space to the Stiefel manifold, which can be achieved by R W l (−ηΩ W l ) = Q(W l − ηΩ W l ), where the operator Q denotes the Q factor of QR matrix decomposition. To obtain the tangent vector Ω W l , one may project the gradient ∂L ∂W l in the embedding space R n l ×n l−1 (or its momentum version [49]) onto the tangent space T W l M l by P W l ∂L ∂W l , where P W l defines a projection operator according to the local geometry of W l ∈ M l . Convergence analysis for such a scheme to obtain the tangent vector is presented in [43]. In Appendix E, we present the algorithmic details for optimization of weight matrices on the Stiefel manifolds.

Achieving near orthogonality via Singular Value Bounding
Constraining solutions of weight matrices of a DNN on their Stiefel manifolds is an interesting direction of research. It also supports analysis of theoretical properties as in section 3 and the related works [43], [46]. However, it arguably has the following shortcomings concerned with computation, empirical performance, and also compatibility with existing deep learning methods, which motivate us to address these shortcomings by developing new algorithms of approximate OrthDNNs.
• Strict constraining of weight matrices on the Stiefel manifolds requires expensive computations -in particular, the operations of projecting the Euclidean gradient onto the tangent space and retraction onto the Stiefel manifold (Steps 2 and 4 in Appendix E) dominate the costs in each iteration. If we allow the solutions slightly away from the manifolds, the expensive projection and retraction operations are not necessary to be performed in each iteration. Instead, similar pulling-back operations can be performed less frequently, e.g., in every a certain number of iterations, and consequently such a burden of pulling back is amortized.
• Theorem 3.2 gives a bound GE(f Sm ) of the expected error R(f Sm ) w.r.t the training error R m (f Sm ). To achieve good performance on practical problems, both R m (f Sm ) and GE(f Sm ) should be small. However, optimization of DNNs is characterized by proliferation of local optima/critical points [12], [29]. When we are motivated to optimize weight matrices on their Stiefel manifolds, obtaining W l ∈ M l , l ∈ {1, . . . , L}, with Ω W l = 0, it is very likely that for a W l , there exists a better local optimum in the embedding Euclidean space that is slightly away from the manifold (e.g., Ω W l = 0 while ∂L ∂W l = 0, or the Euclidean gradient is in the complement null space of the current tangent space), and has a smaller R m (f Sm ). If we allow the optimization to step away from, but still pivot around, the manifold, better solutions could be obtained by escaping from local optima on the manifold.
• Successful training of modern DNNs depends heavily on BN [27], a technique that can greatly improve training convergence and empirical results. However, as analyzed shortly in section 4.4, BN would change the spectrum of singular values of each layer transform (i.e., the combined linear transform of each layer achieved by weight mapping and BN, as specified in (16)). Consequently, the efforts spending on enforcing strict orthogonality of weight matrices become ineffectual. Algorithms of approximate Or-thDNNs seem more compatible with BN transform.
To develop an algorithm of approximate OrthDNNs, we propose a simple yet effective network training method called Singular Value Bounding (SVB). SVB is a sort of projected SGD method and can be summarized as follows: SVB simply bounds, after every T svb iterations of SGD training, all the singular values of each W l , for l = 1, . . . , L, in a narrow band [1/(1 + ), (1 + )] around the value of 1, where ≥ 0 is a specified small constant. Appendix F presents the algorithmic details.
After each bounding step, optimization of SVB in fact proceeds in the embedding Euclidean space, to search for potentially better solutions, before next bounding step that pulls the solutions back onto ( = 0) or near ( > 0) the Stiefel manifolds. With annealed learning rate schedules, we observe empirical convergence of SVB. Compared with manifold optimization in section 4.1, SVB is more efficient since the dominating computation of SVD is invoked only every a certain number of iterations. Experiments of image classification in section 5 show that SVB sometimes outperforms the algorithm of strict OrthDNNs in section 4.1, both of which outperform the commonly used SGD based methods, and in many cases with a large margin.

Alternative algorithms for approximate OrthDNNs
To achieve approximate OrthDNNs, one may alternatively penalize the main objective L ({x i , y i } m i=1 ; Θ) with an augmented term that encourages orthonormality of columns or rows of weight matrices, resulting in the following unconstrained optimization problem where · F denotes Frobenius norm, λ is the penalty parameter, and we have assumed n l ≥ n l−1 for a certain layer l. By using increasingly larger values of λ, the problem (14) approaches to achieve strict OrthDNNs. One can use SGD based methods to solve (14), where the additional computation cost incurred by the regularizer is marginal. However, it is not straightforward to set a proper value of λ in order to strike a good balance between the two terms of (14). Note that soft regularization algorithm of the type (14) is used in the related works [11], [55]. In this work, we compare them with our proposed SVB. Experiments in section 5 show that SVB performs better than this alternative on modern architectures.

Compatibility with Batch Normalization
We start this section by showing that the original design of Batch Normalization [27] is incompatible with our proposed OrthDNNs. Technically, for a network layer that computes, before the nonlinear activation, h = W x ∈ R n , BN inserts a normalization denoted as BN(h) = BN(W x), where we have ignored the bias term for simplicity. BN in fact applies the following linear transformation to h where each entry of µ ∈ R n is the output mean at each of the n neurons of the layer, the diagonal matrix Φ ∈ R n×n contains entries {1/φ i } n i=1 that is the inverse of the neuronwise output standard deviation φ i (obtained by adding a small constant to the variance for numerical stability), Υ ∈ R n×n is a diagonal matrix containing trainable scalar parameters {υ i } n i=1 , and β ∈ R n is a trainable bias term. Note that during training, µ and φ for each neuron are computed using mini-batch examples, and during inference they are fixed representing the statistics of all the training population, which are usually obtained by running average. Thus the computation (15) for each example is deterministic after network training.
Inserting h = W x into (15) we get which is simply a standard layer with change of variables.
The following lemma suggests that BN is incompatible with OrthDNNs: even though W is enforced to have orthonormal rows or columns, BN would change the conditioning of layer transform, i.e., spectrum of singular values of W , by learning Υ and Φ whose product does not necessarily contain diagonal entries {υ i /φ i } n i=1 of equal value. To make BN compatible with strict OrthDNNs, we propose Degenerate Batch Normalization (DBN) that learns layerwiseῡ andφ instead of neuron-wise {υ i } n i=1 and {φ i } n i=1 , so that the learned Υ and Φ respectively contain diagonal entries of equal value. Such a Υ is learned simply by using a single trainable parameterῡ shared by all n neurons of the layer. To learn such a Φ, DBN still computes neuronwise {φ i } n i=1 in each iteration of training, but uses running average to learnφ = 1 n n i=1 φ i that would be shared by the n neurons. The proposed DBN enjoys the benefit of neuronwise normalization in BN, and when training converges, it learns a layer transform W = ΥΦW =ῡ/φW whose conditioning is the same as that of W , thus achieving compatibility with strict OrthDNNs. Section 4.2 discusses the potential advantages of approximate OrthDNNs over the strict ones. To make BN compatible with approximate OrthDNNs, especially with our proposed SVB method, we propose Bounded Batch Normalization (BBN) that controls the variations among {υ i /φ i } n i=1 , so that the conditioning of layer transform is not severely affected by ΥΦ. More specifically, BBN computes α = 1 n n i=1 υ i /φ i in each iteration of training, and bounds each of in a narrow band [1/(1 +˜ ), (1 +˜ )] around the value of 1, where˜ ≥ 0 is a scalar parameter. Appendix H presents algorithmic details of the proposed BBN.
The introduction of DBN makes it possible to empirically compare strict and approximate OrthDNNs in the context of modern architectures, which is presented in section 5. Experiments in section 5 also show that performance is improved when using BBN instead of BN, confirming the benefit by resolving BN's compatibility with OrthDNNs.

Orthogonal Convolutional Neural Networks
In previous sections, we present theories and algorithms of OrthDNNs by writing their layer-wise weights in matrix forms. When applying DNNs to image data, one is actually using networks with convolutional layers. For an l th convolutional layer with weight tensor of the size n l × n l−1 × n h × n w , where n h and n w denote the height and width of the convolutional kernel, we choose to convert the tensor as a matrix of the size n l × n l−1 n h n w based on the following rational. Natural images are usually modeled by first learning filters from (densely overlapped) local patches, and then applying the thus learned filters to images to aggregate the corresponding local statistics. The convolutional layer in fact linearly transforms the n l−1 input feature maps in the same way, by applying each of n l filters of the size n l−1 × n h × n w to n l−1 patches of the size n h × n w in a sliding window fashion, resulting in local responses that are arrayed in the form of n l feature maps, which have the same size as that of input feature maps when padding the boundaries. In other words, the convolutional layer applies linear transform, using n l filters, to n l−1 n h n w -dimensional instances that are collected from local patches of input feature maps. Correspondingly, we choose to convert the weight tensor containing the n l filters of dimension n l−1 n h n w to its matrix form, and apply to it specific algorithms of OrthDNNs.

EXPERIMENTS
We present in this section extensive experiments of image classification to verify the efficacy of OrthDNNs. We are particularly interested in how algorithms of strict or approximate OrthDNNs provide regularization to various architectures of modern DNNs, such as ConvNets [36], [47], ResNets [20], [21], DenseNet [25], and ResNeXt [56]. We use the benchmark datasets of CIFAR10, CIFAR100 [34], and ImageNet [45] for these experiments. We compare empirical performance and efficiency among different algorithms of strict and approximate OrthDNNs. For some of these comparisons, we also investigate behaviours of OrthDNNs under regimes of both small and large sizes of training samples, and robustness of OrthDNNs against corruptions that are commonly encountered in natural images, in order to better understand the empirical strength of OrthDNNs. For network training, we use SGD with momentum and initialize networks using orthogonal weight matrices, where the momentum is set as 0.9 with a weight decay of 0.0001. When our proposed SVB is turned on, we apply it to weight matrices of all layers after every epoch of training.

Comparative studies on algorithms of strict and approximate OrthDNNs
In this section, we use architectures of ConvNet and ResNet on CIFAR10 to study the behaviours of strict and approximate OrthDNNs. The CIFAR10 dataset consists of 60, 000 32 × 32 color images of 10 object categories (50, 000 training and 10, 000 testing ones). We use raw images without preprocessing. Data augmentation follows the standard manner in [38]: during training, we zero-pad 4 pixels along each image side, and sample a 32 × 32 region crop from the padded image or its horizontal flip; during testing, we use the original non-padded image. Our ConvNet architectures follow [20], [47]. Each network starts with a conv layer of 16 3 × 3 filters, and then sequentially stacks three types of 2X conv layers of 3 × 3 filters, each of which has the feature map sizes of 32, 16, and 8, and filter numbers of 16, 32, and 64, respectively; spatial sub-sampling of feature maps is achieved by conv layers of stride 2; the network ends with a global average pooling and a fully-connected layer. The ResNet construction is based on the ConvNets presented above, where we use an "identify shortcut" to connect every two conv layers of 3 × 3 filters and use a "projection shortcut" when sub-sampling of feature maps is needed; we adopt the pre-activation version [21]. Thus, for both types of networks, we have 6X + 2 weight layers in Comparison of strict and approximate OrthDNNs on the CIFAR10 dataset [34], using ConvNet and ResNet architectures respectively of 20 weight layers (referring to the main text for their specifics). Degenerate batch normalization is used due to its compatibility with strict OrthDNNs. total. We set X = 3 for experiments in this section, giving networks of 20 weight layers. The DBN proposed in section 4.4 is designed to be compatible with strict OrthDNNs. We use DBN to enable training and comparison of strict and approximate OrthDNNs on the networks constructed above. We implement strict OrthDNNs as the algorithm presented in section 4.1. We use our proposed SVB and those in [11], [55] (i.e., the problem (14)) for approximate OrthDNNs. The learning rates start at 0.1 and end at 0.001, and decay every two epochs until the end of 160 epochs of training, where we set the mini-batch size as 128. We fix the parameter of SVB as 0.05, while λ of soft regularization in (14) is optimally tuned as 0.1. Table 1 gives the results with the curves of training convergence plotted in fig. 3. Table 1 shows that on both of the two networks, algorithms of strict and approximate OrthDNNs outperform standard SGD based method, confirming the improved generalization by their regularization of network training. Moreover, approximate OrthDNNs via either SVB or soft regularization perform as well as strict ones, but at a much lower computational cost 2 , suggesting their advantage in practical use. Due to the prohibitive computation of strict OrthDNNs on modern architectures of larger sizes, we choose to use approximate OrthDNNs in subsequent experiments, and correspondingly use BN or our proposed BBN to replace DBN.

Comparison of hard and soft regularization for approximate OrthDNNs
In this section, we study algorithms of approximate Or-thDNNs by comparing our proposed SVB with soft regu-2. In table 1, the respective dominating computations of QR decomposition for manifold optimization and singular value decomposition for SVB are based on CUDA implementation.  larization [11], [55]. The experiments are conducted on CI-FAR10 using a pre-activation version of ResNet constructed in the same way as in section 5.1. Setting X = 11 gives a total of 68 weight layers. To train the network, we use learning rates that start at 0.5 and end at 0.001, and decay every two epochs until the end of 160 epochs of training, where we set the mini-batch size as 128. Approximate OrthDNNs via SVB and soft regularization are compared with the baseline of standard SGD with momentum. We also switch BBN on or off to verify its effectiveness. We fix of SVB as 0.5, while the penalty λ of soft regularization is optimally tuned as 0.005. We fix˜ of BBN as 0.2. We run each setting of experiments for five times, and report results in the format of best (mean ± standard deviation). Table 2 shows that approximate OrthDNNs via both SVB and soft regularization provide effective regularization to network training, and SVB outperforms soft regularization with a noticeable margin. Compared with BN, our proposed BBN can better regularize training and give slightly improved performance. Note that algorithmic design of BBN may not be compatible with soft regularization, which explains the degraded performance when using them together.

Experiments with Modern Architectures
In this section, we investigate how our proposed SVB and BBN methods provide regularization to modern architectures of ResNet [21], Wide ResNet [59], DenseNet [25], and ResNeXt [56]. We use CIFAR10, CIFAR100 [34], and ImageNet [45] for these experiments. The CIFAR100 dataset has the same number of 32 × 32 color images as CIFAR10 does, but it has 100 object categories where each category contains one-tenth images of those of CIFAR10. We use data augmentation in the same way as for CIFAR10. The ImageNet dataset contains 1.28 million images of 1, 000 categories for training, and 50, 000 images for validation. We use data augmentation as in [56].
For experiments on CIFAR10 and CIFAR100, we use the following specific architectures. ResNet is constructed in the same way as in section 5.1; we set X = 9 here giving a total of 56 weight layers. Wide ResNet is the same as "WRN-28-10" in [59]. ResNeXt is the same as "ResNeXt-29 (16 × 64d)" in [56], i.e., the depth L = 29, cardinality C = 16, and the feature width in each cardinal branch d = 64. We use consistent hyper-parameters to train these architectures. The learning rates start at 0.5 and end at 0.001, and decay every two epochs until the end of 300 epochs of training, where we set the mini-batch size as 64 -note that this schedule with more training epochs and smaller mini-batch size usually gives better empirical performance than the training schedule used in section 5.2 does. We fix and˜ of SVB and BBN as 0.5 and 0.2 respectively. Table  3 confirms that SVB and BBN improve generalization of various architectures. We also observe that improvements on CIFAR100 are generally greater than those on CIFAR10, which may be due to the problem nature of smaller sample size for CIFAR100. We will investigate how our methods perform with varying sample sizes more thoroughly in the subsequent section.
For experiments on ImageNet, we use top-performing models of the following architectures: "ResNet-152" of [21], "DenseNet-264" of [25], and "ResNeXt-101 (64×4d)" of [56]. To train these models, we use the same hyper-parameters as respectively reported in these methods. The parameters and˜ of SVB and BBN are fixed as 0.5 and 0.5 respectively. Results in table 4 confirm that approximate OrthDNNs via our proposed methods improve generalization by providing effective regularization to large-scale learning.  Error rates (%) on the validation set of ImageNet [45] when applying our proposed SVB and BBN to various modern architectures (referring to the main text for their specifics). Results are based on single-crop testing of the size 320 × 320.

Effects of Varying Sample Sizes
We are also interested in the efficacy of SVB and BBN for problems with varying sizes of training samples. To this end, we respectively sample 1/10, 1/5, 1/2, or all of training images per category from ImageNet [45], which constitute our ImageNet training subsets of varying sizes. We train the "ResNeXt-101 (64×4d)" model of [56] in the same way as in section 5.3 for this investigation. Fig. 4 shows that SVB and BBN consistently improve classification across the regimes from small to large sizes of training samples, and the improvements are more significant for the smaller ones, suggesting the great potential of approximate OrthDNNs for machine learning problems of small sample sizes.

Robustness against Common Corruptions
We have shown in previous experiments that OrthDNNs, particularly our proposed SVB and BBN, have better generalization to testing samples that are drawn from the same distributions of training ones. In this section, we investigate the robustness of OrthDNNs when testing samples are corrupted such that they are getting away from the distributions of training ones. We focus on corruptions that are frequently encountered in natural images, e.g., the "common" corruptions of noise, blur, weather, or digitization [22]. Existing research suggests fragility of deep learning models to corruptions of such kinds [14], and that finetuning on specific corruption types would help, but cannot well generalize to other types of corruptions [16], [50].
To this end, we use the ImageNet-C dataset [22] that is produced by applying 15 corruption types of 5 severity   fig. 5 demonstrates better robustness of our proposed methods against common corruptions, and the robustness stands gracefully with the increase of severity levels.

CONCLUSION
In this paper, we present theoretical analysis to connect with the recent interest of spectrally regularized deep learning methods. Technically, we prove a new generalization error bound for DNNs, which is both scale-and range-sensitive to singular value spectrum of each of networks' weight matrices. The bound is established by first proving that DNNs are of local isometry on data distributions of practical interest, and then introducing the local isometry property of DNNs into a PAC based generalization analysis. We further prove that the optimal bound w.r.t. the degree of isometry is attained when each weight matrix has a spectrum of equal singular values -OrthDNNs with weight matrices of orthonormal rows or columns are thus the most straightforward choice. Based on such analysis, we present algorithms of strict and approximate OrthDNNs, and propose a simple yet effective algorithm called Singular Value Bounding. We also propose Bounded Batch Normalization to make compatible use of batch normalization with OrthDNNs. Experiments on benchmark image classification show the efficacy and robustness of OrthDNNs and our proposed SVB and BNN methods, especially for problems where the sizes of training samples are relatively small. We expect OrthDNNs provide effective regularization to other deep learning problems of smaller sample sizes.

APPENDIX A PROOF OF LEMMA 3.1.
To prove lemma 3.1, we begin with the following lemma regarding matrix pseudo-inverse.
where W † is the pseudo-inverse of W , given as W † = V Σ † U T when W has the singular value decomposition W = U ΣV T , and Σ † is the matrix obtained by first taking the transpose of Σ, and then the inverse of its non-zero elements.
Proof. Let R be the index set such that Σ rr = 0, ∀r ∈ R.
Given any x ∈ X − N (W ), x can be represented as where entries of α have α r = 0 when r ∈ R. Then Now, we are going to prove lemma 3.1.
Proof. For any W i with i ∈ {1, . . . , L}, performing singular value decomposition (SVD) upon it, we have We show that if ∆ i−1,k = 0, then Σ i,kk = 0, which implies that ∆ i−1 lies in the subspace spanned by right singular vectors of W i that have nonzero singular values, where ∆ i−1,k and Σ i,kk are respectively the k th element of ∆ i−1 and k th diagonal element of Σ i . Suppose otherwise, for a set K, ∆ i−1,k = 0 and Σ i,kk = 0, ∀k ∈ K. Let v k denote the k th column of V i , we reparameterize ∆ i−1 as where the terms having ∆ i−1,k = 0 are omitted.
Denote W = i−1 j=1 W j . Since ∆ ∈ X − N (T ), we have ∆ ∈ X − N (W ). By lemma A.1, we have W † ∆ i−1 = W † W ∆ = ∆ ∈ X − N (T ). Since By assumption we have ∆ i−1,k = 0; we also have W † v k = 0 -otherwise ∆ i−1,k would be zero, and W † v k ⊥ W † v k , for k = k ; Considering that ∆ ∈ X − N (T ), we have W † v k ∈ X − N (T ). However, we assume Σ i,kk = 0 for k ∈ K, and have which implies W † v k ∈ N (T ) and leads to a contradiction. We thus prove that if ∆ i−1,k = 0, then Σ i,kk = 0.
With the above result, we can constrain the sample variation more precisely through singular values of weight matrices. To be specific, for any pair of x i−1 ∈ X i−1 and where the second and last equalities use the fact that an orthogonal matrix does not change the norm of operated vectors, and the two inequalities are derived based on our result that for x i−1 − x i−1 ∈ ∆ i−1 , it lies in the subspace spanned by the right singular vectors of W i whose corresponding singular values are great than or equal to the nonzero σ i min . Similarly, we have Cascading on all layers, we have Thus We conclude the proof by showing T is of

APPENDIX B PROOF OF LEMMA 3.2.
Proof. We proceed by induction on layer l. For l = 1, each row in W l corresponds to a hyperplane in X . Thus, W l imposes a hyperplane arrangement A = {W l,i } i=1,...,n l on X , and is associated with an index set T (A, τ ) of the region set R(A), where W l,i denotes the i th row of W l . Denote a li the neuron corresponding to a hyperplane W l,i ∈ A, we have i.e., a li is linear over regions of R that are active on the i th neuron of the layer. We then have Q li = {q ∈ R|π i τ (q) = 1} as the support of a li .
To present the effect in the form of W l , we have W q l = diag(τ (q))W l , which is a linear map over each q ∈ Q l = n l i=1 Q li . The case l = 1 is proved, with τ 1 = τ .
Assume now for all the neurons a lj , j = 1, . . . , n l , of layer l, a lj is a linear functional over its support Q lj , and is 0-valued otherwise. We proceed by building a new set of regions Q (l+1)i , i = 1, . . . , n l+1 , for layer l + 1, whose neurons are linear functionals over regions of Q (l+1)i .
We separately discuss the cases of g with or without max pooling.
First, when g does not include max pooling, for a neuron a (l+1)i of layer l + 1, it is a functional of the form a (l+1)i = g n l j=1 W l+1,ij a lj = g pre(a (l+1)i ), where W l+1,ij is the (i, j)-entry of W l+1 . Since ∀q ∈ n l j=1 Q lj , a lj , j = 1, . . . , n l , is a linear functional over q, so is the linear combination pre(a (l+1)i ) of them.
Consequently, neurons a l+1 further divide the region q into sub-regions, where for each region, a l+1 is a linear map. We say the boundaries of the new set of regions as the hyperplane arrangement induced by neurons, with a alight abuse of terminology. For the newly created set of regions Q l+1 , we define the labeling function τ l+1 for layer l + 1 such that for x ∈ q ∈ Q l+1 , we have Since for each q ∈ Q l+1 , it is a sub-region of q ∈ Q l , q ∈ Q l holds true as well. In this case, the new set of regions are built.
For the case that g includes max pooling, the neuron is of the form Similarly, for each k ∈ K, pre(a (l+1)i ) k may split q into two sub-regions. Denote the boundary as H k if it indeed splits. {H k } k∈K , together with the boundary of q, forms a hyperplane arrangement A * (l+1)i induced by neurons within q. Denote the set of regions in this new arrangement as Q * (l+1)i , for each q * ∈ Q * (l+1)i , consider the set of hyperplanes A (l+1)i = { pre(a (l+1)i ) k − pre(a (l+1)i ) k } k<k ,k,k ∈K . With a similar argument, they will create another hyperplane arrangement A (l+1)i within q * . For q ∈ Q (l+1)i , H ∈ A (l+1)i does not have discontinuity in derivative -does not suddenly switch from constant function 0 to non-zero linear function. Now within each q , we impose an order on A * (l+1)i , if pre(a (l+1)i ) k − pre(a (l+1)i ) k ≥ 0, we say H k ≥ H k . Given that K is a finite totally ordered set, the maximum element w.r.t. the defined order exists, and we denote its index as k max . Thus, for each q , a (l+1)i = pre(a (l+1)i ) kmax , which we have proved to be a linear function over q in the ReLU case, so is a li . Thus similar to the ReLU only case, Max Pooling with ReLU divides q into a new set of regions as well.
Since for each q ∈ Q l+1 , it is a sub-region of q ∈ Q l , q ∈ Q l holds true as well. In this case, the new set of regions are built. The max pooling case is proved.
The same with the case l = 1, to present the effect in the form of W l , denoting τ l+1 (x) = τ relu l+1 for ReLU only case, and τ l+1 (x) = τ max l+1 τ relu l+1 for ReLU with Max Pooling case, we have W q l+1 = diag(τ l+1 (q))W l+1 , which is a linear map over q ∈ Q l+1 = n l+1 i=1 Q (l+1)i . By induction, ∀i ≤ l + 1, W q i is a linear map. Cascading the result, we have the neural network T as a linear map over q ∈ Q (l+1) We now finish the induction and prove that for 0 < l ≤ L, there exists a set Q l such that ∀q ∈ Q l , T l q is linear over q. Given that we are interested in T L q in this paper, we drop the upper index, and denote it as and the corresponding region set Q L is the set of regions over which T is linear. We also drop the index, and denote it as Q.

APPENDIX C PROOF OF LEMMA 3.3.
Proof. By lemma 3.2, a set of regions Q exists such that for q ∈ Q, T q is a linear mapping induced by T over q.
For any given x ∈ X , denote the region it belongs to as q x . Let d min = min x∈S (x) m min x ∈∂qx ρ(x, x ), the shortest distance from x to the boundary of q x , denoted as ∂q x , among all training samples. Denote by a lk the neuron that defines the hyperplane corresponding to the boundary that produces the shortest distance d min . Note that a lk may exist in the intermediate network layers, i.e., 1 ≤ l ≤ L, and the specific value of l depends on the training set S m and the learned T . By Corollary 3.1., we have ∀{x ∈ X | x + (x − x) ∈ q x }, a lk|x (x) − a lk|x (x ) = a lk|x (x − x ). Thus Given x ∈ ∂q x , it implies a lk (x ) = 0, and we have a lower bound on d min as where σ i max |x is the maximum singular value of W q i for i = 1, . . . , l.
Since W q i is a submatrix of W i , by Cauchy interlacing law by rows deletion [10], we have σ i max |x ≤ σ i max , and Denote o(S m , T )/2 = |a lk|x (x)| to stress the fact that it is a fixed value once S m and T are given, we have a lower bound where we have explicitly write l(S m , T ) to emphasize the dependence of l on S m and T . In addition, q x ∈ Q, we have T is linear over q x . Consequently, a covering set of X with radius r is found, such that within each covering ball, T is linear. Then for any given x ∈ X and {x ∈ X | ||x − x || ≤ r}, x ∈ q x , thus T x − T x = T |x (x − x ). The diameter γ of the covering ball is 2r.