{'fraction_non_alphanumeric': 0.0377394953343272, 'fraction_numerical': 0.011484733035171994, 'mean_word_length': 4.222462514417532, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A jamming transition from under- to over-parametrization affects loss landscape and generalization

We argue that in fully-connected networks a phase transition delimits the over- and under-parametrized regimes where fitting can or cannot be achieved. Under some general conditions, we show that this transition is sharp for the hinge loss. In the whole over-parametrized regime, poor minima of the loss are not encountered during training since the number of constraints to satisfy is too small to hamper minimization. Our findings support a link between this transition and the generalization properties of the network: as we increase the number of parameters of a given model, starting from an under-parametrized network, we observe that the generalization error displays three phases: (i) initial decay, (ii) increase until the transition point --- where it displays a cusp --- and (iii) slow decay toward a constant for the rest of the over-parametrized regime. Thereby we identify the region where the classical phenomenon of over-fitting takes place, and the region where the model keeps improving, in line with previous empirical observations for modern neural networks.


Introduction
Despite the remarkable progress in designing [1,2] and training [3] 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 [4]. An analogy between deep networks and glasses has been proposed [5,6], 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 arXiv:1810.09665v5 [cs.LG] 18 Jun 2019 numerical and rigorous works [7,8,9,10,11] 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 [12,13,14] and of the learning dynamics [15,16] support that the landscape is characterized by an abundance of flat directions, even near its bottom, at odds with traditional glasses.

Jamming transition and supervised learning
In a previous article [17] 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 [18], 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 [19,20] 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 [19]. 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 [17]. 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 (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.

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 [21,22,23,24]. Yet another body of works focus on the effects of sheer size of a neural network [25,26,27].
Among previous studies, [28] stands out as a natural predecessor of our work. In [28] 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 * = 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 [29] for random forest models and simple neural networks. Further theoretical studies on regression showed a precise mathematical description of the cusp behaviour in [30,31,32], albeit on models that are practically somewhat further away from modern neural networks. Finally, in [33], 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 → ∞ 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 [17,18]. This will set the stage for the following thorough analysis of the phase transition and its role on generalization.

Set-up
We consider a binary classification problem, with a set of P distinct training data denoted {(x µ , y µ )} P µ=1 . The vector x µ is the input, which lives in a d-dimensional space, and y µ = ±1 is its label. We denote by f (x; W) the output of a fully-connected network corresponding to an input x, parametrized by W. We represent the network as in Fig. 2, and the output function is written recursively as where a (i) α are the preactivations. In our notation the set of parameters W includes, with a slight abuse of notation, both the weights W (i) α,β and the biases B (i) α . ρ(z) is the non-linear activation function, e.g. the ReLU ρ(z) = zθ(z) or the hyperbolic tangent ρ(z) = tanh(z). The parameters are learned by minimizing the quadratic hinge loss: where ∆ µ ≡ 1 − y µ f (x µ ; W) and m is the set of patterns with ∆ µ > 0 and contains N ∆ elements. These patterns describe unsatisfied constraints: they are either incorrectly classified or classified with an insufficient margin (whereas patterns with ∆ µ < 0 are learned with margin 1). We adopt this loss function since it makes the jamming transition simpler to analyze ‡, but this choice does not influence the performance of the network, as we have reported in [17].
We are interested in the transition between an over-parametrized phase where the network can satisfy all the constraints (L = 0) and an under-parametrized phase where some constraints remain unsatisfied (L > 0).

A note on the effective number of parameters
In our discussion the notion of effective number of degrees of freedom is important. In the space of functions going from the neighborhoods of the training set to real numbers, consider the manifold of functions f (x; W) obtained by varying W. We denote by N eff (W) the dimension of the tangent space of this manifold at W. In general we have N eff (W) ≤ N . Several reasons can make N eff (W) strictly smaller than N , including: • The signal does not propagate in the network, i.e. f (x; W) = C 1 for all x in the neighborhood of the training points x µ . In that case, the manifold is of dimension unity and N eff (W) = 1. This situation will occur for a poor initialization of the weights, for example if all biases are too negative on the neurons of one layer for ReLU activation function (see for instance [34]). It can also occur if the data x µ are chosen in an adversarial manner for a given choice of initial weights. For example, one can choose input patterns so as to not activate the first layer of neurons (which is possible if the number of such neurons is not too large). Poor transmission will be enhanced (and adversarial choices of data will be made simpler) if the architecture presents some bottlenecks. In the situation where N eff (W) = 1, it is very simple to obtain local minima of the loss at finite loss values, even when the model has many parameters. • The activation function is linear. Then, the output is an affine function of the input, leading to N eff ≤ d + 1. Dimension-dependent bounds will also exist if the activation function is polynomial (because the output function then is also restricted to be polynomial). • 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. • Some neurons are never active e.g. in the ReLU case, their associated weights do not contribute to N eff .
Thus there are N − N eff directions in parameter space that do not affect the function. These directions will lead to zero modes in the Hessian at any minimum of the loss. In what follows we consider stability with respect to the relevant N eff directions, which do affect the output function. Our results on the impossibility to get stuck in bad minima are expressed in terms of N eff . However, as reported in Appendix A.1, we find empirically that for a proper initialization of the weights and constantwidth fully connected networks, N eff ≈ N (the difference is small and equal to the number of hidden neurons, and only results from the symmetry associated with each ReLU neuron). Henceforth to simplify notations we will use the symbol N to represent the number of effective parameters.

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 L = 0, whereas if N is small after training L > 0. Imagine increasing N starting from a small value: at some N * the loss obtained after training approaches zero §, i.e. lim N →N * L = 0. We refer to this point as the jamming transition. A vanishing training loss implies that ∆ µ → 0 for each pattern µ = 1, . . . , P . As argued in [35], for each µ ∈ m the constraint ∆ µ ≈ 0 defines a manifold of dimension N − 1 . Satisfying N ∆ such equations thus generically leads to a manifold of solutions of dimension N − N ∆ ¶. 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 [37]), 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 ∇ is the gradient operator and ⊗ stands for tensor product). The first term H 0 is positive semi-definite: it is the sum of N ∆ rank-one matrices, thus rk Let us denote by E − the negative eigenspace + of H p and call N − its dimension. Stability then imposes that ker Hence, minima with positive training loss (and therefore the minimum found right above jamming) can only occur for: (The first inequality trivially follows from the fact that the N ∆ patterns belong to the training set of size P ). As reported in [17], we observe empirically that the spectrum of 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 Appendix A.2 we provide a non-rigorous argument supporting that in the case of ReLU activation functions and random data the spectrum of H p is indeed symmetric with lim N →∞ 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 → ∞ (for any fixed ratio P/N ) has a finite fraction C 0 = N − /N of negative eigenvalues for generic architectures and datasets. In [17] we observed C 0 = 1/2 for § 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. Related arguments were recently made for a quadratic loss [11]. In that case, we expect the landscape to be related to that of floppy spring networks, whose spectra are predicted in [36].
¶ Note that this argument implicitly assumes that the N ∆ constraints are independent. In disordered systems this assumption is generally correct, but it may break down if symmetries are present. + The negative eigenspace is the subspace spanned by the eigenvectors associated with negative eigenvalues.
the ReLU activation function as expected, for tanh activation functions at jamming and at the end of training we found C 0 ≈ 0.43. Thus, C 0 is not universal. Finally we assume that the spectrum of H p does not display a finite density of zero eigenvalues (once restricted to the space of parameters that affect the output, of dimension N eff , supposed here to be equal to N ). Note that this assumption breaks down if data can be identical with different labels, a case we exclude here * . Under this assumption we obtain from Equation (7) that local minima with positive training loss cannot be encountered if N > P/C 0 , implying in particular that: Non-smooth activation functions: With ReLU activation functions, the output function f (x; 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 ≡ βN presenting cusps near jamming, stability implies N ∆ > N − − N c and: implying in turn that: Numerically, we find that at jamming the fraction of directions along which there is a cusp is N c /N * ≡ β ∈ (0.21, 0.25) both for random data and images as reported in the Appendix A.3. Using C 0 = 1/2 for Relu, we obtain the bounds: N * ≤ 4P.
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 ) ≤ 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 ∆ /N of unsatisfied constraints per degree of freedom jumps discontinuously to a finite value satisfying C 0 ≤ N ∆ /N ≤ 1.
The complete list of results, including consequences of this analysis on the Hessian, is included in [17]. 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. * Indeed even in the over-parametrized case, if x i = x j but y i = y j then an exact cancellation of terms occurs in the sum defining Hp, which can then be zero while L > 0.

Location of the jamming transition
Here we present the numerical results on random data (uniformly distributed on a hypersphere and with random labels y µ = ±1) and on the MNIST dataset (partitioned into two groups according to the parity of the digits and with labels y µ = ±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 Appendix 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 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 ∆ /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 [20]. All this scatter is however gone when plotting N ∆ /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 Appendix 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 (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 [22] should diminish in the large-size limit. Beyond the jamming point, the accuracy keeps steadily improving as the number of parameters increases [25,26,27], although it does so quite slowly. We have provided a quantitative explanation for this phenomenon in [33]. In  . Appendix 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 [38,39]), 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 [20]. 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 (L = 0) from a region where the number of unsatisfied constraints is a fraction of the number of parameters, so 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 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 ≥ P/C 0 , implying in particular that N * ≤ 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 N * ; then (ii) it increases and culminates in a cusp at N ≈ 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. 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 ∇ W f (x µ ; W) as µ 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 x µ , where the network's function and its gradient can be expanded in the pattern space: Varying the pattern µ and the point x in the neighborhood of x µ , we can build a family M of vectors: We then define the effective dimension N eff as the dimension of M . Because of the linear structure of M , it is sufficient to consider, for each µ, only d + 1 values for x, e.g.
x − x µ = 0,ê 1 , . . . ,ê d , whereê 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 α ≡ (µ, n). The index n ranges from 0 to d, andê 0 ≡ 0.
In Fig. A1 we show the effective number of parameters N eff versus the total number of parameters N , in the case of a network with L = 3 layers trained on the first 10 PCA components of the MNIST dataset. There is no noticeable difference between the two quantities: the only reduction is due to the symmetries induced by the ReLU functions (there is one such symmetry per neuron. Indeed the ReLU function ρ(z) = zΘ(z) satisfies Λρ(z/Λ) ≡ ρ(z).) We observed the same results for random data.

Appendix A.2. sp(H p ) is symmetric for ReLu activation functions and random data
We consider H p = − µ y µ ρ (∆ µ )Ĥ µ , whereĤ µ is the Hessian of the network function f (x µ ; W) and ρ is the Relu function. We want to argue that the spectrum of H p is symmetric in the limit of large N .
We do two main hypothesis: First, the trace of any finite power of H p is selfaveraging (concentrates) with respect to the average over the random data:  Figure A1. Results with the MNIST dataset, keeping the first 10 PCA components. d = 10 and L = 3, varying P and h. Effective N eff vs total number of parameters N . N eff is always smaller than N because there is a symmetry per each ReLU-neuron in the network.

Second,
1 N µ1,···,µn y µ1 ρ(∆ µ1 ) · · · y µn ρ(∆ µn )tr(Ĥ µ1 · · ·Ĥ µn ) = The first hypothesis is natural sinceĤ 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 tr(Ĥ µ1 · · ·Ĥ µn ) = 0 (A. 8) for all n odds, one can conclude that all odds traces ofĤ p are zero. This implies that the spectrum ofĤ p is symmetric, more precisely that the fractions of negative and positive eigenvalues are equal.
In order to show that the statement (A.8) above holds, let us argue first that tr(Ĥ n µ ) = 0 for any odd n.
where the indices i 1 , . . . , i n stand for synapses connecting a pair of neurons (i.e. each index is associated with a synaptic weight W where we denoted with a the inputs in the nodes of the network. Our argument is based on a symmetry of the problem with random data: changing the sign of the weight of the last layer W (L+1) −→ −W (L+1) and changing the labels y µ −→ −y µ leaves the loss unchanged. We will show that this symmetry implies that tr(Ĥ n µ ) averaged over the random labels is zero for odd n.
In fact, note that the sum in Equation (A.10) 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) −→ −W (L+1) , then the sign of all these Hessian elements is inverted as well.
If in the argument of the sum in Equation (A.9) there is no index belonging to the last layer, then the whole term changes sign under the transformation W (L+1) −→ −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 0, as one can see in Equation (A.10)). For each index equal to L + 1 (the last layer), there are exactly two terms:Ĥ µ j,L+1Ĥ µ 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 (A.9) 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Ĥ µ , such asĤ µĤµ Ĥ µ , must also have a symmetric spectrum, concluding our argument.  Figure A2. Density of the pre-activations for each layers with L = 5 and random data, averaged over all the runs just above the jamming transition with that architecture. Black: distribution obtained over the training set. Blue: previously unseen random data (the two curves are on top of each other except for the delta in zero). The values indicate the mass of the peak in zero, which is only present when the training set is considered.

Appendix 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. A2 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 ≈ 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 ≈ 0.19N . By taking L = 2 and random data we find N c ≈ 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 ≈ 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 √ d, x µ ∈ S d , with random label y µ = ±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 [−σ, σ], where σ 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 · 10 5 steps and 10 −3 for 5 · 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. B1 we check that t = 10 6 is enough to converge. When using the hinge loss, we use an orthogonal initialization [42], no batch normalization and t = 2 · 10 6 steps of ADAM [43] 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 ∆ of unsatisfied constraints at the transition (Fig. 3B and inset), we consider three architectures, both with N ≈ 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). Fig 3A hinge We took d = h and trained for 2M steps. For some values of P ∈ (500, 60k), start at large h where we reach N ∆ = 0 and decrease h until N ∆ > 0.1N .

Details about Fig 3B
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 ∆ = 0, 74 have N ∆ > 0.4N . Among the 19 remaining ones, 14 of them have N ∆ 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 ∆ < 0.4N has been stoped too early. The same observation apply for the other depths.

Appendix B.2. Real data
The images in the MNIST dataset are gathered into two groups, with even and odd numbers and with labels y µ = ±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 × 28 = 784 (each picture is 28 × 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 ).