A jamming transition from under- to over-parametrization affects loss landscape and generalization

Introduction

Despite the remarkable progress in designing and training neural networks, there is still no general theory explaining their success, and their understanding remains mostly empirical. Central questions need to be clarified, such as what conditions need to be met in order to fit data properly, why the dynamics does not get stuck in spurious local minima, and how the depth of the network affects its loss landscape.

Complex physical systems with non-convex energy landscapes featuring an exponentially large number of local minima are called glasses . An analogy between deep networks and glasses has been proposed , in which the learning dynamics is expected to slow down and to get stuck in the highest minima of the loss. Yet, in the regime where the number of parameters is large (often considered in practice), several numerical and rigorous works 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 and of the learning dynamics support that the landscape is characterized by an abundance of flat directions, even near its bottom, at odds with traditional glasses.

In a previous article we have introduced an analogy between supervised learning with deep neural networks and a class of glassy systems, namely random dense packings of repulsive particles. It generalized a previous seminal analogy established between the loss landscape of the perceptron (the simplest network without hidden neurons) and the energy landscape of spherical particles , and specified the universality class to which deep learning corresponds to. The critical behavior of these granular systems, although very general, is of easier understanding when we consider particles that interact only within a finite range: upon increasing their density, such systems undergo a critical jamming transition when there is no longer space to accommodate all the particles without them touching one another. Before the transition the energy is zero, and after it increases with the density. The inclusion of longer-range interactions blurs the transition but its effects are still affecting the energy landscape . Deep networks behave similarly when we look at the training loss, and, again, a clear criticality emerges when considering a “finite-range” loss function — the hinge loss: when the number $P$ of training points is small enough, the network is able to learn the whole training set and reaches zero training loss, and upon increasing the dataset size we find a critical “jamming” point where perfect training does not occur and learning gets stuck in a positive minimum of the the training loss.

For the full analogy we point to the aforementioned paper . In the present work we first review the arguments that show that the existence of the jamming transition, studied in the $(N,P)$ plane where $N$ is the number of degrees of freedom of the network (informally, its size) and $P$ is the size of the training set. As it turns out, there is a critical line $N^{\star}(P)$ (whose exact location can depend on the chosen dynamics) delimiting two phases, one where the learning reaches zero training loss, and one where it gets stuck in a minimum with finite loss — see Fig. 1. We present some numerical results that characterize the different phases, both for random data and for the MNIST dataset, using fully-connected networks with ReLU activation functions. Then, we show novel data that illustrate that this transition affects the most crucial aspect of learning, namely the generalization error. We observe that generalization properties are strongly affected by the proximity to the jamming transition: for a gradient descent dynamics, in the under-parameterized phase before jamming (large $P$ or small $N$) the generalization error is increasing; at the transition it displays a cusp; after jamming, in the over-parametrized phase, the error decreases monotonically. If early stopping is used, the cusp disappears, implying that the jamming transition is precisely the point where over-fitting is very strong.

2 Generalization versus over-fitting

The puzzle regarding the good generalization properties of neural networks despite their large size has been the topic of study for several other works. Some of them focus on the effects of various ways of regularizing the network, thereby effectively reducing the dimension . Yet another body of works focus on the effects of sheer size of a neural network .

Among previous studies, stands out as a natural predecessor of our work. In the authors present one of the first empirical observations of the cusp in test error of a non-linear model, a behaviour that is reminiscent of the perceptron. There, the training dynamics is run on a two-layer student network whose training data is provided by a teacher network with a similar architecture. The authors have observed (i) a cusp in generalization error, and (ii) monotonic decay in the test error when early stopping is used.

In this work, we show that the cusp in generalization corresponds to a phase transition where the number of unsatisfied constraints suddenly drops to zero as $N$ increases. We quantify how the location of the transition $N^{*}(P)$ depends on $P$ for both random data and natural images, and find that the data structure significantly affects $N^{*}(P)$. Our analysis makes it clear that $N^{*}\neq P$, an assumption sometimes made in previous studies. Overall, it relates the cusp in generalization to a well-known body of literature in physics associated with the “jamming” transition.

Since the initial preparation of the present work, the field has progressed quickly within a matter of months. The described cusp in the generalization error has been observed empirically in for random forest models and simple neural networks. Further theoretical studies on regression showed a precise mathematical description of the cusp behaviour in , albeit on models that are practically somewhat further away from modern neural networks. Finally, in , our subsequent work, we develop a quantitative theory for (i) the cusp which is associated to the divergence of the norm of the output function at the critical point of the phase transition and (ii) the asymptotic behaviour of the generalization error as $N\rightarrow\infty$ which is associated with the reduced fluctuations of the output function. That work also shows that after ensemble averaging several networks, performance is optimal near the jamming the threshold, emphasizing the practical importance of this transition.

Theoretical framework

In this section we recall in detail the analogy between jamming and supervised learning for deep neural networks . This will set the stage for the following thorough analysis of the phase transition and its role on generalization.

We consider a binary classification problem, with a set of $P$ distinct training data denoted $\{(\mathbf{x}_{\mu},y_{\mu})\}_{\mu=1}^{P}$. The vector $\mathbf{x}_{\mu}$ is the input, which lives in a $d$-dimensional space, and $y_{\mu}=\pm 1$ is its label. We denote by $f(\mathbf{x};\mathbf{W})$ the output of a fully-connected network corresponding to an input $\mathbf{x}$, parametrized by $\mathbf{W}$. We represent the network as in Fig. 2, and the output function is written recursively as

where $a^{(i)}_{\alpha}$ are the preactivations. In our notation the set of parameters $\mathbf{W}$ includes, with a slight abuse of notation, both the weights $W^{(i)}_{\alpha,\beta}$ and the biases $B^{(i)}_{\alpha}$. $\rho(z)$ is the non-linear activation function, e.g. the ReLU $\rho(z)=z\theta(z)$ or the hyperbolic tangent $\rho(z)=\tanh(z)$. The parameters are learned by minimizing the quadratic hinge loss:

where $\Delta_{\mu}\equiv 1-y_{\mu}f(\mathbf{x}_{\mu};\mathbf{W})$ and $m$ is the set of patterns with $\Delta_{\mu}>0$ and contains $N_{\Delta}$ elements. These patterns describe unsatisfied constraints: they are either incorrectly classified or classified with an insufficient margin (whereas patterns with $\Delta_{\mu}<0$ are learned with margin 1). We adopt this loss function since it makes the jamming transition simpler to analyzeThe often used cross-entropy loss function also displays a transition where all data are well-fitted. However, in the over-parametrized regime the dynamic never stops, as the total loss vanishes only if the output and therefore the weights diverge. Imposing a time cut-off is done in practice, but it blurs the criticality near jamming, as exemplified below with the early stopping procedure., but this choice does not influence the performance of the network, as we have reported in .

We are interested in the transition between an over-parametrized phase where the network can satisfy all the constraints ($\mathcal{L}=0$) and an under-parametrized phase where some constraints remain unsatisfied ($\mathcal{L}>0$).

2 A note on the effective number of parameters

Symmetries are present in the network, e.g. the scale symmetry in ReLU networks: since the ReLU function is homogeneous, multiplying the weights of a layer by some factor and dividing the weights in the next layer by the same factor leaves the output function invariant. It will reduce one degrees of freedom per node.

3 Constraints on the stability of minima

In this section we show that the existence of a minimum at a vanishingly small training loss (i.e. approaching jamming) is enough to derive an upper bound for the transition in the $(N,P)$ plane.

Let us suppose (and justify later) that, for a fixed number of data $P$ and with proper initialization of weights, if $N$ is large enough then gradient descent leads to ${\cal L}=0$, whereas if $N$ is small after training ${\cal L}>0$. Imagine increasing $N$ starting from a small value: at some $N^{*}$ the loss obtained after training approaches 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 an instance of finite size effects. , i.e. $\lim_{N\rightarrow N^{*}}{\cal L}=0$. We refer to this point as the jamming transition. A vanishing training loss implies that $\Delta_{\mu}\rightarrow 0$ for each pattern $\mu=1,\dots,P$. As argued in , for each $\mu\in m$ the constraint $\Delta_{\mu}\approx 0$ defines a manifold of dimension $N-1$Related arguments were recently made for a quadratic loss . In that case, we expect the landscape to be related to that of floppy spring networks, whose spectra are predicted in .. Satisfying $N_{\Delta}$ such equations thus generically leads to a manifold of solutions of dimension $N-N_{\Delta}$Note that this argument implicitly assumes that the $N_{\Delta}$ constraints are independent. In disordered systems this assumption is generally correct, but it may break down if symmetries are present.. Imposing that a solution exists implies that at jamming:

Smooth activation function: An opposite bound can be obtained by considerations of stability (as was done for the jamming of repulsive spheres in ), by imposing that in a stable minimum the Hessian must be positive definite if the output function is smooth, as it must be the case if the activation function is smooth (see below for the situation where the function displays cusps, as occurs for ReLU neurons). The Hessian matrix, that is the matrix of second derivatives, is

(Here $\nabla$ is the gradient operator and $\otimes$ stands for tensor product). The first term ${\cal H}_{0}$ is positive semi-definite: it is the sum of $N_{\Delta}$ rank-one matrices, 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}$.

(The first inequality trivially follows from the fact that the $N_{\Delta}$ patterns belong to the training set of size $P$). As reported in , we observe empirically that the spectrum of $\mathcal{H}_{p}$ is statistically symmetric in the cases that we consider in the present work, i.e. for ReLU activation function, both for MNIST and random data, both at initialization and at the end of training. In A.2 we provide a non-rigorous argument supporting that in the case of ReLU activation functions and random data the spectrum of $\mathcal{H}_{p}$ is indeed symmetric with $\lim_{N\rightarrow\infty}N_{-}/N_{+}=1$ independently of depth, where $N_{+}$ is the number of positive eigenvalue. We conjecture that in general the limiting spectrum of $H_{p}$ as $N,P\rightarrow\infty$ (for any fixed ratio $P/N$) has a finite fraction $C_{0}=N_{-}/N$ of negative eigenvalues for generic architectures and datasets. In we observed $C_{0}=1/2$ for the ReLU activation function as expected, for tanh activation functions at jamming and at the end of training we found $C_{0}\approx 0.43$. Thus, $C_{0}$ is not universal.

Non-smooth activation functions: With ReLU activation functions, the output function $f(\mathbf{x};\mathbf{W})$ is not smooth and presents cusps, so that the Hessian needs not be positive definite for stability. A minimum can lie on a point where the second derivative is not defined along some directions (because of the cusp), and we say that the cusp stabilizes those directions. Equation (7) needs to be modified accordingly: introducing the number of directions $N_{c}\equiv\beta N$ presenting cusps near jamming, stability implies $N_{\Delta}>N_{-}-N_{c}$ and:

Numerically, we find that at jamming the fraction of directions along which there is a cusp is $N_{c}/N^{*}\equiv\beta\in(0.21,0.25)$ both for random data and images as reported in the A.3. Using $C_{0}=1/2$ for Relu, we obtain the bounds:

Main results: Overall, our analysis supports that for smooth activation functions there exists a constant $C_{0}$ such that:

there is a transition for $N^{*}(P)\leq P/C_{0}$ below which the training loss converges to some non-zero value (under-parametrized phase) and above which it becomes null (over-parametrized phase).

At the transition, the fraction $N_{\Delta}/N$ of unsatisfied constraints per degree of freedom jumps discontinuously to a finite value satisfying $C_{0}\leq N_{\Delta}/N\leq 1$.

The complete list of results, including consequences of this analysis on the Hessian, is included in .

For $Relu$ activation function, $C_{0}=1/2$ but the analysis is complicated by the presence of cusps. The jamming transition is still sharp, i.e. characterized by a discontinuous jump in constraints as specified by Eq.11.

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

Location of the jamming transition

Here we present the numerical results on random data (uniformly distributed on a hypersphere and with random labels $y_{\mu}=\pm 1$) and on the MNIST dataset (partitioned into two groups according to the parity of the digits and with labels $y_{\mu}=\pm 1$). With MNIST, in order not to have most of the weights in the first layer, we reduce the actual input size by retaining only the first $d=10$ principal components that carry the most variance (this hardly diminishes the performance for such a task). Further description of the protocols is in B.

In Fig. 3A,C we show the location of boundary $N^{*}$ versus the number of samples $P$. $N^{*}$ is estimated numerically for each $P$ by starting from a large value of $N$ and progressively decreasing it until $L>0$ at the end of training. Varying input dimension, depth and loss function (cross entropy or hinge) has little effect on the transition. This result indicates that in the present setup the ability of fully-connected networks to fit random data does not depend crucially on depth. Fig. 3C shows also a comparison of random data with MNIST. A difference between random data and images is that the minimum number of parameters $N^{*}$ needed to fit the real data is significantly smaller and grows less fast as $P$ increases — for $P\gg 1$, $N^{*}(P)$ could be sub-linear or even tend to a finite asymptote: how the data structure affects $N^{*}(P)$ is an important questions for future studies.

From the analysis of Section 2, the number of constraints per parameter $N_{\Delta}/N$ is expected to jump discontinuously at the transition. This is shown in the insets of Fig. 3B,D. The scatter in these plots presumably reflects finite size effects known to occur near the jamming transition of particles . All this scatter is however gone when plotting $N_{\Delta}/N$ as a function of the loss itself, as shown in the main panels of Fig. 3B,D.

Generalization at and beyond jamming

In Fig. 4A we show the evolution of the generalization error for networks at four different locations in the $(N,P)$ plane. The networks are trained on MNIST at fixed $P=50k$, and at different values $N$, both above, at and below jamming. Training is run for a fixed number of steps of vanilla gradient descent (the simulation details are in B). The profile of these curves is typical of most learning problems (if one does not recur to early stopping): notice that the point of minimum generalization error happens before the end of training. The increase of test error at late times is referred to “over-fitting” in the field. Very interestingly, it is clear from this figure that at small and large $N$, over-fitting is a weak effect, which however becomes very significant at intermediate $N$.

To study this effect, we systematically vary $N$ at fixed $P$. In Fig. 4B the solid curve shows the generalization error against the network size $N$ for three different values of $P$ (we sampled subsets of MNIST). The dashed curve represents the value of the smallest error obtained during training, at prior time-steps (extracted from the profiles shown in Fig. 4A). The former displays a cusp at the transition point, as one can see clearly after rescaling the $N$-axis of each curve by the corresponding value of $N^{\star}(P)$. Strong over-fitting, corresponding to the difference between the solid and dashed lines, takes place only in the vicinity of the critical jamming transition (Fig. 4B-C). We thus posit that at fixed $P$, the benefit of early stopping should diminish in the large-size limit. Beyond the jamming point, the accuracy keeps steadily improving as the number of parameters increases , although it does so quite slowly. We have provided a quantitative explanation for this phenomenon in . In B.2 we have verified that the overall trends showed in Fig. 4 qualitatively hold also for other depths.

Notice that although the cusp has been found also in shallow networks (in particular the perceptron ), their behavior is at odds with what we observe: for the perceptron, test error asymptotically increases with $N$.

Conclusions

Understanding the effect of over-parametrization on the behavior of deep neural networks is a central problem in machine learning. In this work, by focusing on the hinge loss, we recast the minimization of the loss function as a constraint-satisfaction problem with continuous degrees of freedom. A similar approach was used in the field of interacting particles, which display a sharp jamming transition affecting the landscape if the interaction is chosen to be finite range . Following the analogy we were able to predict a sharp transition as the number of network parameters is varied, separating a region in the $(P,N)$ plane where a global minimum can be found ($\mathcal{L}=0$) from a region where the number of unsatisfied constraints is a fraction of the number of parameters, so $\mathcal{L}>0$. These results also shed light on several aspects of deep learning:

Not getting stuck in local minima: In the over-parametrized regime, the dynamics does not get stuck in local minima at finite loss value because the number of constraints to satisfy is too small to hamper minimization. It follows from our assumptions on the negative eigenspace of the matrix $\mathcal{H}_{p}$ that in this regime the landscape is flat and local minima do not exist (assuming that the number of effective parameters that affect the output function is $N$). For a smooth activation function we predict that one cannot get stuck in a bad minimum for $N\geq P/C_{0}$, implying in particular that $N^{*}\leq P/C_{0}$ where $C_{0}$ is a constant. We obtain a less demanding bound for ReLU activation functions due to the presence of cusps in the landscape, a situation for which we expect $C_{0}=1/2$. In practice, for random data $N^{*}(P)$ scale linearly with $P$ (in this sense, the bound is tight). By contrast, for structured data $N^{*}(P)$ appears to scale sub-linearly with $P$. Predicting the curve $N^{*}(P)$ remains a challenge for the future.

Reference point for fitting and generalization: There exists a critical curve $N^{*}(P)$ on the $N$-$P$ plane above which the global minima of the landscape become accessible. The curve also appears to be linked to the generalization potential of the model. We show that in the cases that we considered, (i) the generalization error decreases when $N\ll N^{*}$; then (ii) it increases and culminates in a cusp at $N\approx N^{*}$ that is erased by early stopping, most useful in this region; finally, (iii) in the over-parametrized phase, it monotonically decreases, although very slowly.

We thank Marco Baity-Jesi, Carolina Brito, Chiara Cammarota, Taco S. Cohen, Silvio Franz, Yann LeCun, Florent Krzakala, Riccardo Ravasio, Andrew Saxe, Pierfrancesco Urbani and Lenka Zdeborova for helpful discussions. This work was partially supported by the grant from the Simons Foundation (#454935 Giulio Biroli, #454953 Matthieu Wyart). M.W. thanks the Swiss National Science Foundation for support under Grant No. 200021-165509. The manuscript , which appeared at the same time as ours, shows that the critical properties of the jamming transition found for the non-convex perceptron hold more generally in some shallow networks. This universality is an intriguing result. Understanding the connection with our findingsis certainly worth future studies.

References

Appendix A Network properties

In the following, we analyze numerically the networks properties that were used in the previous analysis. This provides a numerical confirmation of our arguments, and an in depth characterization of the networks.

Due to several effects discussed above, 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:

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 Equation (22) 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 Equation (21) 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 Equation (22)). 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 Equation (21) 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.

A.3 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. 6 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 (7) one should use $N/2-N_{c}\approx 0.25N$.

Appendix B Parameters used in simulations

The dataset is composed of $P$ points 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 networks are fully connected, and have an input layer of size $d$ and $L$ layers with $h$ neurons each, culminating in a final layer of size $1$. To find the transition we proceed as follows: we build a network with a number of parameters $N$ large enough for it to be able to fit the whole dataset without errors. Next, we decrease 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^{*}$. As initial conditions 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.

When using the cross entropy, the system 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 ($10^{6}$ steps in total); the batch size is set to $\min(P/2,1024)$, and batch normalization is used. We do not use any explicit regularization in training the networks. In Fig. 7 we check that $t=10^{6}$ is enough to converge.

When using the hinge loss, we use an orthogonal initialization , no batch normalization and $t=2\cdot 10^{6}$ steps of ADAM with batch size $P$ and a learning rate starting at $10^{-4}$. In the experiments of section 3 (not for the experiments of section 4), we progressively divided the learning rate by $10$ every 250k steps. Also in this case we do not use any explicit regularization in training the networks.

To observe the discontinuous jump in the number $N_{\Delta}$ of unsatisfied constraints at the transition (Fig. 3B and inset), we consider three 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 and minimizing for $10^{7}$ steps (a better minimization is needed to improve the precision close to the transition).

We trained networks of depth 2,3,5 with $d=h=$ 62, 51, 40 respectively for 10M steps. For $L=3$ ($d=51$, $h=51$) we ran 128 training varying $P$ from 21991 to 25918. For the value of $N$ we take $7854$ that correspond to the number of parameters minus the number of neurons, per neuron there is a degree of freedom lost in a symmetry induced by the homogeneity of the ReLU function. 37 of the runs have $N_{\Delta}=0$, 74 have $N_{\Delta}>0.4N$. Among the 19 remaining ones, 14 of them have $N_{\Delta}$ between 1 and 4, we think that these runs encounter numerical precision issues, we observed that using 32 bit precision accentuate this issue. We think that the 5 left with $4<N_{\Delta}<0.4N$ has been stoped too early. The same observation apply for the other depths.

B.2 Real data

The images in the MNIST dataset are gathered into two groups, with even and odd numbers and 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 each, that in the end result in a single scalar output. The loss function used is always the hinge loss.

If we kept the original input size of $28\times 28=784$ (each picture is $28\times 28$ pixels) 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 opt for a reduction of the input size. We perform a principal component analysis (PCA) on the whole dataset and we identify the 10 dimensions that carry the most variance on the whole dataset; then we use 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 (which we find to be larger than $90\%$ when using all the data and large $N$).

We trained a network of $L=5$, $d=10$, $h=30$ for 3M steps. With $P$ varying from 31k to 68k (using trainset and testset of MNIST).