Year: 2021

Abstract:
The noise in stochastic gradient descent (SGD), caused by minibatch sampling, is poorly understood despite its practical importance in deep learning. This work presents the first systematic study of the SGD noise and fluctuations close to a local minimum. We first analyze the SGD noise in linear regression in detail and then derive a general formula for approximating SGD noise in different types of minima. For application, our results (1) provide insight into the stability of training a neural network, (2) suggest that a large learning rate can help generalization by introducing an implicit regularization, (3) explain why the linear learning ratebatchsize scaling law fails at a large learning rate or at a small batchsize and (4) can provide an understanding of how discrete-time nature of SGD affects the recently discovered power-law phenomenon of SGD.

heading: 1 INTRODUCTION
text: Stochastic gradient descent (SGD) is the simple and efficient optimization algorithm behind the success of deep learning (Allen-Zhu et al., 2019; Xing et al., 2018; Zhang et al., 2018; Wang et al., 2020; He and Tao, 2020; Liu et al., 2021; Simsekli et al., 2019; Wu et al., 2020). Minibatch noise, also known as the SGD noise, is the primary type of noise in the learning dynamics of neural networks. Practically, minibatch noise is unavoidable because a modern computer’s memory is limited while the size of the datasets we use is large; this demands the dataset to be split into “minibatches" for training. At the same time, using minibatch is also a recommended practice because using a smaller batch size often leads to better generalization performance (Hoffer et al., 2017). Therefore, understanding minibatch noise in SGD has been one of the primary topics in deep learning theory. Dominantly many theoretical studies take two approximations: (1) the continuous-time approximation, which takes the infinitesimal step-size limit; (2) the Hessian approximation, which assumes that the covariance matrix of the SGD noise is equal to the Hessian H . While these approximations have been shown to provide some qualitative understanding, the limitation of these approximations is not well understood. For example, it is still unsure when such approximations are valid, which hinders our capability to assess the correctness of the results obtained by approximations.
In this work, we fill this gap by deriving analytical formulae for discrete-time SGD with arbitrary learning rates and exact minibatch noise covariance. In summary, the main contributions are: (1) we derive the strength and the shape of the minibatch SGD noise in the cases where the noise for discrete-time SGD is analytically solvable; (2) we show that the SGD noise takes a different form in different kinds of minima and propose general and more accurate approximations. This work is organized as follows: Sec. 2 introduces the background. Sec. 3 discusses the related works. Sec. 4 outlines our theoretical results. Sec. 5 derives new approximation formulae for SGD noises. In Sec. 6, we show how our results can provide practical and theoretical insights to problems relevant to contemporary machine learning research. For reference, the relationship of this work to the previous works is shown in Table 1.

heading: 2 BACKGROUND
text: In this section, we introduce the minibatch SGD algorithm. Let {xi, yi}Ni=1 be a training set. We can define the gradient descent (GD) algorithm for a differentiable loss function L as wt =wt−1 − λ∇wL(w,{x,y}), where λ is the learning rate and w ∈ RD is the weights of the model. We consider an additive loss function for applying the minibatch SGD. Definition 1. A loss function L({xi, yi}Ni=1,w) is additive if L({xi, yi}Ni=1,w) = 1 N ∑ N i=1 `(xi, yi,w) for some differentiable, non-negative function `(⋅).
This definition is quite general. Most commonly studied and used loss functions are additive, e.g., the mean-square error (MSE) and cross-entropy loss. For an additive loss, the minibatch SGD with momentum algorithm can be defined. Definition 2. The minibatch SGD with momentum algorithm by sampling with replacement computes the update to the parameter w with the following set of equations:
⎧⎪⎪⎪⎪ ⎨ ⎪⎪⎪⎪⎩ ĝt = 1 S ∑i∈Bt ∇`(xi, yi,wt−1); mt = µmt−1 + ĝt; wt =wt−1 − λmt.
(1)
where µ ∈ [0,1) is the momentum hyperparameter, S ∶= ∣Bt∣ is the minibatch size, and the set Bt = {i1, ...iS} are S i.i.d. random integers sampled uniformly from [1,N].
One can decompose the gradient into a deterministic plus a stochastic term. Note that EB[ĝt] = ∇L is equal to the gradient for the GD algorithm. We use EB(⋅) to denote the expectation over batches, and use Ew(⋅) to denote the expectation over the stationary distribution of the model parameters. Therefore, we can write ĝt = EB[ĝt] + ηt, where ηt ∶= 1S ∑i∈Bt ∇`(xi, yi,wt−1) − EB[ĝt] is the noise term; the noise covariance is C(wt) ∶= cov(ηt, ηt). Of central importance to us is the averaged asymptotic noise covariance C ∶= limt→∞Ewt[C(wt)]. Also, we consider the asymptotic model fluctuation Σ ∶= limt→∞ cov(wt,wt). Σ gives the strength and shape of the fluctuation of w around a local minimum and is another quantity of central importance to this work. Throughout this work, C is called the “noise" and Σ the “fluctuation".

heading: 3 RELATED WORKS
text: Noise and Fluctuation in SGD. Deep learning models are trained with SGD and its variants. To understand the parameter distribution in deep learning, one needs to understand the stationary distribution of SGD (Mandt et al., 2017). Sato and Nakagawa (2014) describes the stationary distribution of stochastic gradient Langevin dynamics using discrete-time Fokker-Planck equation. Yaida (2019) connects the covariance of parameter Σ to that of the noise C through the fluctuationdissipation theorem. When Σ is known, one may obtain by Laplace approximation the stationary distribution of the model parameter around a local minimum w∗ as N (w∗,Σ). Therefore, knowing Σ can be of great practical use. For example, it has been used to estimate the local minimum escape efficiency (Zhu et al., 2019; Liu et al., 2021) and argue that SGD prefers a flatter minimum; it can also be used to assess parameter uncertainty and prediction uncertainty when a Bayesian prior is specified (Mandt et al., 2017; Gal and Ghahramani, 2016; Pearce et al., 2020). Empirically, both the fluctuation and the noise are known to crucially affect the generalization of a deep neural network. Wu et al. (2020) shows that the strength and shape of the Σ due to the minibatch noise lead to better generalization of neural networks in comparison to an artificially constructed noise.
Hessian Approximation of the Minibatch Noise. However, it is not yet known what form C and Σ actually take for SGD in a realistic learning setting. Early attempts assume isotropic noise in the continuous-time limit (Sato and Nakagawa, 2014; Mandt et al., 2017). In this setting, the noise is an isotropic Gaussian with C ∼ ID, and Σ is known to be proportional to the inverse Hessian H−1. More recently, the importance of noise structure was realized (Hoffer et al., 2017; Jastrzebski et al., 2018; Zhu et al., 2019; HaoChen et al., 2020). “Hessian approximation", which assumes C ≈ c0H for some constant c0, has often been adopted for understanding SGD (see Table 1); this assumption is often motivated by the fact that C = Jw ≈H , where Jw is the Fisher information matrix (FIM) (Zhu et al., 2019); the fluctuation can be solved to be isotropic: Σ ∼ ID. However, it is not known under what conditions the Hessian approximation is valid, while previous works have argued that it can be very inaccurate (Martens, 2014; Liu et al., 2021; Thomas et al., 2020; Kunstner et al., 2019). However, Martens (2014) and Kunstner et al. (2019) only focuses on the natural gradient descent (NGD) setting; Thomas et al. (2020) is closest to ours, but it does not apply to the case with momentum, a matrix learning rate, or regularization.
Discrete-time SGD with a Large Learning Rate. Recently, it has been realized that networks trained at a large learning rate have a dramatically better performance than networks trained with a vanishing learning rate (lazy training) (Chizat and Bach, 2018). Lewkowycz et al. (2020) shows that there is a qualitative difference between the lazy training regime and the large learning rate regime; the performance features two plateaus in testing accuracy in the two regimes, with the large learning rate regime performing much better. However, the theory regarding discrete-time SGD at a large learning rate is almost non-existent, and it is also not known what Σ may be when the learning rate is non-vanishing. Our work also sheds light on the behavior of SGD at a large learning rate. Some other works also consider discrete-time SGD in a similar setting (Fontaine et al., 2021; Dieuleveut et al., 2020; Toulis et al., 2017), but the focus is not on deriving analytical formulae or does not deal with the stationary distribution.

heading: 4 SGD NOISE AND FLUCTUATION IN LINEAR REGRESSION
text: This section derives the shape and strength of SGD noise and fluctuation for linear regression; concurrent to our work, Kunin et al. (2021) also studies the same problem but with continuous-time approximation; our result is thus more general. To emphasize the message, we discuss the label noise case in more detail. The other situations also deserve detailed analysis; we delay such discussion to the appendix due to space constraints. Notation: S denotes the minibatch size. w ∈ RD is the model parameter viewed in a vectorized form; λ ∈ R+ denotes a scalar learning rate; when the learning rate takes the form of a preconditioning matrix, we use Λ ∈ RD×D. A ∈ RD×D denotes the covariance matrix of the input data. When a matrix X is positive semi-definite, we write X ≥ 0; throughout, we require Λ ≥ 0. γ ∈ R denotes the weight decay hyperparameter; when the weight decay hyperparameter is a matrix, we write Γ ∈ RD×D. µ is the momentum hyperparameter in SGD. For two matrices X,Y , the commutator is defined as [X,Y ] ∶= XY − Y X . Other notations are introduced in the context.1 The results of this section are numerically verified in Sec. A.

heading: 4.1 KEY PREVIOUS RESULTS
text: When N ≫ S, the following proposition is well-known and gives the exact noise due to minibatch sampling. See Appendix E.1 for derivation. Proposition 1. The noise covariance of SGD as defined in Definition 2 is
C(w) = 1
SN
N
∑ i=1
∇`i(w)∇`i(w) T −
1 S ∇L(w)∇L(w)T, (2)
where the notations `i(w) ∶= l(xi, yi,w) and L(w) ∶= L({xi, yi}Ni=1,w) are used.
This gradient covariance matrix C is crucial to understand the minibatch noise. The standard literature often assumes C(w) ≈H(w); however, the following well-known proposition shows that this approximation can easily break down. Proposition 2. Let w∗ be the solution such that L(w∗) = 0, then C(w∗) = 0.
Proof. Because `i is non-negative for all i, L(w∗) = 0 implies that `i(w∗) = 0. The differentiability in turn implies that each ∇`i(w∗) = 0; therefore, C = 0. ◻
This proposition implies that there is no noise if our model can achieve zero training loss (which is achievable for an overparametrized model). This already suggests that the Hessian approximation C ∼ H is wrong since the Hessian is unlikely to vanish in any minimum. The fact that the noise strength vanishes at L = 0 suggests that the SGD noise might at least be proportional to L(w), which we will show to be true for many cases. The following theorem relates C and Σ of the discrete-time SGD algorithm with momentum for a matrix learning rate. Theorem 1. (Liu et al., 2021) Consider running SGD on a quadratic loss function with Hessian H , learning rate matrix Λ, momentum µ. Assuming ergodicity, then
(1 − µ)(ΛHΣ +ΣHΛ) − 1 + µ2
1 − µ2 ΛHΣHΛ +
µ
1 − µ2 (ΛHΛHΣ +ΣHΛHΛ) = ΛCΛ. (3)
Propostion 1 and Theorem 1 allow one to solve C and Σ. Equation (3) can be seen as a general form of the Lyapunov equation (Lyapunov, 1992) and is hard to solve in general (Hammarling, 1982; Ye et al., 1998; Simoncini, 2016). Solving this analytical equation in settings of machine learning relevance is one of the main technical contributions of this work.
1We use the word global minimum to refer to the global minimum of the loss function, i.e., where L = 0 and a local minimum refers to a minimum that has a non-negative loss, i.e., L ≥ 0.

heading: 4.2 RANDOM NOISE IN THE LABEL
text: We first consider the case when the labels contain noise. The loss function takes the form
L(w) = 1
2N
N
∑ i=1
(wTxi − yi) 2, (4)
where xi ∈ RD are drawn from a zero-mean Gaussian distribution with feature covariance A ∶= EB[xxT], and yi = uTxi + i, for some fixed u and i ∈ R is drawn from a distribution with zero mean and finite second momentum σ2. We redefine w − u→w and let N →∞ with D held fixed. The following lemma finds C as a function of Σ.
Lemma 1. (Covariance matrix for SGD noise in the label) Let the model be updated according to Eq. (1) with random noise in the label while N →∞. Then,
C = 1
S (AΣA +Tr[AΣ]A + σ2A). (5)
The model fluctuation can be obtained using this lemma.
Theorem 2. (Fluctuation of model parameters with random noise in the label) Let the assumptions be the same as in Proposition 1 and [Λ,A] = 0. Then,
Σ = σ2
S (1 +
κµ
S )ΛG−1µ , (6)
where κµ ∶= Tr[ΛAG−1µ ]
1− 1STr[ΛAG−1µ ] with Gµ ∶= 2(1 − µ)ID − (
1−µ 1+µ + 1 S )ΛA.
Remark. This result is numerically validated in Appendix A. The subscript µ refers to momentum. To obtain results for vanilla SGD, one can set µ = 0, which has the effect of reducing Gµ → G = 2ID − (1 + 1 S )ΛA. From now on, we focus on the case when µ = 0 for notational simplicity, but we note that the results for momentum can be likewise studied. The condition [Λ,A] = 0 is not peculiar because this condition holds for a scalar learning rate and common second-order methods such as Newton’s method.
If σ2 = 0, then Σ = 0. This means that when there is no label noise, the model parameter has a vanishing stationary fluctuation, which corroborates Proposition 2. When a scalar learning rate λ≪ 1 and 1 ≪ S, we have
Σ ≈ λσ2
2S ID, (7)
which is the result one would expect from the continuous-time theory with the Hessian approximation (Liu et al., 2021; Xie et al., 2021; Zhu et al., 2019), except for a correction factor of σ2. Therefore, a Hessian approximation fails to account for the randomness in the data of strength σ2. We provide a systematic and detailed comparison with the Hessian approximation in Table 2 of Sec. B.
Moreover, it is worth comparing the exact result in Theorem 2 with Eq. (7) in the regime of nonvanishing learning rate and small batch size. One notices two differences: (1) an anisotropic enhancement, appearing in the matrix Gµ and taking the form −λ(1 + 1/S)A; compared with the result in Liu et al. (2021), this term is due to the compound effect of using a large learning rate and a small batchsize; (2) an isotropic enhancement term κ, which causes the overall magnitude of fluctuations to increase; this term does not appear in the previous works that are based on the Hessian approximation and is due to the minibatch sampling process alone. As the numerical example in Sec. A shows, at large batch size, the discrete-time nature of SGD is the leading source of fluctuation; at small batch size, the isotropic enhancement becomes the dominant source of fluctuation. Therefore, the minibatch sampling process causes two different kinds of enhancement to the fluctuation, potentially increasing the exploration power of SGD at initialization but reducing the convergence speed.
Now, combining Theorem 2 and Lemma 1, one can obtain an explicit form of the noise covariance.
Theorem 3. The noise covariance matrix of minibatch SGD with random noise in the label is
C = σ2
S A +
σ2 S2 (1 + κµ S )(ΛAG−1µ +Tr[ΛAG −1 µ ]ID)A. (8)
By definition, C = J is the FIM. The Hessian approximation, in sharp contrast, can only account for the term in orange. A significant modification containing both anisotropic and isotropic (up to Hessian) is required to fully understand SGD noise, even in this simple example. Additionally, comparing this result with the training loss (127), one can find that the noise covariance contains one term that is proportional to the training loss. In fact, we will derive in Sec. 5 that containing a term proportional to training loss is a general feature of the SGD noise. We also study the case when the input is contaminated with noise. Interestingly, the result is the same with the label noise case with σ2 replaced by a more complicated term of the form Tr[AK−1BU]. We thus omit this part from the main text. A detailed discussion can be found in Appendix E.3.1. In the next section, we study the effect of regularization on SGD noise and fluctuation.

heading: 4.3 LEARNING WITH REGULARIZATION
text: Now, we show that regularization also causes a unique SGD noise. The loss function for Γ − L2 regularized linear regression is
LΓ(w) = 1
2N
N
∑ i=1
[(w − u)Txi] 2 +
1 2 wTΓw = 1 2 (w − u)TA(w − u) + 1 2 wTΓw, (9)
where Γ is a symmetric matrix; conventionally, one set Γ = γID with a scalar γ > 0. For conciseness, we assume that there is no noise in the label, namely yi = uTxi with a constant vector u. One important quantity in this case will be uuT ∶= U . The noise for this form of regularization can be calculated but takes a complicated form. Proposition 3. (Noise covariance matrix for learning with L2 regularization) Let the algorithm be updated according to Eq. (1) with L2 regularization with N →∞ and [A,Γ] = 0. Then,
C = 1 S (AΣA +Tr[AΣ]A +Tr[Γ′TAΓ′U]A + ΓA′UA′Γ) , (10)
where A′ ∶=K−1A, Γ′ ∶=K−1Γ with K ∶= A + Γ.
Notice that the last term ΓA′UA′Γ in C is unique to the regularization-based noise: it is rank-1 because U is rank-1. This term is due to the mismatch between the regularization and the minimum of the original loss. Also, note that the term Tr[AΣ] is proportional to the training loss. Define the test loss to be Ltest ∶= limt→∞Ewt[ 12(wt − u)
TA(wt − u)], we can prove the following theorem. We will show that one intriguing feature of discrete-time SGD is that the weight decay can be negative. Theorem 4. (Test loss and model fluctuation for L2 regularization) Let the assumptions be the same as in Proposition 3. Then
Ltest = λ 2S (Tr[AK−2Γ2U]κ + r) + 1 2 Tr[AK−2Γ2U], (11)
where κ ∶= Tr[A 2K−1G−1] 1− λSTr[A2K−1G−1] , r ∶= Tr[A 3K−3Γ2G−1U] 1− λSTr[A2K−1G−1] , with G ∶= 2ID − λ (K + 1SK −1A2). Moreover, let [Γ, U] = 0, then
Σ = λ S Tr[AK−2Γ2U] (1 + λκ S )AK−1G−1 + λ S (A2K−2Γ2U + λr S A)K−1G−1. (12)
This result is numerically validated in Appendix A. The test loss (11) has an interesting consequence. One can show that there exist situations where the optimal Γ is negative.2 When discussing the test loss, we make the convention that if wt diverges, then Ltest =∞. Corollary 1. Let γ∗ = arg minγ Ltest. There exist a, λ and S such that γ∗ < 0.
The proof shows that when the learning rate is sufficiently large, only negative weight decay is allowed. This agrees with the argument in Liu et al. (2021) that discrete-time SGD introduces an implicit L2 regularization that favors small norm solutions. A too-large learning rate requires a negative weight decay because a large learning rate already over-regularizes the model and one needs
2Some readers might argue that discussing test loss is meaningless when N →∞; however, this criticism does not apply because the size of the training set is not the only factor that affects generalization. In fact, this section’s crucial message is that using a large learning rate affects the generalization by regularizing the model and, if one over-regularizes, one needs to offset this effect.
to introduce an explicit negative weight decay to offset this over-regularization effect of SGD. This is a piece of direct evidence that using a large learning rate can help regularize the models. This result relates to the implicit regularization effect of SGD. It has been hypothesized that the dynamics of SGD implicitly regularizes neural networks such that the training favors simpler solutions (Kalimeris et al., 2019). Our result suggests one new mechanism for such a regularization.

heading: 5 NOISE STRUCTURE FOR GENERIC SETTINGS
text: The results in the previous sections suggest that (1) the SGD noises differ for different kinds of situations, and (2) SGD noise contains a term proportional to the training loss in general. These two facts motivate us to derive the noise covariance differently for different kinds of minima. Let f(w, x) denote the output of the model for a given input x ∈ RD. Here, we consider a more general case; f(w, x) may be any differentiable function, e.g., a non-linear deep neural network. The number of parameters in the model is denoted by P , and hence w ∈ RP . For a training dataset {xi, yi}i=1,2,...,N , the loss function with a L2 regularization is given by
LΓ(w) = L0(w) + 1
2 wTΓw, (13)
where L0(w) = 1N ∑ N i=1 `(f(w, xi), yi) is the loss function without regularization, and H0 is the Hessian of L0. We focus on the MSE loss `(f(w, xi), yi) = [f(w, xi) − yi]2/2. Our result crucially relies on the following two assumptions, which relate to the conditions of different kinds of local minima. Assumption 1. (Fluctuation decays with batch size) Σ is proportional to S−1, i.e. Σ = O(S−1).
This is justified by the results in all the related works (Liu et al., 2021; Xie et al., 2021; Meng et al., 2020; Mori et al., 2021), where Σ is found to be O(S−1). Assumption 2. (Weak homogeneity) ∣L − `i∣ is small; in particular, it is of order o(L).
This assumption amounts to assuming that the current training loss L reflects the actual level of approximation for each data point well. In fact, since L ≥ 0, one can easily show that ∣L− `i∣ = O(L). Here, we require a slightly stronger condition for a more clean expression, when ∣L − `i∣ = O(L) we can still get a similar expression but with some constant that hinders the clarity. Relaxing this condition can be an important and interesting future work. The above two conditions allow us to state our general theorem formally. Theorem 5. Let the training loss be LΓ = L0 + 12w
TΓw and the models be optimized with SGD in the neighborhood of a local minimum w∗. Then,
C(w) = 2L0(w)
S H0(w) −
1 S ∇LΓ(w)∇LΓ(w) T + o(L0). (14)
The noise takes different forms for different kinds of local minima. Corollary 2. Omitting the terms of order o(L0), when Γ ≠ 0,
C = 2L0(w ∗) S H0(w ∗ ) − 1 S Γw∗w∗TΓ +O(S−2) +O(∣w −w∗∣2). (15)
When Γ = 0 and L0(w∗) ≠ 0,
C = 2L0(w ∗) S H0(w ∗ ) +O(S−2) +O(∣w −w∗∣2). (16)
When Γ = 0 and L0(w∗) = 0,
C = 1
S (Tr[H0(w
∗ )Σ]ID −H0(w ∗ )Σ)H0(w ∗ ) +O(S−2) +O(∣w −w∗∣2). (17)
Remark. Assumption 2 can be replaced by a weaker but more technical assumption called the “decoupling assumption", which has been used in recent works to derive the continuous-time distribution of SGD (Mori et al., 2021; Wojtowytsch, 2021). The Hessian approximation was invoked in most of the literature without considering the conditions of its applicability (Jastrzebski et al., 2018; Zhu et al., 2019; Liu et al., 2021; Wu et al., 2020; Xie et al., 2021). Our result does provide such
conditions for applicability. As indicated by the two assumptions, this theorem is applicable when the batch size is not too small and when the local minimum has a loss close to 0. The reason for the failure of the Hessian approximation is that, while the FIM is equal to the expected Hessian J = E[H], there is no reason to expect the expected Hessian to be close to the actual Hessian of the minimum.
The proof is given in Appendix C. Two crucial messages this corollary delivers are (1) the SGD noise is different in strength and shape in different kinds of local minima and that they need to be analyzed differently; (2) the SGD noise contains a term that is proportional to the training loss L0 in general. Recently, it has been experimentally demonstrated that the SGD noise is indeed proportional to the training loss in realistic deep neural network settings, both when the loss function is MSE and cross-entropy (Mori et al., 2021); our result offers a theoretical justification. The previous works all treat all the minima as if the noise is similar (Jastrzebski et al., 2018; Zhu et al., 2019; Liu et al., 2021; Wu et al., 2020; Xie et al., 2021), which can lead to inaccurate or even incorrect understanding. For example, Theorem 3.2 in Xie et al. (2021) predicts a high escape probability from a sharp local or global minimum. However, this is incorrect because a model at a global minimum has zero probability of escaping due to a vanishing gradient. In contrast, the escape rate results derived in Mori et al. (2021) correctly differentiate the local and global minima. We also note that these general formulae are consistent with the exact solutions we obtained in the previous section than the Hessian approximation. For example, the dependence of the noise strength on the training loss in Theorem 2, and the rank-1 noise of regularization are all well-reflected in these formulae. In contrast, the simple Hessian approximation misses these crucial distinctions. Lastly, combining Theorem 5 with Theorem 1, one can also find the fluctuation. Corollary 3. Let the noise be as in Theorem 5, and omit the terms of order O(S−2) and O(∣w −w∗∣2). Then, when Γ ≠ 0 and when Λ, H0(w∗) and Γ commute with each other, Pr′Σ = 1 S Λ 1−µ(2L0H0 − Γw ∗w∗TΓ)(H0 + Γ)+ [2ID − Λ1+µ(H0 + Γ)] −1 . When Γ = 0 and L0(w∗) ≠ 0,
PrΣ = 2L0 S(1−µ)PrΛ (2ID − Λ 1+µH0) −1
. When Γ = 0 and L0(w∗) = 0, PrΣ = 0. Here the superscript + is the Moore-Penrose pseudo inverse, Pr ∶= diag(1, . . . ,1,0, . . . ,0) is the projection operator with r non-zero entries, r ≤ D is the rank of the Hessian H0, and r′ ≤ D is the rank of H0 + Γ. For the null space H0, Σ can be arbitrary.

heading: 6 APPLICATIONS
text: One major advantage of analytical solutions is that they can be applied in a simple “plug-in" manner by the practitioners or theorists to analyze new problems they encounter. In this section, we briefly outline a few examples where the proposed theories can be relevant.

heading: 6.1 HIGH-DIMENSIONAL REGRESSION
text: We first apply our result to the high-dimensional regression problem and show how over-andunderparametrization might play a role in determining the minibatch noise. Here, we take N,D →∞ with the ratio α ∶= N/D held fixed. The loss function is L(w) = 1
2N ∑ N i=1 (w Txi − yi) 2 . As in the
standard literature (Hastie et al., 2019), we assume the existence of label noise: yi = uTxi + i, with Var[ i] = σ2. A key difference between our setting and the standard high-dimensional setting is that, in the standard setting (Hastie et al., 2019), one uses the GD algorithm with vanishing learning rate λ instead of the minibatch SGD algorithm with a non-vanishing learning rate. Tackling the high-dimensional regression problem with non-vanishing λ and a minibatch noise is another main technical contribution of this work. In this setting, we can obtain the following result on the noise covariance matrix. Proposition 4. Let Â = 1
N ∑ N i xix T i and suppose assumptions 1 and 2 hold. With fixed S, λ, then
C = 1 S (Tr[ÂΣ]ID − ÂΣ) Â +max{0,
σ2 S (1 − 1 α )} Â.
We note that this proposition follows from Theorem 5, showing an important theoretical application of our general theory. An interesting observation is that one Σ-independent term proportional to σ2 emerges in the underparametrized regime (α > 1). However, for the overparametrized regime, the noise is completely dependent on Σ, which is a sign that the stationary solution has no fluctuation. This shows that the degree of underparametrization also plays a distinctive role in the fluctuation. In fact, one can prove the following theorem, which is verified in Sec. A.2.
Theorem 6. When a stationary solution exists for w, we have Tr[ÂΣ] = max{0, λσ 2
S (1 − 1 α ) κ̂},
where κ̂ ∶= Tr[Ĝ −1Â]
1− λSTr[Ĝ−1Â] with Ĝ ∶= 2ID − λ (1 − 1S ) Â.

heading: 6.2 IMPLICATION FOR NEURAL NETWORK TRAINING
text: It is commonly believed that the high-dimensional linear regression problem can be a minimal model for deep learning. Taking this stance, Theorem 6 suggests a technique for training neural networks. For SGD to converge, a positive semi-definite Σ must exist; however, Σ ≥ 0 if and only if κ̂ ≥ 0. From κ̂ > 0, we have ∑Di=1 1 2/λai−1+1/S < S, where ai are the eigenvalues of Â. This means that each summand should have the order of D/S. Thus the upper bound of λ should have the order of 2S/aD, where a is the typical value of ai’s. One implication of the dependence on the dimension is that the stability of a neural network trained with SGD may strongly depend on its width d, and one may rescale the learning rate according to the width to stabilize neural network training. See Figure 1-Left and Middle. We train a two-layer tanh neural network on MNIST and plot the variance of its training loss in the first epoch with fixed λ = 0.5. We see that, when d ≥ 200, the training starts to destabilize, and the training loss begins to fluctuate dramatically. When rescaling the learning rate by 1/d, we see that the variance of the training loss is successfully kept roughly constant across all d. This suggests a training technique worth being explored by practitioners in the field. In Figure 1-Middle, we also use Adam for training the same network and find a similar stabilizing trick to work for Adam.

heading: 6.3 A NATURAL LEARNING EXAMPLE WITH NEGATIVE WEIGHT DECAY
text: Sec. 4.3 shows that a too-large learning rate introduces an effective L2 regularization that can only be corrected by setting the weight decay to be negative. This effect can be observed in more realistic learning settings. We train a logistic regressor on the MNIST dataset with a large learning rate (of order O(1)). Figure 1-Right confirms that, at a large learning rate, the optimal weight decay can indeed be optimal. This agrees with our argument that using a large learning rate can effectively regularize the training.

heading: 6.4 SECOND-ORDER METHODS
text: Understanding stochastic second-order methods (including the adaptive gradient methods) is also important for deep learning (Agarwal et al., 2017; Zhiyi and Ziyin, 2021; Martens, 2014; Kunstner et al., 2019). In this section, we apply our theory to two standard second-order methods: damped Newton’s method (DNM) and natural gradient descent (NGD). We provide more accurate results than those derived in Liu et al. (2021). The derivations are given in Appendix D.2. For DNM, the preconditioning learning rate matrix is defined as Λ ∶= λA−1. The model fluctuation is shown to be proportional to the inverse of the Hessian: Σ = λσ 2
gS−λDA −1, where g ∶= 2(1 − µ) − ( 1−µ 1+µ + 1 S )λ.
The main difference with the previous results is that the fluctuation now depends explicitly on the dimension D, and implies a stability condition: S ≥ λD/g, corroborating the stability condition we derived above. For NGD, the preconditioning matrix is defined by the inverse of the Fisher information that Λ ∶= λ
S J(w)−1 = λ S C−1. We show that Σ = λ 2 ( 1 1+D 1 1+µ + 1 1−µ 1 S )A−1 is one
solution when σ = 0, which also contains a correction related to D compared to the result in Liu et al. (2021) which is Σ = λ
2 ( 1 1+µ + 1 1−µ 1 S )A−1. A consequence is that J ∼ Σ−1. The surprising
fact is that the stability of both NGD and DNM now crucially depends on D; combining with the results in Sec. 6.1, this suggests that the dimension of the problem may crucially affect the stability and performance of the minibatch-based algorithms. This result also implies that some features we derived are shared across many algorithms that depend on minibatch noise and that our results may be relevant to a broad class of optimization algorithms other than SGD.
6.5 FAILURE OF THE λ − S SCALING LAW One well-known technique in deep learning training is that one can scale λ linearly as one increases the batch size S to achieve high-efficiency training without hindering the generalization performance; however, it is known that this scaling law fails when the learning rate is too large, or the batch size is too small (Goyal et al., 2017). In Hoffer et al. (2017), this scaling law is established on the ground that Σ ∼ λ/S. However, our result in Theorem 2 suggests the reason for the failure even for the simple setting of linear regression. Recall that the exact Σ takes the form:
Σ = λσ2
S (1 +
κµ
S )G−1µ
for a scalar λ. One notices that the leading term is indeed proportional to λ/S. However, the discrete-time SGD results in a second-order correction in S, and the term proportional to 1/S2 does not contain a corresponding λ; this explains the failure of the scaling law in small S, where the second-order contribution of S becomes significant. To understand the failure at large λ, we need to look at the term Gµ:
Gµ = 2(1 − µ)ID − (λ 1 − µ 1 + µ + λ S )A.
One notices that the second term contains a part that only depends on λ but not on S. This part is negligible compared to the first term when λ is small; however, it becomes significant as the second term approaches the first term. Therefore, increasing λ changes this part of the fluctuation, and thus the scaling law no more holds if λ is large.
6.6 POWER LAW TAIL IN DISCRETE-TIME SGD
It has recently been discovered that the SGD noise causes a heavytail distribution (Simsekli et al., 2019; 2020), with a tail decaying like a power law with tail index β (Hodgkinson and Mahoney, 2020). In continuous-time, the stationary distribution has been found to obey a Student’s t-like distribution, p(w) ∼ L−(1+β)/2 ∼ (σ2 + aw2) −(1+β)/2 (Meng et al., 2020; Mori et al., 2021; Wojtowytsch, 2021). However, this result is only established for continuous-time approximations to SGD and one does not know what affects the exponent β for discrete-time SGD. Our result in Theorem 2 can serve as a tool to find the discrete-time correction to the tail index of the stationary distribution. In Appendix D.3, we show that the tail index of discrete-time SGD in 1d can be estimated as β(λ,S) = 2S
aλ − S. A clear discrete-time contribution is −(S + 1)
which depends only on the batch size, while 2S aλ + 1 is the tail index in the continuous-time limit (Mori et al., 2021). See Figure 2; the proposed formula agrees with the experiment. Knowing the tail index β is important for understanding the SGD dynamics because β is equal to the smallest moment of w that diverges. For example, when β ≤ 4, then the kurtosis of w diverges, and one expects to see outliers of w very often during training; when β ≤ 2, then the second moment of w diverges, and one does not expect w to converge in the minimum under consideration. Our result suggests that the discrete-time dynamics always leads to a heavier tail than the continuous-time theory expects, and therefore is more unstable.

heading: 7 OUTLOOK
text: In this work, we have presented a systematic analysis with a focus on exactly solvable results to promote our fundamental understanding of SGD, and our results are shown to provide a novel and accurate understanding of how SGD works for problems relevant to deep learning. One major limitation is that we have only focused on studying the asymptotic behavior of SGD, and one major future step is to investigate the dynamical, non-asymptotic noise and fluctuations of SGD.

heading: A EXPERIMENTS

heading: A.1 LABEL NOISE AND REGULARIZATION
text: Theorem 2 can be verified empirically. We run 1d experiment in Figure 4(a) and high dimensional experiments in Figures 4(b)-(c), where we choose D = 2 for visualization. We see that the continuous Hessian approximation fails badly for both large and small batch sizes. When the batch size is large, both the discrete-time Hessian approximation and our solution give a accurate estimate of the shape and the spread of the distribution. This suggests that when the batch size is large, discreteness is the determining factor of the fluctuation. When the batch size is small, the discrete Hessian approximation severely underestimates the strength of the noise. This reflects the fact that the isotropic noise enhancement is dominant at a small batch size.
In Figure 3-Left, we run a 1d experiment with λ = 1, N = 10000 and σ2 = 0.25. Comparing the predicted Σ, we see that the proposed theory agrees with the experiment across all ranges of S. The continuous theory with the Hessian approximation fails almost everywhere, while the recently proposed discrete theory with the Hessian approximation underestimates the fluctuation when S is small. In Figure 3-Right, we plot a standard case where the optimal regularization strength γ is vanishing.
Now, we validate the existence of the optimal negative weight decay as predicted by our formula. For illustration, we plot in Figure 5 the test loss (11) for a 1d example while varying either S or λ. The orange vertical lines show the place where the theory predicts a divergence in the test loss. We also
plot a standard case where the optimal γ is close to 0 in Appendix A. Also, we note that the proposed theory agrees better with the experiment.

heading: A.2 HIGH-DIMENSIONAL REGRESSION
text: See Figure 6-Left. We vary N with D = 1000 held fixed. We set λ = 0.01 and S = 32. We see that the agreement between the theory and experiment is good, even for this modest dimension number D. The vertical line shows where the over-to-underparametrization transition takes place. As expected, there is no fluctuation when α < 1, and the fluctuation gradually increases as α →∞. On the other hand, the Hessian approximation gives a wrong picture, predicting fluctuation to rise when there is no fluctuation and predicting a constant fluctuation just when the fluctuation starts to rise.

heading: B COMPARISON WITH CONVENTIONAL HESSIAN APPROXIMATION
text: We compare our results for the three cases with the results obtained with the conventional Hessian approximation of the noise covariance, i.e., C ≈H , where H is the Hessian of the loss function. We summarize the analytical results for a special case in Table 2.

heading: B.1 LABEL NOISE
text: We first consider discrete-time dynamics with the Hessian approximation. The matrix equation is
ΣA +AΣ − λAΣA = λ
S A. (18)
Compared with the exact result (3), it is a large-S limit up to the constant σ2. This constant factor is ignored during the approximation that J(w) ∶= EB[∇l∇lT] ≈ EB[∇∇Tl] ∶=H(w), which is exact only when l({xi},w) is a negative log likelihood function of w. Solving the matrix equation yields
Σ = λ
S (2ID − λA)
−1. (19)
The training loss and the test loss are
Ltrain = λ
2S Tr[A(2ID − λA)
−1 ] +
1 2 σ2, (20)
Ltest = λ
2S Tr[A(2ID − λA)
−1 ]. (21)
On the other hand, by taking the large-S limit directly from the exact equation (3), the factor σ2 is present:
ΣA +AΣ − λAΣA = λ
S σ2A. (22)
For the continuous-time limit with the Hessian approximation, the matrix equation is
ΣA +AΣ = λ
S A, (23)
which is the small-λ limit up to the factor σ2. The variance is
Σ = λ
2S ID. (24)
The training and the test error are
Ltrain = λ
4S Tr[A] +
1 2 σ2, (25)
Ltest = λ
4S Tr[A]. (26)
Again, taking the small-λ limit directly from the exact result (3) shows the presence of the factor σ2 on the right hand side of the matrix equation.
B.2 INPUT NOISE
The case with the input noise is similar to the label noise. This can be understood if we replace A by K and σ2 by Tr[AK−1BU]. The model parameter variance resulting from the discrete-time dynamics under the Hessian approximation is
Σ = λ
S (2ID − λK)
−1. (27)
The training and the test error are
Ltrain = λ
2S Tr[K(2ID − λK)
−1 ] +
1 2 Tr[AK−1BU], (28)
Ltest = λ
2S Tr[A(2ID − λK)
−1 ] +
1 2 Tr[B′TAB′U]. (29)
The large-S limit from the exact matrix equation (144) results in a prefactor Tr[AK−1BU] in the fluctuation:
Σ = λ S Tr[AK−1BU](2ID − λK)−1. (30)
For the continuous-time limit, we take λ→ 0. The Hessian approximation gives
Σ = λ
2S ID, (31)
Ltrain = λ
4S Tr[K] +
1 2 Tr[AK−1BU], (32)
Ltest = λ
4S Tr[A] +
1 2 Tr[B′TAB′U]. (33)
The large-S limit again produces a prefactor Tr[AK−1BU].

heading: B.3 L2 REGULARIZATION
text: For learning with regularization, there is a more difference between the Hessian approximation and the limit taken directly from the exact theory. We first adopt the Hessian approximation for the discrete-time dynamics. The matrix equation is
ΣK +KΣ − λKΣK = λ
S K, (34)
which is similar to the previous subsection. However, it is different from the large-S limit of the exact matrix equation (154):
ΣK +KΣ − λKΣK = λ S (Tr[AK−2Γ2U]A +AK−1ΓUΓK−1A) . (35)
This significant difference suggests that the conventional Fisher-to-Hessian approximation J ≈ H fails badly. The fluctuation, the training loss, and the test loss with the Hessian approximation are
Σ = λ
S (2ID − λK)
−1, (36)
Ltrain = λ
2S Tr[K(2ID − λK)
−1 ] +
1 2 Tr[AK−1ΓU], (37)
Ltest = λ
2S Tr[A(2ID − λK)
−1 ] +
1 2 Tr[AK−2Γ2U], (38)
while the large-S limit of the exact theory yields
Σ = λ S Tr[AK−2Γ2U]AK−1(2ID − λK)−1 + λ S A2K−3Γ2(2ID − λK)−1U, (39)
Ltrain = λ 2S Tr[AK−2Γ2U]Tr[A(2ID − λK)−1] + λ 2S Tr[A2K−2Γ2(2ID − λK)−1U]
+ 1 2 Tr[AK−1ΓU], (40)
Ltest = λ 2S Tr[AK−2Γ2U]Tr[A(2ID − λK)−1] + λ 2S Tr[A3K−3Γ2(2ID − λK)−1U]
+ 1 2 Tr[AK−2Γ2U]. (41)
The continuous-time results are obtained by taking the small-λ limit on Eqs. (36)-(38) for the Hessian approximation and on Eqs. (39)-(41) for the limiting cases of the exact theory. Specifically, for the Hessian approximation, we have
Σ = λ
2S ID, (42)
Ltrain = λ
4S Tr[K] +
1 2 Tr[AK−1ΓU], (43)
Ltest = λ
4S Tr[A] +
1 2 Tr[AK−2Γ2U]. (44)
The small-λ limit of the exact theory yields
Σ = λ 2S Tr[AK−2Γ2U]AK−1 + λ 2S A2K−3Γ2U, (45)
Ltrain = λ 4S Tr[AK−2Γ2U]Tr[A] + λ 4S Tr[A2K−2Γ2U] + 1 2 Tr[AK−1ΓU], (46)
Ltest = λ 4S Tr[AK−2Γ2U]Tr[A] + λ 4S Tr[A3K−3Γ2U] + 1 2 Tr[AK−2Γ2U]. (47)

heading: C PROOF OF THE GENERAL FORMULA

heading: C.1 PROOF OF THEOREM 5 AND COROLLARY 2
text: We restate the theorem. Theorem 7. Let the training loss be LΓ = L0 + 12w
TΓw and the models be optimized with SGD in the neighborhood of a local minimum w∗. When Γ ≠ 0, the noise covariance is given by
C = 2L0(w ∗) S H0(w ∗ ) − 1 S Γw∗w∗TΓ +O(S−2) +O(∣w −w∗∣2). (48)
When Γ = 0 and L0(w∗) ≠ 0,
C = 2L0(w ∗) S H0(w ∗ ) +O(S−2) +O(∣w −w∗∣2). (49)
When Γ = 0 and L0(w∗) = 0,
C = 1
S (Tr[H0(w
∗ )Σ]ID −H0(w ∗ )Σ)H0(w ∗ ) +O(S−2) +O(∣w −w∗∣2). (50)
Proof. We will use the following shorthand notations: `i ∶= `(f(w, xi), yi), `′i ∶= ∂`i ∂f , `′′i ∶= ∂2`i ∂f2
. The Hessian of the loss function without regularization H0(w) = ∇∇TL0(w) is given by
H0(w) = 1
N
N ∑ i=1 `′′i ∇f(w, xi)∇f(w, xi) T + 1 N N ∑ i=1 `′i∇∇ Tf(w, xi). (51)
The last term of Eq. (51) can be ignored when L0 ≪ 1, since
∥ 1
N
N ∑ i=1 `′i∇∇ Tf(w, xi)∥ F ≤ ( 1 N N ∑ i=1 (`′i) 2 )
1/2 ( 1
N
N
∑ i=1
∥∇∇ Tf(w, xi)∥ 2 F)
1/2
= ⟨`′2⟩1/2 ( 1
N
N
∑ i=1
∥∇∇ Tf(w, xi)∥ 2 F)
1 2
,
= √ 2L0(w)( 1
N
N
∑ i=1
∥∇∇ Tf(w, xi)∥ 2 F)
1 2
,
where ∥ ⋅ ∥F stands for the Frobenius norm3, and we have defined the variable ⟨`′2⟩ ∶= 1N ∑ N i=1(` ′ i) 2. Since `′′i = 1 for the mean-square error, we obtain
H0(w) = 1
N
N
∑ i=1
∇f(w, xi)∇f(w, xi) T +O ( √ L0) (52)
near a minimum. The Hessian with regularization HΓ(w) = ∇∇TLΓ(w) is just given by H0(w)+Γ.
On the other hand, the SGD noise covariance C(w) is given by Eq. (2). By assumption 2, the SGD noise covariance is directly related to the Hessian:
C(w) = ⟨`′2⟩ SN N ∑ i=1 ∇f(w, xi)∇f(w, xi) T − 1 S ∇LΓ(w)∇LΓ(w) T
+ 2
SN
N
∑ i=1
(`i −L0)∇f(w, xi)∇f(w, xi) T
= 2L0(w)
S H0(w) −
1 S ∇LΓ(w)∇LΓ(w) T + o(L0). (53)
This finishes the proof. ◻
3In the linear regression problem, the last term of Eq. (51) does not exist since ∇∇Tf(w, xi) = 0.
Now we prove Corollary 2.
Proof. Near a minimum w∗ of the full loss LΓ(w), we have
∇LΓ(w) =H0(w ∗ )(w −w∗) + Γw∗ +O(∣w −w∗∣2), (54)
within the approximation LΓ(w) = LΓ(w∗)+(1/2)(w−w∗)THΓ(w∗)(w−w∗)T+O(∣w−w∗∣2). Equations (14) and (54) give the SGD noise covariance near a minimum of LΓ(w).
Now it is worth discussing two different cases separately: (1) with regularization and (2) without regularization. We first discuss the case when regularization is present. In this case, the regularization Γ is not small enough, and the SGD noise covariance is not proportional to the Hessian. Near a local or global minimum w ≈w∗, the first term of the right-hand side of Eq. (54) is negligible, and hence we obtain
Ew[C(w)] = 2L0(w ∗) S H0(w ∗ ) − 1 S Γw∗w∗TΓ
+Ew [ 1
S H0(w
∗ )(w −w∗)(w −w∗)TH0(w∗)] +O(∣w −w∗∣2)
= 2L0(w ∗) S H0(w ∗ ) − 1 S Γw∗w∗TΓ +O(S−2) +O(∣w −w∗∣2). (55)
where we have used the fact that E[w] = w∗. The SGD noise does not vanish even at a global minimum of LΓ(w). Note that this also agrees with the exact result derived in Sec. 4.3: together with an anisotropic noise that is proportional to the Hessian, a rank-1 noise proportional to the strength of the regularization appears. This rank-1 noise is a signature of regularization.
On the other hand, as we will see below, the SGD noise covariance is proportional to the Hessian near a minimum when there is no regularization, i.e., Γ = 0. We have
C(w) = 2L0(w)
S H0(w) −
1 S H0(w ∗ )(w −w∗)(w −w∗)TH0(w∗) +O(∣w −w∗∣2). (56)
For this case, we need to differentiate between a local minimum and a global minimum. When L0(w ∗) is not small enough (e.g. at a local but not global minimum),
C(w) = 2L0(w ∗) S H0(w) + (w −w∗)TH0(w∗)(w −w∗) S H0(w)
− 1
S H0(w
∗ )(w −w∗)(w −w∗)TH0(w∗) +O(∣w −w∗∣2)
= 2L0(w ∗) S H0(w) +O(S −2 ) +O(∣w −w∗∣2)
= 2L0(w ∗) S H0(w ∗ ) +O(S−2) +O(∣w −w∗∣2), (57)
and so, to leading order,
C = 2L0(w ∗) S H0(w ∗ ), (58)
which is proportional to the Hessian but also proportional to the achievable approximation error.
On the other hand, when L0(w∗) is vanishingly small (e.g. at a global minimum), we have 2L0(w) ≈ (w −w∗)TH0(w∗)(w −w∗), and thus obtain
C(w) = 1 S [(w −w∗)TH0(w∗)(w −w∗)H0(w∗) −H0(w∗)(w −w∗)(w −w∗)TH0(w∗)]
+O(S−2) +O(∣w −w∗∣2), (59)
i.e.,
E[C] = 1
S (Tr[H0Σ]ID −H0Σ)H0 +O(S
−2 ) +O(∣w −w∗∣2). (60)
This completes the proof.
Remark. It should be noted that the second term on the right-hand side of Eq. (59) would typically be much smaller than the first term for large D. For example, when H0(w∗) = aID with a > 0, the first and the second terms are respectively given by (a2/S)∥w−w∗∥2ID and −(a2/S)(w−w∗)(w−w∗)T. The Frobenius norm of the former is given by (Da2/S)∥w −w∗∥2, while that of the latter is given by (a2/S)∥w −w∗∥2, which indicates that in Eq. (59), the first term is dominant over the second term for large D. Therefore the second term of Eq. (59) can be dropped for large D, and Eq. (59) is simplified as
⎧⎪⎪ ⎨ ⎪⎪⎩ C(w) ≈ (w−w ∗)TH0(w∗)(w−w∗) S H0(w ∗); E[C] ≈ Tr[H0Σ] S H0.
(61)
Again, the SGD noise covariance is proportional to the Hessian.
In conclusion, as long as the regularization is small enough, that the SGD noise covariance near a minimum is proportional to the Hessian is a good approximation. This implies that the noise is multiplicative, which is known to lead to a heavy tail distribution (Clauset et al., 2009; Levy and Solomon, 1996). Thus, we have studied the nature of the minibatch SGD noise in three different situations. As an example, we have demonstrated the power of this general formulation by applying it to the high-dimensional linear regression problem in Sec. 6.1.

heading: C.2 PROOF OF COROLLARY 3
text: Proof. We prove the case where Γ = 0 and L(w∗) ≠ 0 as an example. Substituting Theorem 5 into Theorem 1yields
[2ID − 1
1 + µ ΛH0]ΛH0Σ = 2L0 S(1 − µ) Λ2H0, (62)
where we have assumed necessary commutation relations. Suppose that the Hessian H0 is of rank-r with r ≤D. The singular-value decomposition and its Moore-Penrose pseudo inverse are given by

heading: H0 = USV
text: T and H+0 = V S +UT, respectively, where U and V are unitary, S is a rank-r diagonal matrix with elements being singular values of H0, and S+ is obtained by inverting every non-zero entry of S. Multiplying H+0 to both sides of the above equation, we have
PrΣ = 2L0
S(1 − µ) PrΛ(2ID −
Λ
1 + µ H0)
−1 , (63)
where Pr = diag(1,1, . . . ,1,0, . . . ,0) is the projection operator with r non-zero entries. When the Hessian is full-rank, i.e., r =D, the Moore-Penrose pseudo inverse is nothing but the usual inverse. The other cases can be calculated similarly.

heading: D APPLICATIONS
text: D.1 INFINITE-DIMENSIONAL LIMIT OF THE LINEAR REGRESSION PROBLEM
Now we apply the general theory in Sec. 5 to linear regressions in the high-dimensional limit, namely N,D →∞ with α ∶= N/D held fixed.

heading: D.1.1 PROOF OF PROPOSITION 4
text: The loss function
L(w) = 1
2N
N
∑ i=1
(wTxi − yi) 2
(64)
with yi = uT + i can be written as
L(w) = 1 2 (w − u − Â+v) T Â (w − u − Â+v) − 1 2 vTÂ+v + 1 2N
N ∑ i=1 2i , (65)
where Â ∶= 1 N ∑ N i=1 xix T i is an empirical covariance for the training data and v ∶= 1 N ∑ N i=1 xi i. The symbol (⋅)+ denotes the Moore-Penrose pseudoinverse. We also introduce the the averaged traing loss: Ltrain ∶= Ew[L(w)]
The minimum of the loss function is given by
w∗ = u + Â+v +Πr, (66)
where r ∈ RD is an arbitrary vector and Π is the projection onto the null space of Â. Since 1 −Π = Â+Â, w∗ is also expressed as
w∗ = Â+(Âu + v) +Πr. (67)
In an underparameterized regime α > 1, Π = 0 almost surely holds as long as the minimum eigenvalue of A (not Â) is positive (Hastie et al., 2019). In this case, Â+ = Â−1 and we obtain
w∗ = u + Â−1v for α > 1. (68)
On the other hand, in an overparameterized regime α > 1, Π ≠ 0 and there are infinitely many global minima. In the ridgeless regression, we consider the global minimum that has the minimum norm ∥w∗∥, which corresponds to
w∗ = Â+(Âu + v) = (1 −Π)u + Â+v for ridgeless regression with α < 1. (69)
In both cases, the loss function is expressed as
L(w) = 1 2 (w −w∗)T Â (w −w∗) − 1 2 vTÂv + 1 2N
N ∑ i=1 2i . (70)
Asymptotically, wt converges to a stationary point w∗ with fluctuation Σ obeying the following equation (Theorem 1:
λÂΣ + λΣÂ − λ2ÂΣÂ = λ2C. (71)
The SGD noise covariance C is given by Eq. (14). In the present case, the Hessian is given by H = Â and we also have
1
N
N
∑ i=1
(`′i) 2 =
1
N
N
∑ i=1
(wTxi − yi) 2 =
2
N
N ∑ i=1 `i = 2L(w). (72)

heading: On the other hand, ∇L(w)∇L(w)T = Â(w−w∗)(w−w∗)TÂ, and hence Ew[∇L(w)∇L(w)T] =
text: ÂΣÂ. Therefore we obtain
C = Ew[C(w)] = 2Ltrain S Â − 1 S ÂΣÂ. (73)
Now, we find Ltrain. First, we define X ∈ RN×D as Xik = (xi)k, and ⃗ ∈ RN as ⃗i = i. Then w∗ = (1 −Π)u + Â+v = (1 −Π)u + (XTX)+XT⃗.
With this notation, we have Â =XTX/N , and the loss function is expressed as
L(w) = 1 2 (w −w∗)TÂ(w −w∗) − 1 2N ⃗TX(XTX)+XT⃗ + 1 2N
N ∑ i=1 2i . (74)
We therefore obtain
Ltrain = 1
2 Tr[ÂΣ] −
1
2N E[⃗TX(XTX)+XT⃗] +
σ2
2 . (75)
Here, E[⃗TX(XTX)+XT⃗] = σ2Tr[(XTX)(XTX)+] = σ2Tr(1 −Π). (76)
We can prove that the following identity is almost surely satisfied (Hastie et al., 2019) as long as the smallest eigenvalue of A (not Â) is positive:
Tr(1 −Π) = min{D,N}. (77)
We therefore obtain
Ltrain = 1
2 Tr[ÂΣ] −
σ2 2N min{D,N} + σ2 2 = ⎧⎪⎪⎪⎪⎪ ⎨ ⎪⎪⎪⎪⎪⎩ 1 2 Tr[ÂΣ] + 1 2 (1 − 1 α )σ2 for α > 1, 1 2 Tr[ÂΣ] for α ≤ 1
(78)
By substituting Eq. (78) into Eq. (73), we obtain the following SGD noise covariance:
C = ⎧⎪⎪⎪⎪⎪ ⎨ ⎪⎪⎪⎪⎪⎩ 1 S (Tr[ÂΣ] − ÂΣ) Â + σ2 S (1 − 1 α ) Â for α > 1, 1 S (Tr[ÂΣ] − ÂΣ) Â for α ≤ 1.
(79)
This finishes the proof. ◻

heading: D.1.2 PROOF OF THEOREM 6
text: Proof. We have to solve this equation:
ÂΣ +ΣÂ − λÂΣÂ = λC, (80)
where C is given in Proposition 4. Using the similar trick of multiplying by Ĝ ∶= 2ID − λ (1 − 1S ) Â as in Appendix E.2.2, one obtains
Tr[ÂΣ] = { λσ2 S (1 − 1 α ) κ̂ for α > 1;
0 for α ≤ 1, (81)
where κ̂ ∶= Tr[Ĝ −1Â]
1− λSTr[Ĝ−1Â] with Ĝ ∶= 2ID − λ (1 − 1S ) Â.
Substituting the above trace into the matrix equation, we have
Σ = { λσ2 S (1 − 1 α ) (1 + λ S κ̂) Ĝ−1 for α > 1;
0 for α ≤ 1. (82)

heading: D.2 SECOND-ORDER METHODS
text: Proposition 5. Suppose that we run DNM with Λ ∶= λA−1 with random noise in the label. The model fluctuation is
Σ = λσ2
gS − λD A−1, (83)
where g ∶= 2(1 − µ) − ( 1−µ 1+µ + 1 S )λ.
Proof. Substituting Λ = λA−1 into Eqs. (3) and (5) yields
gΣ = λ S (Tr[AΣ] + σ2)A−1, (84)
where g ∶== 2(1 − µ) − ( 1−µ 1+µ + 1 S )λ. Multiplying A and taking trace on both sides, we have
Tr[AΣ] = λDσ2
gS − λD . (85)
Therefore, the model fluctuation is
Σ = λσ2
gS − λD A−1. (86)
Proposition 6. Suppose that we run NGD with Λ ∶= λ S J(w)−1 ≈ λ S C−1 with random noise in the label. The model fluctuation is
Σ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ λ 4 g − 1 2 σ2 1 +D + 1 4 ¿ Á ÁÀλ2g2 + 4λ(g − 2 1 +D 1 1 + µ ) σ2 1 +D + 4( σ2 1 +D ) 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ A−1, (87)
where g ∶= 1 1+D 1 1+µ + 1 1−µ 1 S .
Proof. Similarly to the previous case, the matrix equation satisfied by Σ is
(1−µ)(C−1AΣ+ΣAC−1)− 1 + µ 2 1 − µ2 λ S C−1AΣAC−1 + µ 1 − µ2 λ S (C−1AC−1AΣ+ΣAC−1AC−1) = λ S C−1.
(88)
Although it is not obvious how to directly solve this equation, it is possible to guess one solution according to the hope that Σ be proportional to J−1, in turn, A−1 (Amari, 1998; Liu et al., 2021). We assume that Σ = xA−1 and substitute it into the above equation to solve for x. This yields one solution without claiming its uniqueness. By simple algebra, this x is solved to be
x = λ
4 g −
1
2
σ2
1 +D +
1
4
¿ Á ÁÀλ2g2 + 4λ(g − 2
1 +D
1
1 + µ )
σ2
1 +D + 4(
σ2
1 +D )
2
. (89)
Let σ = 0. We obtain the result in Sec. 6.4.

heading: D.3 ESTIMATION OF TAIL INDEX
text: In Mori et al. (2021); Meng et al. (2020), it is shown that the (1d) discrete-time SGD results in a distribution that is similar to a Student’s t-distribution:
p(w) ∼ (σ2 + aw2)− 1+β 2 , (90)
where σ2 is the degree of noise in the label, and a is the local curvature of the minimum. For large w, this distribution is a power-law distribution with tail index:
p(∣w∣) ∼ ∣w∣−(1+β), (91)
and it is not hard to check that β also equal to the smallest moment of w that diverges: E[wβ] =∞. Therefore, estimating β can be of great use both empirically and theoretically.
In continuous-time, it is found that βcts = 2Saλ + 1 (Mori et al., 2021). For discrete-time SGD, we hypothesize that the discrete-time nature causes a change in the tail index β = βcts + , and we are interested in finding . We propose a “semi-continuous" approximation to give the formula to estimate the tail index. Notice that Theorem 2 gives the variance of the discrete-time SGD, while Eq. (90) can be integrated to give another value of the variance, and the two expressions must be equal for consistency. This gives us an equation that β must satisfy:
∫ p(w;β)(w −E[w])2 = Var[w], (92)
this procedure gives the following formula:
β(λ,S) = 2S
aλ − S = βcts + , (93)
and one immediately recognizes that −(S + 1) is the discrete-time contribution to the tail index. See Figure 7 for additional experiments. We see that the proposed formula agrees with the experimentally measured value of the tail index for all ranges of the learning rate, while the result of Mori et al. (2021) is only correct when λ→ 0+. Hodgkinson and Mahoney (2020) also studies the tail exponent of discrete-time SGD; however, their conclusion is only that the “index decreases with the learning rate and increases with the batch size". In contrast, our result give the functional form of the tail index directly. In fact, this is the first work that gives any functional form for the tail index of discrete-time SGD fluctuation to the best of our knowledge.
The following proposition gives the intermediate steps in the calculation. Proposition 7. (Tail index estimation for discrete-time SGD) Let the parameter distribution be p(w) ∼ (σ2 + aw2) − 1+β2 , and Var[w] be given by Theorem 2. Then
β(λ,S) = 2S
aλ − S. (94)
Proof. The normalization factor for the distribution exists if β > 0:
N =
√ a
π
σβΓ( 1+β 2 )
Γ(β 2 )
. (95)
If β > 2, the variance exists and the value is
Var[w] = σ2
a(β − 2) . (96)
By equating Eq. (96) with the exact variance (6), we are able to solve for an expression of the tail index as
β(λ,S) = 2S
aλ − S. (97)

heading: E PROOFS AND ADDITIONAL THEORETICAL CONSIDERATIONS

heading: E.1 PROOF OF PROPOSITION 1
text: The with-replacement sampling is defined in Definition 2. Let us here define the without-replacement sampling. Definition 3. A minibatch SGD without replacement computes the update to the parameter w with the following set of equations:
{ ĝt =
1 S ∑i∈Bt ∇`(xi, yi,wt−1);
wt =wt−1 − λĝt, (98)
where S ∶= ∣Bt∣ ≤ N is the minibatch size, and the set Bt is an element uniformly-randomly drawn from the set of all S-size subsets of {1, ...,N}.
From the definition of the update rule for sampling with or without replacement, the covariance matrix of the SGD noise can be exactly derived. Proposition 8. The covariance matrices of noise in SGD due to minibatch sampling as defined in Definitions 2 and 3 with an arbitrary N are
C(w) = { 1 S [ 1 N ∑ N i=1∇`i∇` T i −∇L(w)∇L(w) T] , (with replacement) N−S
S(N−1) [ 1 N ∑ N i=1∇`i∇` T i −∇L(w)∇L(w)
T] , (without replacement) (99)
where the shorthand notation `i(w) ∶= l(xi, yi,w) is used.
In the limit of S = 1 or N ≫ S, two cases coincide. In the N ≫ S limit, both methods of sampling have the same noise covariance as stated in Proposition 1:
C(w) = 1
SN
N
∑ i=1
∇`i∇` T i −
1 S ∇L(w)∇L(w)T. (100)
Remark. We also note that a different way of defining minibatch noise exists in Hoffer et al. (2017). The difference is that our definition requires the size of each minibatch to be exactly S, while Hoffer et al. (2017) treats the batch size also as a random variable and is only expected to be S. In comparison, our definition agrees better with the common practice.
Now we prove Proposition 8.
Proof. We derive the noise covariance matrices for sampling with and without replacement. We first derive the case with replacement. According to the definition, the stochastic gradient for sampling with replacement can be rewritten as
ĝ = 1
S
N
∑ n=1
gnsn, (101)
where gn ∶= ∇`n and
sn = l, if l −multiple n ′s are sampled in S, with 0 ≤ l ≤ S. (102)
The probability of sn assuming value l is given by the multinomial distribution
P (sn = l) = ( S
l )(
1
N )
l
(1 − 1
N )
S−l . (103)
Therefore, the expectation value of sn is given by
EB[sn] = S
∑ l=0 lP (sn = l) =
S N , (104)
which gives
EB[ĝ] = g ∶= 1
N
N
∑ n=1
gn = ∇L(w). (105)
For the covariance, we first calculate the covariance between sn and sn′ . Due to the properties of the covariance of multinomial distribution, we have for n ≠ n′
EB[snsn′] = cov[sn, sn′] +E[sn]2
= − S N2 + S2 N2
= S(S − 1)
N2 ; (106)
and for n = n′
EB[snsn] = Var[sn] +E[sn]2
= S
N
N − 1 N + S2 N2
= SN + S(S − 1)
N2 . (107)
Substituting these results into the definition of the noise covariance yields
C(w) = EB[ĝĝT] −EB[ĝ]EB[ĝ]T
= 1
S2
N
∑ n=1
N
∑ n′=1
gng T n′EB[snsn′] − gg T
= 1
S2
N
∑ n,n′=1
gng T n′ S(S − 1)
N2 +
1
S2
N
∑ n=1
gng T n [
SN + S(S − 1) N2 − S(S − 1) N2 ] − ggT
= 1
NS
N
∑ n=1
gng T n −
1 S ggT
= 1
S [
1
N
N
∑ i=1
∇`i∇` T i −∇L(w)∇L(w) T ] . (108)
Then, we derive the noise covariance for sampling without replacement. Similarly, according to the definition, the stochastic gradient for sampling without replacement can be rewritten as
ĝ = 1
S
N
∑ n=1
gnsn, (109)
where
sn = { 0, if n ∉ S;
1, if n ∈ S. (110)
The probability of n that is sampled in S from N is given by
P (sn = 1) = ( N−1 S−1)
( N S )
= S
N . (111)
The expectation value of sn is then given by
EB[sn] = P (sn = 1) = S
N , (112)
which gives
EB[ĝ] = g ∶= 1
N
N
∑ n=1
gn = ∇L(w). (113)
For the covariance, we first calculate the covariance between sn and sn′ . By definition, we have for n ≠ n′
EB[snsn′] = P (sn = 1, s′n = 1) = P (sn = 1∣s ′ n = 1)P (s ′ n = 1)
= ( N−2 S−2)
( N−1 S−1)
( N−1 S−1)
( N S )
= S(S − 1)
N(N − 1) ; (114)
and for n = n′
EB[snsn] = P (sn = l) = S
N . (115)
Substituting these results into the definition of the noise covariance yields
C(w) = EB[ĝĝT] −EB[ĝ]EB[ĝ]T
= 1
S2
N
∑ n=1
N
∑ n′=1
gng T n′EB[snsn′] − gg T
= 1
S2
N
∑ n,n′=1
gng T n′ S(S − 1)
N(N − 1) +
1
S2
N
∑ n=1
gng T n [
S N − S(S − 1) N(N − 1) ] − ggT
= 1
NS
N − S
N − 1
N
∑ n=1
gng T n −
N − S
S(N − 1) ggT
= N − S
S(N − 1) [
1
N
N
∑ i=1
∇`i∇` T i −∇L(w)∇L(w) T ] . (116)

heading: E.2 PROOFS IN SEC. 4.2

heading: E.2.1 PROOF OF LEMMA 1
text: Proof. From the definition of noise covariance (2), the covariance matrix for the noise in the label is
C(w) = 1 NS
N ∑ i=1 ∇li(wt−1)∇li(wt−1)T − 1 S ∇L(wt−1)∇L(wt−1)T
= 1 S 1 N
N ∑ i (wTxi − i)xixTi (wTxi − i)T − 1 S [ 1 N N ∑ i (wTxi − i)xi] ⎡⎢⎢⎢⎣ 1 N N ∑ j xTj (wTxj − j)T ⎤⎥⎥⎥⎦
= 1 S 1 N
N ∑ i (wTxixixTi xTi w + 2ixixTi ) − 1 S [ 1 N N ∑ i (wTxixi)] ⎡⎢⎢⎢⎣ 1 N N ∑ j (xTi xTi w) ⎤⎥⎥⎥⎦
(117)
= 1 S (AwwTA +Tr[AwwT]A + σ2A), (118)
where we have invoked the law of large numbers and the expectation value of the product of four Gaussian random variables in the third line is evaluated as follows.
Because N is large, we invoke the law of large numbers to obtain the (j, k)-th component of the matrix as
lim N→∞
1
N
N
∑ i=1
(wTxixix T i x T i w)jk = EB[w TxxxTxTw]jk = EB [ D
∑ i
wixixjxk D
∑ i′ xi′wi′] . (119)
Because the average is taken with respect to x and each x is a Gaussian random variable, we apply the expression for the product of four Gaussian random variables E[x1x2x3x4] = E[x1x2]E[x3x4] + E[x1x3]E[x2x4]+E[x1x4]E[x2x3]−2E[x1]E[x2]E[x3]E[x4] (Janssen and Stoica, 1988) to obtain
EB [ D
∑ i
wixixjxk D
∑ i′ xi′wi′]
= EB [ D
∑ i
wixixj]EB [xk D
∑ i′ xi′wi′] +EB [
D
∑ i
wixixk]EB [xj D
∑ i′ xi′wi′]
+EB [ D
∑ i
wixi D
∑ i′ xi′wi′]EB [xjxk]
= 2(AwwTA)jk +Tr[Aww T ]Ajk. (120)
Writing Σ ∶= Ew[wwT], we obtain
Ew[C(w)] =∶ C = 1
S (AΣA +Tr[AΣ]A + σ2A). (121)
This method has been utilized repeatedly in this work.

heading: E.2.2 PROOF OF THEOREM 2
text: Proof. We substitute Eq. (5) into Eq. (3) which is a general solution obtained in a recent work (Liu et al., 2021):
(1 − µ)(ΛAΣ +ΣAΛ) − 1 + µ2
1 − µ2 ΛAΣAΛ +
µ
1 − µ2 (ΛAΛAΣ +ΣAΛAΛ) = ΛCΛ. (122)
To solve it, we assume the commutation relation that [Λ,A] ∶= ΛA −AΛ = 0. Therefore, the above equation can be alternatively rewritten as
[(1 − µ)ID − 1
2 (
1 − µ 1 + µ + 1 S )ΛA]ΣAΛ +ΛAΣ [(1 − µ)ID − 1 2 ( 1 − µ 1 + µ + 1 S )ΛA]
− 1
S Tr[AΣ]ΛA =
1 S σ2ΛA. (123)
To solve this equation, we first need to solve for Tr[AΣ]. Multiplying Eq. (123) by G−1µ ∶= [2(1 − µ)ID − ( 1−µ 1+µ + 1 S )ΛA] −1 and taking trace, we obtain
Tr[AΣ] − 1 S Tr[AΣ]Tr[ΛAG−1µ ] = 1 S σ2Tr[ΛAG−1µ ], (124)
which solves to give
Tr[AΣ] = σ2
S
Tr[ΛAG−1µ ]
1 − 1 S Tr[ΛAG−1µ ] ∶= σ2 S κµ. (125)
Therefore, Σ is
Σ = σ2
S (1 +
κ S )Λ [2(1 − µ)ID − ( 1 − µ 1 + µ + 1 S )ΛA]
−1 . (126)

heading: E.2.3 TRAINING ERROR AND TEST ERROR FOR LABEL NOISE
text: In the following theorem, we calculate the expected training and test loss for random noise in the label. Theorem 8. (Approximation error and test loss for SGD noise in the label) The expected approximation error, or the training loss, is defined as Ltrain ∶= Ew[L(w)]; the expected test loss is defined as Ltest ∶= 12EwEB [(w
Tx)2]. For SGD with noise in the label given by Eq. (5), the expected approximation error and test loss are
Ltrain = σ2
2 (1 +
λκ
S ) , (127)
Ltest = λσ2
2S κ. (128)
Remark. Notably, the training loss decomposes into two additive terms. The term that is proportional to 1 is the bias, caused by insufficient model expressivity to perfectly fit all the data points, while the second term that is proportional to λκ/S is the variance in the model parameter, induced by the randomness of minibatch noise. Remark. When the learning rate λ is vanishingly small, the expected test loss diminishes whereas the training error remains finite as long as label noise exists.
Proof. We first calculate the approximation error. By definition,
Ltrain ∶= Ew[L(w)]
= 1
2 Tr[AΣ] +
1 2 σ2 = 1 2
λσ2
S κ +
1 2 σ2
= σ2
2 (1 +
λκ
S ) . (129)
The test loss is
Ltest = 1
2 Ew [wTAw] =
1 2 Tr[AΣ]
= λσ2
2S κ. (130)

heading: E.3 MINIBATCH NOISE FOR RANDOM NOISE IN THE INPUT

heading: E.3.1 NOISE STRUCTURE
text: Similar to label noise, noise in the input data can also cause fluctuation. We assume that the training data points x̃i = xi+ηi can be decomposed into a signal part and a random part. As before, we assume Gaussian distributions, xi ∼ N (0,A) and ηi ∼ N (0,B). The problem remains analytically solvable if we replace the Gaussian assumption by the weaker assumption that the fourth-order moment exists and takes some matrix form. For conciseness, we assume that there is no noise in the label, namely yi = u
Txi with a constant vector u. One important quantity in this case will be uuT ∶= U . Notice that the trick w − u = w no more works, and so we write the difference explicitly here. The loss function then takes the form
L(w) = 1
2N
N
∑ i=1
[(w − u)Txi +w Tηi] 2 =
1 2 (w − u)TA(w − u) + 1 2 wTBw. (131)
The gradient ∇L = (A +B)(w − u) +Bu vanishes at w∗ ∶= (A +B)−1Au, which is the minimum of the loss function and the expectation of the parameter at convergence. It can be seen that, even at the minimum w∗, the loss function remains finite unless u = 0, which reflects the fact that in the presence of input noise, the network is not expressive enough to memorize all the information of the data. The SGD noise covariance for this type of noise is calculated in the following proposition. Proposition 9. (Covariance matrix for SGD noise in the input) Let the algorithm be updated according to Eq. (1) or (98) with random noise in the input while the limit N →∞ is taken with D held fixed. Then the noise covariance is
C = 1 S {KΣK +Tr[KΣ]K +Tr[AK−1BU]K} , (132)
where K ∶= A +B, and Σ ∶= Ew [(w −w∗)(w −w∗)T]. Remark. It can be seen that the form of the covariance (132) of input noise is similar to that of label noise (5) with replacing A by K and σ2 by Tr[AK−1BU], suggesting that these two types of noise share a similar nature.
Defining the test loss as Ltest ∶= 12EwEB [(w Tx − uTx)2], Proposition 9 can then be used to calculate the test loss and the model fluctuation. Theorem 9. (Training error, test loss and model fluctuation for noise in the input) The expected training loss is defined as Ltrain ∶= Ew[L(w)], and the expected test loss is defined as Ltest ∶= 1 2 EwEB [(wTx − uTx)2]. For SGD with noise in the input given in Proposition (9), the expected approximation error and test loss are
Ltrain = 1 2 Tr[AK−1BU] (1 + λ S κ′) , (133)
Ltest = λ 2S Tr[AK−1BU]κ′ + 1 2 Tr[B′TAB′U], (134)
where κ′ ∶= Tr[KG ′−1] 1−λ 1STr[KG′−1] withG′ ∶= 2ID−λ (1 + 1S )K, andB ′ ∶=K−1B. Moreover, let [K,U] = 0. Then the covariance matrix of model parameters is
Σ = λTr[AK−1BU]
S (1 +
λκ′
S )[2ID − λ(1 +
1 S )K]
−1 . (135)
Remark. Note that [K,U] = 0 is necessary only for an analytical expression of Σ. It can be obtained by solving Eq. (144) even without invoking [K,U] = 0. In general, the condition that [K,U] = 0 does not hold. Therefore, only the training and test error can be calculated exactly. Remark. The test loss is always smaller than or equal to the training loss because all matrices involved here are positive semidefinite.

heading: E.3.2 PROOF OF PROPOSITION 9
text: Proof. We define Σ ∶= Ew [(w −w∗)(w −w∗)T]. Then,
Ew[wwT] = Σ + (A +B)−1AUA(A +B)−1 ∶= Σ +A′UA′T ∶= ΣA, (136) Ew [(w − u)(w − u)T] = Σ +B′UB′T ∶= ΣB , (137)
where we use the shorthand notations A′ ∶= (A+B)−1A, B′ ∶= (A+B)−1B and ΣA ∶= Σ+A′UA′T, ΣB ∶= Σ +B
′UB′T. We remark that the covariance matrix Σ here still satisfies the matrix equation (3) with the Hessian being K ∶= A +B.
The noise covariance is
C(w) = 1
S
1
N
N
∑ i
(wTx̃i − u Txi)x̃ix̃ T i (w Tx̃i − u Txi) T −
1 S ∇L(w)∇L(w)T
= 1
S {A(w − u)(w − u)TA +BwwTB +A(w − u)wTB +Bw(w − u)TA
+Tr[A(w − u)(w − u)T]K +Tr[BwwT]K}. (138)
In Eq. (138), there are four terms without trace and two terms with trace. We first calculate the traceless terms. For the latter two terms, we have
Ew[(w − u)wT] = Σ −A′UB′T, (139) Ew[w(w − u)T] = Σ −B′UA′T. (140)
Because A′ +B′ = ID, after simple algebra the four traceless terms result in 2(A +B)Σ(A +B).
The two traceful terms add to Tr[AΣB +BΣA]K. With the relation AB′ = BA′, what inside the trace is
AΣB +BΣA =KΣ +AK −1BU. (141)
Therefore, the asymptotic noise is
C ∶= Ew[C(w)]
= 1
S {KΣK +Tr[AΣB +BΣA]K} (142)
= 1 S {KΣK +Tr[KΣ]K +Tr[AK−1BU]K} . (143)

heading: E.3.3 PROOF OF THEOREM 9
text: Proof. The matrix equation satisfied by Σ is
ΣK +KΣ − λ(1 + 1
S )KΣK =
λ S (Tr[KΣ]K +Tr[AK−1BU]K) . (144)
By using a similar technique as in Appendix E.2.2, the trace Tr[KΣ] can be calculated to give
Tr[KΣ] = λTr[AK−1BU]
S κ′, (145)
where κ′ ∶= Tr[KG ′−1] 1−λ 1STr[KG′−1] with G′ ∶= 2ID − λ (1 + 1S )K.
With Eq. (145), the training error and the test error can be calculated. The approximation error is
Ltrain = Ew[L(w)] = 1
2 Tr[AΣB +BΣA] =
1 2 Tr[AK−1BU] (1 + λ S κ′) . (146)
The test loss takes the form of a bias-variance tradeoff:
Ltest = 1
2 EwEB [(wTx − uTx)2] =
1 2 Ew [(w − u)TA(w − u)] = 1 2 Tr[AΣB]
= λ 2S Tr[AK−1BU] (1 +
λκ′
S )Tr[AG′−1] +
1 2 Tr[B′TAB′U]
= λ 2S Tr[AK−1BU]κ′ + 1 2 Tr[B′TAB′U]. (147)
Let [K,U] = 0. Then Σ can be explicitly solved because it is a function of K and U . Specifically,
Σ = λTr[AK−1BU]
S (1 +
λκ′
S )[2ID − λ(1 +
1 S )K]
−1 . (148)

heading: E.4 PROOFS IN SEC. 4.3

heading: E.4.1 PROOF OF PROPOSITION 3
text: Proof. The covariance matrix of the noise is
C(w) = 1
S
1
N
N
∑ i
[(w − u)Txixi + Γw] [x T i x T i (w − u) +w
TΓ] − 1
S ∇LΓ(w)∇LΓ(w)
T
= 1
S {A(w − u)(w − u)TA +Tr[A(w − u)(w − u)T]A}. (149)
Using a similar trick as in Appendix E.3.2, the asymptotic noise is
C = 1 S (AΣA +Tr[AΣ]A +Tr[Γ′TAΓ′U]A + ΓA′UA′Γ) . (150)

heading: E.4.2 PROOF OF THEOREM 4
text: Besides the test loss and the model fluctuation, we derive the approximation error here as well.
Theorem. (Training error, test loss and model fluctuation for learning with L2 regularization) The expected training loss is defined as Ltrain ∶= Ew[L(w)], and the expected test loss is defined as Ltest ∶= 1 2 EwEB [(wTx − uTx)2]. For noise induced by L2 regularization given in Proposition 3, let [A,Γ] = 0. Then the expected approximation error and test loss are
Ltrain = λ 2S Tr[AK−2Γ2U]Tr[AG−1] (1 + λκ S ) + λ 2S (Tr[A2K−2Γ2G−1U] + λr S Tr[AG−1])
+ 1 2 Tr[AK−1ΓU], (151)
Ltest = λ 2S (Tr[AK−2Γ2U]κ + r) + 1 2 Tr[AK−2Γ2U], (152)
where κ ∶= Tr[A 2K−1G−1] 1− λSTr[A2K−1G−1] , r ∶= Tr[A 3K−3Γ2G−1U] 1− λSTr[A2K−1G−1] , with G ∶= 2ID −λ (K + 1SK −1A2). Moreover, if A, Γ and U commute with each other, the model fluctuation is
Σ = λ S Tr[AK−2Γ2U] (1 + λκ S )AK−1G−1 + λ S (A2K−2Γ2U + λr S A)K−1G−1. (153)
Remark. Because Γ may not be positive semidefinite, the test loss can be larger than the training loss, which is different from the input noise case.
Proof. The matrix equation obeyed by Σ is
ΣK +KΣ − λKΣK − λ
S AΣA =
λ S Tr[AΣ]A + λ S (Tr[AK−2Γ2U]A +AK−1ΓUΓK−1A) ,
(154)
where we use the shorthand notation K ∶= A + Γ. Let [A,Γ] = 0. Using the trick in Appendix E.2.2, the trace term Tr[AΣ] is calculated as
Tr[AΣ] = λ S (Tr[AK−2Γ2U]κ + r) , (155)
where κ ∶= Tr[A 2K−1G−1] 1− λSTr[A2K−1G−1] , r ∶= Tr[A 3K−3Γ2G−1U] 1− λSTr[A2K−1G−1] , and G ∶= 2ID − λ (K + 1SK −1A2).
The training error is
Ltrain = 1
2 Tr[AΣΓ + ΓΣA]
= 1
2 Tr[KΣ] +
1 2 Tr[AK−1ΓU]
= λ 2S Tr[AK−2Γ2U]Tr[AG−1] (1 + λκ S ) + λ 2S (Tr[A2K−2Γ2G−1U] + λr S Tr[AG−1])
+ 1 2 Tr[AK−1ΓU]. (156)
The test loss is
Ltest = 1
2 EwEB [(wTx − uTx)2]
= 1
2 Ew [(w − u)TA(w − u)] =
1 2 Tr[AΣΓ]
= λ 2S (Tr[AK−2Γ2U]κ + r) + 1 2 Tr[AK−2Γ2U]. (157)
Let A, Γ and U commute with each other. Then,
Σ = λ S Tr[AK−2Γ2U] (1 + λκ S )AK−1G−1 + λ S (A2K−2Γ2U + λr S A)K−1G−1. (158)

heading: E.4.3 PROOF OF COROLLARY 1
text: For a 1d example, the training loss and the test loss have a simple form. We use lowercase letters for 1d cases. Corollary 4. For a 1d SGD with L2 regularization, the training loss and the test loss are
Ltrain = aγ
2(a + γ)
2(a + γ) − λ [(a + γ)2 + 2 S a(a − γ)]
2(a + γ) − λ [(a + γ)2 + 2 S a2]
u2, (159)
Ltest = aγ2
2(a + γ)
2 − λ(a + γ)
2(a + γ) − λ [(a + γ)2 + 2 S a2]
u2. (160)
Proof. The training error and the test loss for 1d cases can be easily obtained from Theorem E.4.2.
Now we prove Corollary 1.
Proof. The condition for convergence is 1 − λ S
Tr[A2K−1G−1] > 0. Specifically,
λ(a + γ)2 − 2(a + γ) + λ 2
S a2 < 0. (161)
For a given γ, the learning rate needs to satisfy
λ < 2(a + γ)
(a + γ)2 + 2 S a2 . (162)
For a given λ, γ needs to satisfy
1 − aλ − √
1 − 2 S a2λ2
λ < γ <
1 − aλ + √
1 − 2 S a2λ2
λ , (163)
which indicates a constraint on λ:
aλ <
√ S
2 . (164)
If γ is allowed to be non-negative, the optimal value can only be 0 due to the convergence condition. Therefore, a negative optimal γ requires an upper bound on it being negative, namely
1 − aλ + √
1 − 2 S a2λ2
λ < 0. (165)
Solving it, we have
aλ > 2
1 + 2 S
. (166)
By combining with Eq. (164), a necessary condition for the existence of a negative optimal γ is
2
1 + 2 S
<
√ S
2 → (S − 2)2 > 0→ S ≠ 2. (167)
Hence, a negative optimal γ exists, if and only if
2
1 + 2 S
< aλ <
√ S
2 , and S ≠ 2. (168)
For higher dimension with Γ = γID, it is possible to calculate the optimal γ for minimizing the test loss (11) as well. Specifically, the condition is given by
d
dγ Ltest ∶=
1
2
d
dγ
f(γ) g(γ) = 0, (169)
where
f(γ) ∶= γ2Tr [AK−2 (ID + λ S A2K−1G−1)U] , (170)
g(γ) ∶= 1 − λ S Tr[A2K−1G−1]. (171)
Although it is impossible to solve the equation analytically, it can be solved numerically.

Review summaries:
This paper investigates the importance of noise in mini-batch SGD. The main contribution of this paper is to derive an analytic solution to the shape and strength of SGD minibatch noise for linear regression with random noise in the label, linear regression with additional L2 regularization and non-linear regression given some assumption about the model fluctuation and loss homogeneity.  Those analysis reveal that 1) the SGD noise is proportional to the loss level 2) the shape of the noise differs for different loss minima. This is in contrast with previous work looking at this problem using various approximation which neglects those effect.  Authors finally highlight various applications for those findings. In particular, the analytical formulation they derive for the noise explains why the linear scaling rule of the learning rate does not work for small batch or high learning rate. 

This paper studies the properties of the gradient noise in mini-batch SGD using discrete-time analysis for a fixed learning rate. Their analysis is more general than prior work and considers SGD with momentum, a learning rate matrix (that subsumes preconditioning methods), and regularization. They provide closed form expressions for the noise covariance for various machine learning problems including linear regression with label noise (Theorem 3), with regularization (Proposition 3), and more general models under reasonable assumptions (Theorem 5). Their theory matches the empirical observations about the effect of mini-batch size and learning rate on the gradient noise for small batch sizes and large learning rates where there was previously a gap in theory. Section 6 illustrates the significance of the results as it is able to explain various observations in machine learning theory and deep learning practice as well as providing ideas for practitioners.

This paper studies the minibatch noise for discrete time SGD and discusses its approximation on different applications. The novelty of the paper stands on the derivation of minibatch noise covatiance of discrete time SGD. For special cases with label noise and L2 regularization, this work derives the exact solution of covariance. For a more general setting, the author gives the covariance replying on two assumptions. In application section, a few experiments are presented to show the strength of the theory.