Resurrecting the sigmoid in deep learning through dynamical isometry: theory and practice

Introduction

Deep learning has achieved state-of-the-art performance in many domains, including computer vision , machine translation , human games , education , and neurobiological modeling . A major determinant of success in training deep networks lies in appropriately choosing the initial weights. Indeed the very genesis of deep learning rested upon the initial observation that unsupervised pre-training provides a good set of initial weights for subsequent fine-tuning through backpropagation . Moreover, seminal work in deep learning suggested that appropriately-scaled Gaussian weights can prevent gradients from exploding or vanishing exponentially , a condition that has been found to be necessary to achieve reasonable learning speeds .

These random weight initializations were primarily driven by the principle that the mean squared singular value of a deep network’s Jacobian from input to output should remain close to $1$. This condition implies that on average, a randomly chosen error vector will preserve its norm under backpropagation; however, it provides no guarantees on the worst case growth or shrinkage of an error vector. A stronger requirement one might demand is that every Jacobian singular value remain close to $1$. Under this stronger requirement, every single error vector will approximately preserve its norm, and moreover all angles between different error vectors will be preserved. Since error information backpropagates faithfully and isometrically through the network, this stronger requirement is called dynamical isometry .

A theoretical analysis of exact solutions to the nonlinear dynamics of learning in deep linear networks revealed that weight initializations satisfying dynamical isometry yield a dramatic increase in learning speed compared to initializations that do not. For such linear networks, orthogonal weight initializations achieve dynamical isometry, and, remarkably, their learning time, measured in number of learning epochs, becomes independent of depth. In contrast, random Gaussian initializations do not achieve dynamical isometry, nor do they achieve depth-independent training times.

It remains unclear, however, how these results carry over to deep nonlinear networks. Indeed, empirically, a simple change from Gaussian to orthogonal initializations in nonlinear networks has yielded mixed results , raising important theoretical and practical questions. First, how does the entire distribution of singular values of a deep network’s input-output Jacobian depend upon the depth, the statistics of random initial weights, and the shape of the nonlinearity? Second, what combinations of these ingredients can achieve dynamical isometry? And third, among the nonlinear networks that have neither vanishing nor exploding gradients, do those that in addition achieve dynamical isometry also achieve much faster learning compared to those that do not? Here we answer these three questions, and we provide a detailed summary of our results in the discussion.

Theoretical Results

In this section we derive expressions for the entire singular value density of the input-output Jacobian for a variety of nonlinear networks in the large-width limit. We compute the mean squared singular value of $\mathbf{J}$ (or, equivalently, the mean eigenvalue of $\mathbf{J}\mathbf{J}^{T}$), and deduce a rescaling that sets it equal to $1$. We then examine two metrics that help quantify the conditioning of the Jacobian: $s_{\text{max}}$, the maximum singular value of $\mathbf{J}$ (or, equivalently, $\lambda_{\text{max}}$, the maximum eigenvalue of $\mathbf{J}\mathbf{J}^{T}$); and $\sigma^{2}_{JJ^{T}}$, the variance of the eigenvalue distribution of $\mathbf{J}\mathbf{J}^{T}$. If $\lambda_{\text{max}}\gg 1$ and $\sigma^{2}_{JJ^{T}}\gg 1$ then the Jacobian is ill-conditioned and we expect the learning dynamics to be slow.

Here $\mathbf{D}^{l}$ is a diagonal matrix with entries ${D}^{l}_{ij}=\phi^{\prime}({h}^{l}_{i})\,\delta_{ij}$. The input-output Jacobian $\mathbf{J}$ is closely related to the backpropagation operator mapping output errors to weight matrices at a given layer; if the former is well conditioned, then the latter tends to be well-conditioned for all weight layers. We therefore wish to understand the entire singular value spectrum of $\mathbf{J}$ for deep networks with randomly initialized weights and biases.

In particular, we will take the biases ${b}^{l}_{i}$ to be drawn i.i.d. from a zero mean Gaussian with standard deviation $\sigma_{b}$. For the weights, we will consider two random matrix ensembles: (1) random Gaussian weights in which each ${W}^{l}_{ij}$ is drawn i.i.d from a Gaussian with variance $\sigma_{w}^{2}/N$, and (2) random orthogonal weights, drawn from a uniform distribution over scaled orthogonal matrices obeying $(\mathbf{W}^{l})^{T}\mathbf{W}^{l}=\sigma_{w}^{2}\,\mathbf{I}$.

2 Review of signal propagation

The random matrices $\mathbf{D}^{l}$ in eqn. (2) depend on the empirical distribution of pre-activations $\mathbf{h}^{l}$ entering the nonlinearity $\phi$ in eqn. (1). The propagation of this empirical distribution through different layers $l$ was studied in . There, it was shown that in the large-$N$ limit this empirical distribution converges to a Gaussian with zero mean and variance $q^{l}$, where $q^{l}$ obeys a recursion relation induced by the dynamics in eqn. (1),

with initial condition $q^{0}=\frac{1}{N}\sum_{i=1}^{N}({h}^{0}_{i})^{2}$, and where $\mathcal{D}h=\frac{dh}{\sqrt{2\pi}}\,\exp{(-\frac{h^{2}}{2})}$ denotes the standard Gaussian measure. This recursion has a fixed point obeying,

If the input $\mathbf{h}^{0}$ is chosen so that $q^{0}=q^{*}$, then we start at the fixed point, and the distribution of $\mathbf{D}^{l}$ becomes independent of $l$. Also, if we do not start at the fixed point, in many scenarios we rapidly approach it in a few layers (see ), so for large $L$, assuming $q^{l}=q^{*}$ at all depths $l$ is a good approximation in computing the spectrum of $\mathbf{J}$.

Another important quantity governing signal propagation through deep networks is

where $\phi^{\prime}$ is the derivative of $\phi$. Here $\chi$ is the mean of the distribution of squared singular values of the matrix $\mathbf{DW}$, when the pre-activations are at their fixed point distribution with variance $q^{*}$. As shown in and Fig. 1, $\chi(\sigma_{w},\sigma_{b})$ separates the $(\sigma_{w},\sigma_{b})$ plane into two phases, chaotic and ordered, in which gradients exponentially explode or vanish respectively. Indeed, the mean squared singular value of $\mathbf{J}$ was shown simply to be $\chi^{L}$ in , so $\chi=1$ is a critical line of initializations with neither vanishing nor exploding gradients.

3 Free probability, random matrix theory and deep networks.

While the previous section revealed that the mean squared singular value of $\mathbf{J}$ is $\chi^{L}$, we would like to obtain more detailed information about the entire singular value distribution of $\mathbf{J}$, especially when $\chi=1$. Since eqn. (2) consists of a product of random matrices, free probability becomes relevant to deep learning as a powerful tool to compute the spectrum of $\mathbf{J}$, as we now review.

In general, given a random matrix $\mathbf{X}$, its limiting spectral density is defined as

where $\langle\cdot\rangle_{X}$ denotes the mean with respect to the distribution of the random matrix $\mathbf{X}$. Also,

is the definition of the Stieltjes transform of $\rho_{X}$, which can be inverted using,

The Stieltjes transform $G_{X}$ is related to the moment generating function $M_{X}$,

where the $m_{k}$ is the $k$th moment of the distribution $\rho_{X}$, $m_{k}=\int d\lambda\;\rho_{X}(\lambda)\lambda^{k}=\frac{1}{N}\langle\text{tr}\mathbf{X}^{k}\rangle_{X}\,.$ In turn, we denote the functional inverse of $M_{X}$ by $M_{X}^{-1}$, which by definition satisfies $M_{X}(M_{X}^{-1}(z))=M_{X}^{-1}(M_{X}(z))=z$. Finally, the S-transform is defined as,

The utility of the S-transform arises from its behavior under multiplication. Specifically, if $\mathbf{A}$ and $\mathbf{B}$ are two freely-independent random matrices, then the S-transform of the product random matrix ensemble $\mathbf{A}\mathbf{B}$ is simply the product of their S-transforms,

Our first main result will be to use eqn. (11) to write down an implicit definition of the spectral density of $\mathbf{J}\mathbf{J}^{T}$. To do this we first note that (see Result 1 of the supplementary material),

where we have used the identical distribution of the weights to define $S_{WW^{T}}=S_{W_{l}W_{l}^{T}}$ for all $l$, and we have also used the fact the pre-activations are distributed independently of depth as $h_{l}\sim\mathcal{N}(0,q^{*})$, which implies that $S_{D^{2}_{l}}=S_{D^{2}}$ for all $l$.

Eqn. (12) provides a method to compute the spectrum $\rho_{JJ^{T}}(\lambda)$. Starting from $\rho_{W^{T}W}(\lambda)$ and $\rho_{D^{2}}(\lambda)$, we compute their respective S-transforms through the sequence of equations eqns. (7), (9), and (10), take the product in eqn. (12), and then reverse the sequence of steps to go from $S_{JJ^{T}}$ to $\rho_{JJ^{T}}(\lambda)$ through the inverses of eqns. (10), (9), and (8). Thus we must calculate the S-transforms of $\mathbf{W}\mathbf{W}^{T}$ and $\mathbf{D}^{2}$, which we attack next for specific nonlinearities and weight ensembles in the following sections. In principle, this procedure can be carried out numerically for an arbitrary choice of nonlinearity, but we postpone this investigation to future work.

4 Linear networks

As a warm-up, we first consider a linear network in which $\mathbf{J}=\prod_{l=1}^{L}\mathbf{W}^{l}$. Since criticality ($\chi=1$ in eqn. (5)) implies $\sigma_{w}^{2}=1$ and eqn. (4) reduces to $q^{*}=\sigma_{w}^{2}q^{*}+\sigma_{b}^{2}$, the only critical point is $(\sigma_{w},\sigma_{b})=(1,0)$. The case of orthogonal weights is simple: $\mathbf{J}$ is also orthogonal, and all its singular values are $1$, thereby achieving perfect dynamic isometry. Gaussian weights behave very differently. The squared singular values $s^{2}_{i}$ of $\bm{J}$ equal the eigenvalues $\lambda_{i}$ of $\mathbf{J}\mathbf{J}^{T}$, which is a product Wishart matrix, whose spectral density was recently computed in . The resulting singular value density of $\mathbf{J}$ is given by,

Fig. 2(a) demonstrates a match between this theoretical density and the empirical density obtained from numerical simulations of random linear networks. As the depth increases, this density becomes highly anisotropic, both concentrating about zero and developing an extended tail.

Note that $\phi=\pi/(L+1)$ corresponds to the minimum singular value $s_{\text{min}}=0$, while $\phi=0$ corresponds to the maximum eigenvalue, $\lambda_{\text{max}}=s_{\text{max}}^{2}=L^{-L}(L+1)^{L+1}$, which, for large $L$ scales as $\lambda_{\text{max}}\sim eL$. Both eqn. (13) and the methods of Section 2.5 yield the variance of the eigenvalue distribution of $\mathbf{J}\mathbf{J}^{T}$ to be $\sigma^{2}_{JJ^{T}}=L$. Thus for linear Gaussian networks, both $s_{\text{max}}$ and $\sigma^{2}_{JJ^{T}}$ grow linearly with depth, signaling poor conditioning and the breakdown of dynamical isometry.

5 ReLU and hard-tanh networks

We first discuss the criticality conditions (finite $q^{*}$ in eqn. (4) and $\chi=1$ in eqn. (5)) in these two nonlinear networks. For both networks, since the slope of the nonlinearity $\phi^{\prime}(h)$ only takes the values and $1$, $\chi$ in eqn. (5) reduces to $\chi=\sigma_{w}^{2}p(q^{*})$ where $p(q^{*})$ is the probability that a given neuron is in the linear regime with $\phi^{\prime}(h)=1$. As discussed above, we take the large-width limit in which the distribution of the pre-activations $h$ is a zero mean Gaussian with variance $q^{*}$. We therefore find that for ReLU, $p(q^{*})=\frac{1}{2}$ is independent of $q^{*}$, whereas for hard-tanh, $p(q^{*})=\int_{-1}^{1}dh\frac{e^{-h^{2}/2q^{*}}}{\sqrt{2\pi q^{*}}}=\operatorname{erf}(1/\sqrt{2q^{*}})$ depends on $q^{*}$. In particular, it approaches $1$ as $q^{*}\to 0$.

Thus for ReLU, $\chi=1$ if and only if $\sigma_{w}^{2}=2$, in which case eqn. (4) reduces to $q^{*}=\frac{1}{2}\sigma_{w}^{2}q^{*}+\sigma_{b}^{2}$, implying that the only critical point is $(\sigma_{w},\sigma_{b})=(2,0)$. For hard-tanh, in contrast, $\chi=\sigma_{w}^{2}p(q^{*})$, where $p(q^{*})$ itself depends on $\sigma_{w}$ and $\sigma_{b}$ through eqn. (4), and so the criticality condition $\chi=1$ yields a curve in the $(\sigma_{w},\sigma_{b})$ plane similar to that shown for the tanh network in Fig. 1. As one moves along this curve in the direction of decreasing $\sigma_{w}$, the curve approaches the point $(\sigma_{w},\sigma_{b})=(1,0)$ with $q^{*}$ monotonically decreasing towards , i.e. $q^{*}\rightarrow 0$ as $\sigma_{w}\rightarrow 1$.

The critical ReLU network and the one parameter family of critical hard-tanh networks have neither vanishing nor exploding gradients, due to $\chi=1$. Nevertheless, the entire singular value spectrum of $\mathbf{J}$ of these networks can behave very differently. From eqn. (12), this spectrum depends on the non-linearity $\phi(h)$ through $S_{D^{2}}$ in eqn. (10), which in turn only depends on the distribution of eigenvalues of $\mathbf{D}^{2}$, or equivalently, the distribution of squared derivatives $\phi^{\prime}(h)^{2}$. As we have seen, this distribution is a Bernoulli distribution with parameter $p(q^{*})$: $\rho_{D^{2}}(z)=(1-p(q^{*}))\,\delta(z)+p(q^{*})\,\delta(z-1)$. Inserting this distribution into the sequence eqn. (7), eqn. (9), eqn. (10) then yields

To complete the calculation of $S_{JJ^{T}}$ in eqn. (12), we must also compute $S_{WW^{T}}$. We do this for Gaussian and orthogonal weights in the next two subsections.

We re-derive the well-known expression for the $S$-transform of products of random Gaussian matrices with variance $\sigma_{w}^{2}$ in Example 3 of the supplementary material. The result is $S_{WW^{T}}=\sigma_{w}^{-2}(1+z)^{-1}$, which, when combined with eqn. (14) for $S_{D^{2}}$, eqn. (12) for $S_{JJ^{T}}$, and eqn. (10) for $M^{-1}_{X}(z)$, yields

Using eqn. (15) and eqn. (9), we can define a polynomial that the Stieltjes transform $G$ satisfies,

The correct root of this equation is the one for which $G\sim 1/z$ as $z\to\infty$ . From eqn. (8), the spectral density is obtained from the imaginary part of $G(\lambda+i\epsilon)$ as $\epsilon\to 0^{+}$.

The positions of the spectral edges, namely locations of the minimum and maximum eigenvalues of $\mathbf{J}\mathbf{J}^{T}$, can be deduced from the values of $z$ for which the imaginary part of the root of eqn. (16) vanishes, i.e. when the discriminant of the polynomial in eqn. (16) vanishes. After a detailed but unenlightening calculation, we find, for large $L$,

Recalling that $\chi=\sigma_{w}^{2}p(q^{*})$, we find exponential growth in $\lambda_{\text{max}}$ if $\chi>1$ and exponential decay if $\chi<1$. Moreover, even at criticality when $\chi=1$, $\lambda_{\max}$ still grows linearly with depth.

Next, we obtain the variance $\sigma_{JJ^{T}}^{2}$ of the eigenvalue density of $\mathbf{J}\mathbf{J}^{T}$ by computing its first two moments $m_{1}$ and $m_{2}$. We employ the Lagrange inversion theorem ,

which relates the expansions of the moment generating function $M_{JJ^{T}}(z)$ and its functional inverse $M_{JJ^{T}}^{-1}(z)$. Substituting this expansion for $M_{JJ^{T}}^{-1}(z)$ into eqn. (15), expanding the right hand side, and equating the coefficients of $z$, we find,

Both moments generically either exponentially grow or vanish. However even at criticality, when $\chi=\sigma_{w}^{2}p(q^{*})=1$, the variance $\sigma_{JJ^{T}}^{2}=m_{2}-m_{1}^{2}=\frac{L}{p(q^{*})}$ still exhibits linear growth with depth.

Note that $p(q^{*})$ is the fraction of neurons operating in the linear regime, which is always less than $1$. Thus for both ReLU and hard-tanh networks, no choice of Gaussian initialization can ever prevent this linear growth, both in $\sigma_{JJ^{T}}^{2}$ and $\lambda_{\text{max}}$, implying that even critical Gaussian initializations will always lead to a failure of dynamical isometry at large depth for these networks.

5.2 Orthogonal weights

For orthogonal $\mathbf{W}$, we have $\mathbf{W}\mathbf{W}^{T}=\bm{I}$, and the S-transform is $S_{I}=1$ (see Example 2 of the supplementary material). After scaling by $\sigma_{w}$, we have $S_{WW^{T}}=S_{\sigma_{w}^{2}I}=\sigma_{w}^{-2}S_{I}=\sigma_{w}^{-2}$. Combining this with eqn. (14) and eqn. (12) yields $S_{JJ^{T}}(z)$ and, through eqn. (10), yields $M^{-1}_{JJ^{T}}$:

Now, combining eqn. (20) and eqn. (9), we obtain a polynomial that the Stieltjes transform $G$ satisfies:

From this we can extract the eigenvalue and singular value density of $\mathbf{J}\mathbf{J}^{T}$ and $\mathbf{J}$, respectively, through eqn. (8). Figs. 2(b) and 2(c) demonstrate an excellent match between our theoretical predictions and numerical simulations of random networks. We find that at modest depths, the singular values are peaked near $\lambda_{\text{max}}$, but at larger depths, the distribution both accumulates mass at and spreads out, developing a growing tail. Thus at fixed critical values of $\sigma_{w}$ and $\sigma_{b}$, both deep ReLU and hard-tanh networks have ill-conditioned Jacobians, even with orthogonal weight matrices.

As above, we can obtain the maximum eigenvalue of $\mathbf{J}\mathbf{J}^{T}$ by determining the values of $z$ for which the discriminant of the polynomial in eqn. (21) vanishes. This calculation yields,

For large $L$, $\lambda_{\text{max}}$ either exponentially explodes or decays, except at criticality when $\chi=\sigma_{w}^{2}p(q^{*})=1$, where it behaves as $\lambda_{\text{max}}=\frac{1-p(q^{*})}{p(q^{*})}\left(eL-\frac{e}{2}\right)+\mathcal{O}(L^{-1})$. Also, as above, we can compute the variance $\sigma^{2}_{JJ^{T}}$ by expanding $M^{-1}_{JJ^{T}}$ in eqn. (20) and applying eqn. (18). At criticality, we find $\sigma_{JJ^{T}}^{2}=\frac{1-p(q^{*})}{p(q^{*})}L$ for large $L$. Now the large $L$ asymptotic behavior of both $\lambda_{\text{max}}$ and $\sigma^{2}_{JJ^{T}}$ depends crucially on $p(q^{*})$, the fraction of neurons in the linear regime.

For ReLU networks, $p(q^{*})=1/2$, and we see that $\lambda_{\text{max}}$ and $\sigma_{JJ^{T}}^{2}$ grow linearly with depth and dynamical isometry is unachievable in ReLU networks, even for critical orthogonal weights. In contrast, for hard tanh networks, $p(q^{*})=\operatorname{erf}(1/\sqrt{2q^{*}})$. Therefore, one can always move along the critical line in the $(\sigma_{w},\sigma_{b})$ plane towards the point $(1,0)$, thereby reducing $q^{*}$, increasing $p(q^{*})$, and decreasing, to an arbitrarily small value, the prefactor $\frac{1-p(q^{*})}{p(q^{*})}$ controlling the linear growth of both $\lambda_{\text{max}}$ and $\sigma_{JJ^{T}}^{2}$. So unlike either ReLU networks, or Gaussian networks, one can achieve dynamical isometry up to depth $L$ by choosing $q^{*}$ small enough so that $p(q^{*})\approx 1-\frac{1}{L}$. In essence, this strategy increases the fraction of neurons operating in the linear regime, enabling orthogonal hard-tanh nets to mimic the successful dynamical isometry achieved by orthogonal linear nets. However, this strategy is unavailable for orthogonal ReLU networks. A demonstration of these results is shown in Fig. 3.

Experiments

Having established a theory of the entire singular value distribution of $\mathbf{J}$, and in particular of when dynamical isometry is present or not, we now provide empirical evidence that the presence or absence of this isometry can have a large impact on training speed. In our first experiment, summarized in Fig. 4, we compare three different classes of critical neural networks: (1) $\tanh$ with small $\sigma_{w}^{2}=1.05$ and $\sigma_{b}^{2}=2.01\times 10^{-5}$; (2) $\tanh$ with large $\sigma_{w}^{2}=2$ and $\sigma_{b}^{2}=0.104$; and (3) ReLU with $\sigma_{w}^{2}=2$ and $\sigma_{b}^{2}=2.01\times 10^{-5}$. In each case $\sigma_{b}$ is chosen appropriately to achieve critical initial conditions at the boundary between order and chaos , with $\chi=1$. All three of these networks have a mean squared singular value of $1$ with neither vanishing nor exploding gradients in the infinite width limit. These experiments therefore probe the specific effect of dynamical isometry, or the entire shape of the spectrum of $\mathbf{J}$, on learning. We also explore the degree to which more sophisticated optimizers can overcome poor initializations. We compare SGD, Momentum, RMSProp , and ADAM .

We train networks of depth $L=200$ and width $N=400$ for $10^{5}$ steps with a batch size of $10^{3}$. We additionally average our results over 30 different instantiations of the network to reduce noise. For each nonlinearity, initialization, and optimizer, we obtain the optimal learning rate through grid search. For SGD and SGD+Momentum we consider logarithmically spaced rates between $[10^{-4},10^{-1}]$ in steps $10^{0.1}$; for ADAM and RMSProp we explore the range $[10^{-7},10^{-4}]$ at the same step size. To choose the optimal learning rate we select a threshold accuracy $p$ and measure the first step when performance exceeds $p$. Our qualitative conclusions are fairly independent of $p$. Here we report results on a version of CIFAR-10We use the standard CIFAR-10 dataset augmented with random flips and crops, and random saturation, brightness, and contrast perturbations.

Based on our theory, we expect the performance advantage of orthogonal over Gaussian initializations to be significant in case (1) and somewhat negligible in cases (2) and (3). This prediction is verified in Fig. 4 (blue solid and dashed learning curves are well-separated, compared to red and black cases). Furthermore, the extent of dynamical isometry at initialization strongly predicts the speed of learning. The effect is large, with the most isometric case (orthogonal $\tanh$ with small $\sigma_{w}^{2}$) learning faster than the least isometric case (ReLU networks) by several orders of magnitude. Moreover, these conclusions robustly persist across all optimizers. Intriguingly, in the case where dynamical isometry helps the most ($\tanh$ with small $\sigma_{w}^{2}$), the effect of initialization (orthogonal versus Gaussian) has a much larger impact on learning speed than the choice of optimizer.

These insights suggest a more quantitative analysis of the relation between dynamical isometry and learning speed for orthogonal $\tanh$ networks, summarized in Fig. 5. We focus on SGD, given the lack of a strong dependence on optimizer. Intriguingly, Fig. 5(a) demonstrates the optimal training time is $O(\sqrt{L})$ and so grows sublinearly with depth $L$. Also Fig. 5(b) reveals that increased dynamical isometry enables faster training by making available larger (i.e. faster) learning rates. Finally, Fig. 5(c,d) and their similarity to Fig. 3(b,d) suggest a strong positive correlation between training time and max singular value of $\mathbf{J}$. Overall, these results suggest that dynamical isometry is correlated with learning speed, and controlling the entire distribution of Jacobian singular values may be an important design consideration in deep learning.

In Fig. 6, we explore the relationship between dynamical isometry and performance going beyond initialization by studying the evolution of singular values throughout training. We find that if dynamical isometry is present at initialization, it persists for some time into training. Intriguingly, perfect dynamical isometry at initialization ($q^{*}=0$) is not the best choice for preserving isometry throughout training; instead, some small but nonzero value of $q^{*}$ appears optimal. Moreover, both learning speed and generalization accuracy peak at this nonzero value. These results bolster the relationship between dynamical isometry and performance beyond simply the initialization.

Discussion

In summary, we have employed free probability theory to analytically compute the entire distribution of Jacobian singular values as a function of depth, random initialization, and nonlinearity shape. This analytic computation yielded several insights into which combinations of these ingredients enable nonlinear deep networks to achieve dynamical isometry. In particular, deep linear Gaussian networks cannot; the maximum Jacobian singular value grows linearly with depth even if the second moment remains $1$. The same is true for both orthogonal and Gaussian ReLU networks. Thus the ReLU nonlinearity destroys the dynamical isometry of orthogonal linear networks. In contrast, orthogonal, but not Gaussian, sigmoidal networks can achieve dynamical isometry; as the depth increases, the max singular value can remain $O(1)$ in the former case but grows linearly in the latter. Thus orthogonal sigmoidal networks rescue the failure of dynamical isometry in ReLU networks.

Correspondingly, we demonstrate, on CIFAR-10, that orthogonal sigmoidal networks can learn orders of magnitude faster than ReLU networks. This performance advantage is robust to the choice of a variety of optimizers, including SGD, momentum, RMSProp and ADAM. Orthogonal sigmoidal networks moreover have sublinear learning times with depth. While not as fast as orthogonal linear networks, which have depth independent training times , orthogonal sigmoidal networks have training times growing as the square root of depth. Finally, dynamical isometry, if present at initialization, persists for a large amount of time during training. Moreover, isometric initializations with longer persistence times yield both faster learning and better generalization.

Overall, these results yield the insight that the shape of the entire distribution of a deep network’s Jacobian singular values can have a dramatic effect on learning speed; only controlling the second moment, to avoid exponentially vanishing and exploding gradients, can leave significant performance advantages on the table. Moreover, by pursuing the design principle of tightly concentrating the entire distribution around $1$, we reveal that very deep feedfoward networks, with sigmoidal nonlinearities, can actually outperform ReLU networks, the most popular type of nonlinear deep network used today.

In future work, it would be interesting to extend our methods to other types of networks, including for example skip connections, or convolutional architectures. More generally, the performance advantage in learning that accompanies dynamical isometry suggests it may be interesting to explicitly optimize this property in reinforcement learning based searches over architectures .

S.G. thanks the Simons, McKnight, James S. McDonnell, and Burroughs Wellcome Foundations and the Office of Naval Research for support.

References

Theoretical results

The $S$-transform for $JJ^{T}$ is given by,

First notice that, by eqn. (9), $M(z)$ and thus $S(z)$ depend only on the moments of the distribution. The moments, in turn, can be defined in terms of traces, which are invariant to cyclic permutations, i.e.,

where the last line follows since each weight matrix is identically distributed. ∎

Products of Gaussian random matrices with variance $\sigma_{w}^{2}$ have the $S$ transform,

It is well-known (see, e.g. ) that the moments of a Wishart are proportional to the Catalan numbers, i.e.,

It is straightforward to invert this function,

The $S$-transform of the identity is given by $S_{I}=1$.

The moments of the identity are all equal to one, so we have,