Year: 2021

Abstract:
In this work we explore the limiting dynamics of deep neural networks trained with stochastic gradient descent (SGD). As observed previously, long after performance has converged, networks continue to move through parameter space by a process of anomalous diffusion in which distance travelled grows as a power law in the number of gradient updates with a nontrivial exponent. We reveal an intricate interaction between the hyperparameters of optimization, the structure in the gradient noise, and the Hessian matrix at the end of training that explains this anomalous diffusion. To build this understanding, we first derive a continuoustime model for SGD with finite learning rates and batch sizes as an underdamped Langevin equation. We study this equation in the setting of linear regression, where we can derive exact, analytic expressions for the phase space dynamics of the parameters and their instantaneous velocities from initialization to stationarity. Using the Fokker-Planck equation, we show that the key ingredient driving these dynamics is not the original training loss, but rather the combination of a modified loss, which implicitly regularizes the velocity, and probability currents, which cause oscillations in phase space. We identify qualitative and quantitative predictions of this theory in the dynamics of a ResNet-18 model trained on ImageNet. Through the lens of statistical physics, we uncover a mechanistic origin for the anomalous limiting dynamics of deep neural networks trained with SGD.

heading: 1 INTRODUCTION
text: Deep neural networks have demonstrated remarkable generalization across a variety of datasets and tasks. Essential to their success has been a collection of good practices on how to train these models with stochastic gradient descent (SGD). Yet, despite their importance, these practices are mainly based on heuristic arguments and trial and error search. Without a general theory connecting the hyperparameters of optimization, the architecture of the network, and the geometry of the dataset, theory-driven design of deep learning systems is impossible. Existing theoretical works studying this interaction have leveraged the random structure of neural networks at initialization [1, 2, 3] and in their infinite width limits in order to study their dynamics [4, 5, 6, 7, 8]. Here we take a different approach and study the training dynamics of pre-trained networks that are ready to be used for inference. By leveraging the mathematical structures found at the end of training, we uncover an intricate interaction between the hyperparameters of optimization, the structure in the gradient noise, and the Hessian matrix that corroborates previously identified empirical behavior such as anomalous limiting dynamics. Not only is understanding the limiting dynamics of SGD a critical stepping stone to building a complete theory for the learning dynamics of neural networks, but recently there have been a series of works demonstrating that the performance of pre-trained networks can be improved through averaging and ensembling [9, 10, 11]. Combining empirical exploration and theoretical tools from statistical physics, we identify and uncover a mechanistic explanation for the limiting dynamics of neural networks trained with SGD.

heading: 2 DIFFUSIVE BEHAVIOR IN THE LIMITING DYNAMICS OF SGD
text: A network that has converged in performance will continue to move through parameter space [12, 13, 14, 15]. To demonstrate this behavior, we resume training of pre-trained convolutional networks while tracking the network trajectory through parameter space. Let θ∗ ∈ Rm be the parameter vector for a pre-trained network and θk ∈ Rm be the parameter vector after k steps of resumed training.
We track two metrics of the training trajectory, namely the local parameter displacement δk between consecutive steps, and the global displacement ∆k after k steps from the pre-trained initialization:
δk = θk − θk−1, ∆k = θk − θ∗. (1) As shown in Fig. 1, neither of these differences converge to zero across a variety of architectures, indicating that despite performance convergence, the networks continue to move through parameter space, both locally and globally. The squared norm of the local displacement ‖δk‖22 remains near a constant value, indicating the network is essentially moving at a constant instantaneous speed. This observation is quite similar to the “equilibrium" phenomenon or “constant angular update" observed in Li et al. [17] and Wan et al. [13] respectively. However, these works only studied the displacement for parameters immediately preceding a normalization layer. The constant instantaneous speed behavior we observe is for all parameters in the model and is even present in models without normalization layers.
While the squared norm of the local displacement is essentially constant, the squared norm of the global displacement ‖∆k‖22 is monotonically growing for all networks, implying even once trained, the network continues to diverge from where it has been. Indeed Fig. 1 indicates a power law relationship between global displacement and number of steps, given by ‖∆k‖22 ∝ kc. As we’ll see in section 8, this relationship is indicative of anomalous diffusion where c corresponds to the anomalous diffusion exponent. Standard Brownian motion corresponds to c = 1. Similar observation were made by Baity-Jesi et al. [14] who noticed distinct phases of the training trajectory evident in the dynamics of the global displacement and Chen et al. [15] who found that the exponent of diffusion changes through the course of training. A parallel observation is given by Hoffer et al. [18] for the beginning of training, where they measure the global dis-
placement from the initialization of an untrained network and observe a rate ∝ log(k), a form of ultra-slow diffusion. These empirical observations raise the natural questions, where is the network moving to and why? To answer these questions we will build a diffusion based theory of SGD, study these dynamics in the setting of linear regression, and use lessons learned in this fundamental setting to understand the limiting dynamics of neural networks.

heading: 3 RELATED WORK
text: There is a long line of literature studying both theoretically and empirically the learning dynamics of deep neural networks trained with SGD. Our analysis and experiments build upon this literature.
Continuous models for SGD. Many works consider how to improve the classic gradient flow model for SGD to more realistically reflect momentum [19], discretization due to finite learning rates [20, 21], and stochasticity due to random batches [22, 23]. One line of work has studied the dynamics of networks in their infinite width limits through dynamical mean field theory [24, 25, 26, 27], while a different approach has used stochastic differential equations (SDEs) to model SGD directly, the approach we take in this work. However, recently, the validity of this approach has been questioned. The main argument, as nicely explained in Yaida [28], is that most SDE approximations simultaneously assume that ∆t → 0+, while maintaining that the learning rate η = ∆t is finite. The works Simsekli et al. [29] and Li et al. [30] have questioned the correctness of the using the central limit theorem (CLT) to model the gradient noise as Gaussian, arguing respectively that the heavy-tailed structure in the gradient noise and the weak dependence between batches leads the CLT to break down. In our work, we maintain the CLT assumption holds, which we discuss fur-
ther in appendix A, but importantly we avoid the pitfalls of many previous SDE approximations by simultaneously modeling the effect of finite learning rates and stochasticity.
Limiting dynamics. A series of works have applied SDE models of SGD to study the limiting dynamics of neural networks. In the seminal work by Mandt et al. [31], the limiting dynamics were modeled with a multivariate Ornstein-Uhlenbeck process by combining a first-order SDE model for SGD with assumptions on the geometry of the loss and covariance matrix for the gradient noise. This analysis was extended by Jastrzębski et al. [12] through additional assumptions on the covariance matrix to gain tractable insights and applied by Ali et al. [32] to the simpler setting of linear regression, which has a quadratic loss. A different approach was taken by Chaudhari and Soatto [33], which did not formulate the dynamics as an OU process, nor assume directly a structure on the loss or gradient noise. Rather, this analysis studied the same first-order SDE via the Fokker-Planck equation to propose the existence of a modified loss and probability currents driving the limiting dynamics, but did not provide explicit expressions. Our analysis deepens and combines ideas from all these works, where our key insight is to lift the dynamics into phase space. By studying the dynamics of the parameters and their velocities, and by applying the analysis first in the setting of linear regression where assumptions are provably true, we are able to identify analytic expressions and explicit insights which lead to concrete predictions and testable hypothesis.
Stationary dynamics. A different line of work avoids modeling the limiting dynamics of SGD with an SDE and instead chooses to leverage the property of stationarity. These works [28, 34, 35, 36] assume that eventually the probability distribution governing the model parameters reaches stationarity such that the discrete SGD process is simply sampling from this distribution. Yaida [28] used this approach to derive fluctuation-dissipation relations that link measurable quantities of the parameters and hyperparameters of SGD. Liu et al. [35] used this approach to derive properties for the stationary distribution of SGD with a quadratic loss. Similar to our analysis, this work identifies that the stationary distribution for the parameters reflects a modified loss function dependent on the relationship between the covariance matrix of the gradient noise and the Hessian matrix for the original loss.
Empirical exploration. Another set of works analyzing the limiting dynamics of SGD has taken a purely empirical approach. Building on the intuition that flat minima generalize better than sharp minima, Keskar et al. [37] demonstrated empirically that the hyperparameters of optimization influence the eigenvalue spectrum of the Hessian matrix at the end of training. Many subsequent works have studied the Hessian eigenspectrum during and at the end of training. Jastrzębski et al. [38], Cohen et al. [39] studied the dynamics of the top eigenvalues during training. Sagun et al. [40], Papyan [41], Ghorbani et al. [42] demonstrated the spectrum has a bulk of values near zero plus a small number of larger outliers. Gur-Ari et al. [43] demonstrated that the learning dynamics are constrained to the subspace spanned by the top eigenvectors, but found no special properties of the dynamics within this subspace. In our work we also determine that the top eigensubspace of the Hessian plays a crucial role in the limiting dynamics and by projecting the dynamics into this subspace in phase space, we see that the motion is not random, but consists of incoherent oscillations leading to anomalous diffusion.

heading: 4 MODELING SGD AS AN UNDERDAMPED LANGEVIN EQUATION
text: Following the route of previous works [31, 12, 33] studying the limiting dynamics of neural networks, we first seek to model SGD as a continuous stochastic process. We consider a network parameterized by θ ∈ Rm, a training dataset {x1, . . . , xN} of size N , and a training loss L(θ) = 1N ∑N i=1 `(θ, xi) with corresponding gradient g(θ) = ∂L ∂θ . The state of the network at the kth step of training is defined by the position vector θk and velocity vector vk of the same dimension. The gradient descent update with learning rate η, momentum β, and weight decayλ is given by vk+1 = βvk − g(θk)− λθk, θk+1 = θk + ηvk+1, (2) where we initialize the network such that v0 = 0 and θ0 is the parameter initialization. In order to understand the dynamics of the network through position and velocity space, which we will refer to as phase space, we express these discrete recursive equations as the discretization of some unknown ordinary differential equation (ODE), sometimes referred to as a modified equation as in [44, 20]. While this ODE models the gradient descent process even at finite learning rates, it fails to account for the stochasticity introduced by choosing a random batch B of size S drawn uniformly from the set of N training points. This sampling yields the stochastic gradient gB(θ) = 1S ∑ i∈B∇`(θ, xi). To model this effect, we make the following assumption:
Assumption 1 (CLT). We assume the batch gradient is a noisy version of the true gradient such that gB(θ)− g(θ) is a Gaussian random variable with mean 0 and covariance 1SΣ(θ).
Incorporating this model of stochastic gradients into the previous finite difference equation and applying the stochastic counterparts to Euler discretizations, results in the standard drift-diffusion stochastic differential equation (SDE), referred to as an underdamped Langevin equation,
d [ θ v ] = [ v − 2η(1+β) (g(θ) + λθ + (1− β)v) ] dt+ [ 0 0 0 2√ ηS(1+β) √ Σ(θ) ] dWt, (3)
where Wt is a standard Wiener process. This is the continuous model we will study in this work: Assumption 2 (SDE). We assume the underdamped Langevin equation (3) accurately models the trajectory of the network driven by SGD through phase space such that θ(ηk) ≈ θk and v(ηk) ≈ vk.
See appendix A for further discussion on the nuances of modeling SGD with an SDE.

heading: 5 LINEAR REGRESSION WITH SGD IS AN ORNSTEIN-UHLENBECK PROCESS
text: Equipped with a model for SGD, we seek to understand its dynamics in the fundamental setting of linear regression, one of the few cases where we have a complete model for the interaction of the dataset, architecture, and optimizer. Let X ∈ RN×d be the input data, Y ∈ RN be the output labels, and θ ∈ Rd be our vector of regression coefficients. The least squares loss is the convex quadratic lossL(θ) = 12N ‖Y −Xθ‖2 with gradient g(θ) = Hθ−b, whereH = X ᵀX N and b = XᵀY N . Plugging this expression for the gradient into the underdamped Langevin equation (3), and rearranging terms, results in the multivariate Ornstein-Uhlenbeck (OU) process,
d [ θt vt ] = − [ 0 −I
2 η(1+β) (H + λI) 2(1−β) η(1+β)I ] ︸ ︷︷ ︸
A
([ θt vt ] − [ µ 0 ]) dt+ √ 2κ−1 √√√√√ [ 0 0 0 2(1−β)η(1+β)Σ(θ) ] ︸ ︷︷ ︸
D
dWt,
(4)
where A and D are the drift and diffusion matrices respectively, κ = S(1− β2) is an inverse temperature constant, and µ = (H + λI)−1b is the ridge regression solution. The solution to an OU process is a Gaussian process. By solving for the temporal dynamics of the first and second moments of the process, we can obtain an analytic expression for the trajectory at any time t. In particular, we can decompose the trajectory as the sum of a deterministic and stochastic component defined by the first and second moments respectively.
Deterministic component. Using the form of A we can decompose the expectation as a sum of harmonic oscillators in the eigenbasis {q1, . . . , qm} of the Hessian,
E [[ θt vt ]] = [ µ 0 ] + m∑ i=1 ( ai(t) [ qi 0 ] + bi(t) [ 0 qi ]) . (5)
Here the coefficients ai(t) and bi(t) depend on the optimization hyperparameters η, β, λ, S and the respective eigenvalue of the Hessian ρi as further explained in appendix F. We verify this expression nearly perfectly matches empirics on complex datasets under various hyperparameter settings as shown in Fig. 2.
Stochastic component. The cross-covariance of the process between two points in time t ≤ s, is
Cov ([ θt vt ] , [ θs vs ]) =κ−1 ( B−e−AtBe−Aᵀt ) eA ᵀ(t−s),
(6)
where B solves the Lyapunov equation AB +BAᵀ = 2D. In order to gain analytic expressions for B in terms of the optimization hyperparameters, eigendecomposition of the Hessian, and covariance of the gradient noise, we must introduce the following assumption: Assumption 3 (Simultaneously Diagonalizable). We assume the covariance of the gradient noise is spatially independent Σ(θ) = Σ and commutes with the Hessian HΣ = ΣH , therefore sharing a common eigenbasis.

heading: 6 UNDERSTANDING STATIONARITY VIA THE FOKKER-PLANCK EQUATION
text: The OU process is unique in that it is one of the few SDEs which we can solve exactly. As shown in section 5, we were able to derive exact expressions for the dynamics of linear regression trained with SGD from initialization to stationarity by simply solving for the first and second moments. While the expression for the first moment provides an understanding of the intricate oscillatory relationship in the deterministic component of the process, the second moment, driving the stochastic component, is much more opaque. An alternative route to solving the OU process that potentially provides more insight is the Fokker-Planck equation.
The Fokker-Planck (FP) equation is a PDE describing the time evolution for the probability distribution of a particle governed by Langevin dynamics. For an arbitrary potential Φ and diffusion matrix D, the Fokker-Planck equation (under an Itô integration prescription) is
∂tp = ∇ · ( ∇Φp+∇ · ( κ−1Dp ))︸ ︷︷ ︸ −J , (7)
where p represents the time-dependent probability distribution, and J is a vector field commonly referred to as the probability current. The FP equation is especially useful for explicitly solving for the stationary solution, assuming one exists, of the Langevin dynamics. The stationary solution pss by definition obeys ∂tpss = 0 or equivalently ∇ · Jss = 0. From this second definition we see that there are two distinct settings of stationarity: detailed balance when Jss = 0 everywhere, or broken detailed balance when ∇ · Jss = 0 and Jss 6= 0. For a general OU process, the potential is a convex quadratic function Φ(x) = xᵀAx defined by the drift matrix A. When the diffusion matrix is isotropic (D ∝ I) and spatially independent (∇ · D = 0) the resulting stationary solution is a Gibbs distribution pss(x) ∝ e−κΦ(x) determined by the original loss Φ(x) and is in detailed balance. Lesser known properties of the OU process arise when the diffusion matrix is anisotropic or spatially dependent [45, 46]. In this setting the solution is still a Gaussian process, but the stationary solution, if it exists, is no longer defined by the Gibbs distribution of the original loss Φ(x), but actually a modified loss Ψ(x). Furthermore, the stationary solution may be in broken detailed balance leading to a non-zero probability current Jss(x). Depending on the relationship between the drift matrix A and the diffusion matrix D the resulting dynamics of the OU process can have very nontrivial behavior.
In the setting of linear regression, anisotropy in the data distribution will lead to anisotropy in the gradient noise and thus an anisotropic diffusion matrix. This implies that for most datasets we should expect that the SGD trajectory is not driven by the original least squares loss, but by a modified loss and converges to a stationary solution with broken detailed balance, as predicted by Chaudhari and Soatto [33]. Using the explicit expressions for the drift A and diffusion D matrices we can compute analytically the modified loss and stationary probability current,
Ψ(θ, v) = ([ θ v ] − [ µ 0 ])ᵀ( U
2
)([ θ v ] − [ µ 0 ]) , Jss(θ, v) = −QU ([ θ v ] − [ µ 0 ]) pss, (8)
where Q is a skew-symmetric matrix and U is a positive definite matrix defined as,
Q =
[ 0 −Σ(θ)
Σ(θ) 0
] , U = [ 2 η(1+β)Σ(θ) −1 (H + λI) 0
0 Σ(θ)−1
] . (9)
These new fundamental matrices, Q and U , relate to the original drift A and diffusion D matrices through the unique decomposition A = (D + Q)U , introduced by Ao [47] and Kwon et al. [48]. Using this decomposition we can easily show that B = U−1 solves the Lyapunov equation and indeed the stationary solution pss is the Gibbs distribution defined by the modified loss Ψ(θ, v) in equation (8). Further, the stationary cross-covariance solved in section 5 reflects the oscillatory dynamics introduced by the stationary probability currents Jss(θ, v) in equation (8). Taken together, we gain the intuition that the limiting dynamics of SGD in linear regression are driven by a modified loss subject to oscillatory probability currents.

heading: 7 EVIDENCE OF A MODIFIED LOSS AND OSCILLATIONS IN DEEP LEARNING
text: Does the theory derived in the linear regression setting (sections 5, 6) help explain the empirical phenomena observed in the non-linear setting of deep neural networks (section 2)? In order for the theory built in the previous sections to apply to the limiting dynamics of neural networks, we must introduce simplifying assumptions on the loss landscape and gradient noise at the end of training:
Assumption 4 (Quadratic Loss). We assume that at the end of training the loss for a neural network can be approximated by the quadratic loss L(θ) = (θ − µ)ᵀ ( H 2 ) (θ − µ), where H 0 is the training loss Hessian and µ is some unknown mean vector, corresponding to a local minimum.
Assumption 5 (Covariance Structure). We assume the covariance of the gradient noise is proportional to the Hessian of the quadratic loss Σ(θ) = σ2H where σ ∈ R+ is some unknown scalar.
Under these simplifications, then the expressions derived in the linear regression setting would apply to the limiting dynamics of deep neural networks and depend only on quantities that we can easily estimate empirically. Of course, these simplifications are quite strong, but without arguing their theoretical validity, we can empirically test their qualitative implications: (1) a modified isotropic loss driving the limiting dynamics through parameter space, (2) implicit regularization of the velocity trajectory, and (3) oscillatory phase space dynamics determined by the Hessian eigen-structure.
Modified loss. As discussed in section 6, due to the anisotropy of the diffusion matrix, the loss landscape driving the dynamics at the end of training is not the original training loss L(θ), but a modified loss Ψ(θ, v) in phase space. As shown in equation (8), the modified loss decouples into a term Ψθ that only depends on the parameters θ and a term Ψv that only depends on the velocities v. Under assumption 5, the parameter dependent component is proportional to the convex quadratic,
Ψθ ∝ (θ − µ)ᵀ ( H−1(H + λI)
η(1 + β)
) (θ − µ) . (10)
This quadratic function has the same mean µ as the training loss, but a different curvature. Using this expression, notice that when λ ≈ 0, the modified loss is isotropic in the column space of H , regardless of what the nonzero eigenspectrum of H is. This striking prediction suggests that no matter how anisotropic the original training loss – as reflected by poor conditioning of the Hessian eigenspectrum – the training trajectory of the network will behave isotropically, since it is driven not by the original anisotropic loss, but a modified isotropic loss.
We test this prediction by studying the limiting dynamics of a pre-trained ResNet-18 model with batch normalization that we continue to train on ImageNet according to the last setting of its hyperparameters [49]. Let θ∗ represent the initial pre-trained parameters of the network, depicted with the white dot in figures 3 and 4.
We estimate1 the top thirty eigenvectors q1, . . . , q30 of the Hessian matrix H∗ evaluated at θ∗ and project the limiting trajectory for the parameters onto the plane spanned by the top q1 and bottom q30 eigenvectors to maximize the illustrated anisotropy with our estimates. We sample the train and test loss in this subspace for a region around the projected trajectory. Additionally, using the hyperparameters of the optimization, the eigenvalues ρ1 and ρ30, and the estimate for the mean µ = θ∗−H−1∗ g∗ (g∗ is the gradient evaluated at θ∗), we also sample from the modified loss equation (10) in the same region. Figure 3 shows the projected parameter trajectory on the sampled train, test and modified losses. Contour lines of both the train and test loss exhibit anisotropic structure, with sharper curvature along eigenvector q1 compared to eigenvector q30, as expected. However, as predicted, the trajectory appears to cover both directions equally. This striking isotropy of the trajectory within a highly anisotropic slice of the loss landscape indicates qualitatively that the trajectory evolves in a modified isotropic loss landscape.
Implicit velocity regularization. A second qualitative prediction of the theory is that the velocity is regulated by the inverse Hessian of the training loss. Of course there are no explicit terms in either the train or test losses that depend on the velocity. Yet, the modified loss contains a component, Ψv ∝ vᵀH−1v, that only depends on the velocities This additional term can be understood as a form of implicit regularization on the velocity trajectory. Indeed, when we project the velocity trajectory onto the plane spanned by the q1 and q30 eigenvectors, as shown in Fig. 4, we see that the trajectory closely resembles the curvature of the inverse Hessian H−1. The modified loss is effectively penalizing SGD for moving in eigenvectors of the Hessian with small eigenvalues. A similar qualitative effect was recently proposed by Barrett and Dherin [21] as a consequence of the discretization error due to finite learning rates.
Phase space oscillations. A final implication of the theory is that at stationarity the network is in broken detailed balance leading to non-zero probability currents flowing through phase space:
Jss(θ, v) =
[ v − 2η(1+β) (H + λI) (θ − µ) ] pss. (11)
These probability currents encourage oscillatory dynamics in the phase space planes characterized by the eigenvectors of the Hessian, at rates proportional to their eigenvalues. We consider the same projected trajectory of the ResNet-18 model visualized in figures 3 and 4, but plot the trajectory in phase space for the two eigenvectors q1 and q30 separately. Shown in Fig. 5, we see that both trajectories look like noisy clockwise rotations. Qualitatively, the trajectories for the different eigenvectors appear to be rotating at different rates.
The integral curves of the stationary probability current are one-dimensional paths confined to level sets of the modified loss. These paths might cross themselves, in which case they are limit cycles, or they could cover the entire surface of the level sets, in which case they are space-filling curves. This distinction depends on the relative frequencies of the oscillations, as determined by the pairwise
1To estimate the eigenvectors of H∗ we use subspace iteration, and limit ourselves to 30 eigenvectors to constrain computation time. See appendix H for details.
ratios of the eigenvalues of the Hessian. For real-world datasets, with a large spectrum of incommensurate frequencies, we expect to be in the latter setting, thus contradicting the suggestion that SGD in deep networks converges to limit cycles, as claimed in Chaudhari and Soatto [33].

heading: 8 UNDERSTANDING THE DIFFUSIVE BEHAVIOUR OF THE LIMITING DYNAMICS
text: Taken together the empirical results shown in section 7 indicate that many of the same qualitative behaviors of SGD identified theoretically for linear regression are evident in the limiting dynamics of neural networks. Can this theory quantitatively explain the results we identified in section 2?
Constant instantaneous speed. As noted in section 2, we observed that at the end of training, across various architectures, the squared norm of the local displacement ‖δt‖22 remains essentially constant. Assuming the limiting dynamics are described by the stationary solution the expectation of the local displacement is
Ess [ ‖δt‖2 ] =
η2
S(1− β2)σ 2tr (H) , (12)
as derived in appendix G. We cannot test this prediction directly as we do not know σ2 and computing tr(H) is computationally prohibitive. However, we can estimate σ2tr(H) by resuming training for a model, measuring the average ‖δt‖2, and then inverting equation (12). Using this single estimate, we find that for a sweep of models with varying hyperparameters, equation (12) accurately predicts their instantaneous speed. Indeed, Fig. 6 shows an exact match between the empirics and theory, which strongly suggests that despite changing hyperparameters at the end of training, the model remains in the same quadratic basin.
Exponent of anomalous diffusion. The expected value for the global displacement under the stationary solution can also be analytically expressed in terms of the optimization hyperparameters and the eigendecomposition of the Hessian as,
Ess [ ‖∆t‖2 ] =
η2
S(1− β2)σ 2
( tr (H) t+ 2t
t∑ k=1 ( 1− k t ) m∑ l=1 ρlCl(k) ) , (13)
where Cl(k) is a trigonometric function describing the velocity of a harmonic oscillator with damping ratio ζl = (1 − β)/ √ 2η(1 + β) (pl + λ), see appendix G for details. As shown empirically in section 2, the squared norm ‖∆t‖2 monotonically increases as a power law in the number of steps, suggesting its expectation is proportional to tc for some unknown, constant c. The exponent c determines the regime of diffusion for the process. When c = 1, the process corresponds to standard Brownian diffusion. For c > 1 or c < 1 the process corresponds to anomalous super-diffusion or sub-diffusion respectively. Unfortunately, it is not immediately clear how to extract the explicit
exponent c from equation (13). However, by exploring the functional form of Cl(k) and its relationship to the hyperparameters of optimization through the damping ratio ζl, we can determine overall trends in the diffusion exponent c.
Akin to how the exponent c determines the regime of diffusion, the damping ratio ζl determines the regime for the harmonic oscillator describing the stationary velocity-velocity correlation in the lth eigenvector of the Hessian. When ζl = 1, the oscillator is critically damped implying the velocity correlations converge to zero as quickly as possible. In the extreme setting of Cl(k) = 0 for all l, k, then equation (13) simplifies to standard Brownian diffusion, Ess [ ‖∆t‖2 ] ∝ t. When ζl > 1, the oscillator is overdamped implying the velocity correlations dampen slowly and remain positive even over long temporal lags. Such long lasting temporal correlations in velocity lead to faster global displacement. Indeed, in the extreme setting of Cl(k) = 1 for all l, k, then equation (13) simplifies to a form of anomalous super-diffusion, Ess [ ‖∆t‖2 ] ∝ t2. When ζl < 1, the oscillator is underdamped implying the velocity correlations will oscillate quickly between positive and negative values. Indeed, the only way equation (13) could describe anomalous sub-diffusion is if Cl(k) took on negative values for certain l, k.
Using the same sweep of models described previously, we can empirically confirm that the optimization hyperparameters each influence the diffusion exponent c. As shown in Fig. 6, the learning rate, batch size, and momentum can each independently drive the exponent c into different regimes of anomalous diffusion. Notice how the influence of the learning rate and momentum on the diffusion exponent c closely resembles their respective influences on the damping ratio ζl. Interestingly, a larger learning rate leads to underdamped oscillations, and the resultant temporal velocities’ anti-correlations reduce the exponent of anomalous diffusion. Thus contrary to intuition, a larger learning rate actually leads to slower global transport in parameter space. The batch size on the other hand, has no influence on the damping ratio, but leads to an interesting, non-monotonic influence on the diffusion exponent. Overall, the hyperparameters of optimization and eigenspectrum of the Hessian all conspire to govern the degree of anomalous diffusion at the end of training.

heading: 9 DISCUSSION
text: Through combined empirics and theory based on statistical physics, we uncovered an intricate interplay between the optimization hyperparameters, structure in the gradient noise, and the Hessian matrix at the end of training.
Significance. The significance of our work lies in (1) the identification/verification of multiple empirical phenomena (constant instantaneous speed, anomalous diffusion in global displacement, isotropic parameter exploration despite anisotopic loss, velocity regularization, and slower global parameter exploration with faster learning rates) present in the limiting dynamics of deep neural networks, (2) the emphasis on studying the dynamics in velocity space in addition to parameter space, and (3) concrete quantitative as well as qualitative predictions of an SDE based theory that we empirically verified in deep networks trained on large scale datasets (indeed some of the above nontrivial phenomena were predictions of this theory). Of course, these contributions directly build upon a series of related works studying the immensely complex process of deep learning. To this end, we further clarify the originality of our contributions with respect to some relevant works.
Originality. The empirical phenomena we present provide novel insight with respect to the works of Wan et al. [13], Hoffer et al. [18], and Chen et al. [15]. We observe that all parameters in the network (not just those with scale symmetry) move at a constant instantaneous speed at the end of training and diffuse anomalously at rates determined by the hyperparameters of optimization. In contrast to the work by Liu et al. [35], we modeled the entire SGD process as an OU process which allows us to provide insight into the transient dynamics and identify oscillations in parameter and velocity space. We build on the theoretical framework used by Chaudhari and Soatto [33] and provide explicit expressions for the limiting dynamics in the simplified linear regression setting and conclude that the oscillations present in the limiting dynamics are more likely to be space-filling curves (and not limit cycles) in deep learning due to many incommensurate oscillations.
Overall, by identifying key phenomena, explaining them in a simpler setting, deriving predictions of new phenomena, and providing evidence for these predictions at scale, we are furthering the scientific study of deep learning. We hope our newly derived understanding of the limiting dynamics of SGD, and its dependence on various important hyperparameters like batch size, learning rate, and momentum, can serve as a basis for future work that can turn these insights into algorithmic gains.

heading: A MODELING SGD WITH AN SDE
text: As explained in section 4, in order to understand the dynamics of stochastic gradient descent we build a continuous Langevin equation in phase space modeling the effect of discrete updates and stochastic batches simultaneously.

heading: A.1 MODELING DISCRETIZATION
text: To model the discretization effect we assume that the system of update equations (2) is actually a discretization of some unknown ordinary differential equation. To uncover this ODE, we combine the two update equations in (2), by incorporating a previous time step θk−1, and rearrange into the form of a finite difference discretization, as shown in equation (??). Like all discretizations, the Euler discretizations introduce error terms proportional to the step size, which in this case is the learning rate η. Taylor expanding θk+1 and θk−1 around θk, its easy to show that both Euler discretizations introduce a second-order error term proportional to η2 θ̈.
θk+1 − θk η = θ̇ + η 2 θ̈ +O(η2),
θk − θk−1 η = θ̇ − η 2 θ̈ +O(η2).
Notice how the momentum coefficient β ∈ [0, 1] regulates the amount of backward Euler incorporated into the discretization. When β = 0, we remove all backward Euler discretization leaving just the forward Euler discretization. When β = 1, we have equal amounts of backward Euler as forward Euler resulting in a central second-order discretization2 as noticed in [19].

heading: A.2 MODELING STOCHASTICITY
text: In order to model the effect of stochastic batches, we first model a batch gradient with the following assumption:
Assumption 1 (CLT). We assume the batch gradient is a noisy version of the true gradient such that gB(θ)− g(θ) is a Gaussian random variable with mean 0 and covariance 1SΣ(θ). The two conditions needed for the CLT to hold are not exactly met in the setting of SGD. Independent and identically distributed. Generally we perform SGD by making a complete pass through the entire dataset before using a sample again which introduces a weak dependence between samples. While the covariance matrix without replacement more accurately models the dependence between samples within a batch, it fails to account for the dependence between batches. Finite variance. A different line of work has questioned the Gaussian assumption entirely because of the need for finite variance random variables. This work instead suggests using the generalized central limit theorem implying the noise would be a heavy-tailed α-stable random variable [29]. Thus, the previous assumption is implicitly assuming the i.i.d. and finite variance conditions apply for large enough datasets and small enough batches.
Under the CLT assumption, we must also replace the Euler discretizations with Euler–Maruyama discretizations. For a general stochastic process, dXt = µdt+ σdWt, the Euler–Maruyama method extends the Euler method for ODEs to SDEs, resulting in the update equation Xk+1 = Xk + ∆tµ+√
∆tσξ, where ξ ∼ N (0, 1). Notice, the key difference is that if the temporal step size is ∆t = η, then the noise is scaled by the square root √ η. In fact, the main argument against modeling SGD with an SDE, as nicely explained in Yaida [28], is that most SDE approximations simultaneously assume that ∆t → 0+, while maintaining that the square root of the learning rate √η is finite. However, by modeling the discretization and stochastic effect simultaneously we can avoid this argument, bringing us to our second assumption:
Assumption 2 (SDE). We assume the underdamped Langevin equation (3) accurately models the trajectory of the network driven by SGD through phase space such that θ(ηk) ≈ θk and v(ηk) ≈ vk. This approach of modeling discretization and stochasticity simultaneously is called stochastic modified equations, as further explained in Li et al. [22].
2The difference between a forward Euler and backward Euler discretization is a second-order central discretization, ( θk+1−θk
η
) − ( θk−θk−1
η
) = η ( θk+1−2θk+θk−1
η2
) = ηθ̈ +O(η2).

heading: B STRUCTURE IN THE COVARIANCE OF THE GRADIENT NOISE
text: As we’ve mentioned before, SGD introduces highly structured noise into an optimization process, often assumed to be an essential ingredient for its ability to avoid local minima.
Assumption 5 (Covariance Structure). We assume the covariance of the gradient noise is proportional to the Hessian of the quadratic loss Σ(θ) = σ2H where σ ∈ R+ is some unknown scalar. In the setting of linear regression, this is a very natural assumption. If we assume the classic generative model for linear regression data yi = x ᵀ i θ̄+σ where, θ̄ ∈ Rd is the true model and ∼ N (0, 1), then provably Σ(θ) ≈ σ2H . Proof. We can estimate the covariance as Σ(θ) ≈ 1N ∑N i=1 gig ᵀ i − ggᵀ. Near stationarity ggᵀ 1 N ∑N i=1 gig ᵀ i , and thus,
Σ(θ) ≈ 1 N N∑ i=1 gig ᵀ i .
Under the generative model yi = x ᵀ i θ̄ + σ where ∼ N (0, 1) and σ ∈ R+, then the gradient gi is gi = (x ᵀ i (θ − θ̄)− σ )xi, and the matrix gig ᵀ i is
gig ᵀ i = (x ᵀ i (θ − θ̄)− σ )2(xixᵀi ).
Assuming θ ≈ θ̄ at stationarity, then (xᵀi (θ − θ̄)− σ )2 ≈ σ2. Thus,
Σ(θ) ≈ σ 2
N N∑ i=1 xix ᵀ i = σ2 N XᵀX = σ2H
Also notice that weight decay is independent of the data or batch and thus simply shifts the gradient distribution, but leaves the covariance of the gradient noise unchanged.
While the above analysis is in the linear regression setting, for deep neural networks it is reasonable to make the same assumption. See the appendix of Jastrzębski et al. [12] for a discussion on this assumption in the non-linear setting.
Recent work by Ali et al. [32] also studies the dynamics of SGD (without momentum) in the setting of linear regression. This work, while studying the classic first-order stochastic differential equation, made a point to not introduce an assumption on the diffusion matrix. In particular, they make the point that even in the setting of linear regression, a constant covariance matrix will fail to capture the actual dynamics. To illustrate this point they consider the univariate responseless least squares problem,
minimize θ∈R
1
2n n∑ i=1 (xiθ) 2.
As they explain, the SGD update for this problem would be
θk+1 = θk − η
S (∑ i∈B xi ) θk = k∏ i=1 (1− η( 1S ∑ i∈B xi))θ0,
from which they conclude for a small enough learning rate η, then with probability one θk → 0. They contrast this with the Ornstein-Uhlenbeck process given by a constant covariance matrix where while the mean for θk converges to zero its variance converges to a positive constant. So is this discrepancy evidence that an Ornstein-Uhlenbeck process with a constant covariance matrix fails to capture the updates of SGD? In many ways this problem is not a simple example, rather a pathological edge case. Consider the generative model that would give rise to this problem,
y = 0x+ 0ξ = 0.
In otherwords, the true model θ̄ = 0 and the standard deviation for the noise σ = 0. This would imply by the assumption used in our paper that there would be zero diffusion and the resulting SDE would simplify to a deterministic ODE that exponentially converges to zero.

heading: C A QUADRATIC LOSS AT THE END OF TRAINING
text: Assumption 4 (Quadratic Loss). We assume that at the end of training the loss for a neural network can be approximated by the quadratic loss L(θ) = (θ − µ)ᵀ ( H 2 ) (θ − µ), where H 0 is the training loss Hessian and µ is some unknown mean vector, corresponding to a local minimum.
This assumption has been amply used in previous works such as Mandt et al. [31], Jastrzębski et al. [12], and Poggio et al. [50]. Particularly, Mandt et al. [31] discuss how this assumption makes sense for smooth loss functions for which the stationary solution to the stochastic process reaches a deep local minimum from which it is difficult to escape.
It is a well-studied fact, both empirically and theoretically, that the Hessian is low-rank near local minima as noted by Sagun et al. [51], and Kunin et al. [20]. This degeneracy results in flat directions of equal loss. Kunin et al. [20] discuss how differentiable symmetries, architectural features that keep the loss constant under certain weight transformations, give rise to these flat directions. Importantly, the Hessian and the covariance matrix share the same null space, and thus we can always restrict ourselves to the image space of the Hessian, where the drift and diffusion matrix will be full rank. Further discussion on the relationship between the Hessian and the covariance matrix can be found in Thomas et al. [52].
It is also a well known empirical fact that even at the end of training the Hessian can have negative eigenvalues [41]. This empirical observation is at odds with our assumption that the Hessian is positive semi-definite H 0. Further analysis is needed to alleviate this inconsistency.

heading: D SOLVING AN ORNSTEIN-UHLENBECK PROCESS WITH ANISOTROPIC NOISE
text: We will study the multivariate Ornstein-Uhlenbeck process described by the stochastic differential equation
dXt = A(µ−Xt)dt+ √
2κ−1DdWt X0 = x0, (14) whereA ∈ Sm++ is a positive definite drift matrix, µ ∈ Rm is a mean vector, κ ∈ R+ is some positive constant, and D ∈ Sm++ is a positive definite diffusion matrix. This OU process is unique in that it is one of the few SDEs we can solve explicitly. We can derive an expression for XT as,
XT = e −ATx0 + ( I − e−AT ) µ+ ∫ T 0 eA(t−T ) √ 2κ−1DdWt. (15)
Proof. Consider the function f(t, x) = eAtx where eA is a matrix exponential. Then by Itô’s Lemma3 we can evaluate the derivative of f(t,Xt) as
df(t,Xt) = ( AeAtXt + e AtA(µ−Xt) ) dt+ eAt √ 2κ−1DdWt
= AeAtµdt+ eAt √
2κ−1DdWt Integrating this expression from t = 0 to t = T gives
f(T,XT )− f(0, X0) = ∫ T
0
AeAtµdt+ ∫ T 0 eAt √ 2κ−1DdWt
eATXT − x0 = ( eAT − I ) µ+ ∫ T 0 eAt √ 2κ−1DdWt
which rearranged gives the expression for XT .
From this expression it is clear that XT is a Gaussian process. The mean of the process is E [XT ] = e −ATx0 + ( I − e−AT ) µ, (16)
and the covariance and cross-covariance of the process are Var(XT ) = κ −1 ∫ T
0
eA(t−T )2DeA ᵀ(t−T )dt, (17)
Cov(XT , XS) = κ −1 ∫ min(T,S)
0
eA(t−T )2DeA ᵀ(t−S)dt. (18)
These last two expressions are derived by Itô Isometry4.

heading: D.1 THE LYAPUNOV EQUATION
text: We can explicitly solve the integral expressions for the covariance and cross-covariance exactly by solving for the unique matrix B ∈ Sm++ that solves the Lyapunov equation, AB +BAᵀ = 2D. (19) If B solves the Lyapunov equation, notice
d
dt
( eA(t−T )BeA ᵀ(t−S) ) = eA(t−T )ABeA ᵀ(t−S) + eA(t−T )BAᵀeA ᵀ(t−S)
= eA(t−T )2DeA ᵀ(t−S)
Using this derivative, the integral expressions for the covariance and cross-covariance simplify as, Var(XT ) = κ −1 ( B − e−ATBe−AᵀT ) , (20)
Cov(XT , XS) = κ −1 ( B − e−ATBe−AᵀT ) eA ᵀ(T−S), (21)
where we implicitly assume T ≤ S. 3Itô’s Lemma states that for any Itô drift-diffusion process dXt = µtdt + σtdWt and twice differentiable
scalar function f(t, x), then df(t,Xt) = ( ft + µtfx + σ2t 2 fxx ) dt+ σtfxdWt.
4Itô Isometry states for any standard Itô process Xt, then E [(∫ t
0 XtdWt
)2] = E [∫ t 0 X2t dt ] .

heading: D.2 DECOMPOSING THE DRIFT MATRIX
text: While the Lyapunov equation simplifies the expressions for the covariance and cross-covariance, it does not explain how to actually solve for the unknown matrix B. Following a method proposed by Kwon et al. [48], we will show how to solve for B explicitly in terms of the drift A and diffusion D.
The drift matrix A can be uniquely decomposed as, A = (D +Q)U (22) whereD is our symmetric diffusion matrix,Q is a skew-symmetric matrix (i.e. Q = −Qᵀ), and U is a positive definite matrix. Using this decomposition, then B = U−1, solves the Lyapunov equation.
Proof. Plug B = U−1 into the left-hand side of equation (19), AU−1 + U−1Aᵀ = (D +Q)UU−1 + U−1U(D −Q)
= (D +Q) + (D −Q) = 2D
Here we used the symmetry of A,D,U and the skew-symmetry of Q.
All that is left is to do is solve for the unknown matricesQ and U . First notice the following identity, AD −DA = QA+AQ (23)
Proof. Multiplying A = (D +Q)U on the right by (D −Q) gives, A(D −Q) = (D +Q)U(D −Q)
= (D +Q)Aᵀ,
which rearranged and using A = Aᵀ gives the desired equation.
Let V ΛV ᵀ be the eigendecomposition of A and define the matrices D̃ = V ᵀDV and Q̃ = V ᵀQV . These matrices observe the following relationship,
Q̃ij = λi − λj ρi + λj D̃ij . (24)
Proof. Replace A in the previous equality with its eigendecompsoition, V ΛV ᵀD −DV ΛV ᵀ = QV ΛV ᵀ + V ΛV ᵀQ.
Multiply this equation on the right by V and on the left by V ᵀ, ΛD̃ − D̃Λ = Q̃Λ + ΛQ̃.
Looking at this equality element-wise and using the fact that Λ is diagonal gives the scalar equality for any i, j, (λi − λj)D̃ij = (λi + λj)Q̃ij , which rearranged gives the desired expression.
Thus, Q and U are given by, Q = V Q̃V ᵀ, U = (D +Q)−1A. (25)
This decomposition always holds uniquely when A,D 0, as λi−λjλi+λj exists and (D +Q) is invertible. See [48] for a discussion on the singularities of this decomposition.

heading: D.3 STATIONARY SOLUTION
text: Using the Lyapunov equation and the drift decomposition, then XT ∼ pT , where pT = N ( e−ATx0 + ( I − e−AT ) µ, κ−1 ( U−1 − e−ATU−1e−AᵀT )) . (26)
In the limit as T →∞, then e−AT → 0 and pT → pss where pss = N ( µ, κ−1U−1 ) . (27)
Similarly, the cross-covariance converges to the stationary cross-covariance, Covss(XT , XS) = κ −1BeA ᵀ(T−S). (28)

heading: E A VARIATIONAL FORMULATION OF THE OU PROCESS WITH ANISOTROPIC NOISE
text: In this section we will describe an alternative, variational, route towards solving the dynamics of the OU process studied in appendix D.
Let Φ : Rn → R be an arbitrary, non-negative potential and consider the stochastic differential equation describing the Langevin dynamics of a particle in this potential field,
dXt = −∇Φ(Xt)dt+ √ 2κ−1D(Xt)dWt, X0 = x0, (29)
where D(Xt) is an arbitrary, spatially-dependent, diffusion matrix, κ is a temperature constant, and x0 ∈ Rm is the particle’s initial position. The Fokker-Planck equation describes the time evolution for the probability distribution p of the particle’s position such that p(x, t) = P(Xt = x). The FP equation is the partial differential equation5,
∂tp = ∇ · ( ∇Φ(Xt)p+ κ−1∇ · (D(Xt)p) ) , p(x, 0) = δ(x0), (30)
where ∇· denotes the divergence and δ(x0) is a dirac delta distribution centered at the initialization x0. To assist in the exploration of the FP equation we define the vector field,
J(x, t) = −∇Φ(Xt)p−∇ · (D(Xt)p) , (31) which is commonly referred to as the probability current. Notice, that this gives an alternative expression for the FP equation, ∂tp = −∇·J , demonstrating that J(x, t) defines the flow of probability mass through space and time. This interpretation is especially useful for solving for the stationary solution pss, which is the unique distribution that satisfies,
∂tpss = −∇ · Jss = 0, (32) where Jss is the probability current for pss. The stationary condition can be obtained in two distinct ways:
1. Detailed balance. This is when Jss(x) = 0 for all x ∈ Ω. This is analogous to reversibility for discrete Markov chains, which implies that the probability mass flowing from a state i to any state j is the same as the probability mass flowing from state j to state i.
2. Broken detailed balance. This is when ∇ · Jss(x) = 0 but Jss(x) 6= 0 for all x ∈ Ω. This is analogous to irreversibility for discrete Markov chains, which only implies that the total probability mass flowing out of state i equals to the total probability mass flowing into state i.
The distinction between these two cases is critical for understanding the limiting dynamics of the process.

heading: E.1 THE VARIATIONAL FORMULATION OF THE FOKKER-PLANCK EQUATION WITH ISOTROPIC DIFFUSION
text: We will now consider the restricted setting of standard, isotropic diffusion (D = I). It is easy enough to check that in this setting the stationary solution is
pss(x) = e−κΦ(x)
Z , Z = ∫ Ω e−κΦ(x)dx, (33)
where pss is called a Gibbs distribution and Z is the partition function. Under this distribution, the stationary probability current is zero (Jss(x) = 0) and thus the process is in detailed balance. Interestingly, the Gibbs distribution pss has another interpretation as the unique minimizer of the the Gibbs free energy functional, F (p) = E [Φ]− κ−1H(p), (34) where E [Φ] is the expectation of the potential Φ under the distribution p and H(p) = − ∫
Ω p(x)log(p(x))dx is the Shannon entropy of p.
5This PDE is also known as the Forward Kolmogorov equation.
Proof. To prove that indeed pss is the unique minimizer of the Gibbs free energy functional, consider the following equivalent expression
F (p) = ∫ Ω p(x)Φ(x)dx+ κ−1 ∫ Ω p(x)log(p(x))dx
= κ−1 ∫
Ω
p(x) (log(p(x))− log(pss(x))) dx− κ−1 ∫
Ω
log(Z)
= κ−1DKL(p ‖ pss)− κ−1log(Z) From this expressions, it is clear that the Kullback–Leibler divergence is uniquely minimized when p = pss.
In other words, with isotropic diffusion the stationary solution pss can be thought of as the limiting distribution given by the Fokker-Planck equation or the unique minimizer of an energetic-entropic functional.
Seminal work by Jordan et al. [53] deepened this connection between the Fokker-Planck equation and the Gibbs free energy functional. In particular, their work demonstrates that the solution p(x, t) to the Fokker-Planck equation is the Wasserstein gradient flow trajectory on the Gibbs free energy functional.
Steepest descent is always defined with respect to a distance metric. For example, the update equation, xk+1 = xk − η∇Φ(xk), for classic gradient descent on a potential Φ(x), can be formulated as the solution to the minimization problem xk+1 = argminxηΦ(x) + 1 2d(x, xk)
2 where d(x, xk) = ‖x− xk‖ is the Euclidean distance metric. Gradient flow is the continuous-time limit of gradient descent where we take η → 0+. Similarly, Wasserstein gradient flow is the continuous-time limit of steepest descent optimization defined by the Wasserstein metric. The Wasserstein metric is a distance metric between probability measures defined as,
W 22 (µ1, µ2) = inf p∈Π(µ1,µ2) ∫ Rn×Rn |x− y|2p(dx, dy), (35)
where µ1 and µ2 are two probability measures on Rn with finite second moments and Π(µ1, µ2) defines the set of joint probability measures with marginals µ1 and µ2. Thus, given an initial distribution and learning rate η, we can use the Wasserstein metric to derive a sequence of distributions minimizing some functional in the sense of steepest descent. In the continuous-time limit as η → 0+ this sequence defines a continuous trajectory of probability distributions minimizing the functional. Jordan et al. [54] proved, through the following theorem, that this process applied to the Gibbs free energy functional converges to the solution to the Fokker-Planck equation with the same initialization: Theorem 1 (JKO). Given an initial condition p0 with finite second moment and an η > 0, define the iterative scheme pη with iterates defined by
pk = argminpη ( E [Φ]− κ−1H(p) ) +W 22 (p, p k−1).
As η → 0+, then pη → p weakly in L1 where p is the solution to the Fokker-Planck equation with the same initial condition.
See [54] for further explanation and [53] for a complete derivation.

heading: E.2 EXTENDING THE VARIATIONAL FORMULATION TO THE SETTING OF ANISOTROPIC DIFFUSION
text: While the JKO theorem provides a very powerful lens through which to view solutions to the FokkerPlanck equation, and thus distributions for particles governed by Langevin dynamics, it only applies in the very restricted setting of isotropic diffusion. In this section we will review work by Chaudhari and Soatto [33] extending the variational interpretation to the setting of anisotropic diffusion.
Consider when D(Xt) is an anisotropic, spatially-dependent diffusion matrix. In this setting, the original Gibbs distribution given in equation (33) does not necessarily satisfy the stationarity condition equation (32). In fact, it is not immediately clear what the stationary solution is or if the dynamics even have one. Thus, Chaudhari and Soatto [33] make the following assumption:
Stationary Assumption. Assume there exists a unique distribution pss that is the stationary solution to the Fokker-Planck equation irregardless of initial conditions.
Under this assumption we can implicitly define the potential Ψ(x) = −κ−1log(pss(x)). Using this modified potential we can express the stationary solution as a Gibbs distribution,
pss(x) ∝ e−κΨ(x). (36) Under this implicit definition we can define the stationary probability current as Jss(x) = j(x)pss(x) where j(x) = −∇Φ(x)− κ−1∇ ·D(x) +D(x)∇Ψ(x). (37) The vector field j(x) reflects the discrepancy between the original potential Φ and the modified potential Ψ according to the diffusion D(x). Notice that in the isotropic case, when D(x) = I , then Φ = Ψ and j(x) = 0. Chaudhari and Soatto [33] introduce another property of j(x) through assumption,
Conservative Assumption. Assume that the force j(x) is conservative (i.e. ∇ · j(x) = 0). Using this assumption, Chaudhari and Soatto [33] extends the variational formulation provided by the JKO theorem to the anisotropic setting, Theorem 2 (CS). Given an initial condition p0 with finite second moment, then the energeticentropic functional, F (p) = Ep [Ψ(x)]− κ−1H(p) monotonically decreases throughout the trajectory given by the solution to the Fokker-Planck equation with the given initial condition.
In other words, the Fokker-Plank equation (30) with anisotropic diffusion can be interpreted as minimizing the expectation of a modified loss Ψ, while being implicitly regularized towards distributions that maximize entropy. The derivation requires we assume a stationary solution pss exists and that the force j(x) implicitly defined by pss is conservative. However, rather than implicitly define Ψ(x) and j(x) through assumption, if we can explicitly construct a modified loss Ψ(x) such that the resulting j(x) satisfies certain conditions, then the stationary solution exists and the variational formulation will apply as well. We formalize this statement with the following theorem, Theorem 3 (Explicit Construction). If there exists a potential Ψ(x) such that either j(x) = 0 or ∇ · j(x) = 0 and ∇Ψ(x) ⊥ j(x), then pss is the Gibbs distribution ∝ e−κΨ(x) and the variational formulation given in Theorem 2 applies.

heading: E.3 APPLYING THE VARIATIONAL FORMULATION TO THE OU PROCESS
text: Through explicit construction we now seek to find analytic expressions for the modified loss Ψ(x) and force j(x) hypothesised by Chaudhari and Soatto [33] in the fundamental setting of an OU process with anisotropic diffusion, as described in section D. We assume the diffusion matrix is anisotropic, but spatially independent, ∇ · D(x) = 0. For the OU process the original potential generating the drift is Φ(x) = (x− µ)ᵀA2 (x− µ). (38) Recall, that in order to extend the variational formulation we must construct some potential Ψ(x) such that∇ · j(x) = 0 and∇Ψ ⊥ j(x). It is possible to construct Ψ(x) using the unique decomposition of the drift matrix A = (D +Q)U discussed in appendix D. Define the modified potential,
Ψ(x) = (x− µ)ᵀ U2 (x− µ). (39) Using this potential, the force j(x) is
j(x) = −A(x− µ) +DU(x− µ) = −QU(x− µ). (40) Notice that j(x) is conservative, ∇ · j(x) = ∇ · −QU (x− µ) = 0 because Q is skew-symmetric. Additionally, j(x) is orthogonal, j(x)ᵀ∇Ψ(x) = (x− µ)ᵀ UᵀQU (x− µ) = 0, again because Q is skew-symmetric. Thus, we have determined a modified potential Ψ(x) that results in a conservative orthogonal force j(x) satisfying the conditions for Theorem 3. Indeed the stationary Gibbs distribution given by Theorem 3 agrees with equation (27) derived via the first and second moments in appendix D,
e−κΨ(x) ∝ N ( µ, κ−1U−1 )
In addition to the variational formulation, this interpretation further details explicitly the stationary probability current, Jss(x) = j(x)pss, and whether or not the the stationary solution is in broken detailed balance.

heading: F EXPLICIT EXPRESSIONS FOR THE OU PROCESS GENERATED BY SGD
text: We will now consider the specific OU process generated by SGD with linear regression. Here we repeat the setup as explained in section 5.
Let X ∈ RN×d, Y ∈ RN be the input data, output labels respectively and θ ∈ Rd be our vector of regression coefficients. The least squares loss is the convex quadratic loss L(θ) = 12N ‖Y −Xθ‖2 with gradient g(θ) = Hθ − b, where H = XᵀXN and b = X ᵀY N . Plugging this expression for the gradient into the underdamped Langevin equation (3), and rearranging terms, results in the multivariate Ornstein-Uhlenbeck (OU) process,
d [ θt vt ] = A ([ µ 0 ] − [ θt vt ]) dt+ √ 2κ−1DdWt, (41)
where A and D are the drift and diffusion matrices respectively,
A =
[ 0 −I
2 η(1+β) (H + λI) 2(1−β) η(1+β)I
] , D = [ 0 0
0 2(1−β)η(1+β)Σ(θ)
] , (42)
κ = S(1− β2) is a temperature constant, and µ = (H + λI)−1b is the ridge regression solution.

heading: F.1 SOLVING FOR THE MODIFIED LOSS AND CONSERVATIVE FORCE
text: In order to apply the expressions derived for a general OU process in appendix D and E, we must first decompose the drift as A = (D + Q)U . Under the simplification Σ(θ) = σ2H discussed in appendix B, then the matrices Q and U , as defined below, achieve this,
Q =
[ 0 −σ2H
σ2H 0
] , U = [ 2 η(1+β)σ2H −1 (H + λI) 0
0 1σ2H −1
] . (43)
Using these matrices we can now derive explicit expressions for the modified loss Ψ(θ, v) and conservative force j(θ, v). First notice that the least squares loss with L2 regularization is proportional to the convex quadratic, Φ(θ) = (θ − µ)ᵀ(H + λI)(θ − µ). (44) The modified loss Ψ is composed of two terms, one that only depends on the position,
Ψθ(θ) = (θ − µ)ᵀ ( H−1(H + λI)
η(1 + β)σ2
) (θ − µ) , (45)
and another that only depends on the velocity,
Ψv(v) = v ᵀ
( H−1
σ2
) v. (46)
The conservative force j(θ, v) is
j(θ, v) =
[ v − 2η(1+β) (H + λI) (θ − µ) ] , (47)
and thus the stationary probability current is Jss(θ, v) = j(θ, v)pss.

heading: F.2 DECOMPOSING THE TRAJECTORY INTO THE EIGENBASIS OF THE HESSIAN
text: As shown in appendix D, the temporal distribution for the OU process at some time T ≥ 0 is,
pT ([ θ v ]) = N ( e−AT [ θ0 v0 ] + ( I − e−AT ) [µ 0 ] , κ−1 ( U−1 − e−ATU−1e−AᵀT )) .
Here we will now use the eigenbasis {q1, . . . , qm} of the Hessian with eigenvalues {ρ1, . . . , ρm} to derive explicit expressions for the mean and covariance of the process through time.
Deterministic component. We can rearrange the expectation as
E [[ θ v ]] = [ µ 0 ] + e−AT [ θ0 − µ v0 ] .
Notice that the second, time-dependent term is actually the solution to the system of ODEs ˙[θ v ] = −A [ θ v ] with initial condition [θ0 − µ v0]ᵀ. This system of ODEs can be block diagonalized by factorizing A = OSOᵀ where O is orthogonal and S is block diagonal defined as
O =  q1 0 . . . qm 0
. . . 0 q1 . . . 0 qm  S =  0 −1 2 η(1+β) (ρ1 + λ) 2(1−β) η(1+β) . . . . . . . . .
0 −1 2
η(1+β) (ρm + λ) 2(1−β) η(1+β)  In otherwords in the plane spanned by [qi 0] ᵀ and [0 qi] ᵀ the system of ODEs decouples into the 2D system ˙[ai bi ] = [ 0 1 − 2η(1+β) (ρi + λ) − 2(1−β) η(1+β) ] [ ai bi
] This system has a simple physical interpretation as a damped harmonic oscillator. If we let bi = ȧi, then we can unravel this system into the second order ODE
äi + 2 1− β
η(1 + β) ȧi +
2
η(1 + β) (ρi + λ)ai = 0
which is in standard form (i.e. ẍ + 2γẋ + ω2x = 0) for γ = 1−βη(1+β) and ωi = √ 2 η(1+β) (ρi + λ). Let ai(0) = 〈θ0 − µ, qi〉 and bi(0) = 〈v0, qi〉, then the solution in terms of γ and ωi is
ai(t) =  e−γt ( ai(0) cosh (√ γ2 − ω2i t ) + γai(0)+bi(0)√ γ2−ω2i sinh (√ γ2 − ω2i t )) γ > ωi e−γt(ai(0) + (γai(0) + bi(0))t) γ = ωi e−γt ( ai(0) cos (√ ω2i − γ2t ) + γai(0)+bi(0)√ ω2i−γ2 sin (√ ω2i − γ2t )) γ < ωi
Differentiating these equations gives us solutions for bi(t)
bi(t) =  e−γt ( bi(0) cosh (√ γ2 − ω2i t ) − ω 2 i ai(0)+γbi(0)√ γ2−ω2i sinh (√ γ2 − ω2i t )) γ > ωi e−γt ( bi(0)− ( ω2i ai(0) + γbi(0) ) t ) γ = ωi e−γt ( bi(0) cos (√ ω2i − γ2t ) − ω 2 i ai(0)+γbi(0)√ ω2i−γ2 sin (√ ω2i − γ2t )) γ < ωi
Combining all these results, we can now analytically decompose the expectation as the sum,
E [[ θ v ]] = [ µ 0 ] + m∑ i=1 ( ai(t) [ qi 0 ] + bi(t) [ 0 qi ]) .
Intuitively, this equation describes a damped rotation (spiral) around the OLS solution in the planes defined by the the eigenvectors of the Hessian at a rate proportional to the respective eigenvalue.
Stochastic component. Using the previous block diagonal decomposition A = OSOᵀ we can simplify the variance as
Var ([ θ v ]) = κ−1 ( U−1 − e−ATU−1e−AᵀT ) = κ−1 ( U−1 − e−OSOᵀTU−1e−OSᵀOᵀT
) = κ−1O ( OᵀU−1O − e−ST (OᵀU−1O)e−ST ᵀ ) Oᵀ
Interestingly, the matrix OᵀU−1O is also block diagonal,
OᵀU−1O = Oᵀ [ η(1+β)σ2 2 (H + λI) −1 H 0
0 σ2H
] O =  η(1+β)σ2 2 ρ1 ρ1+λ 0 0 σ2ρ1 . . . . . . . . . η(1+β)σ2
2 ρm ρm+λ 0
0 σ2ρm  Thus, similar to the mean, we can simply consider the variance in each of the planes spanned by [qi 0] ᵀ and [0 qi] ᵀ. If we define the block matrices,
Di =
[ ησ2
2S(1−β) ρi ρi+λ 0
0 σ 2
S(1−β2)ρi
] Si = [ 0 1
− 2η(1+β) (ρi + λ) − 2(1−β) η(1+β) ] then the projected variance matrix in this plane simplifies as
Var ([ qᵀi θ qᵀi v ]) = Di − e−SiTDie−SiT ᵀ
Using the solution to a damped harmonic osccilator discussed previously, we can express the matrix exponential e−SiT explicitly in terms of γ = 1−βη(1+β) and ωi =
√ 2
η(1+β) (ρi + λ). If we let αi =√ |γ2 − ω2i |, then the matrix exponential is
e−Sit =  e−γt [ cosh (αit) + γ αi sinh (αit) 1 αi sinh (αit) −ω 2 i αi sinh (αit) cosh (αit)− γαi sinh (αit) ] γ > ωi e−γt [ 1 + γt t −ω2i t 1− γt ] γ = ωi e−γt [ cos (αit) + γ αi sin (αit) 1 αi sin (αit)
−ω 2 i
αi sin (αit) cos (αit)− γαi sin (αit)
] γ < ωi

heading: G ANALYZING PROPERTIES OF THE STATIONARY SOLUTION
text: Assuming the stationary solution is given by equation (??) we can solve for the expected value of the norm of the local displacement and gain some intuition for the expected value of the norm of global displacement.
G.1 INSTANTANEOUS SPEED
Ess [ ‖δk‖2 ] = Ess [ ‖θk+1 − θk‖2 ] = η2Ess [ ‖vk+1‖2
] = η2tr ( Ess [ vk+1v ᵀ k+1
]) = η2tr (Varss (vk+1) + Ess [vk+1] Ess [vk+1] ᵀ )
= η2tr ( κ−1U−1 ) = η2
S(1− β2) tr ( σ2H ) Note that this follows directly from the definition of δk in equation (1) and the mean and variance of the stationary solution in equation ( ??), as well as the follow-up derivation in appendix F.

heading: G.2 ANOMALOUS DIFFUSION
text: Notice, that the global movement ∆t = θt−θ0 can be broken up into the sum of the local movements ∆t = ∑t i=1 δi, where δi = θi − θi−1. Applying this decomposition,
Ess [ ‖∆t‖2 ] = Ess ∣∣∣∣∣ ∣∣∣∣∣ t∑ i=1 δi ∣∣∣∣∣ ∣∣∣∣∣ 2 
= t∑ i=1 Ess [ ‖δi‖2 ] + t∑ i 6=j Ess [〈δi, δj〉]
As we solved for previously,
Ess [ ‖δi‖2 ] = η2Ess [ ‖vi‖2 ] = η2tr (Varss(vi)) =
η2 S(1− β2) tr ( σ2H ) .
By a similar simplification, we can express the second term in terms of the stationary crosscovariance, Ess [〈δi, δj〉] = η2Ess [〈vi, vj〉] = η2tr (Covss(vi, vj)) . Thus, to simplify this expression we just need to consider the velocity-velocity covariance Covss(vi, vj). At stationarity, the cross-covariance for the system in phase space, zi = [θi vi] is
Covss(zi, zj) = κ −1U−1e−A ᵀ|i−j|
where κ = S(1− β2), and
U =
[ 2 η(1+β)σ2H −1 (H + λI) 0
0 1σ2H −1
] A = [ 0 −I
2 η(1+β) (H + λI) 2(1−β) η(1+β)I ] As discussed when solving for the mean of the OU trajectory, the drift matrix A can be block diagonalized as A = OSOᵀ where O is orthogonal and S is block diagonal defined as
O =  q1 0 . . . qm 0
. . . 0 q1 . . . 0 qm  , S =  0 −1 2 η(1+β) (ρ1 + λ) 2(1−β) η(1+β) . . . . . . . . .
0 −1 2
η(1+β) (ρm + λ) 2(1−β) η(1+β)
 .
Notice also that O diagonalizes U−1 such that,
Λ = OᵀU−1O =  η(1+β)σ2 2 ρ1 ρ1+λ 0 0 σ2ρ1 . . . . . . . . . η(1+β)σ2
2 ρm ρm+λ 0
0 σ2ρm
 .
Applying these decompositions, properties of matrix exponentials, and the cyclic invariance of the trace, allows us to express the trace of the cross-covariance as
tr (Covss(zi, zj)) = κ −1tr ( U−1e−A ᵀ|i−j| )
= κ−1tr ( U−1Oe−S ᵀ|i−j|Oᵀ )
= κ−1tr ( Λe−S ᵀ|i−j| ) = κ−1
n∑ k=1 tr ( Λke −Sᵀk |i−j| ) where Λk and Sk are the blocks associated with each eigenvector of H . As solved for previously in the variance of the OU process, we can express the matrix exponential e−Sk|i−j| explicitly in terms of γ = 1−βη(1+β) and ωk = √ 2 η(1+β) (ρk + λ). If we let τ = |i − j| and αk = √ |γ2 − ω2k|, then the matrix exponential is
e−Sk|i−j| =  e−γτ [ cosh (αkτ) + γ αk sinh (αkτ) 1 αk sinh (αkτ) −ω 2 k αk sinh (αkτ) cosh (αkτ)− γαk sinh (αkτ) ] γ > ωk e−γτ [ 1 + γτ τ −ω2kτ 1− γτ ] γ = ωk e−γτ [ cos (αkτ) + γ αk sin (αkτ) 1 αk sin (αkτ)
−ω 2 k
αk sin (αkτ) cos (αkτ)− γαk sin (αkτ)
] γ < ωk
Plugging in these expressions into previous expression and restricting to just the kth velocity component, we see
Covss(vi,k, vj,k) = κ −1 [ Λke −Sᵀk |i−j| ]
1,1 =  κ−1σ2ρke −γτ ( cosh (αkτ)− γαk sinh (αkτ) ) γ > ωk κ−1σ2ρke −γτ (1− γτ) γ = ωk
κ−1σ2ρke −γτ ( cos (αkτ)− γαk sin (αkτ) ) γ < ωk
Pulling it all together,
Ess [ ‖∆t‖2 ] = η2σ2
S(1− β2)
( tr (H) t+ 2t
t∑ k=1 ( 1− k t ) m∑ l=1 ρlCl(k) )
where Cl(k) is defined as
Cl(k) =  e−γk ( cosh (αlk)− γαl sinh (αlk) ) γ > ωl e−γk (1− γk) γ = ωl e−γk ( cos (αlk)− γαl sin (αlk) ) γ < ωl
for γ = 1−βη(1+β) , ωl = √ 2 η(1+β) (ρl + λ), and αl = √ |γ2 − ω2l |.

heading: G.3 TRAINING LOSS AND EQUIPARTITION THEOREM
text: In addition to solving for the expected values of the local and global displacements, we can consider the expected training loss and find an interesting relationship to the equipartition theorem from classical statistical mechanics.
The regularized training loss is Lλ(θ) = 12 (θ − µ)ᵀH(θ − µ) + λ2 ‖θ‖2, where H is the Hessian matrix and µ is the true mean. Taking the expectation with respect to the stationary distribution,
Ess [Lλ(θ)] = 1
2 tr ((H + λI)Ess [θθ
ᵀ])− µᵀHEss [θ] + 1
2 µᵀHµ
The first and second moments of the stationary solution are
Ess [θ] = µ Ess [θθ ᵀ] =
η S(1− β) σ2 2 (H + λI)−1H + µµᵀ
Plugging these expressions in and canceling terms we get
Ess [Lλ(θ)] = η 4S(1− β) tr ( σ2H ) + λ 2 ‖µ‖2
Define the kinetic energy of the network as K(v) = 12m‖v‖2, where m = η 2 (1 + β) is the perparameter “mass" of the network according to our previously derived Langevin dynamics. At stationarity,
Ess [K(v)] = η(1 + β)
4 tr (Ess [vv
ᵀ]) = η 4S(1− β) tr ( σ2H ) where we used the fact that Ess [vvᵀ] = 1S(1−β2)σ 2H . In otherwords, at stationarity,
Ess [Lλ(θ)] = Ess [K(v)] + λ
2 ‖µ‖2.
This relationship between the expected potential and kinetic energy can be understood as a form of the equipartition theorem.

heading: H EXPERIMENTAL DETAILS

heading: H.1 COMPUTING THE HESSIAN EIGENDECOMPOSITION
text: Computing the full Hessian of the loss with respect to the parameters is computationally intractable for large models. However, equipped with an autograd engine, we can compute Hessian-vector products. We use the subspace iteration on Hessian-vector products computed on a variety of datasets. For Cifar-10 we use the entire train dataset to compute the Hessian-vector products. For Imagenet, we use a subset 40,000 images sampled from the train dataset to keep the computation within reasonable limits.
For experiments on linear regression, the Hessian is independent of the model (it only depends on the data) and can be computed using any model checkpoint. For all other experiments, the Hessian eigenvectors were computed using the model at its initial pre-trained state.

heading: H.2 FIGURE 1
text: We resumed training for a variety of ImageNet pre-trained models from Torchvision [16] for 10 epochs, with the hyperparameters used at the end of training, shown in table 1.
We kept track of the norms of the local and global displacement, ‖δk‖22 and ‖∆k‖22, every 250 steps in the training process, to keep the length of the trajectories within reasonable limits. ‖δk‖22 is visualized directly, along with its 15 step moving average. We then fit a power law of the form αkc to the ‖∆k‖22 trajectories for each model, using the last 2/3 of the saved trajectories. We visualize the ‖∆k‖22 trajectories along with their fits on a log-log plot.

heading: H.3 FIGURE 2
text: For this figure we trained a linear regression model on Cifar-10, using MSE loss on the one-hot encoded labels. The hyperparameters used during training are shown in table 2.
At every step in training the full set of model weights and velocities were stored. The top 30 eigenvectors of the hessian were computed as described in appendix H.1, using 10 subspace iterations.
The saved weight and velocity trajectories were then projected onto the top eigenvector of the hessian and were visualized in black. Using the initial weights and velocities, the red trajectories were computed according to equation (5).

heading: H.4 FIGURE 3
text: For this figure we constructed an arbitrary Ornstein-Uhlenbeck process with anisotropic noise which would help contrast the original and modified potentials. We sampled from a 2 dimensional OU process of the form
d [ x1 x2 ] = A ( b− [ x1 x2 ]) dt+ 10−4 √ DdWt.
where we set b = [−0.1, 0.05]ᵀ and arbitrarily construct A such that it’s eigenvectors are aligned with q1 = [−1, 1]ᵀ and q2 = [1, 1]ᵀ and it’s eigenvalues are 4 and 1 as follows:
D = ( 4 0 0 1 ) , V = ( −1 1 1 1 ) , A = V −1DV
The background for the left panel was computed from the convex quadratic potential Φ(x) = 1 2x
ᵀAx − bx. The background for the right panel was computed from the modified quadratic Ψ(x) = 12x
ᵀUx − ux, with U = (D + Q)−1A and u = UA−1b (see equation (25)). Both were sampled in a regular 40× 40 grid in [−0.1, 0.1]× [−0.1, 0.1].

heading: H.5 FIGURES 4, 5, AND 6
text: Starting with the ImageNet pre-trained ResNet-18 from Torchvision [16], we resumed training for 5 epochs with the hyperparameters used at the end of training, shown in table3. The top 30 Hessian eigenvectors were computed as described in appendix H.1, using 10 subspace iterations. During training, we tracked the projection of the weights and velocities onto eigenvectors q1 and q30.
For Figure 3, we show the projection of the position trajectory onto eigenvectors q1 and q30 in 2D in black. The background for the left and center panels was computed taking the ImageNet pre-trained ResNet-18 from Torchvision [16] and perturbing its weights in the q1 and q30 directions in a region close to the projected trajectory. The training and test loss were computed for a grid of 20 × 20 perturbed models. The background for the right panel was computed according to equation (10).

heading: H.6 FIGURE 7
text: We resumed training for an ImageNet pre-trained ResNet-18 from Torchvision [16] for 2 epochs, using the sweeps of hyperparameters shown in table 4. We indicate a sweep in a particular hyperparameter by [A,B], which denotes 20 evenly spaced numbers between A and B, inclusive.
We kept track of the norms of the local and global displacement, ‖δk‖22 and ‖∆k‖22, every step in the training process. The value for ‖δk‖22 at the end of the two epochs is shown in the top row of the figure. We fitted power law of the form αkc to the ‖∆k‖22 trajectories for each model on the full trajectories. The fitted exponent c for each model is plotted in the bottom row of the figure.
H.7 INCREASING RATE OF ANOMALOUS DIFFUSION
Upon further experimentation with the fitting procedure for the rate of anomalous diffusion explained in Figure 1, we observed an interesting phenomenon. The fitted exponent c for the power law relationship ‖∆k‖22 ∝ kc increases as a function of the length of the trajectory we fit to. As can
be seen in Figure 7, c increases at a diminishing rate with the length of the trajectory. This could be indicative of ‖∆k‖22 being governed by a sum of power laws where the leading term becomes dominant for longer trajectories.

Review summaries:
The paper studies the long-time dynamics of deep neural networks. The author(s) (1) show some empirical findings related to the mean square displacement, (2) model SGD as an underdamped Langevin Equation, relate it to an Ornstein Uhlenbeck process in a linear regression setting, and use it to study the limiting dynamics of SGD, (3) use the Fokker-Planck formalism to show that the steady state weight distribution obeys a modified loss which is isotropic in the absence of L2-like regularization, (4) provide empirical evidence that their findings are relevant also beyond the context of linear regression. According to the authors, these are the novelties provided by the paper:  (a) "We find empirically that long after performance has converged, networks continue to move through parameter space by a process of anomalous diffusion in which distance travelled grows as a power law in the number of gradient updates with a nontrivial exponent."  (b) "We reveal an intricate interaction between the hyperparameters of optimization, the structure in the gradient noise, and the Hessian matrix at the end of training that explains this anomalous diffusion."  (c) "we show that the key ingredient driving these dynamics is not the original training loss, but rather the combination of a modified loss, which implicitly regularizes the velocity, and probability currents, which cause oscillations in phase space"  (d) "We identify qualitative and quantitative predictions of this theory in the dynamics of a ResNet-18 model trained on ImageNet."  (e) "We uncover a mechanistic origin for the anomalous limiting dynamics of deep neural networks trained with SGD"  (f) "one of the most significant results of our analysis is that, depending on the relationship between the gradient noise covariance and the Hessian, the stationary distribution in phase space will generically violate detailed balance"  (g) 'Natural intuitions, such as “the network converges in parameter space” or “the network stays within a local region”, are wrong' 

The paper uses a stochastic differential equation approximation of stochastic gradient descent momentum to model the neural networks dynamics. At the theoretical level the study focuses on linear regression. They showed that there is agreement between the results in the simpler model and deep neural networks.  The most relevant fact is that the networks can have a super or sub diffusive behavior at the end of training. Unfortunately, their theory cannot explain this finding nor the practical consequences on the test loss.

The starting point of the present paper is the observation that the *parameters* of a neural network continue to change when training the network with stochastic gradient descent (SGD) even after the performance of the network has stabilised at a final value. The authors study these "limiting dynamics" in linear regression, which they model as an underdamped Langevin equation resulting in an Ornstein-Uhlenbeck (OU) process and find an oscillatory dynamics in the (position, velocity)-phase space. The authors discuss how to connect their findings in this simple model with the dynamics of Resnet18 (see below). Finally, they attempt to predict features of the diffusive behaviour of the ResNet18 from their model (Sec. 8)