{'fraction_non_alphanumeric': 0.035174704148931195, 'fraction_numerical': 0.01969222133869685, 'mean_word_length': 4.5117298272750705, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Subaging in underparametrized Deep Neural Networks

We consider a simple classification problem to show that the dynamics of finite-width Deep Neural Networks in the underparametrized regime gives rise to effects similar to those associated with glassy systems, namely a slow evolution of the loss function and aging. Remarkably, the aging is sublinear in the waiting time (subaging) and the power-law exponent characterizing it is robust to different architectures under the constraint of a constant total number of parameters. Our results are maintained in the more complex scenario of the MNIST database. We find that for this database there is a unique exponent ruling the subaging behavior in the whole phase.


Introduction
Understanding the (supervised) learning of Deep Neural Networks (DNNs) continues to be one of the biggest challenges from a theoretical point of view [Carleo et al., 2019]. Despite the existence of specific cases in which a significant understanding has been reached, such as the perceptron [Seung et al., 1992, Franz andParisi, 2016], the single hidden layer in the mean field limit [Mei et al., 2018] and the infinite-width limit [Jacot et al., 2018, Geiger et al., 2021 (see also [Carleo et al., 2019]), in the more realistic case of finite-width and a finite number of nodes, a deeper understanding is still demanded. Steps toward this direction require the analysis of the interplay between the landscape of a generic Deep Neural Network and the (learning) dynamics.
The landscape refers to the loss (hyper)-surface generated by the parameters of a network (weights and biases) given a dataset, and its generic properties are crucial in database together with the relevant information regarding the landscape structure; moreover, we discuss the dynamical functions to be considered. In section 3 we provide our main results and discuss them. Finally, in section 4 we draw the conclusions of our analysis.

Landscape structure and dynamics
We consider a simple classification problem as originally introduced in [Ballard et al., 2016]. This refers to the database LJAT19 [Verpoort et al., 2020a], which consists of data with three coordinates (input) that describe a triatomic system and a label (output) that corresponds to one out of four possible stable configurations that are reached using standard molecular dynamics. In [Verpoort et al., 2020b], the difficulty of the classification problem was increased by reducing to two the input coordinates, with the implication that the maximum possible accuracy during training, defined as the ratio of correctly classified samples over the total number of samples, is around 0.8.
In [Verpoort et al., 2020b], the landscape of small fully-connected DNNs is studied as the depth (number of hidden layers, H = 1, 2, 3) and dataset size (N ) is varied, while keeping the number of parameters low and approximately constant (see details in appendix 5.1). The idea behind this procedure is to keep the dimension of the parameter space approximately equal and, in this way, make a meaningful comparison of architectures.
The loss function to be considered is the cross entropy with softmax outputs and a L2 regularization, normalized by the number of data. This is: where D refers to the dataset, θ = (θ 1 , . . . , θ P ) is the vector of all the weights (including biases), a (k) n refers to the softmax output for the true k-th class associated with the n-th data and λ is the regularization constant.
Specified the dataset D and the architecture of the network, the set of local minima {θ * } of the function L D (θ) can be extracted with a basinhopping algorithm [Wales and Bogdan, 2006], and the transition states (saddles of index-1) connecting those minima are searched with Discrete Path Sampling [Wales and Bogdan, 2006]. For the LJAT19 database, this characterization was done in reference [Verpoort et al., 2020b] and the results unveil how the structure of the landscape is database size dependent. In particular they show that for a small number of data N , the multilayer architectures present many competing low-lying minima with similar loss values, whereas the single hidden layer landscape is funneled with a small number of local minima; on the other hand, in the regime of high N , all the landscapes turn out to be funneled with a small number of local (shallow) minima (see [Wales et al., 1998] for a visualization of funneled landscapes).
In this manuscript, unless otherwise explicitly stated, we show results for the small dataset regime N = 1000. To help the reader to have an idea of the complexity of the landscape, we report here the order of the number of local minima found (n min ) for the considered dataset (cf. Table 1 in [Verpoort et al., 2020b]): single-hidden layer, n min ∼ 10, two-hidden layer, n min ∼ 10 3 , three-hidden layer, n min ∼ 10 4 . In all the cases, the loss value of the deepest found minimum is around 0.5.
Once the landscape is defined, the next step is to perform the optimization dynamics that allows its exploration. We implement the standard Stochastic Gradient Descent (SGD) with a fixed learning rate and fixed batch size (details in appendix 5.1. We remark that the validity of our results is independent from the chosen values). Then we keep track of the evolution of the loss function during training together with the mean square displacement (MSD). This is defined between times t w and t w + t as [Baity-Jesi et al., 2019] with t w the so-called "waiting time" and P the total number of parameters. The function (2) characterizes the overall motion between two times in the parameter space.
In systems with glassy dynamics, a way of characterizing the slowness of the evolution and in particular the aging is via the scaling of a characteristic time t * with t w [Berthier, 2000, Arceri et al., 2020. As a matter of fact, t * may be extracted by fixing a value of ∆ = ∆ * and then solving for t as a function of the waiting time. In particular, the scaling t * ∼ t µ w with µ a positive number is quite generic in glassy models [Bouchaud, 1999, Rinn et al., 2000, Berthier, 2000. Case µ = 1 is denoted as simple aging, whereas subaging and superaging denote the cases µ < 1 and µ > 1, respectively. We obtain numerically t * (t w ) for the architectures/datasets discussed in this work and show that a subaging scenario is compatible with the data.
Finally, we characterize the local geometry of the landscape along the evolution by looking at the distribution of the Hessian eigenvalues of the loss function (as proposed in [Sagun et al., 2017] As the stochastic evolution is driven towards local minima, the number of negative eigenvalues is highly reduced. However, because of the stochastic evolution, the dynamics do not get trapped in the local minima and negative curvatures prevail to some extent [Wei and Schwab, 2019]. In the next section, we illustrate that the change in the distribution helps to gain intuition about the training.

Results
We trained networks with one, two and three hidden layers with approximately P ≈ 70 each using the database LJAT19 (details in appendix 5.1). We show the train loss function and accuracy against the number of steps in figure 1; a step refers to one iteration across a single minibatch. As a general observation, we see that for the three The solid line corresponds to one hidden layer, the discontinuous line to two hidden layers, and the dotted line to three hidden layers. We consider 15000 epochs with a batch size of 100 that leads to 150000 steps.
architectures both the loss and the accuracy reach similar values and they enter into the regime where the lowest local minima lie (L ≈ 0.5) at t ≈ 10 3 . Regarding the specific form of the loss evolution, which exhibits an intermediate plateau especially for the multilayer cases, we point out that the corresponding landscapes exhibit many local minima at values of L ≈ 0.8 + , which explains this transient behavior. According to Figure 1 and the information on the landscapes, we separate the evolution of the loss into two phases. In the first one, learning is relatively fast, despite the eventual existence of intermediate plateaus. In the second phase, learning is slow, and the system explores the bottom of the landscape composed of multiple local minima at around the same depth. This division is consistent with other descriptions, as for instance references [Frankle et al., 2020, Feng andTu, 2021]. Following [Feng and Tu, 2021], in the first phase, learning is fast, and the difference between consecutive gradient updates is aligned; whereas the second phase is an exploration phase, in which the difference between updates is distributed in all directions (see also [Kunin et al., 2021]).
We characterize further the slow learning phase with the mean square displacement. In this phase, there is a non-translational invariant time dependence of the MSD, as Fig. 2 reveals. According to this figure we can say that the overall motion becomes slower with the "age" t w of the system as it happens for structural glasses [Scalliet andBerthier, 2019, Arceri et al., 2020]. This is the situation regardless of the architecture. To substantiate this, we compute numerically the function t * (t w ) obtained by fixing ∆ * = 10 −2 and the result is shown in Fig. 2d. Remarkably, the asymptotic scaling is independent of the architecture and reveals a subaging behavior. (2)) for the same data shown in Fig. 1 evaluated at waiting times t w > 10 3 , where the system has entered (or is close to enter) in the slow learning phase. Figure (a) corresponds to the architecture with one hidden layer, (b) to two hidden layers and (c) to three hidden layers. Figure (d) shows the characteristic time t * for each t w ; for this plot, each curve is an average over 100 independent runs. The scaling with t w shown as a dot-dashed line is set by hand with the purpose of illustration.
Few comments on this result are in order. Firstly, even though the value for ∆ * was chosen arbitrarily, Figure 2 (a)-(c) shows that there exists an interval ∆ * ∈ [∆ min , ∆ max ] such that the scaling exponent t * ∼ t µ w is independent or weakly dependent from ∆ * . This can be read as the vertical range in which the curves for different waiting times remain parallel. Secondly, we see that for long relative times t, the different MSD curves tend to approach together; this indicates a change in the nature of aging for long times. This regime has not analogue in standard models of glassy systems. Its full characterization is left for future work.
We also comment on the crossover around t w ≈ 10 3.6 in Figure 2d. This time is related to the overcoming of the regions with high loss values that lead to intermediate plateaus in the evolution of the loss function (Fig. 1a) and is characteristic of the LJAT19 database.
To show that our observations for the dynamics are quite generic in the underparametrized regime, we repeat the characterization for a minimal architecture (a chain) that contains three hidden layers and one neuron per layer. The landscape turns out to be rather simple, with a single pair of minima found with the basin-hopping algorithm * . Nevertheless, the training dynamics exhibit similar features to the higherdimensional parameter space as shown in Figure 3. Additionally, we consider a more complex scenario and repeat the previous analysis using the MNIST database. We fix the ratio α = P/N ≈ 0.07 as investigated for the LJAT19 and use this as a constraint for different architectures with H = 1, 2, 3 hidden layers. As dataset we use the complete database, so that N = 60000. Therefore, P ≈ 4200 (details in appendix 5.1). Although for this scenario we do not have a * We remark that the search is not exhaustive and new minima may appear with a longer run. We ran the GMNIN software for a time period of one day. characterization of the structure of the landscapes, we observe the same generic features in the evolution of the loss function, the MSD and the subaging behavior as for the LJAT19 database (see figures 4 and 5).
As a further step to characterize the underparametrized phase for the MNIST database, we study the scaling of the characteristic time with the waiting time for different values of α for a fixed depth H = 3. . The data together with the associated exponent µ(α) is shown in Figure 6. We observe that the data is consistent with a unique exponent across the phase, that we found to be µ = 0.709 ± 0.017.
Finally, and in order to gain a better intuition on the interplay between landscape and dynamics, we explore the distribution of eigenvalues of the Hessian of the loss function during the training. This gives a notion of the local geometry of the landscape around the point where they are calculated [Sagun et al., 2017, Wei andSchwab, 2019].
In practice, we consider that a network belongs to the underparametrized phase as long as α 1 and the accuracy reached after t = 10 6 steps is less than 0.99. The width of the networks is varied from 7 to 17 neurons, in steps of 2.    Figure (a) corresponds to the mean t * over an ensemble of 100 independent networks for each α. Figure (b) shows the corresponding exponent t * ∼ t µ w extracted via least squares. The horizontal dashed line corresponds to the mean value, µ = 0.709.
In figure 7 we show boxplots of the distribution for different training steps for a single architecture and the two different databases. At the start of the training process, we see that the eigenvalues are distributed in a wide range of values. In contrast, in the end, the majority of the eigenvalues seem to be spread around zero, with some outliers being mainly positive. The number of negative eigenvalues decreases as training advances, but they do not disappear entirely. which is naturally a consequence of the stochasticity of the equations of motion (SGD). † †

Conclusions
Following the original observation in , we show that the interplay between the structure of the landscape and the stochastic gradient descent in the underparametrized regime of a Deep Neural Network leads to glassy features. This glassiness is manifested in the slow evolution (after a transient) of the loss function and the non time-translation invariant mean square displacement. In fact, for the cases considered we presented evidence of subaging, in the sense that the characteristic time grows more slowly than the age of the system. Furthermore, we show that the depth of the network plays a marginal role by comparing architectures with roughly the same number of parameters. In addition, along the lines of [Feng andTu, 2021, Frankle et al., 2020] we identify two phases of learning, and associate the dynamics of the last (slow) phase with the exploration of the † † We remark that we use a hyperbolic tangent as the activation function for the hidden layers (see methods in appendix 5.1). Therefore we do not expect cusps in the associated landscapes, and the negative eigenvalues can be attributed entirely to the dynamics. This would not be the case with other activation functions, such as ReLU (for a discussion we refer the reader to ).
bottom of the landscape where the deepest local minima lie. Additionally, we state that the slow phase may be separated into two regimes according to the behavior of the MSD. The first regime exhibits subaging as characterized by a timescale that grows sublinearly with the waiting time. The second regime shows also a watiting-time dependence of the MSD with a longer associated timescale. It is left for a future work to characterize the transition between regimes by computing, for instance, a generalized diffusion coefficient associated with each phase.
As stated in references , Geiger et al., 2021, the transition from the underparametrized to the overparametrized regime is analogous to the solid-fluid change happening in amorphous materials with a variation of density (analogue of 1/α) and known as the jamming transition [Altieri, 2019]. We characterize  Figure 7: Boxplots for the eigenvalues of the hessian at different training steps for the architecture with two hidden layers for the LJAT19 dataset (top) and the MNIST dataset (bottom). The boxes extend from the first to the third quartile, while the whiskers cover 1.5 times the interquartile range. For scalling, some outliers are not shown. It is evident there is a greater amount of eigenvalues that are close to zero or positive at the end of the training process. the underparametrized (jammed) phase for the MNIST, and our results point out to a unique subaging exponent across the whole phase. This is consistent with observations for amorphous solids obtained via athermal simulations of viscous soft particles (see for instance [Head, 2009]). Our result is certainly intriguing and calls attention to a theoretical study that elaborates, firstly, on the origin of subaging in the underparametrized phase, secondly, on the nature of the scaling exponent, and thirdly, on the connection between the subaging behavior and other relevant properties of learning, such as the generalization error .
Finally, considering that our results are robust to different landscapes, it is worth exploring the robustness to different types of noises. In this direction, an interesting option is to consider the persistent-SGD [Mignacco and Urbani, 2021], which has the appealing property of being amenable to a dynamical mean field theory treatment [Mignacco et al., 2020]. Here it is easier to visualize there is a decrease in the amount of negative eigenvalues at the end of the training process.

Methods
We worked with the architectures described in reference [Verpoort et al., 2020b], that is, three different configurations with one, two and three hidden layers, each with 10, 5 and 4 neurons per layer respectively. The total number of parameters is 74, 69 and 72 respectively. As mentioned in the main text, we also consider the same dataset LJAT19 [Verpoort et al., 2020a]. As input, we take two of the three coordinates that describe the triatomic system, and as output, one of the four possible stable fixed points (see the appendix of [Verpoort et al., 2020b] for more details). We used Tensorflow (v2.4) and the Keras API (v2.4) for training. For the loss function, we used cross entropy with soft-max outputs and tanh (x) as activation for the hidden layers, as well as L2 regularization with a constant of λ = 10 −4 . The initialization for each layer is Keras default initialization (Xavier uniform). We used SGD with a learning rate of 0.05, a batch size of 100 and the training is run through 150000 steps.
For the minimum architecture, we put only one neuron in each hidden layer, resulting in 11, 13 and 15 parameters respectively. All other hyperparameters stay the same. For brevity, we show only the results for the case of three hidden layers.
For the MNIST dataset, 5 neurons are used in each hidden layer, resulting in architectures with 3985, 4015, 4045 parameters for one, two and three hidden layers respectively. The batch size is changed to 50 and all other hyperparameters are kept the same. The training is run through 300000 steps.  [Verpoort et al., 2020b] against the generalization error at the end of the training, t = 150000. Each marker indicates a different dataset size (described in the labels). There are three points for each marker type, corresponding to the architectures of one, two and three hidden layers.

Eigenvalues of the Hessian
Following the results presented by [Sagun et al., 2017], we show the distribution of the hessian eigenvalues in the first and last training steps in Fig. A.1, ordered from least to greatest. This is a different way to visualize the change in distribution from the start to the end of the training process, where we can observe a reduction in the number of negative eigenvalues.
The trace of the hessian has been used as a way to estimate the generalization error of neural networks [Jastrzębski et al., 2017]. This error is defined as the difference between the test and training loss [Bahri et al., 2020]. With the intention of confirming the correlation found by [Verpoort et al., 2020b] between the local curvatures and the generalizability of minima, in figure A.2 we plot the generalization error against the normalized (by the dataset size) trace of the hessian for different dataset sizes at the last training step. The positive correlation supports the findings in the aforementioned reference.

Gradient
Using the same architectures discussed by [Verpoort et al., 2020b], we calculate the magnitude of the gradient of the loss function in each training step ( figure A.3). In general, we observe a fast decrease to a value close to zero, around which it fluctuates. The noise in this part is dependent on the depth of the network, being more prominent for the network with three hidden layers and less so for the one with one hidden layer. The root-mean-square of the points after the step t = 10 4 is 0.02, 0.07 and 0.12 for one, two and three hidden layers respectively. This result can be rationalized from the perspective of the roughness of the landscape. As the network becomes deeper, the number of (shallow) local minima increases, and hence the stochastic dynamics induces higher fluctuations of the gradient during the last phase of the training.