summaries: The paper studies the mean and variance of the gradient norm at each layer for vanilla feedforward, ResNet and DenseNet, respectively, at the initialization step, which is related with Hanin & Ronick 2018 studying the mean and variance of forward activations. They show that ResNet and DenseNet preserve the variance of the layer gradient norm through depths. In comparison, for the vanilla feedforward network, although the mean of the gradient norm is preserved if  is properly initialized, the variance of the layer gradient norm increases over depths, which may explode or decay the gradient at deeper layers.  This paper studies the effects of residual and dense net type connections on the moments of per-layer gradients at random initialization. In particular, using duality, bounds on the variance of the square norm of Jacobian (with respect to the randomness of random initialization) are derived for vanilla networks. By noticing that connections that bypass a specific layer does not affect the expected square norm of the gradients of a particular layer, the paper continue to characterize the bounds for residual networks and densely connected networks that have special skip connections. In particular, with properly chosen initialization scales for each layer, the architectures with skip connections can be initialized so that the gradient norm does not explode with increasing depth.
abstractText: 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.
references.authors: Djork-Arné Clevert Thomas Unterthiner Sepp Hochreiter
references.title: Fast and accurate deep network learning by exponential linear units (elus)
references.venue: arXiv preprint arXiv:1511.07289,
references.authors: Boris Hanin
references.title: Which neural net architectures give rise to exploding and vanishing gradients
references.venue: Advances in Neural Information Processing Systems
references.authors: Boris Hanin David Rolnick
references.title: How to start training: The effect of initialization and architecture
references.venue: Advances in Neural Information Processing Systems
references.authors: Kaiming He Xiangyu Zhang Shaoqing Ren Jian Sun
references.title: Delving deep into rectifiers: Surpassing human-level performance on imagenet classification
references.venue: In Proceedings of the IEEE international conference on computer vision,
references.authors: Kaiming He Xiangyu Zhang Shaoqing Ren Jian Sun
references.title: Deep residual learning for image recognition
references.venue: In Proceedings of the IEEE conference on computer vision and pattern recognition,
references.authors: Gao Huang Zhuang Liu Kilian Q. Weinberger
references.title: Densely connected convolutional networks
references.venue: In arXiv:1608.06993,
references.authors: Günter Klambauer Thomas Unterthiner Andreas Mayr Sepp Hochreiter
references.title: Self-normalizing neural networks
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Jeffrey Pennington Samuel Schoenholz Surya Ganguli
references.title: Resurrecting the sigmoid in deep learning through dynamical isometry: theory and practice
references.venue: Advances in Neural Information Processing Systems
references.authors: Jeffrey Pennington Samuel S. Schoenholz Surya Ganguli
references.title: The emergence of spectral universality in deep networks
references.venue: In AISTATS,
references.authors: Samuel S Schoenholz Justin Gilmer Surya Ganguli Jascha Sohl-Dickstein
references.title: Deep information propagation
references.venue: In ICLR,
references.authors: Wenling Shang Kihyuk Sohn Diogo Almeida Honglak Lee
references.title: Understanding and improving convolutional neural networks via concatenated rectified linear units
references.venue: In ICML,
references.authors: Lechao Xiao Yasaman Bahri Jascha Sohl-Dickstein Samuel Schoenholz Jeffrey Pennington
references.title: Dynamical isometry and a mean field theory of CNNs: How to train 10,000-layer vanilla convolutional neural networks
references.venue: Proceedings of the 35th International Conference on Machine Learning,
references.authors: Greg Yang Sam Schoenholz
references.title: Mean field residual networks: On the edge of chaos
references.authors: Hongyi Zhang Yann N. Dauphin Tengyu Ma
references.title: Fixup initialization: Residual learning without normalization
references.venue: In ICLR,
references.authors: Following Klambauer
references.title: 2017), an averaging operator with a mask size of 5 was applied
sections.text: 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.
sections.heading: 1 INTRODUCTION
sections.text: 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 y0 to an L layer neural network with weights given by W , and intermediate output for any layer l denoted by yl. Denoting by wkij the weights ij of layer k, we are interested in the quantities:
  ‖yl‖2, ‖Jk‖2 = ∑ ij ‖ ∂y L ∂wkij ‖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 ‖yl‖2 and ‖Jk‖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 Jk. 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:
  ‖∇WkL‖ = ‖ ∂L ∂yL Jk‖ ≤ ‖ ∂L ∂yL ‖‖Jk‖ (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 ‖yl‖2 and ‖Jk‖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(‖Jk‖2) and var(‖yl‖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.
sections.heading: 2 RELATED WORK
sections.text: 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 ( ( ∂yL ∂yli )2 ) = E ( ( ∂yL ∂yl+1i )2 ) σ2w ∫ 1√ 2π e− z2 2 φ′( √ q∗z)dz (3)
  where φ′ the derivative of the activation function and σ2w denotes the variance of the weights. By initializing with σ2w = 2 fan−in , and considering that for the ReLU activation function φ
  ′(z) = 1z>0, we have that
  E(( ∂yL
  ∂yli )2) = E((
  ∂yL
  ∂yl+1i )2). (4)
  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:
  JIO = ∂yL
  ∂y0 (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 JIO 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 2fan−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:
  J(k → L) = (∂y L
  ∂yli , i = 1...nk) =
  ∂yL ∂yl (6)
  with nk 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:
  Jkij = ∂yL ∂yk ∂yk ∂wkij = 〈J(k → L), ∂y k ∂wkij 〉. (7)
  Moreover, the analysis presented by Hanin (2018) does not directly carry over to ResNets or DenseNets.
sections.heading: 3 PRELIMINARIES AND NOTATIONS
sections.text: We make use of the following notations: N denotes a neural network architecture with intermediate outputs yl of layer l for a fixed input y0. yli denotes the i’th component of the vector y
  l, and n1...nL denote the width of the corresponding layers, with n the length of the input vector. ‖ · ‖2 is the squared Euclidean norm, and we assume ‖y0‖2 = 1. We denote the weight matrix associated with layer l by W l ∈ Rnl−1×nl , with lower case letters wlij 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((wlij)m) = clm,
  such that ∀ijl, ∀l, cl1 = 0, cl2 = 2αnl−1 . We further identify the following property of a distribution ∆ = 2(cl4 − 3(cl2)2)/nl−1 (for Gaussian distributions, ∆ = 0). The derivative of the output unit yL with respect to wlij is given by J l ij = ∂yL
  ∂wlij = ∂y
  L ∂yl ∂yl
  ∂wlij , with ‖J l‖2 =
  ∑ ij(J l ij) 2. We will
  be concerned in our analysis with the quantities E(‖yl‖2), var(‖yl‖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 y0.
  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
  yl = φ(W l>yl−1), (8) where φ() denotes the ReLU non-linearity.
  For residual architectures, we study models with m layer residual branches of the form
  yl = yl−1 +Rm(yl−1) (9) where: Rm(yl−1) = W lm>...φ(W l2>φ(W l1>yl−1). (10) We denote the k’th intermediate output of the l’th residual branch by
  ỹl−1,k = W lk>...φ(W l2>φ(W l1>yl−1), (11) so that, for example
  Rm(yl−1) = W lm>ỹl−1,m−1 = W lm>φ(W lm−1>ỹ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 clm2 = 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 y0...yl−1, and is of dimensionsW l ∈ Rln×n. For our analysis, we breakW l into l−1 weight matrices W l1...W ll−1 ∈ Rn×n, so that the output of layer l is given by yl = φ( ∑l−1 l′=1W ll′>yl ′ ).
sections.heading: 4 FORWARD-BACKWARD NORM PROPAGATION DUALITY
sections.text: 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 ‖yL‖2, and those of the Jacobian norm ‖Jk‖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 yLt can be decomposed to paths that go through weight matrix W
  k, denoted by yL,kt (neglecting an additional superscript for ResNets and DenseNets), and paths that skip W
  k, denoted by ŷL,kt . Assuming the input ∀j , y0j = 1, we have:
  yLt = y L,k t + ŷ L,k t = ∑ k∈γt |γt|∏ l=1 wlγt,lzγt + ∑ k/∈γt |γt|∏ l=1 wlγt,lzγ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 nl 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 ‖Jkij‖2 is given by:
  ‖Jkij‖2 = ‖ ∂yL,k
  ∂wkij ‖2 = nL∑ t=1 ( ∑ ijk∈γt 1 wkij |γt|∏ l=1 wlγt,lzγt )2 (14)
  That is, the derivative of output unit t with respect to weight wkij is given by the sum of all paths γ that contain wkij divided by w k ij . In order to link J k with yL, we make the following definition:
  Definition 1. Reduced network: the outputs of a reduced network N(k) denoted by y1(k)...y L (k) given input y0 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, yl(k) = y l. The following Lemma links between the moments of yL(k) and those of y L:
  Lemma 1. For any architecture described in Sec. 3, it holds that for m ∈ {2, 4}: E ( ‖yL(k)‖ m ) = E ( ‖yL,k‖m ) (15)
  Note that if network N contains connections bypassing layer k, then for a specific realization of the weights, the equality yL(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 inN(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 ‖Jk‖2 with those of ‖yL(k)‖: Theorem 1. For any architecture described in Sec. 3,the following hold:
  E ( ‖Jk‖2 ) = 1 ck2 E ( ‖yL(k)‖ 2 ) , var(‖yL(k)‖ 2) ck4 ≤ var(‖Jk‖2) ≤ var(‖yL(k)‖ 2) (ck2) 2
  (16)
  where ck2 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 yL(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.
sections.heading: 4.1 VANILLA NETWORKS
sections.text: In order to compare our results for ResNets and DenseNets, we first derive the moments of yl(k) in our setting for vanilla architectures. Note that for vanilla networks, the reduced network matches the original, and so ∀0<l≤L, yl(k) = y l.
  Lemma 2. For a vanilla ReLU neural network, as described in Sec, 3, it holds for α = 1:
  C1 exp [ L−1∑ l=1 5− |∆| 2nl ] − C2 ≤ var ( ‖Jk‖2 ) ≤ C3 exp [ L−1∑ l=1 5 + |∆| nl ] − C4 (17)
  where C1, C2, C3, C4 are constants that do not depend on L.
  By setting α = 1, the variance of the Jacobian squared norm ‖Jk‖2 will increase exponentially with ∑L l=1 5 nl
  , as noted by Hanin & Rolnick (2018). This can be mitigated by linearly increasing the width of each layer with depth, so that nl = nl, resulting in constant asymptotic behavior
  exp [∑L−1 l=1 1 nl ] ∼ L 1n . 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.
sections.heading: 4.2 RESIDUAL NETWORKS
sections.text: 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(‖yl‖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 ‖Jkk′‖2 = ∑ ij ( ∂yL
  ∂wkk ′
  ij )2 for k′ < m. Unlike vanilla networks, the outputs of the reduced network yl(k) for l > k are not equal to those of the full network yl, and are given by:
  yl(k) = { Rm(yl−1) + yl−1 l 6= k Rm(yl−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:
  C1ρ m − exp
  [ 1
  2 Lρm−
  ] − C2µ2 ≤ var ( ‖Jkk ′ ‖2 ) ≤ C3ρm+ exp [ 2Lρ m 2 + ] − C4µ2 (19)
  where ρ± = α2(1 + 5±|∆| n ), C1, C2, C3, C4 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.
sections.heading: 4.3 DENSELY CONNECTED NETWORKS
sections.text: 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 2nl−1 initialization, gives rise to well-behaving gradients in any depth. By the definition of DenseNets given in Sec. 3, the reduced architecture outputs yl(k) are given by:
  yl(k) = { φ( ∑l−1 l′=1W ll′>yl ′ ) l ≤ k
  φ( ∑l−1 l′=kW ll′>yl ′ ) 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) ( C1 exp
  [ −2|∆|L
  n
  ] − C2 ) ≤ var(‖Jk‖2) ≤ L2(α−1) ( C3 exp [ 2|∆|L n ] − C4 ) (21)
  where C1, C2, C3, C4 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 nl
  . 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.
sections.heading: 5 THE CURIOUS CASE OF CONCATENATED RELU
sections.text: 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 wlij ∼ N (0, 1nl−1 ), then it holds that:
  E ( ‖yL‖2 ) = nL n , var ( ‖yL‖2 ) = n2L n2 ( L∏ l=1 (1 + 2 nL )− 1 ) (22)
  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 yl = W l1>φ(yl−1)−W l2>φ(−yl−1). An algebraically equivalent way to write the layer’s computation is yl = W l2>yl−1 + (W l1>+W l2>)φ(yl−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 byW l1 andW 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 yl = W l1>φ(yl−1)−W l2>φ(yl−1) conditioned on yl−1, where wl1ij , w l2 ij ∼ N (0, αnl−1 ) equals the distribution of y = W l>yl−1, where the elements of W l are distributed according to an i.i.d normal distribution wlij ∼ N (0, αnl−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 ( ∂y l
  ∂wlij )2 = ( ∂y
  l
  ∂wl1ij )2 + ( ∂y
  l
  ∂wl2ij )2. A surprising
  property emerges. Lem. 4 shows that similar to vanilla ReLU based networks, CR based ReLU networks have exponential behavior in ∑L l=1 2 nl . Since CR employ twice as many parameters, then
  for the same budget of parameters (having nl = nl√2 for CR) we have ∑L l=1 2 √ 2 nl < ∑L l=1 5 nl
  , 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 yl = W lk1>φ(ỹl,k−1) + W lk2>φ(−ỹl,k−1) and yl = ∑l−1 l′=0W ll′>yl ′
  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.
sections.heading: 6 CONCLUSIONS
sections.text: 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.
sections.heading: 7 PROOFS
sections.text: 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
  Qγt = |γt|∏ l=1 wγt,lzγt = |γt|∏ l=1 qγt,l (23)
  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,lzγt .
  We generalize the above expression the following way:
  Qγt = L∏ l=1 qγt,l (24)
  where:
  qγt,l = { 1 l /∈ γt w̃γt,lzγt,l else
  (25)
  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,l1wγt,l2 for residual layers (wγt,l1 ∈W l1 and wγt,l2 ∈W l2). proposition 1. Given a random vector w = [w1...wn] such that each component is identically and symmetrically distributed iid random variable with moments E(wmi ) = cm, and random binary variable z such that z|w ∈ {0, 1}, z|w = 1− z| − w, then it holds that for n > l:
  E( l∏ i=1 wmi z) =  0 m is odd, l is even (cm) l 2 m is even
  E( ∏l i=1 w m i z) else
  (26)
  Proof. We have:
  clm = ∫ w l∏ i=1 wmi p(w)dw = ∫ w|z=1 l∏ i=1 wmi p(w)dw + ∫ w|z=0 l∏ i=1 wmi p(w)dw (27)
  = ∫ w|z=1 l∏ i=1 wmi p(w)dw + ∫ w|z=1 l∏ i=1 (−wi)mp(w)dw (28)
  = ∫ w l∏ i=1 wmi zp(w)dw + ∫ w l∏ i=1 (−wi)mzp(w)dw (29)
  For even m or even l: it follows that:
  clm = 2 ∫ w l∏ i=1 wmi zp(w)dw (30)
  and so:
  E( l∏ i=1 wmi z) = clm 2
  (31)
  For odd l, since w is symmetrically distributed, we have that clm = 0, and so E( ∏l i=1 w m i z) = clm 2 = 0
  proposition 2. For any set of paths γ1t , γ2t , it holds that:
  E(Qγ1tQγ2t ) = { 0 γ1t 6= γ2t E(Q2
  γ1t ) =
  ∏L l=1 E ( (qγ1t ,l) 2 ) γ1t = γ 2 t
  (32)
  Proof. From proposition 2 we have:
  E(Qγ1tQγ2t ) = E( L∏ l=1 qγ1t ,lqγ2t ,l) (33)
  Since the computations done by all considered architectures form a markov chain, such that the output of layer yl depends only on weights Rl = {W 1...W l}, using conditional expectation,:
  E(Qγ1tQγ2t ) = E ( E(Qγ1tQγ2t ∣∣∣RL−1)) = E( L−1∏ l=1 qγ1t ,lqγ2t ,lE(qγ1t ,Lqγ2t ,L ∣∣∣RL−1)) (34) Note that since both are paths leading to output unit t, the following holds:
  E(qγLt ,lqγ2t ,L ∣∣∣RL−1) =  E(w̃γ1t ,Lw̃γ2t ,Lz L t ) L ∈ γ1t , L ∈ γ2t E(w̃γ1t ,Lz L t ) L ∈ γ1t , L /∈ γ2t E(w̃γ2t ,Lz L t ) L /∈ γ1t , L ∈ γ2t
  1 L /∈ γ1t , L /∈ γ2t
  (35)
  and so it follows from proposition 1:
  E(qγLt ,lqγ2t ,L ∣∣∣RL−1) =  E(w̃2 γ1t ,L zLt ) L ∈ γ1t , L ∈ γ2t , w̃γ1t ,L = w̃γ2t ,L
  1 L /∈ γ1t , L /∈ γ2t 0 else
  = E ( (qγ1t ,L) 2 ) (36)
  Recursively going through l = l − 1...1 completes the proof.
  proposition 3. For any set of paths γ1t , γ2t , γ1s , γ2s , from input to output units t and s, it holds that:
  E(Qγ1tQγ2tQγ1sQγ2s ) = L∏ l=1 E ( ξγ12st ,l ) (37)
  such that:
  E ( ξγ12st ,l ) =  E ( q2 γ1t ,l q2γ1s ,l ) qγ1t ,l = qγ2t ,l, qγ1s ,l = qγ2s ,l E ( q2 γ1t ,l q2 γ2t ,l ) qγ1t ,l = qγ1s ,l, qγ2t ,l = qγ2s ,l
  0 else
  (38)
  Proof. The proof is practically identical to proposition 3, and will be omitted.
  For notation simplicity,
sections.heading: 7.1 PROOF OF LEMMA 1
sections.text: 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 yl depend only weights W 1...W l, then for any particular realization of the weights, yL,k and yLk differ only in their activation patterns in layers l ≥ k due to the removal of bypassing connections in yL(k). Using Eq. 23:
  E ( ‖yL,k‖m ) = nL∑ t=1 E ( ( ∑
  γt|k∈γt
  |γt|∏ l=1 wγt,lzγt) m )
  (39)
  while:
  E ( ‖yL(k)‖ m ) = nL∑ t=1 E ( ( ∑
  γt|k∈γt
  |γt|∏ l=1 wγt,lz ′ γ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 ( ‖yL,k‖2 ) = nL∑ t=1 ∑ γ1t |k∈γt E(Q2γ1t ) = nL∑ t=1 ∑ γ1t |k∈γt L∏ l=1 E ( (qγ1t ,l) 2 )
  (41)
  Similarly for m = 4, from proposition 3 the expectation factorizes:
  E ( ‖yL,k‖4 ) = nL∑ t=1 ∑ k∈γ1t ,γ2t ,γ1s ,γ2s L∏ l=1 E ( ξγ12st ,l ) (42)
  where ξγ12st ,l is defined is Eq. 38. Note that switching the sign of W l for all architectures considered will switch the state of both zlj , 0 < j ≤ nl and z ′l j , 0 < j ≤ nl, and so from proposition 1, the factorized expectations for m = 2 and m = 4 meet the required invariance.
sections.heading: 7.2 PROOF OF THEOREM 1
sections.text: From the definition of the Jacobian for weights W k, we have using lemma 1 and proposition 2:
  E(‖Jkij‖2) = E(‖ ∂yL(k)
  ∂wkij ‖2) = nL∑ t=1 ∑ ijk∈γ1t ,γ2t E ( 1 (wkij) 2 L∏ l=1 Qγ1tQγ2t )
  = nL∑ t=1 ∑ ijk∈γt E ( 1 (wkij) 2 L∏ l=1 Q2γt ) (43)
  Since by definition wkij ∈ γ1t , from proposition 2 the term E (
  1 (wkij) 2
  ∏L l=1Q
  2 γ1t
  ) factorizes:
  E ( 1
  (wkij) 2 L∏ l=1 Q2γt ) = E( q2γt,k (wkij) 2 ) L∏ l 6=k,l=1 E(q2γt,l) = E( (wkij) 2zki (wkij) 2 ) L∏ l 6=k,l=1 E(q2γt,l) (44)
  = E(zki ) L∏
  l 6=k,l=1
  E(q2γt,l) = 1
  2 L∏ l 6=k,l=1 E(q2γt,l) = E 1 ck2 ( L∏ l=1 Q2γt )
  (45)
  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:
  E(‖Jk‖4) = E( ∑ ijuv ‖ ∂yL(k) ∂wkij ‖2‖ ∂yL(k) ∂wkuv ‖2)
  = nL∑ t=1 nL∑ s=1 ∑ ijk∈γ1t ,γ2t ∑ uvk∈γ1s ,γ2s E ( 1 (wkij) 2(wkuv) 2 L∏ l=1 Qγ1tQγ2tQγ1sQγ2s ) (46)
  From proposition 3, for uvk ∈ γ1s , γ2s , ijk ∈ γ1t , γ2t the expectation factorizes:
  E ( 1
  (wkij) 2(wkuv) 2 L∏ l=1 Qγ1tQγ2tQγ1sQγ2s ) = E ( ξγ12st ,k (wkij) 2(wkuv) 2 ) L∏ l 6=k,l=1 E ( ξγ12st ,l ) (47)
  (48)
  where from proposition 1 and proposition 3, for uvk ∈ γ1s , γ2s , ijk ∈ γ1t , γ2t :
  E ( ξγ12st ,k
  (wkij) 2(wkuv) 2
  ) =  E(ξ γ12st ,k ) (ck2 ) 2 {u, v} 6= {i, j} E(ξ γ12st ,k )
  ck4 {u, v} = {i, j}
  (49)
  Note that since ck4 ≥ (ck2)2, it holds that:
  E(ξγ12st ,k) ck4
  ≤ E ( ξγ12st ,k
  (wkij) 2wkuv)
  2 ) ) ≤
  E(ξγ12st ,k) (ck2) 2 (50)
  Plugging into Eq. 46, it follows:
  var(‖yL(k)‖ 2)
  ck4 ≤ var(‖Jk‖2) = E(‖Jk‖4)− E(‖Jk‖2)2 ≤
  var(‖yL(k)‖ 2)
  (ck2) 2
  (51)
  Furthermore, it holds that:
  E(‖Jk‖4) = ∑ st ∑ ijk∈γ1t ,γ2t ∑ uvk∈γ1s ,γ2s E ( ξγ12st ,k (wkij) 2(wkuv) 2 ) L∏ l 6=k,l=1 E ( ξγ12st ,l ) (52)
  = ∑ st ∑ ijk∈γ1t ,γ2t ∑ uvk∈γ1s ,γ2s E ( ξγ12st ,k (wkij) 2(wkuv) 2 )( 1{u,v}={i,j} + 1{u,v}6={i,j} ) L∏ l 6=k,l=1 E ( ξγ12st ,l ) (53)
  Therefor for large width networks:
  E(‖Jk‖4) ∼ ∑ st ∑ ijk∈γ1t ,γ2t ∑ uvk∈γ1s ,γ2s E ( ξγ12st ,k (wkij) 2(wkuv) 2 ) 1{u,v}6={i,j} L∏ l 6=k,l=1 E ( ξγ12st ,l ) (54)
  = 1
  (ck2) 2 ∑ st ∑ ijk∈γ1t ,γ2t ∑ uvk∈γ1s ,γ2s E ( ξγ12st ,k ) 1{u,v}6={i,j} L∏ l 6=k,l=1 E ( ξγ12st ,l ) (55)
  ∼ 1 (ck2) 2 E(‖yL(k)‖ 4) (56)
  and so for large width networks:
  var(‖Jk‖2) ∼ var(‖yL(k)‖ 2)
  (ck2) 2
  (57)
sections.heading: 7.3 PROOF OF LEMMA 2
sections.text: Denoting by Zl the diagonal matrix holding in its diagonal the activation variables zlj for unit j in layer l, and so we have:
  yl = ZlW l>yl−1 (58)
  Conditioning on Rl−1 = {W 1...WL−1} and taking expectation: E ( ‖yL‖2|RL−1 ) = yL−1> E ( WLZLWL> ) yL−1 (59)
  = nL∑ j=1 nL−1∑ i1,i2=1 yl−1i1 y l−1 i2 E ( wLi1,jw L i2,jz L j |RL−1 ) (60)
  From proposition 1, it follows that: E ( ‖yL‖2|RL−1 ) =
  nL nL−1
  E ( ‖yL−1‖2 ) (61)
  Recursing through l = L− 1...1: E ( ‖yL‖2 ) = nL n
  (62)
  Similarly:
  E ( ‖yL‖4|RL−1 ) = E ( (yL−1>WLZLWL>yL−1)2 ) (63)
  = ∑
  j1,j2,i1,i2,i3,i4
  yL−1i1 y L−1 i2 yL−1i3 y L−1 i4 E ( wLi1,j1w L i2,j1w L i3,j2w L i4,j2z L j1z L j2 |R L−1 )
  (64)
  From proposition 3:
  E ( wLi1,j1w L i2,j1w L i3,j2w L i4,j2z L j1z L j2 |R L−1 ) (65)
  = E ( wLi1,j1w L i2,j1w L i3,j2w L i4,j2z L j1z L j2 |R L−1 )( 1j1=j2,i1=i2=i3=i4 + 31j1=j2,i1=i2,i3=i4,i1 6=i3 ... (66)
  +1j1 6=j2,i1=i2,i3=i4
  ) (67)
  and so: E ( ‖yL‖4 ) = nLc L 4
  2 ∑ i=1 E ( (yL−1i ) 4 ) + 6nL n2L−1 ∑ i1 6=i2 E ( (yL−1i1 ) 2(yL−1i2 ) 2 )
  (68)
  + nL(nL − 1)
  n2L−1
  ∑ i1,i2 E ( (yL−1i1 ) 2(yL−1i2 ) 2 )
  (69)
  = nL(c
  L 4 − 3(cLc )2)
  2
  ∑ i E ( (yL−1i ) 4 ) + nL(nL + 5) n2L−1 E ( ‖yL−1‖4 ) (70)
  = nL∆
  n2L−1 ∑ i E ( (yL−1i ) 4 ) + nL(nL + 5) n2L−1 E ( ‖yL−1‖4 ) (71)
  finally:
  nL(nL + 5 + ∆) n2L−1 E ( ‖yL−1‖4 ) ≤ E ( ‖yL‖4 ) (72)
  ≤ nL(nL + 5 + |∆|) n2L−1
  E ( ‖yL−1‖4 ) (73)
  Recursing through L − 1...1, and denoting C1 = nL(nL+5−|∆|)n2 , and C2 = nL(nL+5+|∆|)
  n2 and combining with Eq. 62:
  C1 L−1∏ l=1 (1 + 5− |∆| nl )− n 2 L n2 ≤ var ( ‖yL‖2 ) ≤ C2 L−1∏ l=1 (1 + 5 + |∆| nl )− n 2 L n2 (74)
  Using ∀1≥x≥−1, e0.5x ≤ (1 + x) ≤ ex:
  C1 exp [ L∑ l=1 5− |∆| 2nl ] − C2 ≤ var ( ‖yL‖2 ) ≤ C3 exp [ L∑ i=1 5 + |∆| nl ] − C4 (75)
  For Gaussian distribution we have that ∆ = 0, and so:
  var ( ‖yL‖2 ) = nL(nL + 5)
  n2
  L−1∏ l=1 (1 + 5 nl )− n 2 L n2 (76)
sections.heading: 7.4 PROOF OF THEOREM 2
sections.text: In the following proof, for the sake of notation simplicity, we omit the notation k in yl(k), and assume that yl stands for the reduced network yl(k). We have for layer L:
  E ( ‖yL‖2 ) = E ( ‖yL−1‖2 ) + E ( ỹL−1,m−1>WLmWLm>ỹL−1,m−1 ) = E ( ‖yL−1‖2 ) + αE ( ‖ỹL−1,m−1‖2 ) = E ( ‖yL−1‖2 )( 1 + αm
  ) = E ( ‖yk‖2
  ) L∏ l=k+1 ( 1 + αm ) = E ( ‖yk−1‖2 ) L∏ l=k+1 ( 1 + αm ) αm
  = ∏ l 6=k ( 1 + αm ) αm (77)
  E ( ‖yL‖4 ) = E ( ‖yL−1‖4 ) + E ( ‖ỹL−1,m‖4 ) + 4E ( (ỹL−1,m>yL−1)2 ) + 2E ( ‖ỹL−1,m‖2‖yL−1‖2 ) (78)
  We now handle each term separately: E ( ‖ỹL−1,m‖4 ) = E ( E ( ‖ỹL−1,m‖4|RL−1 )) (79)
  Similarly to the vanilla case with a linear weight layer on top, we have:
  α2m(1 + 5− |∆| n
  )m−1‖yL−1‖4 ≤ E ( ‖ỹL−1,m‖4|RL−1 ) (80)
  ≤ α2m(1 + 5 + |∆| n )m‖yL−1‖4 (81)
  E ( (ỹL−1,m−1>yL−1)2 ) = ∑
  j1j2i1i2
  E ( ỹL−1,m−1i1 ỹ L−1,m−1 i2 yL−1j1 y L−1 j2 wLmi1,j1w Lm i2,j2 ) (82)
  = α n E ( ‖ỹL−1,m−1‖2‖yL−1‖2 ) = αm n E ( ‖yL−1‖4 ) (83)
  and: E ( ‖ỹL−1,m‖2‖yL−1‖2 ) = αmE ( ‖yL−1‖4 ) (84)
  Plugging it all into Eq. 78, and denoting:
  β+l = 1 + 2α m(1 +
  2 n ) + α2m(1 + 5 + |∆| n )m (85)
  β−l = 1 + 2α m(1 +
  2 n ) + α2m(1 + 5− |∆| n )m−1 (86)
  we have after recursing through l = L− 1...1:
  E ( ‖yk‖4 ) L∏ l=k+1 β−l ≤ E ( ‖yL‖4 ) ≤ E ( ‖yk‖4 ) L∏ l=k+1 β+l (87)
  Since for the reduced architecture yL(k), the transformation between y k−1 and yk is given by an m layer fully connected network without any residual connections, we can use the results from the vanilla case:
  L∏ l 6=k β−l α 2m(1 + 5− |∆| n )m−1 ≤ E ( ‖yL‖4 ) ≤ L∏ l 6=k β+l α 2m(1 + 5 + |∆| n )m (88)
  Denoting ρ± = α2(1 + 5±|∆| n ) it follows:
  C1 L∏ l 6=k β−l ρ m − ≤ E ( ‖yL‖4 ) ≤ L∏ l 6=k β+ρm+ (89)
  From the definition of β±l , it holds for ρ+ < 1:
  β+l ≤ ( 1 + ρ m 2 + )2 ≤ exp [ 2ρ m 2 + ] (90)
  and:
  β−l ≥ ( 1 + ρm− ) ≥ exp [ 1
  2 ρm−
  ] (91)
  and so:
  C1ρ m − exp
  [ 1
  2 Lρm−
  ] − C2µ2 ≤ var ( ‖yL‖2 ) ≤ C3ρm+ exp [ 2Lρ m 2 + ] − C3µ2 (92)
  For Gaussian distributions, using a similar derivation and assigning ∆ = 0 we have:
  E ( ‖yL‖4 ) = ρ L∏ l 6=k βl (93)
  where:
  βl = 1 + 2α m(1 +
  2 n ) + α2m(1 + 5 n )m−1(1 + 2 n ) (94)
  and:
  ρ = a2m(1 + 5
  n )m−1(1 +
  2 n ) (95)
sections.heading: 7.5 PROOF OF THEOREM 3
sections.text: In the following proof, for the sake of notation simplicity, we omit the notation k in yl(k), and assume that yl stands for the reduced network yl(k). We have:
  E ( ‖yL‖2 ) = µL = E ( ( L−1∑ l=k yl>WLl)ZL( L−1∑ l=k yl>WLl) ) = aL L L−1∑ l=k µl−1 (96)
  E ( ‖yL‖4 ) = E (( L−1∑ l=k (yl>WLl)ZL L−1∑ l=k (yl>WLl) )2)
  = E (( L−1∑
  l1=k (yl1>WLl1)ZL L−1∑ l2=k (yl2>WLl2) L−1∑ l3=k (yl3>WLl3)ZL L−1∑ l4=k (yl4>WLl4) ) (97)
  Note that the above is a simple Relu block with concatenated y1...yL − 1 as input. Using the results from the vanilla architecture, and denoting Cl,l′ = E ( ‖yl‖2‖yl′‖2 ) , it then follows:
  (n+ 5− |∆|)α2
  nL2
  L−1∑ l1,l2=k Cl1,l2 ≤ CL,L ≤ (n+ 5 + |∆|)α2 nL2 L−1∑ l1,l2=k Cl1,l2 (98)
  From Eq. 98, it holds that:
  n(L− 1)2
  n+ 5 + |∆| CL−1,L−1 ≤ L−2∑ l1,l2=k Cl1,l2 ≤ n(L− 1)2 n+ 5− |∆| CL−1,L−1 (99)
  It then follows:
  E(‖yL‖4) = CL,L ≤ α2
  L2 (1 + 5 + |∆| n ) L−1∑ l1,l2=k Cl1l2
  = α2
  L2 (1 + 5 + |∆| n
  ) ( CL−1,L−1 + L−2∑ l1,l2=k Cl1l2 + 2 L−2∑ l=k CL−1,l )
  = α2
  L2 (1 + 5 + |∆| n
  ) ( CL−1,L−1 +
  (L− 1)2n a2L−1(n+ 5− |∆|) CL−1,L−1 + 2 L−2∑ l=k CL−1,l ) (100)
  (101) The following also holds:
  ∀l1>l2≥k, Cl1,l2 = E ( ( l1−1∑ l=k yl>W l1,lZl1)2‖yl2‖2 ) = α l1 l1−1∑ l=k Cl,l2 (102)
  and so:
  CL,L ≤ (n+ 5 + |∆|)α2
  nL2
  ( CL−1,L−1 +
  (L− 1)2n α2(n+ 5− |∆|) CL−1,L−1... (103)
  + 2α
  L− 1 L−2∑ l1=k L−2∑ l2=k Cl1,l2 ) (104)
  = (n+ 5 + |∆|)α2
  nL2
  ( CL−1,L−1 +
  (L− 1)2n α2(n+ 5− |∆|) CL−1,L−1 (105)
  + 2n(L− 1)
  α(n+ 5− |∆|) CL−1,L−1
  ) (106)
  = (n+ 5 + |∆|)α2
  nL2 CL−1,L−1
  ( 1 +
  (L− 1)2n α2(n+ 5− |∆|) + 2n(L− 1) α(n+ 5− |∆|)
  ) (107)
  = CL−1,L−1 (n+ 5 + |∆|) (n+ 5− |∆|) ((α L + α(L− 1) Lα )2 + α2(5− |∆|) nL2 ) (108)
  CL,L ≤ CL−1,L−1 ( 1 + 2|∆|
  n+ 5− |∆|
  )(( 1 +
  α− 1 L
  )2 + α2(5− |∆|)
  nL2
  ) (109)
  = CL−1,L−1
  ( 1 +
  2|∆| n+ 5− |∆|
  )( 1 +
  α2(5− |∆|) n(L+ α− 1)2
  )( 1 +
  α− 1 L
  )2 (110)
  Telescoping through l = L− 1...1:
  CL,L ≤ ( 1 + 2|∆|
  n+ 5− |∆| )L−k+1 L∏ l 6=k (( 1 + α2(5− |∆|) n(l + α− 1)2 )( 1 + α− 1 l )2) (111)
  ≤ exp [ 2|∆|(L− 1) n+ 5− |∆| ] exp [∑ l 6=k α2(5− |∆|) n(l + α− 1)2) ] exp [∑ l 6=k 2(α− 1) l ] (112)
  ∼ exp [ 2|∆|(L− 1) n+ 5− |∆| ] exp [ α2(5− |∆|) nL ] L2(α−1) (113)
  ∼ C1 exp [ 2|∆|L n ] L2(α−1) (114)
  A similar relation holds for the mean:
  µL = α
  L L−1∑ l=k µl = α L ( µL−1,L−1 + (L− 1) α )µL−1,L−1 ) (115)
  µ2L = µ 2 L−1,L−1
  ( 1 +
  α− 1 L
  )2 ≤ L2(α−1) (116)
  Finally:
  var(‖yL‖2) = CL,L − µ2L (117)
  ≤ L2(α−1) ( C1 exp [ 2|∆|L n ] − 1 )
  (118)
  A lower bound can easily be derived using Eq. 98, 99, and follow the same derivation by flipping the sign of |∆|:
  L(α−1) ( C1 exp
  [ −2|∆|L
  n
  ] − C2 ) ≤ var(‖yL‖2) ≤ L2(α−1) ( C3 exp [ 2|∆|L n ] − C4 ) (119)
  For Gaussian distributions, using a similar derivation and assigning ∆ = 0 we have:
  CL,L = L∏ l 6=k (( 1 +
  5α2 n(l + α− 1)2 )( 1 + α− 1 l )2)(n+ 5 nk2 ) a2 (120)
  µ2L = L∏ l 6=k ( 1 + α− 1 l )2 a2 k2
  (121)
  and so: var(‖yL‖2) = ( a2 (n+ 5 nk2 ) L∏ l 6=k ( 1 +
  5α2 n(l + α− 1)2 ) − a 2 k2 ) L∏ l 6=k ( 1 + α− 1 l )2 (122)
sections.heading: 7.6 PROOF OF LEMMA 3
sections.text: The proof is identical to the Relu case, removing the activation variables Zl, and hence will be omitted.
sections.heading: 7.7 PROOF OF LEMMA 4
sections.text: Proof. We construct the matrix W l as follows:
  wlij =
  { wl1ij y
  l−1 ≥ 0 wl2ij y l−1 < 0 (123)
  And so yl = W l>yl−1. Note that the elements of W l are i.i.d Gaussian random variables wlij ∼ N (0, αnl−1 ).
sections.heading: 8 EMPIRICAL SUPPORT FOR CONCATENATIVE RELU
sections.text: 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.,
  each layer had 512 hidden neurons. This was done since CR employs more weights than the baselines due to the dual pathway 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 τ .
sections.heading: 9 UCI RESULTS
sections.text: In Tab. 1 we provide the full numbers of the UCI experiments in the paper.