The jamming transition as a paradigm to understand the loss landscape of deep neural networks

I Introduction

Deep neural networks are now central tools for a variety of tasks including image classification Krizhevsky et al. (2012); LeCun et al. (2015), speech recognition Hinton et al. (2012) and the development of artificial intelligence that can for example master the game of Go beyond human level Silver et al. (2016, 2017). A neural network represents a (very high-dimensional) function $f$ that depends on a large number of parameters $N$ LeCun et al. (2015). These parameters are learned so as to correctly classify $P$ training data by minimizing some loss function ${\cal L}$, generally via stochastic gradient descent (a kind of noisy version of gradient descent). There is great flexibility in the network architecture, loss function and minimization protocol one can use. These features are ultimately selected to optimize the classification of previously unseen data, or generalization. Although the current progress in designing LeCun et al. (1995); He et al. (2016) and training Ioffe and Szegedy (2015) networks that generalize well is undeniable, it remains mostly empirical. A general theory explaining and fostering this success is lacking, and central questions remain to be clarified. First, since the loss function is generally not convex, why doesn’t the learning dynamics get stuck in poorly performing minima with high loss? In other words, under which conditions can one guarantee that training data are well fitted? Second, what are the benefits of deeper networks? On the one hand it is often argued, and proved in some cases, that the advantage of deep networks stems from their enhanced expressive power, i.e. their ability to build complex functions with a much smaller number of parameters than needed for shallow networks Montufar et al. (2014); Bianchini and Scarselli (2014); Raghu et al. (2017); Eldan and Shamir (2016); Lee et al. (2017). Indeed if deep networks are able to fit data with less parameters, then they are likely to generalize better. On the other hand, one can handcraft neural networks that fit even structure-less, random data with a rather small number of parameters $N\sim P$ Gardner (1988); Monasson and Zecchina (1995); Zhang et al. (2017); Baum (1988). These results for the static capacity of networks appear to be independent of depth Zhang et al. (2017); Baum (1988). Yet, it is unclear whether such parsimonious solutions can be found dynamically in practice simply by descending the loss function, and whether depth can help finding them. More generally, how is the loss landscape affected by depth?

Complex physical systems with non-convex energy landscapes featuring an exponentially large number of local minima are called glassy Berthier and Biroli (2011). Does the landscape of deep learning fall into a known class of glassy systems? Along this line, an analogy between deep networks and mean-field glasses ($p$-spins) has been proposed Choromanska et al. (2015), in which the learning dynamics is expected to get stuck in the highest minima of the loss, which are the most abundant. Yet, several numerical and rigorous works Freeman and Bruna (2017); Hoffer et al. (2017); Soudry and Carmon (2016); Cooper (2018) (the latter focusing on shallow and very overparametrized networks) suggest a different landscape geometry where the loss function is characterized by a connected level set. Furthermore, studies of the Hessian of the loss function Sagun et al. (2017a, b); Ballard et al. (2017) and of the learning dynamics Lipton (2016); Baity-Jesi et al. (2018) support that the landscape is characterized by an abundance of flat directions, even near its bottom, at odds with traditional glassy systems.

In the last decade several works have unveiled an analogy between the physical phenomenon of jamming Wyart (2005); J. Liu et al. (2010) and phase transitions taking place in certain classes of computational optimization and learning problems Krzakala and Kurchan (2007); Zdeborová and Krzakala (2007); Franz et al. (2017), in particular the perceptron Franz and Parisi (2016); Franz et al. (2017) — the simplest neural network performing linear classification. In this work we push this analogy further and show that the geometry of the training loss landscape and the training dynamics of fully connected deep neural networks is affected by a jamming transition similar to that of repulsive ellipses J. Liu et al. (2010). As illustrated in Fig. 2, jamming occurs in packings of particles interacting through a finite-range potential ${\cal U}$, when the particle density $\phi$ reaches some critical value $\phi_{c}$. At that point, particles can no longer be accommodated without touching each other and the system becomes a solid with singular landscape properties, embodied for example in the spectrum of the Hessian of ${\cal U}$ Wyart et al. (2005); Silbert et al. (2005), that at the transition displays many (almost) flat directions. Particles of different shapes, such as spheres and ellipses, can lead to different jamming scenarios Donev et al. (2004); Mailman et al. (2009); Zeravcic et al. (2009); Brito et al. (2018).

Here we show that for two commonly used loss functions (cross-entropy and quadratic hinge), fully-connected deep networks undergo a jamming transition too, below which all data are correctly fitted and above which they are not, both for real data (images) and random data This transition influences the generalization properties of deep networks, too. This has been observed, for instance, in Advani and Saxe (2017); Spigler et al. (2018), and studied by the authors in Geiger et al. (2019) (preprint).. In both cases the transition appears to be solely controlled by the number of parameters of the network $N$, independently of depth. For random data, the transition takes place as the quantity $P/N$ increases toward some critical value $P/N^{*}$. For the hinge loss, using results from the jamming literature we argue that $P/N^{*}\geq C_{0}$ where $C_{0}$ is a constant that we can measure a posteriori once learning took place. To hold, this result requires the network output to remain sensitive to all its weights during training, as we observe empirically in the examples we study. This view supports that the dynamics cannot get stuck in poor minima in the over-parametrized regime where networks tend to operate, because there are not enough constraints to form minima in that regime. We also find that the jamming transition is sharp and the landscape appears to fall in the same universality class independently of depth (as long as at least one hidden layer is present). Differently from the (non-convex) perceptron, that was proven to lie in the same universality class as spherical particles Franz and Parisi (2016); Franz et al. (2017), we show that deep networks instead jam in a manner similar to ellipses. From our analysis we deduce the singular properties of the spectrum of the Hessian of the loss, which indeed must display many flat directions. We find empirically that other key quantities (the fraction of data which are almost correctly or almost incorrectly classified) display power-law behaviours on several decades, with new exponents. In glassy systems, such power-laws reveal properties that cannot be reached by studying the Hessian, in particular the fact that the dynamics occurs via broadly distributed avalanches Wyart (2012); Müller and Wyart (2015); Franz and Spigler (2017), indicative of a hierarchical organization of the landscape Charbonneau et al. (2014a). This observation thus suggests that these properties also characterize deep networks near the transition. Note that in this work we focus on training and the ability of deep neural networks to fit a dataset. The implications and the relations with generalization between jamming and generalization are investigated in Spigler et al. (2018); Geiger et al. (2019).

II Analogy between jamming and deep learning

Understanding the energy landscape — in particular the properties of the Hessian, referred to as vibrational properties in this context — in disordered systems of interacting particles is a long-standing and practically important problem Anderson (1981). It was realized that for purely repulsive, finite-range particles, such properties are singular near the jamming transition where the system becomes a solid Silbert et al. (2005); Wyart et al. (2005), allowing one to develop and test theories for the vibrations of glasses, that turn out to apply in a broader class of systems where the interactions do not necessarily have finite range Wyart (2005). Here we shall follow the same strategy for deep networks, where the role of the “interaction potential” is played by the choice of loss function. Finite-range interactions are mimicked by the hinge loss, for which we predict a sharp transition when going from an overparametrized to an underparametrized regime. At the transition, the Hessian is singular and displays an abundance of low-energy modes. For other types of losses — such as for the commonly used cross-entropy loss defined below — the transition exists but its effects on the spectrum are expected to be less sharp (see discussion below).

As the jamming transition is approached from above (large density $\phi$), ${\cal U}\rightarrow 0$ as sketched in Fig. 2, implying that $\Delta_{\mu}\rightarrow 0$ $\forall\mu\in m$. As argued in Tkachenko and Witten (1999), for each $\mu\in m$ the constraint $\Delta_{\mu}=0$ defines a manifold of dimension $N-1$. Satisfying $N_{\Delta}$ such equations thus generically leads to a manifold of solutions of dimension $N-N_{\Delta}$. Imposing that solutions exist thus implies that, at jamming, one has

Note that this argument implicitly assumes that the $N_{\Delta}$ constraints are independent. In disordered systems this assumption is generally correct in practice, but it may break down if symmetries are present, which is the case e.g. in crystals where Eq. (2) can be violated.

An opposite bound can be obtained for spheres by considerations of stability, by imposing that in a stable minimum the Hessian must be positive definite Wyart et al. (2005). The Hessian is an $N\times N$ matrix which can be written as A similar decomposition has been used in Pennington and Bahri (2017); Martens (2014); Sagun et al. (2017b).

where ${\cal H}_{0}$ and ${\cal H}_{p}$ correspond to the first and second sum, respectively. ${\cal H}_{0}$ is positive semi-definite, since it is the sum of $N_{\Delta}$ matrices of rank unity; thus ${\rm rk}({\cal H}_{0})\leq N_{\Delta}$, implying that the kernel of ${\cal H}_{0}$ is at least of dimension $N-N_{\Delta}$. On the other hand for spheres — but not for ellipses, and this will have major consequences — ${\cal H}_{p}$ is negative definite, which simply stems from the fact that the second-order contribution of the displacement to the distance between two points is always positive - a straightforward application of the Pythagoras theorem. It is easy to show Wyart et al. (2005) that any non-zero vector $|u\rangle$ belonging to the kernel of ${\cal H}_{0}$ must satisfy $\langle u|{\cal H}_{U}|u\rangle=\langle u|{\cal H}_{p}|u\rangle<0$ Once again, this statement is true except for global translation or rotation of the systems, whose number however is fixed in the large $N$ limit and disappears when the ratio $N_{\Delta}/N$ is considered.. Thus stability requires that ${\rm rk}({\cal H}_{0})=N$, implying that $N_{\Delta}\geq N$. Together with Eq. (2) that leads to $N_{\Delta}=N$: as spheres jam the number of degrees of freedom and the number of constraints (stemming from contacts) are equal, as empirically observed O’Hern et al. (2003). This property is often called isostaticity: when it holds, mean-field arguments DeGiuli et al. (2014); Yan et al. (2016); Franz et al. (2015) predict that the density of vibrational modes $D(\sqrt{\lambda})$ displays a plateau up to vanishingly small $\lambda$, as observed numerically Silbert et al. (2005); Wyart et al. (2005) and sketched in Fig. 2G.

However, for ellipses Donev et al. (2004) (and as we shall see, for deep networks), this argument breaks down because ${\cal H}_{p}$ is not negative definite. Whether such a matrix has positive eigenvalues or not plays a role of utmost importance in the selection of the universality class of the jamming transition, and it has major consequences on the singularity of the landscape, as it can be evinced from the spectrum of the Hessian matrix. Indeed, for ellipses stability and jamming can, and generically do, occur at:

a situation that is referred to as hypostatic. The density of vibrational modes at jamming must then display a delta function in zero of magnitude $1-N_{\Delta}/N$, corresponding to the kernel of ${\cal H}_{0}$ (${\cal H}_{p}$ vanishes at jamming since $\Delta_{\mu}\rightarrow 0$ $\forall\mu\in m$). Mean-field arguments applied to hypostatic materials Düring et al. (2013); Brito et al. (2018) predict that at larger $\lambda$, the spectrum presents a gap before becoming continuous again, as sketched in Fig. 2H. Away from jamming the effects of ${\cal H}_{p}$ kick in and broaden the delta function by an amount proportional to the typical value of the overlap $\Delta\sim\sqrt{\cal U}$, as sketched in Fig. 2H.

We now show that even in the hypostatic case, stability can be constraining. Let us denote by $E_{-}$ the vector space spanned by the negative eigenvalues of ${\cal H}_{p}$, whose dimension very close to jamming is denoted $N_{-}$. Stability then imposes that the intersection of the kernel of ${\cal H}_{0}$ and $E_{-}$ is zero, which is possible only if

Finally, another key structural property of the jamming transition is contained in the distribution $P_{+}(\Delta)$ of positive overlaps, sometimes referred to as forces (the force between two particles is $\Delta$ when $\Delta>0$), and the distribution $P_{-}(\Delta)$ of gaps ($\Delta<0$) between particles. It was shown that even if a packing of spherical particles is linearly stable, paths in the phase space that lower the energy are easily found unless both distributions are critical, with $P_{+}(\Delta)\sim\Delta^{\theta}$ and $P_{-}(\Delta)\sim(-\Delta)^{-\gamma}$, with $\gamma\geq(1-\theta)/2$ Wyart (2012); Lerner et al. (2013), as numerically confirmed in Lerner et al. (2012); Charbonneau et al. (2012). For a broad class of dynamics, this bound must be saturated Müller and Wyart (2015), a scenario referred to as marginal stability which implies that the dynamics proceeds via power-law distributed events (called avalanches) in which the set of constraints change. Calculations in infinite dimensions Charbonneau et al. (2014a, b) showed that marginal stability is associated with a hierarchical organization of minima of the energy (a phenomenon referred to as replica symmetry breaking Mézard et al. (1987)), and exponents were found to follow $\gamma=0.41269\ldots$ and $\theta=0.42311\ldots$ which appear accurate even in finite dimensions Lerner et al. (2013); Charbonneau et al. (2015).

II.2 Deep Learning

Performance of the hinge loss and its extension to multi-class problems: This section can be skipped at a first reading. We tested in the context of image classification that the hinge loss performs as well as the cross entropy on a state-of-the-art architecture Gastaldi (2017): we ran the implementation https://github.com/mariogeiger/pytorch_shake_shake for CIFAR-10 and we retrained it by replacing the cross entropy by the hinge loss. To compare the two losses in a standard setting we adapted the hinge loss for multiple classes, although in what follows we only study binary classification. To predict the label of an input $\mathbf{x}_{\mu}$ among $10$ possible labels $c=0,\dots,9$, the network’s last layer returns as output a list of 10 values $f_{\mu,c}$: each $f_{\mu,c}$ can be interpreted as the probability that $c$ is the predicted label. Let $t_{\mu,c}$ be the true target labels: for each $\mu$, $t_{\mu,c}$ is equal to $1$ if $c$ is equal to the label of $\mathbf{x}_{\mu}$ and $-1$ otherwise. Multiclass hinge loss can then be written as

We obtained an error of $3.72\%$ by running their original code (they report on github an error of $3.68\%$) and $3.61\%$, $3.65\%$, $3.82\%$ in three runs with the hinge loss.

Symmetries are present in the network, e.g. the scale symmetry in ReLU networks. It will reduce one degrees of freedom per node.

Constraints on the stability of minima: Let us suppose (and justify later) that for a fixed number of data $P$, if $N$ is sufficiently large then gradient descent with proper weights initialization leads to ${\cal L}=0$, whereas if $N$ is very small after training ${\cal L}>0$. Consider that $N$ is increased from a small value. At some value $N^{*}$ the loss obtained after training will approach zero For finite $P$, $N^{*}$ will present fluctuations induced by differences of initial conditions. The fluctuations of $P/N^{*}$ are however expected to vanish in the limit where $P$ and $N^{*}$ become large. This phenomenon is well-known for the jamming of particles, and is referred to as finite size effects. , i.e. $\lim_{N\rightarrow N^{*}}{\cal L}=0$. In analogy with the behavior of packings of particles, we refer to this point as the jamming transition. At the transition the stability constraint developed in Eq. (5) above also applies if the derivative of $f(\mathbf{x};\mathbf{W})$ is continuous, which holds true if the non-linear function $\rho$ is smooth. Thus we have:

since $P\geq N_{\Delta}$ (the number of unsatisfied patterns is obviously smaller than the total number of patterns).

We shall assume that the fraction $N_{-}/N\equiv C_{0}$ of negative eigenvalues of ${\cal H}_{p}$ does not vanish in the large $N$ limit. In Appendix B we provide an argument supporting this result in the case of a specific non-linear function (ReLU) and random data, that yields $C_{0}=1/2$ independently of depth. It implies that unlike for spheres, but just like ellipses, ${\cal H}_{p}$ is not negative definite: we are therefore in the hypostatic scenario where one expects $N_{\Delta}<N^{*}$ at jamming, a point at which the spectrum must display a fraction of flat directions, as well as stiff ones, as described in Fig. 2H.

Moreover from this assumption and Eq.11, we obtain that stability cannot be obtained for $N\geq P/C_{0}$. For larger $N$, the dynamics cannot get stuck in a bad minimum, because in this over-parametrized regime there are not enough constraints to form them. It implies for the jamming transition that:

Notice that this bound is expected to be valid for any monotonic cost function, as for instance the cross entropy (the Hessian can always be decomposed as in Eq. (3)). However, the spectrum of the Hessian would be different For the cross entropy, when the data become separable true minima exists only for diverging weights, a complication that does not occur with the hinge loss..

Smooth vs non-smooth output function: In our numerical study below, we consider the most common choice for the non-linear function $\rho$, namely the rectified linear unit (ReLU): $\rho(a)=a\,\Theta(a)=\max(0,a)$. In that case, as stated above we expect for random data the spectrum of ${\cal H}_{p}$ to be symmetric (a fact that appears to also hold true for the image dataset we use, see below), thus $N_{-}/N=1/2$. Yet, with the ReLU, $f(\mathbf{x};\mathbf{W})$ is not continuous and presents cusps, so that the Hessian needs not be positive definite for stability and Eq. (11) needs to be modified. Introducing the number of directions $N_{c}$ presenting cusps, stability implies $N_{\Delta}\geq N_{-}-N_{c}$ leading to $C_{0}N^{*}\leq P+N_{c}$. Empirically we find that $N_{c}/N^{*}\in[0.21,0.25]$ both for random data and images as reported in Appendix C, implying that:

Main results: Overall, our analysis supports that

In the case of hinge loss there is a sharp transition for $N^{*}\leq P/C_{0}$ ($N^{*}<4P$ with the ReLU), below which the loss converges to some non-zero value (under-parametrized phase) and above which it becomes null (over-parametrized phase).

At that point the fraction $N_{\Delta}/N$ of unsatisfied constraints per degree of freedom jumps to a finite value, see Fig. 6 (a,b).

Unlike for spheres or the perceptron, isostaticity $N_{\Delta}/N=1$ cannot be guaranteed. Instead one expects generically $N_{\Delta}/N<1$ as for ellipses.

We are thus in the hypostatic universality class, where the scaling properties of the spectrum of the Hessian near jamming are prescribed in Fig.2.

In the next sections, we confirm these predictions in numerical experiments and observe the generalization properties at and beyond the transition point.

III For random data the transition occurs for 𝑵∼𝑷similar-to𝑵𝑷\bm{N\sim P}

We begin the numerical study of the transition between the overparametrized and underparametrized regime in the case of random data, taken to lie on the $d$-dimensional hyper-sphere of radius $\sqrt{d}$, ${\bf x}_{\mu}\in{\cal S}^{d}$ with random label $y_{\mu}=\pm 1$. The source code used to generate the simulations described in this section and the following ones is available at https://github.com/mariogeiger/nn_jamming. We proceed as follows: we build a network with a number of weights $N$ large enough for it to be able to fit the whole dataset without errors. Next, we reduce the number of weights by decreasing the width $h$ while keeping the depth $L$ fixed, until the network cannot correctly classify all the data anymore within the chosen learning time. We denote this transition point $N^{*}$.

We have noticed (data not shown) that the precise location of the transition point $P/N^{*}$ has a mild dependence on the dynamics (ADAM versus regular SGD, choice of batch size, learning rate schedule, etc…): the same holds true for the jamming of repulsive particles, where the choice of the dynamics affects the precise value of the critical density $\phi_{c}$, but not the critical behaviour close to this point.

Cross-entropy loss: We first consider the cross-entropy loss — the results are qualitatively similar to those with the hinge loss. As initial condition for the dynamics we use the default initialization of pytorch Weights and biases are initialized with a uniform distribution on $[-\sigma,\sigma]$, where $\sigma^{2}=1/f_{in}$ and $f_{in}$ is the number of incoming connections.. The system then evolves according to a stochastic gradient descent (SGD) with a learning rate of $10^{-2}$ for $5\cdot 10^{5}$ steps and $10^{-3}$ for $5\cdot 10^{5}$ steps; the batch size is set to $\min(P/2,1024)$; only in this case, with the cross-entropy loss, batch normalization is also used. In Fig. 5 (a) we show how $N^{*}$ depends on the total learning time: the larger is the learning time the more the asymptotic relationship $N^{*}$ vs $P$ is consistent with an asymptotic linear behaviour. Note that for large $P$ and small times, errors are always present and the transition cannot be found.

In Fig. 5 (b) we show $N^{*}$ versus the number of data $P$ after $t=10^{6}$ steps for several depths $L$ and input dimensions $d$ (we checked that $t=10^{6}$ is enough to get convergence to the conjectured asymptotic linear behaviour for all depths investigated). It is noteworthy that (i) the points always lie below the theoretical upper bound $P/N^{*}=1/2-N_{c}/N$, and (ii) the transition does not appear to depend on $L$ and $d$. Surprisingly, this result indicates that in the present setup the ability of fully connected networks to fit random data is independent of the depth. As we shall see, we observe the same independence on depth for the image data studied below.

Hinge loss: In order to test the dependence of our results on the specific choice of the loss function, we performed the same experiment using the hinge loss. In this case we used an orthogonal initialization Saxe et al. (2014), no batch normalization and $t=2\cdot 10^{6}$ steps of ADAM Kingma and Ba (2015) with batch size $=P$ and a learning rate starting at $10^{-4}$, progressively divided by $10$ every 250k steps. The location of the transition is shown in Fig. 5 (c): results are very similar to that of the cross-entropy loss.

Hinge v.s. cross-entropy loss from a conceptual perspective: As shown above, both losses appears to lead to similar performances. As shown in this section, both of them also displays a transition where all data are fitted. Yet, the nature of this transition is harder to investigate for the cross-entropy. Indeed in that case the total loss is never zero, except if the output and therefore the weights diverge. Thus in the over-parametrized phase, the learning dynamic never settles, and the weights slowly drift to infinity. In practice, users stop learning at finite times (which is not needed for the hinge loss where the dynamics really stops in the over-parametrized regime when the loss vanishes). Working at finite time however blurs true critical behavior near jamming, as discussed in Geiger et al. (2019).

IV The transition is hypostatic

From the analysis of Section II, the number of constraints per parameter $N_{\Delta}/N$ is expected to jump discontinuously at the transition. To test this prediction we consider several architectures, both with $N\approx 8000$ and $d=h$ but with different depths $L=2$, $L=3$ and $L=5$. The vicinity of the transition is studied by varying $P$ around the transition value. We used the hinge loss with the same gradient descent dynamics as described above, for a duration of $10^{7}$ steps. Fig. 6 (a) reports the ratio $N_{\Delta}/N$ as a function of the ratio $P/N$ and of the learning time, as detailed in caption. It is clear that in the range where $N_{\Delta}/N$ has reached a stationary value (i.e. for $P/N<2.8$ and $P/N>2.9$), a jump has occurred from 0 to $N_{\Delta}/N\approx 0.75$, a result consistent with the bound of Eq. (5) implying $N_{\Delta}/N\geq(N_{-}-N_{c})/N\gtrapprox 0.25$. For $P/N\in[2.8,2.9]$, the dynamics has not yet converged and the data are somewhat scattered. This observation is presumably the signature of the usual slowing down that occurs near critical points.

Fig. 6 (b) shows the same quantity $N_{\Delta}/N$, now plotted as a function of the loss $\mathcal{L}$. Strikingly, all the scatter is gone, and one observes a clear discontinuous behaviour for ${\cal L}\rightarrow 0$. Interestingly, this state of affairs is very similar to the jamming transition of particles, for which the noise in the data due to finite size effects is quite strong when quantities are expressed in terms of the density $\phi$ (analogous to $P/N$) but very small when quantities are expressed in terms of potential energy ${\cal U}$ (analogous to ${\cal L}$) O’Hern et al. (2003).

For the sake of completeness we also show the number of misclassified data as a function of the loss in Fig. 6 (c). The number of misclassified data increases monotonically — and initially very slowly — with the loss. Indeed, close to the jamming threshold in the underparametrized phase, if $0<\Delta_{\mu}<\epsilon$ the pattern $\mu$ is well classified but the corresponding gap $\Delta_{\mu}$ is positive: unsatisfied constraints do not lead to misclassification right away.

In Fig.6 (d) we show that $N_{\Delta}/N$ vs the loss $\mathcal{L}$ exhibits a sharp transition also for networks with tanh activation functions.

V Spectrum of the Hessian of the loss near Jamming

The Hessian is a key feature of landscapes, as it characterizes its curvature, and it is also a central aspect of the theoretical description above. In this section we systematically analyze the spectra of $\mathcal{H}$, $\mathcal{H}_{0}$ and $\mathcal{H}_{p}$. To test the predictions on the singularity of the Hessian matrix, we need to focus on the underparametrized data points near the transition. These points are contained in the black rectangle on the left side of Fig. 6 (b). The networks that we use are relatively small, but, for reference, it would be possible to compute the spectrum of the Hessian also for large networks, as discussed e.g. in Adams et al. (2018); Ghorbani et al. (2019). The setting is as above: the network uses the hinge loss and is trained with ADAM with full batches (batch size $=P$), orthogonal initialization and no batch normalization.

Relu networks: At the end of each run, we compute the hessian $\cal H$ of the loss $\cal L$, as well as the two terms ${\cal H}_{0}$ and ${\cal H}_{p}$ contributing to it, as defined in Eq. (3). Fig. 7 (a) shows the positive part of the spectrum of ${\cal H}_{p}$ for different values of the loss, illustrating that the dependence on the latter is very significant. In Fig. 7 (b) we confirm that the spectrum of ${\cal H}_{p}$ collapses when the eigenvalues are re-scaled by ${\cal L}^{1/2}$, as expected from Section II. The key observation is that these spectra are symmetric, as argued in Appendix B. We also don’t observe any accumulation of eigenvalues at $\lambda=0$, except for the trivial zero modes stemming from the scaling symmetry of ReLU neurons (whose number is the total number of hidden neurons, much smaller than the number of weights). Fig. 7 (c) shows the spectrum of $\mathcal{H}_{0}$ at the end of training for runs close to the jamming transition. As expected it is semi-positive definite, with a delta peak at $\lambda=0$ corresponding to $N-N_{\Delta}$ modes. It is followed by a gap and a continuous spectrum, as predicted near the jamming transition of particles if $N_{\Delta}<N$ Düring et al. (2013) (which occurs for elliptic particles Brito et al. (2018)). As the loss increases, $N_{\Delta}$ increases and the gap is reduced. Finally in Fig. 7 (d), the spectrum of $\cal H$ is shown. Interestingly the spectrum of the Hessian is not positive definite, but present some unstable modes. This phenomenon stems from our choice of ReLU activation function, which leads to cusps in the landscape as quantified in the C. Such cusps can stabilize directions that would be unstable according to the Hessian.

Overall, as we move from the under-parametrized phase to the over-parametrized one the situation is as follows:

$N$ below $N^{*}$: There are many constraints with respect to the number of variable $N$, $\mathcal{H}_{0}$ is almost full rank and can easily compensate the negative eigenvalues of $\mathcal{H}_{p}$. The spectrum of $\mathcal{H}_{p}$ is symmetric.

$N$ approaching $N^{*}$ from below: The rank of $\mathcal{H}_{0}$ decreases but it does not go below $C_{0}N$ since it has to compensate the vanishingly small negative eigenvalues of $\mathcal{H}_{p}$.

As $N$ is large enough, the dynamics finds a global minimum at $\mathcal{L}=0$ and $\mathcal{H}_{p}$ vanishes.

VI Distribution of gaps reveals new singular behaviour

We now study the distribution of gaps $\Delta<0$ and overlaps $\Delta>0$, which play an important role near jamming. Positive $\Delta$’s are associated with unsatisfied patterns — which increase the loss of the system — whereas negative $\Delta$’s correspond to satisfied patterns — which are correctly classified with a margin $\epsilon$ and do not contribute to the loss. The latter offer an important measure not only at the jamming transition, but also in the overparametrized regime, where they signal how much room is left around a minimum of the loss to fit additional patterns. In Fig. 8 (a,b) we show the two distributions for different depths $L=2,3,5$ (positive $\Delta$’s have been rescaled by $\mathcal{L}^{1/2}$). Remarkably, they behave as power laws for about two decades, $P_{+}(\Delta/\sqrt{\cal L})\sim(\Delta/\sqrt{\cal L})^{\theta}$ and $P_{-}(\Delta)\sim|\Delta|^{-\gamma}$, with novel exponents $\theta\approx 0.3$ and $\gamma\approx 0.2$ that appear to differ from those found for the jamming of particles (which are $\theta\approx 0.42311\ldots$ and $\gamma\approx 0.41269\ldots$ ). For comparison, in Fig. 8 (c,d) we show the distribution of the same variables for tanh-networks, that also display power-law behaviors but with different exponents $\theta\approx 0.2$ and $\gamma\approx 0.16$.

In the case of spheres, the two exponents are related by an inequality that happens to be saturated Wyart (2012); Lerner et al. (2013). The inequality comes from arguments on the stability of jammed packings, and the fact that it is saturated (which can be proven for certain dynamics Müller and Wyart (2015)) implies that such systems are marginally stable: they display an abundance of low-energy excitations and are prone to avalanche dynamics and crackling response when perturbed Müller and Wyart (2015), a property associated with a hierarchical organization of the loss landscape Charbonneau et al. (2014a); Franz and Spigler (2017). The presence of such power laws for deep networks thus suggests they are marginally stable as well, and that the learning dynamics may occur by avalanches where the unsatisfied constraints change by bursts. This will be subject of detailed studies in a future paper.

VII Image data: MNIST

We now consider a dataset called MNIST, which consists of a collection of black and white pictures of $28\times 28$ pixels depicting handwritten digits from 0 to 9. The labels $y_{\mu}$ in principle would be the digits themselves ($y_{\mu}\in\left\{0,\dots,9\right\}$), but to compare more directly with our previous experiments we gathered all the digits into two groups (even and odd numbers) with labels $y_{\mu}=\pm 1$. The architecture of the network is as in the previous sections: the $d$ inputs are fed to a cascade of $L$ fully-connected layers with $h$ neurons, that in the end result in a single scalar output. The loss function used is the hinge loss.

If we kept the original input size of $28\times 28=784$ then the majority of the network’s weights would be necessarily concentrated in the first layer (the width $h$ cannot be too large in order to be able to compute the Hessian). To avoid this issue, we opted for a reduction of the input size. We performed a principal component analysis (PCA) on the whole dataset and we identified the 10 dimensions that carry the most variance; then we used the components of each image along these directions as a new input of dimension $d=10$. This projection hardly diminishes the performance of the network (we find the generalization accuracy to be larger than $90\%$ at the jamming transition in Fig. 10 for $P\geq 10^{4}$).

In Fig. 9 we show that a jamming transition is also found for real data with a discontinuous behavior of $N_{\Delta}/N$. Fig. 9 (a) shows the number of unsatisfied patterns per parameter $N_{\Delta}/N$ increasing $P$ at fixed $N$, and in Fig. 9 (b) the same quantity is plotted against the loss. As for random data, the latter is less noisy. In Fig. 9 (c) we show that the number of misclassified data (i.e. the number of patterns with $y_{\mu}f(\mathbf{x}_{\mu})<0$) grows smoothly with the loss. These plots depict the same scenario as we found for random data, namely the one presented in Fig. 6 (a-c), except for the magnitude of the density of constraints at the transition with $N_{\Delta}/N\approx 0.5$ rather than $N_{\Delta}/N\approx 0.7$ as observed before. Hence, the number of unsatisfied patterns at the transition is not universal.

Also the spectrum of the Hessian matrix is similar to that of random data. In Fig. 9 (e-h) we show the positive part of the spectrum of $\mathcal{H}_{p}$, the total spectrum of $\mathcal{H}_{p}$, the spectrum of $\mathcal{H}_{0}$ and the spectrum of the total Hessian $\mathcal{H}$, respectively. As with random data: the matrix $\mathcal{H}_{p}$ has a symmetric spectrum and the matrix $\mathcal{H}_{0}$ has a finite number of zero modes and a gapped continuous distribution of modes at high energy. The spectrum of the total Hessian is again similar to that of $\mathcal{H}_{0}$, where the delta function in zero has been smeared.

The distribution of gaps (negative $\Delta$’s) is plotted in Fig. 9 (d), suggesting a power law with an exponent $\gamma=0.25$ that is slightly larger than the value found for random data, $\gamma\approx 0.2$. It is unclear whether this difference is significant. We observed that the distribution of overlaps (positive $\Delta$’s) has large sample to sample variations (not shown), and the acquisition of enough statistics to measure it extensively will be done elsewhere.

A key difference between random and structured data however is the location $N^{*}$ of the transition, shown in Fig. 10 versus the number $P$ of patterns. For a fixed number $P$ of MNIST pictures we ran several simulations with networks of different sizes, and found in this way the lowest value $N^{*}$ for which all patterns could still be classified correctly. In the figure we present the results for two network architectures of different depths $L=1,3,5$ (the width $h$ was varied in order to control the network size). Key results are that (i) $N^{*}$ is essentially independent of depth, especially at larger $P$ and (ii) the minimum number of parameters $N^{*}$ to fit the data is significantly smaller than for random data, a difference that seems to increase with $P$. The behavior of $N^{*}$ in the (hypothetical) limit $P\rightarrow\infty$ could be indeed different from the linear scaling of random data: a sub-linear scaling or even a finite asymptotic value are possible alternatives. More generally, how the data structure affects the location of the transition $N^{*}(P)$ is an important question for the future.

VIII Conclusion

By slightly changing the loss function — i.e. by considering the hinge loss rather than the commonly used cross entropy, a change that does not degrade performance — we could recast the problem of minimizing the loss function of deep networks into a constraint satisfaction problem with continuous degrees of freedom. This kind of problem has been abundantly studied in physics, in particular in the context of the jamming of particles, and some theoretical tools developed in that field readily apply to deep networks. In particular from this analogy one predicts a sharp transition as the number of parameters is reduced, separating a region where all constraints can be satisfied (that is, all the data are perfectly fitted) and the loss is zero after learning, and a region where the ratio of the number of unsatisfied constraints to the number of parameters is of order one. This ratio jumps discontinuously at the transition, where it attains a value smaller than one. Near that point, the spectrum of the Hessian is singular, reminiscent of a critical behavior. One key finding is that deep learning falls into the hypostatic universality class, similar to that of ellispes. We also observe a scaling behavior and new exponents characterizing how well constraints are satisfied or not (through the distributions $P_{-}(\Delta)$ and $P_{+}(\Delta)$, respectively). This bears comparison with the known behavior of packings of particles — where such singularities signal marginality and avalanche-type response — and of the perceptron (the simplest, shallow, neural network), that lie in the same universality class. Yet there is no theory so far to explain these exponents for deep networks. These results also shed light on some aspects of deep learning:

Not getting stuck in poor minima of the loss: Our analysis supports that in the overparametrized regime, the dynamics does not get stuck in poor minima because the number of constraints to satisfy (data to fit) P is too small to hamper minimization: the system is in an easy satisfiable phase. In particular assuming that a certain operator (namely the matrix $\mathcal{H}_{p}$) has a fraction of negative eigenvalues (which we could show in the case of the ReLU activation function and random data, and confirm numerically) implies that no poor minima exist if $P/N<P/N^{*}=\mathcal{O}(1)$. Here $N$ is the number of effective degrees of freedom of the network, which in all the cases we studied is essentially equal to the number of parameters. This argument does not rule out the possibility that, with a very poor choice of initial condition, a poor minimum of the loss can be found. This is the case in particular if the network does not propagate the signal (then $N=1$ in our formalism, independently of the number of parameters). Presumably usual tricks used to train deep networks (batch normalization, residual links, proper weight initialization, …) ensure that the sensitivity of the network to its parameters is preserved during training so that $N$ is indeed similar to the number of parameters, a hypothesis that would be useful to test in a broader setting.

In the under-parametrized phase the network gets stuck at a positive loss, either because the ground state is no longer at zero loss or because the system is trapped in an excited local minimum. The fact that the jamming transition itself depends on the dynamics (as is the case for the jamming of particles) suggests that in the underparametrized case the network is in a local minimum.

Role of depth: We observed that depth is not helpful to fit random data in fully connected networks: increasing depth and reducing width so that the total number of weights is fixed does not allow to fit the data with less parameters. We have also observed that this finding continues to hold in a realistic case based on MNIST. This may seem to clash with mathematical results, such as Montufar et al. (2014); Bianchini and Scarselli (2014); Raghu et al. (2017); Lee et al. (2017), which establish that depth enhances expressivity. However, we tackle the question of expressivity for realistic data and learning protocols, which is quite different. Our results, that need confirmation by further studies on a broader range of data, point toward a negative answer for fully connected networks. It may be that the added expressive power of deep networks is only useful for architectures exploiting the symmetry and hierarchy in the data (e.g. as in convolutional networks). Alternatively, depth may play a role in accelerating the learning dynamics Shwartz-Ziv and Tishby (2017).

Reference point for network architectures: key properties of deep networks, including the learning dynamics and the generalization power, are believed to be affected by the landscape geometry. We have argued that there exists a critical line $N^{*}(P)$ where the landscape is singular (with both flat and stiff directions), suggesting that it will be a useful reference point to study dynamics and generalization. Concerning the former, our observations suggest that learning near threshold may occur by avalanches, that is, by abrupt changes in the set of data that are correctly classified. In practice, networks are generally trained in the overparametrized regime $N\gg N^{*}$. It would be interesting to investigate whether the learning dynamics, at intermediate times where many data are not fitted yet, resembles the dynamics near threshold and displays bursts of changes in the constraints. Concerning the latter, we have studied the effect of jamming on generalization since this article was first written, as appears in Spigler et al. (2018); Geiger et al. (2019).

References

Appendix A Effective number of degrees of freedom

Due to several effects discussed in the main text, the function $f(\mathbf{x};\mathbf{W})$ can effectively depend on less variables that the number of parameters, and thus reduce the dimension of the space spanned by the gradients $\nabla_{\mathbf{W}}f(\mathbf{x};\mathbf{W})$ that enters in the theory. For instance, there could be symmetries that reduce the number of effective degrees of freedom (e.g. each ReLU activation function has one of such symmetries, since one can rescale inputs and outputs in such a way that the post-activation is left invariant); another reason could be that a neuron might never activate for all the training data, thus effectively reducing the number of neurons in the network; furthermore, we expect that the network’s true dimension would also be reduced if its architecture presents some bottlenecks, is poorly designed or poorly initialized. For example if all biases are too negative on the neurons of one layer in the Relu case, the network does not transmit any signals, leading to $N=1$ and to the possible absence of unstable directions even if the number of parameters is very large.

It is tempting to define the effective dimension by considering the dimension of the space spanned by $\nabla_{\mathbf{W}}f(\mathbf{x_{\mu}};\mathbf{W})$ as $\mu$ varies. This definition is not practical for small number of samples $P$ however, because this dimension would be bounded by $P$. We can overcome such a problem by considering a neighborhood of each point $\mathbf{x}_{\mu}$, where the network’s function and its gradient can be expanded in the pattern space:

Varying the pattern $\mu$ and the point $\mathbf{x}$ in the neighborhood of $\mathbf{x}_{\mu}$, we can build a family $M$ of vectors:

We then define the effective dimension $N$ as the dimension of $M$. Because of the linear structure of $M$, it is sufficient to consider, for each $\mu$, only $d+1$ values for $x$, e.g. $x-x_{\mu}=0,\hat{\mathbf{e}}_{1},\dots,\hat{\mathbf{e}}_{d}$, where $\hat{\mathbf{e}}_{n}$ is the unit vector along the direction $n$. The effective dimension is therefore

where the elements of the matrix $G$ are defined as

with $\alpha\equiv(\mu,n)$. The index $n$ ranges from to $d$, and $\hat{\mathbf{e}}_{0}\equiv 0$.

We consider $\mathcal{H}_{p}=-\sum_{\mu}y_{\mu}\rho\,(\Delta_{\mu})\,\hat{\mathcal{H}}_{\mu}$, where $\hat{\mathcal{H}}_{\mu}$ is the Hessian of the network function $f(\mathbf{x}_{\mu};\mathbf{W})$ and $\rho$ is the Relu function. We want to argue that the spectrum of $\mathcal{H}_{p}$ is symmetric in the limit of large $N$.

We do two main hypothesis: First, the trace of any finite power of $\mathcal{H}_{p}$ is self-averaging (concentrates) with respect to the average over the random data:

The first hypothesis is natural since $\hat{\mathcal{H}_{p}}$ is a very large random matrix, for which the density of eigenvalues is expected to become a non-fluctuating quantity. The second hypothesis is more tricky: it is natural to assume that the trace concentrates, however one also need to show that the sub-leading corrections to the self-averaging of the trace can be neglected.

Using these two hypothesis and the result, showed below, that

for all $n$ odds, one can conclude that all odds traces of $\hat{\mathcal{H}_{p}}$ are zero. This implies that the spectrum of $\hat{\mathcal{H}_{p}}$ is symmetric, more precisely that the fractions of negative and positive eigenvalues are equal.

where the indices $i_{1},\dots,i_{n}$ stand for synapses connecting a pair of neurons (i.e. each index is associated with a synaptic weight $W^{(j)}_{\alpha,\beta}$: we are not writing all the explicit indexes for the sake of clarity). The term of the hessian obtained when differentiating with respect to weights $W^{(j)}_{\alpha,\beta}$ and $W^{(k)}_{\gamma,\delta}$ reads

In fact, note that the sum in Eq.21 contains a weight per each layer in the network, with the exception of the two layers $j,k$ with respect to which we are deriving. This implies that any element of the hessian matrix where we have not differentiated with respect to the last layer ($j,k<L+1$) is an odd function of the last layer $W^{(L+1)}$, meaning that if $W^{(L+1)}\longrightarrow-W^{(L+1)}$, then the sign of all these Hessian elements is inverted as well.

If in the argument of the sum in Eq. (20) there is no index belonging to the last layer, then the whole term changes sign under the transformation $W^{(L+1)}\longrightarrow-W^{(L+1)}$. Suppose now that, on the contrary, there are $m$ terms with one index belonging to the last layer (we need not consider the case of two indices both belonging to the last layer because the corresponding term in the Hessian would be , as one can see in Eq. (21)). For each index equal to $L+1$ (the last layer), there are exactly two terms: $\hat{\mathcal{H}}^{\mu}_{j,L+1}\hat{\mathcal{H}}^{\mu}_{L+1,k}$ (for some indexes $j,k$). Since $j,k$ cannot be $L+1$ too, this implies that the number $m$ of terms with an index belonging to the last layer is always even. Consequently, when the sign of $W^{(L+1)}$ is reversed, the argument of the sum in Eq. (20) is multiplied by $(-1)^{n-m}$ (once for each term without an index belonging to the last layer), which is equal to $-1$ if $n$ is odd. The same symmetry can be used to show that a matrix made of an odd product of matrices $\hat{\mathcal{H}}_{\mu}$, such as $\hat{\mathcal{H}}_{\mu}\hat{\mathcal{H}}_{\mu^{\prime}}\hat{\mathcal{H}}_{\mu^{\prime\prime}}$, must also have a symmetric spectrum, concluding our argument.

Appendix C Density of pre-activations for ReLU activation functions

The densities of pre-activation (i.e. the value of the neurons before applying the activation function) is shown in Fig. 12 for random data. It contains a delta distribution in zero. The number $N_{c}$ of pre-activations equal to zero when feeding a network $L=5$ all its random dataset is $N_{c}\approx 0.21N$, corresponding to the number of directions in phase space where cusps are present in the loss function. For MNIST data we find $N_{c}\approx 0.19N$. By taking $L=2$ and random data we find $N_{c}\approx 0.25N$. In these directions, stability can be achieved even if the hessian would indicate an instability. For this reason, instead of $N_{-}$ in Equation 11 one should use $N/2-N_{c}\approx 0.25N$.