Title: THE EFFECT OF RESIDUAL ARCHITECTURE ON THE PER-LAYER GRADIENT OF DEEP NETWORKS Anonymous authors Paper under double-blind review

Abstract: A critical part of the training process of neural networks takes place in the very first gradient steps post-initialization. In this work, we study the connection between the network's architecture and initialization parameters, to the statistical properties of the gradient in random fully connected ReLU networks, through the study of the the per layer Jacobian in finite sized networks. We compare three types of architectures: vanilla networks, ResNets and DenseNets. In our analysis, we show that while the variance of Jacobian squared norm is exponential in depth for ResNets, and polynomial for DenseNets, there exists an initialization strategy for both, such that the norm is preserved through arbitrary depths, preventing exploding or decaying gradients in deep networks.We also show that the statistics of the per layer Jacobian norm is a function of the architecture and the layer's size, but surprisingly, not the layer's depth.

1. INTRODUCTION:

Understanding the effect of different architectures on the ability to train deep networks has long been a major topic of research. A recent line of work has focused on the statistical properties of both the activations and the gradients of different architectures at the point of initialization, in order to better understand the effect of different architectural choices on the training dynamics.

A well-known hurdle in the ability to train deep networks, is the effect of depth on the magnitude of the gradients. Specifically, the exploding and vanishing gradient phenomena have been observed in various deep architectures. they are considered a main culprit in the failure of training deep networks successfully. A flurry of architectural tricks and initialization techniques have been proposed, in order to mitigate the problem of exploding and vanishing gradients, ranging from residual architectures, different normalization techniques and orthogonal initializations.

In this work, we perform a rigorous analysis of the norm, of both the activations and the gradients of each layer, in fully connected architectures at the point of initialization. Our setting is as follows: we assume a fixed input y 0 to an L layer neural network with weights given by W , and intermediate output for any layer l denoted by y l . Denoting by w k ij the weights ij of layer k, we are interested in the quantities:

y l 2 , J k 2 = ij ∂y L ∂w k ij 2 (1)

representing the norm of the activations of layer l and Jacobian with respect to layer k. We consider these quantities random variables at initialization, and, therefore, devote our analysis to deriving expressions for the mean and variance of both y l 2 and J k 2 , over the random initializations of the weights. Note that the full gradient requires an additional loss term which is independent of the architecture, and so our attention is focused on J k . While most initialization techniques are designed to maintain the squared length of intermediate activations through the layers, the multiplication of i.i.d random weights along input-to-output paths in increasingly deep networks might lead to unbounded higher order statistics. In the extreme case, while the mean of activation and gradients remain constant in deep layers, it would be improbable to achieve mean values for any specific realization of the weights, rendering training difficult. Specifically, given some loss function L, the norm of the gradient with respect to W k is given by:

∇ W k L = ∂L ∂y L J k ≤ ∂L ∂y L J k (2)

Large fluctuations to the Jacobian norm with respect to small perturbations of the weights will, therefore, cause wildly fluctuating gradients during training, hindering training dynamics.

Our analysis makes use of a duality theorem, which captures an intuitive concept. In residual architectures, there are many alternative paths of varying lengths from input to output. The theorem states that when considering the moments of the gradient at a certain layer, one can change the architecture, so all connections which that skip that layer are removed. This is true, despite the activations of the two networks being very different. See Fig. 1 .

Our contributions are as follows: (i) We rigorously study the effect of different architectures, namely vanilla, residual and densely connected architectures, on the mean and variance of y l 2 and J k 2 in finite depth and width networks. (ii) From our forward-backward norm propagation duality theorem, we derive conditions that prevent exploding gradients in each of these architectures. (iii) We prove that for densely connected networks, the commonly used initialization of He et al. (2015) is enough to ensure the boundness of var( J k 2 ) and var( y l 2 ) for arbitrarily deep networks, shedding new light on the apparent advantage of these architectures in practical applications. (iv)

We prove that replacing the standard conv-ReLU blocks with blocks of the concatenated ReLU (CR) activation presented by Shang et al. (2016) , while maintaining the same parameter budget, significantly reduces the fluctuations of both the activations and gradients at initialization for any ReLU based architecture. We evaluate the applicability of this simple alteration on a large number of UCI datasets, and demonstrate the advantage for CR based networks over their standard counterparts.

2. RELATED WORK:

The effect of architecture and depth on the gradients has been a topic of considerable interest recently, leading to several insightful research directions. A major line of work has been focused on the mean field approach, where the propagation of signals in random networks are analyzed in the regime of infinite width (Yang & Schoenholz (2017 ),Schoenholz et al. (2017) , Xiao et al. (2018) ). In that case, the pre-activations in any layer are approximated by a Gaussian distribution with moments depending on the specific layer. For ReLU networks (without residual connections), controlling the exponential explosion or decay of both the activations and gradients with depth is a simple matter of correctly scaling the initialization parameters. Indeed, it is shown in Schoenholz et al. (2017) that classical feed forward networks exhibit exponential behavior on average, when not initialized properly. For the gradients, a recursion relation is derived:

E ( ∂y L ∂y l i ) 2 = E ( ∂y L ∂y l+1 i ) 2 σ 2 w 1 √ 2π e -z 2 2 φ ( √ q * z)dz (3)

where φ the derivative of the activation function and σ 2 w denotes the variance of the weights. By initializing with σ 2 w = 2 f an-in , and considering that for the ReLU activation function φ (z) = 1 z>0 , we have that E((

EQUATION

Mean field analysis of residual networks is presented in Yang & Schoenholz (2017) , where it is shown that ResNets exhibit sub-exponential, or even polynomial dynamics, depending on the activation function. While mean field theory captures the "average" dynamics of random neural networks, it does not take into account the scale of gradient fluctuations that are crucial to the stability of gradient descent, which is the focus of this work. A more refined analysis presented in Pennington et al. (2017) , considers the full spectrum of the input output Jacobian of infinitely wide networks, given by:

J IO = ∂y L ∂y 0 (5)

Exploding or vanishing gradients are then prevented by analyzing conditions that give rise to dynamical isometry, a state in which the squared eigenvalues of J IO are all concentrated around 1. In (Pennington et al., 2017) it is shown that fully connected ReLU networks are incapable of reaching dynamical isometry, as apposed to sigmoidal networks, while in (Pennington et al., 2018) it is shown that dynamical isometry is achieved in a universal manner for a variety of activation functions. Dynamical isometry is a much stronger condition than the Jacobian norm, since it requires access to the entire spectrum. However, in previous work, this is done using various forms of infinite width approximations. Instead, this work focuses on finite depth corrections to the statistics of various quantities concerning the squared Frobenius norm of the per-layer Jacobian (derivative of the output with respect to the weights), and can be viewed as a complementary approach to analyzing the dynamics at the start of training. Our results cannot be obtained using large width approximations, since we seek to characterise the effect of both depth and width on the Jacobian.

A recent line of work, more related to the analysis presented in this paper, deals with the statistical properties of the norm of activations and gradients in random finite width networks. Hanin & Rolnick (2018) tackle two failure modes that are caused by exponential explosion or decay of the norm of intermediate layers. It is shown that for random fully connected vanilla ReLU networks, the variance of the squared norm of the activations exponentially increases, even when initializing with the 2 f an-in initialization. For ResNets, this failure mode can be overcome by correctly rescaling the residual branches. However, it is not clear how such a rescaling affects the back propagation of gradients. In our work, we prove an equivalence relation between the forward-and back-propagation of norms in most random ReLU networks. In addition, we present an analysis of DenseNets, revealing a surprising norm preservation quality that is independent of depth. Hanin (2018) explores the conditions at initialization that give rise to the exploding or vanishing gradient problem, by analyzing the individual entries of the input-output Jacobian. In the case of a single output network, the layer k to output Jacobian is given by the row vector:

EQUATION

with n k denoting the width of layer k. In our work, we take a more direct approach to analyzing the full expression for the Jacobian of random networks given by:

EQUATION

Moreover, the analysis presented by Hanin (2018) does not directly carry over to ResNets or DenseNets.

3. PRELIMINARIES AND NOTATIONS:

We make use of the following notations: N denotes a neural network architecture with intermediate outputs y l of layer l for a fixed input y 0 . y l i denotes the i'th component of the vector y l , and n 1 ...n L denote the width of the corresponding layers, with n the length of the input vector. • 2 is the squared Euclidean norm, and we assume y 0 2 = 1. We denote the weight matrix associated with layer l by W l ∈ R n l-1 ×n l , with lower case letters w l ij denoting the individual components of W l . Additional superscripts W lk are used, when several weight matrices are associated with layer l, for instance in residual, CR and DenseNet architectures. Throughout the paper, we assume that the weights are sampled i.i.d from some symmetric distribution with moments E((

w l ij ) m ) = c l m , such that ∀ ijl , ∀ l , c l 1 = 0, c l 2 = 2α n l-1 .

We further identify the following property of a distribution ∆ = 2(c l 4 -3(c l 2 ) 2 )/n l-1 (for Gaussian distributions, ∆ = 0). The derivative of the output unit y L with respect to w l ij is given by J

l ij = ∂y L ∂w l ij = ∂y L ∂y l ∂y l ∂w l ij , with J l 2 = ij (J l ij ) 2 .

We will be concerned in our analysis with the quantities E( y l 2 ), var( y l 2 ), E( J l 2 ) and var( J l 2 ), where the expectations are taken across the random sampling of the weights, given the fixed input y 0 .

Multi Pathway Architectures:

We compare a vanilla fully-connected network to ResNets He et al. (2016) and DenseNets Huang et al. (2016) to a network that is built with CR Shang et al. (2016) layers. For the standard vanilla fully-connected network, the output for of layer l is given by y l = φ(W l y l-1 ), (8) where φ() denotes the ReLU non-linearity.

For residual architectures, we study models with m layer residual branches of the form

y l = y l-1 + R m (y l-1 )

(9) where:

R m (y l-1 ) = W lm ...φ(W l2 φ(W l1 y l-1 ).

(10) We denote the k'th intermediate output of the l'th residual branch by ỹl-1,k = W lk ...φ(W l2 φ(W l1 y l-1 ),

(11) so that, for example

R m (y l-1 ) = W lm ỹl-1,m-1 = W lm φ(W lm-1 ỹl-1,m-2 ).

(12) Furthermore, since the last layer of any residual branch W lm does not end with a non-linearity, we assume c lm 2 = 1 n . DenseNets were recently introduced, demonstrating faster training, as well as improved performance on several popular datasets. The main architectual features introduced by DenseNets include the connection of each layer output to all subsequent layers, using concatenation operations instead of summation, such that the weights of layer l multiply the concatenation of the outputs y 0 ...y l-1 , and is of dimensions W l ∈ R ln×n . For our analysis, we break W l into l -1 weight matrices W l1 ...W ll-1 ∈ R n×n , so that the output of layer l is given by y l = φ( l-1 l =1 W ll y l ).

4. FORWARD-BACKWARD NORM PROPAGATION DUALITY:

In this section, we introduce a link between the propagation of the norm of the activations, and the norm of the Jacobian in different layers in random ReLU networks of finite width. This link will then allow us to study the statistical properties of the gradient in general architectures incorporating residual connections and concatenations with relative ease. Specifically, we would like to establish a connection between the first and second moments of the squared norm of the output y L 2 , and those of the Jacobian norm J k 2 . The family of architectures we consider is restricted, so that the ReLU activations are only applied after multiplication with weights. Using a path-based notation, any unit t of the output y L t can be decomposed to paths that go through weight matrix W k , denoted by y L,k t (neglecting an additional superscript for ResNets and DenseNets), and paths that skip W k , denoted by ŷL,k t . Assuming the input ∀ j , y 0 j = 1, we have:

y L t = y L,k t + ŷL,k t = k∈γt |γt| l=1 w l γt,l z γt + k / ∈γt |γt| l=1 w l γt,l z γt (13)

where the summation is over paths from input to unit t of the output, indexed by γ t , with |γ t | denoting the length of the path indexed by γ t . In non-residual networks, we have |γ t | = L and γt sums over L-1 l=0 n l paths. The term |γt| l=1 w l γt,l z γt denotes the product of weights along path γ t , multiplied by a binary variable z γt ∈ [0, 1], indicating whether path γ t is active. Finally, k / ∈γ indicates summation over all possible paths that do not include a weight from layer k. In all relevant architectures, the squared norm of the Jacobian entry J k ij 2 is given by:

J k ij 2 = ∂y L,k ∂w k ij 2 = n L t=1 ijk∈γt 1 w k ij |γt| l=1 w l γt,l z γt 2 (14)

That is, the derivative of output unit t with respect to weight w k ij is given by the sum of all paths γ that contain w k ij divided by w k ij . In order to link J k with y L , we make the following definition: Definition 1. Reduced network: the outputs of a reduced network N (k) denoted by y 1 (k) ...y L (k) given input y 0 are obtained by removing all connections bypassing weight W k from the network N .

Note that for vanilla networks, it holds that N (k) = N , and ∀ 0<l≤L , y l (k) = y l . The following Lemma links between the moments of y L (k) and those of y L : Lemma 1. For any architecture described in Sec. 3, it holds that for m ∈ {2, 4}:

E y L (k) m = E y L,k m (15)

Note that if network N contains connections bypassing layer k, then for a specific realization of the weights, the equality y L (k) = y L,k does not hold in general, due to different activation patterns in N and N (k) induced by additional residual connections. In other words, paths that are open in N (indicated by the variables z γt ) are not necessarily so in N (k) , and vice versa. Lem. 1, however, states that the moments of both are equal in the family of considered ReLU networks. The following theorem relates the moments of J k 2 with those of y L (k) : Theorem 1. For any architecture described in Sec. 3,the following hold:

EQUATION

where c k 2 and c k 4 are the second and fourth moments of the weights in layer k. Thm. 1 indicates that we can compute the moments of the Jacobian norm of layer k, by computing the moments of the output of the reduced network y L (k) . Thm. 1 also reveals a surprising property of the gradients in general Relu networks. That is, when the weights are sampled from the same distribution in each layer, the gradient's magnitude for each layer are equal in expectation, and depend only on the output statistics.

4.1. VANILLA NETWORKS:

In order to compare our results for ResNets and DenseNets, we first derive the moments of y l (k) in our setting for vanilla architectures. Note that for vanilla networks, the reduced network matches the original, and so ∀ 0<l≤L , y l (k) = y l . Lemma 2. For a vanilla ReLU neural network, as described in Sec, 3, it holds for α = 1:

EQUATION

)

where C 1 , C 2 , C 3 , C 4 are constants that do not depend on L.

By setting α = 1, the variance of the Jacobian squared norm J k 2 will increase exponentially with L l=1

5 n l , as noted by Hanin & Rolnick (2018) . This can be mitigated by linearly increasing the width of each layer with depth, so that n l = nl, resulting in constant asymptotic behavior exp L-1 l=1

1 nl ∼ L 1 n

. This, however, would require significant resources, since the weight matrices would grow quadratically with depth. For comparison, in DenseNets the variance is maintained with depth, while the number of weights increases linearly.

4.2. RESIDUAL NETWORKS:

Residual networks have reintroduced the concept of bypass connections, allowing the training of deep and narrow models with relative ease. Hanin & Rolnick (2018) showed shown that the fluctuations of the norm of subsequent layers do not exponentially increase, provided that the weights are initialized, such that the expected norm E( y l 2 ) is bounded for any layer l. Such an initialization is not constant, and is dependant on the depth of the network L. It is also worth mentioning that the same depth dependent initialization for ResNets was found by Zhang et al. (2019) to be crucial for training convolutional ResNets without applying batchnorm, as introduced in the original formulation of ResNets. Recall that in our analysis, we consider residual branches with m layers with weights W kk , 0 < k ≤ L, 0 < k ≤ m,, and so we will analyze the variance of

J kk 2 = ij ∂y L ∂w kk ij 2

for k < m. Unlike vanilla networks, the outputs of the reduced network y l (k) for l > k are not equal to those of the full network y l , and are given by:

y l (k) = { R m (y l-1 ) + y l-1 l = k R m (y l-1 ) l = k (18)

The following theorem states our results for ResNets.

Theorem 2. For a ResNet as in Sec, 3 with m ≥ 2, it holds that for 0

< k ≤ L, 0 < k < m: C 1 ρ m -exp 1 2 Lρ m --C 2 µ 2 ≤ var J kk 2 ≤ C 3 ρ m + exp 2Lρ m 2 + -C 4 µ 2 (19)

where

ρ ± = α 2 (1 + 5±|∆| n ), C 1 , C 2 , C 3 , C 4 are constants that do not depend on L, m, and µ = α m 1 + α m L-1 .

The gradients of weights in residual branches are, therefore, sensitive to the depth of each residual branch (large m). However, unlike their vanilla counterparts, by scaling the residual branches, such that Lρ m + is kept constant, ensures that both E( J 2 ) and var( J 2 ) remain bounded as the overall depth of the network grows (increasing L). Note that a similar initialization was proposed in Zhang et al. (2019) for training ResNets without batchnorm, albeit for a more restricted class of networks, and without an explicit dependency on width. Our result indicates that for narrow and deep networks with residual connections, training without batchnorm requires scaling in a both depth and width dependent manner.

4.3. DENSELY CONNECTED NETWORKS:

DenseNets have shown improvement over ResNets in terms of accuracy and speed of training on several notable image classification tasks, since their inception. As far as we can ascertain, there have been no theoretical studies on the causes of DenseNets success, or the implications on the training process of its architectural novelties, namely dense connections and feature reuse by concatenation of previous layers. In the following, it is shown that the DenseNet architecture, coupled with the standard, depth independent 2 n l-1 initialization, gives rise to well-behaving gradients in any depth. By the definition of DenseNets given in Sec. 3, the reduced architecture outputs y l (k) are given by:

y l (k) = { φ( l-1 l =1 W ll y l ) l ≤ k φ( l-1 l =k W ll y l ) l > k (20)

The following theorem states our results for DenseNets: Theorem 3. For a DenseNet as described in Sec, 3, it holds that for 0 < k ≤ L:

L (α-1) C 1 exp -2|∆|L n -C 2 ≤ var( J k 2 ) ≤ L 2(α-1) C 3 exp 2|∆|L n -C 4 (21)

where C 1 , C 2 , C 3 , C 4 are constants that do not depend on L.

It is worth noting the difference between DenseNets, and both vanilla and ResNet architectures. In vanilla architectures, the mean of the activation norm, as well as the gradient norm are exponential in α. However, the variance of the same quantities is exponential in

L l=1 1 n l .

That is, fixing the mean by correctly setting α = 1 does not prevent the gradients and activations from exploding in deep networks. In ResNets, scaling of the residual branches in order to avoid exploding gradients requires scaling in a depth dependent fashion. However, in DenseNets, the variance of the squared Jacobian norm is in fact polynomial in depth for Gaussian distributions. Constant mean and variance is achieved (at least when ∆ = 0) by setting α = 1, for any depth. Moreover, the constant variance result is achieved with a linear growth rate of the weight complexity, instead of quadratic with vanilla networks.

Empirical support We employed Monte Carlo simulations in order to verify our theoretical findings. As shown in Fig. 2 , there is an excellent match between the simulations and our derivations.

5. THE CURIOUS CASE OF CONCATENATED RELU:

Our results stated above demonstrate how different architectures and initialization schemes affect the training process at the start of training, by analyzing the expected norm and variance of the gradients for each layer, and hold for ReLU based architectures. Since a similar analysis can be performed on linear random networks, it is interesting to see the effect of the ReLU non-linearity on the outcome. As it turns out, the ReLU function has a detrimental effect on random networks, as shown in the following lemma, effectively exacerbating the fluctuations of both the activations and the gradients:

Lemma 3. Given a vanilla linear neural network with w l ij ∼ N (0, 1 n l-1 ), then it holds that:

EQUATION

Linear networks, of course, have limited capacity, considerably limiting their practical use. Fortunately, we can use the ReLU function to construct computation blocks that behave as linear networks at the initialization stage, while not compromising in terms of capacity. This can be achieved by the use of CR units. In concatenated ReLUs blocks, the input signal x is split into two complementary components x + and x -such that the output of layer l is parameterized by two weight matrices W l1 , W l2 , and is given by y l = W l1 φ(y l-1 ) -W l2 φ(-y l-1 ). An algebraically equivalent way to write the layer's computation is y l = W l2 y l-1 + (W l1 + W l2 )φ(y l-1 ), which can be viewed as a variant of a residual connection. However, it is important to note that this equivalence holds for the forward-pass only, since deriving by W l1 and W l2 during the optimization process is not the same as performing an optimization that is based on a derivation by W l1 -W l2 and W l2 . The equivalence at initialization between the CR based architectures and random linear networks is depicted in the following Lemma: Lemma 4. The distribution of the output y l = W l1 φ(y l-1 ) -W l2 φ(y l-1 ) conditioned on y l-1 , where w l1 ij , w l2 ij ∼ N (0, α n l-1 ) equals the distribution of y = W l y l-1 , where the elements of W l are distributed according to an i.i.d normal distribution w l ij ∼ N (0, α n l-1 ).

Using CR, therefore, provides initial conditions identical to random linear networks in terms of the output statistics. However, this equivalence also holds for back propagation. Indeed, from the construction of W l in the lemma, it is easy to see that ( 5 n l , hence more efficient than standard ReLU blocks. We can similarly define residual and densely connected architecture comprised of concatenated ReLU blocks. In that case, the conditions of both architectures in terms of output distributions at initialization will be identical to the linear counterparts of both architectures, given by y l = W lk1 φ(ỹ l,k-1 ) + W lk2 φ(-ỹ l,k-1 ) and y l = l-1 l =0 W ll y l for ResNets and DenseNets, respectively. Further evidence for the practical effectiveness of CR layers is provided in the appendix, in which we present experiments on the UCI datasets.

6. CONCLUSIONS:

The interplay between architecture and initialization determines much of the early dynamics of training deep networks. In this work, we have demonstrated how careful initialization in general residual architectures effects the gradients at the start of training, by a rigorous study of the variance of the per layer gradient norm at the point of initialization. Through the use of the duality principle between forward and backward propagation, we have derived the expression for the variance in three architectures: standard networks, ResNets and DenseNets. We show that DenseNets benefit from a simple variance preserving initialization, which is independent of depth, while a depth dependent initialization exists for ResNet. We have also shown that CR blocks perform better then ReLU blocks, given a fixed budget of weights, when it comes to norm preservation in deep layers. As future work, we would like to extend our analysis to include normalization techniques, such as batchnorm.

7. PROOFS:

We make use of the following propositions and definitions in order to prove lemma 1 and theorem 1. For ResNet, we assume m = 2 for notation simplicity. The results hold for m > 2 as well with no additional arguments.

A path γ t from input to output unit t defines a product of weights along the path denoted by

EQUATION

where |γ t | is the length of path γ t (for vanilla and DenseNets |γ t | = L), z γt is a binary variable z γt ∈ {0, 1} indicating if the path is "active" (i.e, all relevant Relu activations along the path are open), and w γt,l are the weights along path γ t , and q γt,l = w γt,l z γt .

We generalize the above expression the following way:

EQUATION

where:

EQUATION

where l / ∈ γ t indicates if layer l is skipped, z γt,l is the activation variable relevant for weight wγt,l . That is, if wγt,l = w l ij then z γt,l = z l j where z l j ∈ {0, 1} is the activation variable unit j in layer l. wγt,l denotes a generalized weight that can take the value of a weight wγt,l = w γt,l for vanilla layers, or a product of weights wγt,l = w γt,l1 w γt,l2 for residual layers (w γt,l1 ∈ W l1 and w γt,l2 ∈ W l2 ).

proposition 1. Given a random vector w = [w 1 ...w n ] such that each component is identically and symmetrically distributed iid random variable with moments E(w m i ) = c m , and random binary variable z such that z|w ∈ {0, 1}, z|w = 1 -z| -w, then it holds that for n > l:

EQUATION

Proof. We have: For even m or even l: it follows that:

c l m = w l i=1 w m i p(w)dw =

EQUATION

and so:

E( l i=1 w m i z) = c l m 2 (31)

For odd l, since w is symmetrically distributed, we have that c l m = 0, and so E(

l i=1 w m i z) = c l m 2 = 0

proposition 2. For any set of paths γ 1 t , γ 2 t , it holds that:

E(Q γ 1 t Q γ 2 t ) = 0 γ 1 t = γ 2 t E(Q 2 γ 1 t ) = L l=1 E (q γ 1 t ,l ) 2 γ 1 t = γ 2 t (32)

Proof. From proposition 2 we have:

EQUATION

Since the computations done by all considered architectures form a markov chain, such that the output of layer y l depends only on weights R l = {W 1 ...W l }, using conditional expectation,:

E(Q γ 1 t Q γ 2 t ) = E E(Q γ 1 t Q γ 2 t R L-1 ) = E L-1 l=1 q γ 1 t ,l q γ 2 t ,l E(q γ 1 t ,L q γ 2 t ,L R L-1 ) (34)

Note that since both are paths leading to output unit t, the following holds:

EQUATION

)

and so it follows from proposition 1:

E(q γ L t ,l q γ 2 t ,L R L-1 ) =      E( w2 γ 1 t ,L z L t ) L ∈ γ 1 t , L ∈ γ 2 t , wγ 1 t ,L = wγ 2 t ,L 1 L / ∈ γ 1 t , L / ∈ γ 2 t 0 else = E (q γ 1 t ,L ) 2 (36)

Recursively going through l = l -1...1 completes the proof.

proposition 3. For any set of paths γ 1 t , γ 2 t , γ 1 s , γ 2 s , from input to output units t and s, it holds that:

EQUATION

such that:

E ξ γ 12 st ,l =      E q 2 γ 1 t ,l q 2 γ 1 s ,l q γ 1 t ,l = q γ 2 t ,l , q γ 1 s ,l = q γ 2 s ,l E q 2 γ 1 t ,l q 2 γ 2 t ,l q γ 1 t ,l = q γ 1 s ,l , q γ 2 t ,l = q γ 2 s ,l 0 else (38)

Proof. The proof is practically identical to proposition 3, and will be omitted.

For notation simplicity,

7.1. PROOF OF LEMMA 1:

The proof of lemma 1 stems from the following observation: Since the computations done by all considered architectures form a markov chain, such that the output of layer y l depend only weights W 1 ...W l , then for any particular realization of the weights, y L,k and y L k differ only in their activation patterns in layers l ≥ k due to the removal of bypassing connections in y L (k). Using Eq. 23:

EQUATION

while:

E y L (k) m = n L t=1 E ( γt|k∈γt |γt| l=1 w γt,l z γt ) m (40)

We must therefor prove the expectation is invariant to a switch between z γt and z γt for layers l > k.

For m = 2, using proposition 2, the expectation factorizes:

E y L,k 2 = n L t=1 γ 1 t |k∈γt E(Q 2 γ 1 t ) = n L t=1 γ 1 t |k∈γt L l=1 E (q γ 1 t ,l ) 2 (41)

Similarly for m = 4, from proposition 3 the expectation factorizes:

EQUATION

where ξ γ 12 st ,l is defined is Eq. 38. Note that switching the sign of W l for all architectures considered will switch the state of both z l j , 0 < j ≤ n l and z l j , 0 < j ≤ n l , and so from proposition 1, the factorized expectations for m = 2 and m = 4 meet the required invariance.

7.2. PROOF OF THEOREM 1:

From the definition of the Jacobian for weights W k , we have using lemma 1 and proposition 2:

EQUATION

)

Since by definition

EQUATION

)

where the last transition stems from proposition 1. Plugging back into Eq. 43 and summing over ij completes the proof for m = 2. For m = 4, we have:

EQUATION

where from proposition 1 and proposition 3, for uvk

EQUATION

Note that since c k 4 ≥ (c k 2 ) 2 , it holds that:

EQUATION

Plugging into Eq. 46, it follows:

EQUATION

Furthermore, it holds that:

EQUATION

Therefor for large width networks:

EQUATION

EQUATION

EQUATION

and so for large width networks:

EQUATION

7.3. PROOF OF LEMMA 2:

Denoting by Z l the diagonal matrix holding in its diagonal the activation variables z l j for unit j in layer l, and so we have:

EQUATION

Conditioning on R l-1 = {W 1 ...W L-1 } and taking expectation:

EQUATION

From proposition 1, it follows that:

EQUATION

Recursing through l = L -1...1:

EQUATION

Similarly:

EQUATION

From proposition 3:

EQUATION

and so:

EQUATION

EQUATION

finally:

EQUATION

Recursing through L -1...1, and denoting

C 1 = n L (n L +5-|∆|) n 2

, and

C 2 = n L (n L +5+|∆|) n 2

and combining with Eq. 62:

EQUATION

Using ∀ 1≥x≥-1 , e 0.5x ≤ (1 + x) ≤ e x :

EQUATION

For Gaussian distribution we have that ∆ = 0, and so:

EQUATION

7.4 PROOF OF THEOREM 2

In the following proof, for the sake of notation simplicity, we omit the notation k in y l (k) , and assume that y l stands for the reduced network y l (k) . We have for layer L:

E y L 2 = E y L-1 2 + E ỹL-1,m-1 W Lm W Lm ỹL-1,m-1 = E y L-1 2 + αE ỹL-1,m-1 2 = E y L-1 2 1 + α m = E y k 2 L l=k+1 1 + α m = E y k-1 2 L l=k+1 1 + α m α m = l =k 1 + α m α m (77) E y L 4 = E y L-1 4 + E ỹL-1,m 4 + 4E (ỹ L-1,m y L-1 ) 2 + 2E ỹL-1,m 2 y L-1 2 (78)

We now handle each term separately:

EQUATION

Similarly to the vanilla case with a linear weight layer on top, we have:

EQUATION

EQUATION

and:

E ỹL-1,m 2 y L-1 2 = α m E y L-1 4 (84)

Plugging it all into Eq. 78, and denoting:

EQUATION

EQUATION

we have after recursing through l = L -1...1:

E y k 4 L l=k+1 β - l ≤ E y L 4 ≤ E y k 4 L l=k+1 β + l (87)

Since for the reduced architecture y L (k) , the transformation between y k-1 and y k is given by an m layer fully connected network without any residual connections, we can use the results from the vanilla case:

L l =k β - l α 2m (1 + 5 -|∆| n ) m-1 ≤ E y L 4 ≤ L l =k β + l α 2m (1 + 5 + |∆| n ) m (88) Denoting ρ ± = α 2 (1 + 5±|∆| n ) it follows: C 1 L l =k β - l ρ m -≤ E y L 4 ≤ L l =k β + ρ m + (89)

From the definition of β ± l , it holds for ρ + < 1:

EQUATION

)

and:

EQUATION

)

and so:

C 1 ρ m -exp 1 2 Lρ m --C 2 µ 2 ≤ var y L 2 ≤ C 3 ρ m + exp 2Lρ m 2 + -C 3 µ 2 (92)

For Gaussian distributions, using a similar derivation and assigning ∆ = 0 we have:

EQUATION

where:

EQUATION

and:

EQUATION

7.5 PROOF OF THEOREM 3

In the following proof, for the sake of notation simplicity, we omit the notation k in y l (k) , and assume that y l stands for the reduced network y l (k) . We have:

EQUATION

Note that the above is a simple Relu block with concatenated y 1 ...y L -1 as input. Using the results from the vanilla architecture, and denoting C l,l = E y l 2 y l 2 , it then follows:

EQUATION

From Eq. 98, it holds that:

EQUATION

It then follows:

EQUATION

(101) The following also holds:

EQUATION

and so:

EQUATION

EQUATION

EQUATION

Telescoping through l = L -1...1:

EQUATION

A similar relation holds for the mean:

EQUATION

Finally:

EQUATION

A lower bound can easily be derived using Eq. 98, 99, and follow the same derivation by flipping the sign of |∆|:

EQUATION

For Gaussian distributions, using a similar derivation and assigning ∆ = 0 we have:

EQUATION

and so:

var( y L 2 ) = a 2 n + 5 nk 2 L l =k 1 + 5α 2 n(l + α -1) 2 - a 2 k 2 L l =k 1 + α -1 l 2 (122) 7.6 PROOF OF LEMMA 3

The proof is identical to the Relu case, removing the activation variables Z l , and hence will be omitted.

7.7. PROOF OF LEMMA 4:

Proof. We construct the matrix W l as follows:

EQUATION

And so y l = W l y l-1 . Note that the elements of W l are i.i.d Gaussian random variables w l ij ∼ N (0, α n l-1 ).

8. EMPIRICAL SUPPORT FOR CONCATENATIVE RELU:

Our theoretical results are verified in Fig. 3 . As can be seen, the simulations closely match the predictions. The advantage of CR over conventional ReLU networks is also clearly visible.

CR has fallen out of flavor and was never shown to be effective outside CNNs. We, therefore, provide empirical support for its effectiveness, to which our analysis above points. In our experiments, we evaluate the performance of fully connected networks on the 40 UCI datasets with more than 1000 samples. The train/test splits were provided by Klambauer et al. (2017) . In the experiment, we tested three different depths: 8, 16, and 32. Each layer contained n = 256 hidden neurons. Dropout was not employed. A learning rate of either 0.1, 0.01, or 0.001 was used for the first 50 epochs and then a learning rate of one tenth for an additional 50 epochs. All runs were terminated after 100 epochs. Batches were of size N = 50.

Following Klambauer et al. (2017) , an averaging operator with a mask size of 5 was applied to the validation error, and the epoch with the best smoothed validation error was selected. This was done separately for each method, each dataset, and each split. Out of the nine options of depth and learning rate for most methods, the mean validation error across the four splits was used to select the option for which we report the results on the test splits. We compared the following activation functions and architectures: vanilla neural networks (single pathway) with either SELU (Klambauer et al., 2017) or ELU (Clevert et al., 2015) activation functions, ResNets, and CR based architecture. The single pathway ReLU networks were not competitive and the experiments with this architecture were terminated to allow for more experiments with the other architectures. The single pathway experiments were done without batchnorm, as is usually done with ELU and SELU. For ResNets, we report results both with batchnorm and without. Batchnorm was applied along the residual pathway, twice, as is typically practiced. The first and last layers of the ResNet and CR architecture, i.e., the layers that project from the input and to the output layer, are conventional projections that do not employ multiple pathways. For the baseline architectures, we also evaluate using double width, i.e., is the threshold (τ ). Note that since for accuracy scores, higher is better, whereas typical Dolan-More plots show cost, the axis is ordered in reverse. The y-axis depicts, for a given architecture out of the nine, the ratio of datasets in which the obtained accuracy is above the threshold τ times the maximal accuracy obtained by any architecture . In Fig. 4 , we provide the Dolan-More profile of the methods. In this plot, the x-axis is a threshold τ and the y-axis is the ratio, out of all datasets, in which the method has reached an accuracy that is τ times the best accuracy for the dataset. It can be seen that the CR method dominates across all threshold values, except for SELU double width at low values of τ .

9. UCI RESULTS:

In Tab. 1 we provide the full numbers of the UCI experiments in the paper.

:



Figures and Tables:

Figure 1: An illustration of Thm. 1. The activations of the network in (a) are completely different from those of the network in (b), in which all skip connections bypassing layer k = 2 are removed. However, the moments of the gradient norms at layer k = 2 are exactly the same in both (a) and (b).

Figure 2: The variance of the squared norm in log scale of the per layer Jacobian as a function of α (held constant for all the layers) obtained from the simulated results of 200 independent runs, for (a) ResNet with m = 2 and n = 40, and (b) DenseNet with n = 40. The x-axis represents the depth of the network, and y-axis the variance as computed from 200 runs. The dashed lines represent our theoretical predictions for both architectures, and full lines represent the empirical results. All networks were initialized using Gaussian distributions.

Figure 3: Activations norm per layer, as obtained from simulated results of 5000 independent runs. (a) The mean of the squared norm of the activations in each layer (y-axis) as a function of depth (x-axis) for randomly initialized networks. The red plot is for a ReLU network, the blue for a CR based network. (b) the y-axis is the standard deviation of the squared norm of the activation. In addition to the blue and red plots, the black and green plots are the theoretical predictions according to Lem. 2,3.

Lem. 4  shows that similar to vanilla ReLU based networks, CR based ReLU networks have exponential behavior in L l=1 2 n l . Since CR employ twice as many parameters, then for the same budget of parameters (having n l = n l

Figure4: A Dolan-More profile, based on the accuracy obtained for the UCI experiments. The x-axis is the threshold (τ ). Note that since for accuracy scores, higher is better, whereas typical Dolan-More plots show cost, the axis is ordered in reverse. The y-axis depicts, for a given architecture out of the nine, the ratio of datasets in which the obtained accuracy is above the threshold τ times the maximal accuracy obtained by any architecture .

Error rates for each method per dataset. 512 denote the double width settings. 512 hidden neurons. This was done since CR employs more weights than the baselines due to the dual pathway architecture.

Table content: <table><tr><td/><td>CR</td><td>ELU ELU 512 SELU SELU 512 RES RES 512 RES RES 512</td></tr><tr><td/><td/><td>w/o bn w/o bn</td></tr><tr><td>semeion</td><td colspan="2">0.102 0.1357 0.118 0.1543 0.1356 0.1192 0.1215 0.121 0.1308</td></tr><tr><td>oocytes..._2f</td><td colspan="2">0.0941 0.0918 0.0949 0.0938 0.0892 0.0897 0.0921 0.0944 0.1084</td></tr><tr><td>wine-quality-white</td><td colspan="2">0.4532 0.4633 0.4488 0.4597 0.4547 0.4521 0.4552 0.4573 0.4645</td></tr><tr><td>ringnorm</td><td colspan="2">0.0208 0.0253 0.0251 0.026 0.0222 0.0261 0.0220 0.0205 0.0341</td></tr><tr><td>statlog-image</td><td colspan="2">0.0448 0.0483 0.0472 0.042 0.0495 0.0467 0.0639 0.0489 0.0382</td></tr><tr><td colspan="3">cardiotocography-10 0.2090 0.2207 0.212 0.2222 0.2073 0.2147 0.2134 0.2126 0.2141</td></tr><tr><td>abalone</td><td colspan="2">0.3408 0.3439 0.3391 0.3476 0.3492 0.3452 0.3479 0.3464 0.344</td></tr><tr><td>ozone</td><td colspan="2">0.0295 0.0300 0.0287 0.0339 0.0288 0.0287 0.0312 0.0288 0.0321</td></tr><tr><td>thyroid</td><td colspan="2">0.0207 0.0179 0.0183 0.0211 0.0198 0.0192 0.0186 0.018 0.0247</td></tr><tr><td>oocyte..._4d</td><td colspan="2">0.1855 0.1946 0.194 0.2052 0.2076 0.2021 0.2006 0.2047 0.2156</td></tr><tr><td>adult</td><td colspan="2">0.1515 0.1539 0.1530 0.1559 0.1550 0.1515 0.1509 0.1507 0.1538</td></tr><tr><td>chess-krvkp</td><td colspan="2">0.0215 0.0333 0.0257 0.0307 0.0255 0.0205 0.0262 0.0275 0.0178</td></tr><tr><td>waveform</td><td colspan="2">0.1403 0.1419 0.1441 0.1468 0.1455 0.1355 0.1399 0.1377 0.1656</td></tr><tr><td>plant-texture</td><td colspan="2">0.2614 0.3412 0.2991 0.3022 0.2777 0.2814 0.2704 0.2738 0.3296</td></tr><tr><td>statlog-landsat</td><td colspan="2">0.1061 0.1025 0.1011 0.1011 0.0987 0.1053 0.1111 0.1069 0.0960</td></tr><tr><td>steel-plates</td><td colspan="2">0.275 0.2933 0.2815 0.2832 0.2776 0.2961 0.2879 0.2881 0.2926</td></tr><tr><td>statlog-vehicle</td><td colspan="2">0.2188 0.2115 0.2159 0.2213 0.2233 0.2198 0.2285 0.2164 0.2075</td></tr><tr><td>cardiotocography-3</td><td colspan="2">0.0967 0.1057 0.1022 0.1023 0.1008 0.0988 0.1053 0.0994 0.1055</td></tr><tr><td>page-blocks</td><td colspan="2">0.0324 0.0367 0.0345 0.0367 0.0355 0.0356 0.0343 0.0376 0.0367</td></tr><tr><td>wall-following</td><td colspan="2">0.1131 0.1232 0.1272 0.1269 0.1235 0.1204 0.1236 0.1190 0.1192</td></tr><tr><td>miniboone</td><td colspan="2">0.0793 0.0751 0.0750 0.0790 0.0790 0.0849 0.0889 0.0772 0.2481</td></tr><tr><td>statlog-shuttle</td><td colspan="2">0.0016 0.002 0.0017 0.0015 0.0015 0.0025 0.0025 0.0022 0.0017</td></tr><tr><td>mushroom</td><td colspan="2">0.0004 0.0004 0.0003 0.0010 0.0003 0.0010 0.0008 0.0007 0.0001</td></tr><tr><td>waveform-noise</td><td colspan="2">0.1464 0.1501 0.1520 0.1626 0.1576 0.1437 0.1421 0.7090 0.1705</td></tr><tr><td>twonorm</td><td colspan="2">0.0231 0.027 0.0252 0.0261 0.0258 0.0236 0.0243 0.0241 0.0341</td></tr><tr><td>connect-4</td><td colspan="2">0.1442 0.1508 0.1493 0.1483 0.1442 0.1423 0.1431 0.1424 0.1628</td></tr><tr><td>letter</td><td colspan="2">0.0488 0.0700 0.0650 0.0606 0.0583 0.0880 0.0778 0.0818 0.0485</td></tr><tr><td>titanic</td><td colspan="2">0.2197 0.2148 0.2151 0.2170 0.2118 0.2149 0.2181 0.212 0.2149</td></tr><tr><td>statlog-heart</td><td colspan="2">0.1844 0.1939 0.1808 0.2093 0.2113 0.1967 0.2016 0.1794 0.2276</td></tr><tr><td>plant-margin</td><td colspan="2">0.2798 0.3330 0.2835 0.3258 0.2918 0.2772 0.2764 0.2951 0.3749</td></tr><tr><td>contrac</td><td colspan="2">0.4676 0.4757 0.4716 0.4816 0.4852 0.4606 0.4622 0.4693 0.4889</td></tr><tr><td>yeast</td><td colspan="2">0.4239 0.4205 0.4252 0.4333 0.4253 0.4141 0.4199 0.4431 0.4649</td></tr><tr><td>bank</td><td colspan="2">0.1121 0.1114 0.1131 0.1180 0.1170 0.1084 0.1161 0.1099 0.1143</td></tr><tr><td>hill-valley</td><td colspan="2">0.3667 0.3699 0.3774 0.3700 0.3554 0.3899 0.3657 0.3765 0.4062</td></tr><tr><td>optical</td><td colspan="2">0.0299 0.0456 0.0408 0.0416 0.0417 0.0477 0.0475 0.0426 0.0366</td></tr><tr><td>plant-shape</td><td colspan="2">0.4127 0.5439 0.4862 0.4300 0.4220 0.4593 0.4581 0.4509 0.4228</td></tr><tr><td>led-display</td><td colspan="2">0.2906 0.2929 0.2827 0.2911 0.3060 0.2852 0.2856 0.2781 0.3093</td></tr><tr><td>magic</td><td colspan="2">0.1293 0.1310 0.1331 0.1345 0.1344 0.1375 0.1348 0.1329 0.1332</td></tr><tr><td>pendigits</td><td colspan="2">0.0388 0.0446 0.0434 0.0426 0.0432 0.0466 0.0483 0.0467 0.0439</td></tr><tr><td>chess-krvk</td><td colspan="2">0.2452 0.2592 0.2426 0.2399 0.2284 0.2956 0.2848 0.2741 0.3289</td></tr><tr><td>each layer had</td><td/></tr></table>

Peer review (iclr2020):
This paper analyzes the statistics of activation norms and Jacobian norms for randomly-initialized ReLU networks in the presence (and absence) of various types of residual connections. Whereas the variance of the gradient norm grows with depth for vanilla networks, it can be depth-independent for residual networks when using the proper initialization. The main theoretical results stem from a norm propagation duality that allows computations of certain statistics using simplified architectures. As far as I know, this is a novel way of analyzing multi-pathway architectures, and the method and the implied results will be of interest to deep-learning theorists. There are two types of empirical results, one to verify the theoretical predictions, and the other to see if those predictions translate into improved generalization performance on some real-world tasks. Figures 2 and 3 show good agreement with theory. I am less convinced by Fig 4. Unless I missed something, I thought the implication from the theoretical analysis was that multi-pathway models should behave better than single-pathway models at large depth, and for single-pathway models concatenated ReLU may have slightly better large-depth behavior than standard ReLU. Doesn*t this suggest that for large depth DenseNet ~= ResNet > CR > ReLU? Is the good performance of CR supposed to be understood as supporting evidence for the theoretical analysis? Or is the point more like *the theory suggests CR might be interesting, and indeed it performs well, but we leave the details of generalization performance to a future study*? It would be important to specify the desired interpretation here. But I think regardless of the interpretation, the empirical results are not strong enough to stand on their own as support for CR, especially since this type of architecture has been introduced and studied previously (c.f. also the looks-linear setup from [1]). Therefore my evaluation is mostly based on the theoretical contributions, which I think are themselves are probably sufficient for publication. One main area for improvement is in the related work. A significant line of work is omitted that studies higher-order statistics of Jacobian matrices in the large-depth, infinite-width regime [2-8], and while this is at infinite width, it work for any non-linearity, so is an important point for comparison. Other kinds of initialization schemes are also relevant for comparison, e.g. [9-10]. [1] Balduzzi, David, et al. *The shattered gradients problem: If resnets are the answer, then what is the question?.* Proceedings of the 34th International Conference on Machine Learning-Volume 70. JMLR. org, 2017. [2] Pennington, Jeffrey, Samuel Schoenholz, and Surya Ganguli. *Resurrecting the sigmoid in deep learning through dynamical isometry: theory and practice.* Advances in neural information processing systems. 2017. [3] Pennington, Jeffrey, Samuel Schoenholz, and Surya Ganguli. *The emergence of spectral universality in deep networks.* International Conference on Artificial Intelligence and Statistics. 2018. [4] Hayase, Tomohiro. *Almost Surely Asymptotic Freeness for Jacobian Spectrum of Deep Network.* arXiv preprint arXiv:1908.03901 (2019). [5] Tarnowski, Wojciech, et al. *Dynamical Isometry is Achieved in Residual Networks in a Universal Way for any Activation Function.* The 22nd International Conference on Artificial Intelligence and Statistics. 2019. [6] Burkholz, Rebekka, and Alina Dubatovka. *Initialization of ReLUs for Dynamical Isometry.* arXiv preprint arXiv:1806.06362 (2018). [7] Xiao, Lechao, et al. *Dynamical Isometry and a Mean Field Theory of CNNs: How to Train 10,000-Layer Vanilla Convolutional Neural Networks.* International Conference on Machine Learning. 2018. [8] Chen, Minmin, Jeffrey Pennington, and Samuel Schoenholz. *Dynamical Isometry and a Mean Field Theory of RNNs: Gating Enables Signal Propagation in Recurrent Neural Networks.* International Conference on Machine Learning. 2018. [9]Sutskever, Ilya, et al. *On the importance of initialization and momentum in deep learning.* International conference on machine learning. 2013. [10] Le, Quoc V., Navdeep Jaitly, and Geoffrey E. Hinton. *A simple way to initialize recurrent networks of rectified linear units.* arXiv preprint arXiv:1504.00941 (2015).