Convergence Rates of Stochastic Gradient Descent under Infinite Noise Variance
Year: 2022

Abstract:
Recent studies have provided both empirical and theoretical evidence illustrating that heavy tails can emerge in stochastic gradient descent (SGD) in various scenarios. Such heavy tails potentially result in iterates with diverging variance, which hinders the use of conventional convergence analysis techniques that rely on the existence of the second-order moments. In this paper, we provide convergence guarantees for SGD under a state-dependent and heavy-tailed noise with a potentially infinite variance, for a class of strongly convex objectives. In the case where the p-th moment of the noise exists for some p ∈ [1, 2), we first identify a condition on the Hessian, coined ‘p-positive (semi-)definiteness’, that leads to an interesting interpolation between the positive semi-definite cone (p = 2) and the cone of diagonally dominant matrices with non-negative diagonal entries (p = 1). Under this condition, we provide a convergence rate for the distance to the global optimum in L. Furthermore, we provide a generalized central limit theorem, which shows that the properly scaled Polyak-Ruppert averaging converges weakly to a multivariate α-stable random vector. Our results indicate that even under heavy-tailed noise with infinite variance, SGD can converge to the global optimum without necessitating any modification neither to the loss function nor to the algorithm itself, as typically required in robust statistics. We demonstrate the implications of our results over misspecified models, in the presence of heavy-tailed data.

heading: 1 Introduction
text: We consider the unconstrained minimization problem
minimize x∈Rn f(x), (1.1)
using the stochastic gradient descent (SGD) algorithm. Initialized at x0 ∈ Rn, the SGD algorithm is given by the iterations,
xt+1 = xt − γt+1 ( ∇f(xt) + ξt+1(xt) ) , t = 0, 1, 2, ... (1.2)
where {γt}t∈N+ denotes the step-size sequence, and {ξt}t∈N+ is a martingale difference sequence adapted to a filtration {Ft}t∈N, characterizing the noise in the gradient (the sequence {xt}t∈N is also adapted to the same filtration, if we assume x0 is F0-measurable). Our focus is on the case where the noise is state dependent, and its variance is infinite, i.e., E [ ‖ξt‖22 ] =∞.
∗Work partially conducted while affiliated with the Vector Institute.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Many problems in modern statistical learning can be written in the form (1.1), where f(x) typically corresponds to the population risk, that is, f(x) := Ez∼ν [`(x, z)] for a given loss function ` and an unknown data distribution ν. In practice, one observes independent and identically distributed (i.i.d.) samples zi ∼ ν for i ∈ [n], and estimates the population gradient ∇f(x) with a noisy gradient at each iteration, which is based on an empirical average over a subset of the samples {zi}i∈[n]. Due to its simplicity, superior generalization performance, and well-understood theoretical guarantees, SGD has been the method of choice for minimization problems arising in statistical machine learning.
Starting from the pioneering works of Robbins and Monro [1951], Chung [1954], Sacks [1958], Fabian [1968], Ruppert [1988], Shapiro [1989], Polyak and Juditsky [1992], theoretical properties of the SGD algorithm and its variants have been receiving a growing attention under different scenarios. Recent works, for example Tripuraneni et al. [2018], Su and Zhu [2018], Duchi and Ruan [2021], Toulis and Airoldi [2017], Fang et al. [2018], Anastasiou et al. [2019], Yu et al. [2020] established convergence rates for SGD in various settings. By building on the analysis of Polyak and Juditsky [1992] to prove a central limit theorem (CLT) for the Polyak-Ruppert averaging, these works led to novel methodologies to compute confidence intervals using SGD. However, a recurring assumption in this line of work is the finite noise variance, which may be violated frequently in modern frameworks.
Heavy-tailed behavior in statistical methodology may naturally arise from the underlying model, or through the iterative optimization algorithm used during model training. In robust statistics, one often encounters heavy-tailed noise behavior in data, which in conjunction with standard loss functions leads to infinite noise variance in SGD. Very recently, heavy-tailed behavior is shown to emerge from the multiplicative noise in SGD, when the step-size is large and/or the batch-size is small [Hodgkinson and Mahoney, 2021, Gürbüzbalaban et al., 2021]. On the other hand, there is strong empirical evidence in modern machine learning that the gradient noise often exhibits a heavy-tailed behavior, which indicates an infinite variance. For example, this is observed in fully connected and convolutional neural networks [Şimşekli et al., 2019, Gürbüzbalaban and Hu, 2021] as well as recurrent neural networks [Zhang et al., 2020]. Thus, understanding the behavior of SGD under infinite noise variance becomes extremely important for at least two reasons. A computational complexity reason: modern machine learning and robust statistics frameworks lead to heavy-tailed behavior in SGD; thus, understanding the performance of this algorithm in terms of precise convergence rates as well as the required conditions on the step-size sequence as a function of the ‘heaviness’ of the tail become crucial in this setup. A statistical reason: many inference methods that rely on Polyak-Ruppert averaging utilize a CLT which holds under finite noise variance (see e.g. online bootstrap and variance estimation approaches [Fang et al., 2018, Su and Zhu, 2018, Chen et al., 2020]). Using the same methodology in the aforementioned modern frameworks (under heavy-tailed noise) will ultimately result in incorrect confidence intervals, jeopardizing the statistical procedure. Thus, establishing the limit distribution in this setting is of great importance.
In this work, we study the behavior of the SGD algorithm with diminishing step-sizes for a class of strongly convex problems when the noise variance is infinite. We establish the convergence rates of the SGD iterates towards the global minimum, and identify a sufficient condition on the Hessian of f , which interpolates between the positive semi-definite cone and the cone of diagonally dominant matrices (with non-negative diagonal entries). We further study the Polyak-Ruppert averaging of the SGD iterates, and show that the limit distribution is a multivariate α-stable distribution. We illustrate our theory on linear regression and generalized linear models, demonstrating how to verify the conditions of our theorems. Perhaps surprisingly, our results show that even under heavy-tailed noise with infinite variance, SGD with diminishing step-sizes can converge to the global optimum without requiring any modification neither to the loss function nor to the algorithm itself, as opposed to the conventional techniques used in robust statistics [Huber, 2004]. Finally, we argue that our work has potential implications in constructing confidence intervals in the infinite noise variance setting.

heading: 2 Preliminaries and Technical Background
text: Notational Conventions. By N, N+ and R we denote the set of non-negative integers, positive integers, and real numbers, respectively. For m ∈ N+, we define [m] = {1, . . . ,m}. We use italic letters (e.g. x, ξ) to denote scalars and scalar-valued functions, sign(x) to denote the sign of x, bold face italic letters (e.g. x, ξ) to denote vectors and vector-valued functions, and bold face upper case letters (e.g. A) to denote matrices. We use |x| and ‖x‖p to denote the 2-norm and p-norm of a vector x; ‖A‖ and ‖A‖p the operator 2-norm and operator p-norm of a matrix A. The transpose
of a matrix A and a vector x (viewed as a matrix with 1 column) are denoted by AT and xT. If {Ai}i∈N is a sequence of matrices and k > `, the empty product ∏` i=k Ai is understood to be the identity matrix I. For two sequences of real numbers {at}t∈N, {bt}t∈N, we write at = O(bt) if lim supt→∞ |at|/|bt| <∞, at = o(bt) if lim supt→∞ |at|/|bt| = 0, at = Θ(bt) if both at = O(bt) and bt = O(at) hold, and at bt if limt→∞ |at|/|bt| exists and is in (0,∞). If at = O(bttε) for any ε > 0, we say at = Õ(bt). Sufficiently large or sufficiently small positive constants whose values do not matter are written as C,C0, C1, . . ., sometimes without prior introduction. IfX1,X2, . . . is a sequence of random vectors taking value in Rn and µ is a probability measure on Rn, we write Xt
D−−−→ t→∞ µ if {Xt}t∈N+ converges in distribution (also called ‘weak convergence’) to µ.
Stochastic Approximation. In the SGD recursion (1.2), we can replace ∇f with an arbitrary continuous functionR : Rn → Rn, and consider the same iterations that stochastically approximate the zero x∗ ofR,
xt+1 = xt − γt+1 ( R(xt) + ξt+1(xt) ) . (2.1)
This is called the stochastic approximation process [Robbins and Monro, 1951], which is a predecessor of stochastic gradient descent (SGD) and describes a larger family of iterative algorithms (see e.g. [Kushner and Yin, 2003, Chapters 2 and 3]). Theoretical investigation of the recursion (2.1) has been active ever since its invention, especially under finite noise variance assumption: Robbins and Monro [1951] prove that the recursion (2.1) can lead to the L2 convergence limt→∞ E[|xt − x∗|2] = 0; Chung [1954] further calculates an exact convergence rate (see (3.6) in Section 3); Blum [1954] presents an elegant proof that the convergence of xt to x∗ can hold almost surely. The asymptotic distribution of (2.1) can be attributed to [Chung, 1954, Theorem 6], which states that the expression γ −1/2 t (xt − x∗) converges weakly to a normal distribution. In their seminal works, Polyak and Juditsky [1992] and Ruppert [1988] independently introduce the concept of ‘averaging the iterates’,
xt = x0 + . . .+ xt−1
t ,
showing the striking result that √ t(xt−x∗) converges weakly to a fixed normal distribution regardless of the choice of the step-size {γt}t∈N+ as long as it satisfies mild conditions. Recently, optimization algorithms that can handle heavy-tailed noise sequence {ξt}t∈N+ have been proposed [Davis et al., 2019, Nazin et al., 2019, Gorbunov et al., 2020]; however, they still rely on a uniformly bounded variance assumption, hence do not cover our setting.
Compared with the copious collection of theoretical studies on stochastic approximation algorithms with finite variance as mentioned above, papers that study stochastic approximation under infinite noise variance are extremely scarce; we shall summarize only a few papers known to us. Krasulina [1969] is the first to consider such problems, proving almost sure and Lp convergence for the one-dimensional stochastic approximation process without variance. The weak convergence of the iterates (without averaging) t1/α(xt − x∗) is also considered by Krasulina [1969], but only for the fastest-decaying step-size γt = 1/t. Goodsell and Hanson [1976] discuss how xt → x∗ in probability can imply xt → x∗ almost surely, when no finite variance is assumed, and Li [1994] provides a necessary and sufficient condition for almost sure convergence of xt → x∗, stating that faster-decaying step-size γt = o(t−1/p) is required when moments of lower orders E[|ξt|p] are not in place. Anantharam and Borkar [2012] show that although step-size that decays slower than t−1/p cannot yield almost sure convergence, Lp convergence can still hold under what they call the ‘stability assumption’, but their analysis technique provides no convergence rate. Recently, Şimşekli et al. [2019] and Zhang et al. [2020] considered SGD with heavy-tailed noise ξt having uniformly bounded p-th order moments. Besides not being able to handle state-dependent noise due to this uniform moment condition, Şimşekli et al. [2019] imposed further conditions on R = ∇f such as global Hölder continuity for a non-convex f , whereas Zhang et al. [2020] modified SGD with ‘gradient clipping’, in order to be able to compensate the effects of the heavy-tailed noise.
Finally, we shall mention that a class of stochastic recursions similar to (2.1) have been considered in the dynamical systems theory [Mirek, 2011, Buraczewski et al., 2012, 2016], for which generalized central limit theorems with α-stable limits have been established. However, such techniques typically requireR to be (asymptotically) linear and the step-sizes to be constant as they heavily rely on the theory of time-homogeneous Markov processes. Hence, their approach does not readily generalize to the setting of our interest, i.e., non-linearR and diminishing step-sizes, where the latter is crucial for ensuring convergence towards the global optimum.
Stable Distributions. In probability theory, a random variable X is stable if its distribution is non-degenerate and satisfies the following property: Let X1 and X2 be independent copies of X . Then, for any constants a, b > 0, the random variable aX1 + bX2 has the same distribution as cX+d for some constants c > 0 and d (see e.g. [Samorodnitsky and Taqqu, 1994]). The stable distribution is also referred to as the α-stable distribution, first proposed by Lévy [1937], where α ∈ (0, 2] denoting the stability parameter. The case α = 2 corresponds to the normal distribution, and the variance under this distribution is undefined for any α < 2. The multivariate α-stable distribution dates back to Feldheim [1937], which is a multivariate generalization of the univariate α-stable distribution, which is also uniquely characterized by its characteristic function. In particular, an Rn-valued random vector X has a multivariate α-stable distribution, denoted as X ∼ S(α,Λ, δ) if the joint characteristic function ofX is given by
E [ exp ( iuTX )] = exp { − ∫ s∈S2 (|uTs|α + iν(uTs, α))Λ(ds) + iuTδ } ,
for any u ∈ Rn, and 0 < α 6 2. Here, α is the tail-index, Λ is a finite measure on S2 known as the spectral measure, δ ∈ Rn is a shift vector, and ν(y, α) := − sign(y) tan(πα/2)|y|α for α 6= 1 and ν(y, α) := (2/π)y log |y| for α = 1 for any y ∈ R, and S2 denotes the unit sphere in Rn; i.e. S2 = {s ∈ Rn : |s| = 1}. Stable distribution also emerges as the limit distribution in the Generalized Central Limit Theorem (GCLT) [Gnedenko and Kolmogorov, 1954], which states that for a sequence of i.i.d. random variables whose distributions have a power-law tail with index 0 < α < 2, the normalized sum converges to an α-stable distribution as the number of summands go to∞.
Domains of Normal Attraction of Stable Distributions. LetX1,X2, . . . ,Xt be an i.i.d. sequence of random vectors in Rn with a common distribution function F (x). If there exists some constant a > 0 and a sequence bt ∈ Rn such that
X1 + · · ·+Xt at1/α − bt D−−−→ t→∞ µ, (2.2)
then F (x) is said to belong to the domain of normal attraction of the law µ, and α is the characteristic exponent of µ [Gnedenko and Kolmogorov, 1954, page 181]. If µ is an α-stable distribution, then we say F (x) belongs to the domain of normal attraction of an α-stable distribution. For example, the Pareto distribution belongs to the domain of normal attraction of an α-stable law. In Section C in the supplementary document, we provide more details as well as a sufficient and necessary conditions for being in the domain of normal attraction of an α-stable law.

heading: 3 Convergence of SGD under Heavy-tailed Gradient Noise
text: In this section, we identify sufficient conditions for the convergence of SGD under heavy-tailed gradient noise, and derive explicit rate estimates. In the standard setting when the noise variance is finite, it is sufficient to assume Hessian is uniformly positive definite in order to achieve contraction in the subsequent SGD iterations (see for example Polyak and Juditsky [1992], Tripuraneni et al. [2018], Su and Zhu [2018], Duchi and Ruan [2021], Toulis and Airoldi [2017], Fang et al. [2018], Anastasiou et al. [2019]). When the noise variance is infinite with a finite p-th moment for p ∈ [1, 2), a stronger notion of positive definiteness is required in our analysis to achieve such a contraction, which leads to an interesting interpolation between the positive semi-definite cone (as p→ 2), and the cone of diagonally dominant matrices with non-negative diagonal entries (p = 1).
3.1 p-Positive Definiteness
First, we introduce the signed power of vectors which will be used to define a family of matrices.
For a vector v = (v1, . . . , vn)T ∈ Rn and q > 0, the signed power of v is defined as v〈q〉 := ( sign ( v1 )∣∣v1∣∣q, . . . , sign(vn)|vn|q)T. (3.1)
Denoting the n-dimensional `p unit sphere with Sp = {v ∈ Rn : ‖v‖p = 1}, and the set of n× n symmetric matrices with S, we now define the following subset of S. Definition 1 (p-positive definiteness). Let p > 1 and Q be a symmetric matrix. We say that Q is p-positive definite if for all v ∈ Sp, we have vTQv〈p−1〉 > 0. Similarly, we say that Q is p-positive semi-definite if for all v ∈ Sp, we have vTQv〈p−1〉 > 0.
It is not hard to see that the set of p-positive semi-definite matrices (p-PSD) defines a closed pointed cone, which we denote by Sp+, with interior as the set of p-positive definite matrices (p-PD), denoted by Sp++. We are mainly interested in the case 1 6 p < 2. Note that S2+ coincides with the standard PSD cone, and we show in Section A.2 that S1+ is exactly the cone of diagonally dominant matrices with nonnegative diagonal entries, denoted by D+. For any p ∈ [1, 2], these cones satisfy the following
D+ = S1+ ⊆ S p + ⊆ S2+.
Figure 1 is a hypothetical illustration of the inclusion relationship between these cones.
Similar to the uniform PD condition on the Hessian, which is commonly used in classical analysis (i.e. strong convexity), we also define a uniform version of Definition 1. We recall that every operator norm ‖ · ‖p induces the same topology
on the set of n-dimensional matrices, which is just the usual topology on Rn×n. Further, the set of symmetric matrices S, as the set of zeros of the continuous function X 7→ X−XT, is a closed set. Hence for a setM ⊆ S, denoting its topological closure withM, we also haveM ⊆ S. We are interested in the case whereM is a bounded set. Definition 2 (uniform p-PD). Let p > 1 andM⊂ S be a non-empty set of symmetric matrices. We say thatM is uniformly p-PD if for all Q ∈M, we have Q ∈ Sp++.
Notice thatM is uniformly 2-PD if and only if the eigenvalues of the symmetric matrices in the set M are all lower bounded by a positive real number. Notice also that a finite subset of symmetric matrices is uniformly p-PD if and only if each element of the set is p-PD.
p-PSD cone emerges naturally when analyzing SGD algorithm in the heavy-tailed setting, interpolating between the standard PSD cone to the cone of diagonally dominant matrices with non-negative diagonal entries. To the best of our knowledge, we are the first to study such families of matrices and their application in stochastic optimization. For further details about these cones, we refer interested reader to Section A.2 in the supplementary document.
We make the following uniform smoothness and curvature assumptions on the objective function.
Assumption 1. The set of matrices {∇2f(x) : x ∈ Rn} is bounded and uniformly p-PD.

heading: 3.2 Rate of Convergence in Lp
text: We fix a probability space (Ω,F ,P) with filtration {Ft}t∈N, and we let x0 be F0-measurable. We make the following assumption on the gradient noise sequence. Assumption 2. The gradient noise sequence {ξt}t∈N+ admits the following decomposition
ξt+1(xt) = mt+1(xt) + ζt+1, (3.2)
where {ζt}t∈N+ is an i.i.d. sequence with E[ζt] = 0, and E[|ζt|p] <∞ for some p, and {mt}t∈N+ is a martingale difference sequence, and both sequences are adapted to the filtration {Ft}t∈N. Further, the state dependent component of the noise satisfies, for some K > 0,
E [ |mt+1(xt)|2 | Ft ] 6 K ( 1 + |xt|2 ) . (3.3)
We callmt the state-dependent component of the gradient noise, which naturally has a state-dependent conditional second moment. The variance of this component can be arbitrarily large depending on the state xt. The heavy-tailed noise behavior is due to ζt, which may have an infinite variance for p < 2 (i.e., the second moment is undefined). Compared to recent works on SGD with heavy-tailed noise, our noise model in Assumption 2 is significantly more general. For instance, the noise model in the recent work Zhang et al. [2020] assumes E [ |ξt+1(x)|p ] ≤ σp for all x, where σ does not depend on x. Therefore, this noise model cannot handle state-dependent noise, and does not even hold in the
linear regression example given in Section 5 (as the moments of the noise must scale with the norm of x). On the contrary, in many instances of stochastic approximation methods subject to heavy-tailed noise with long-range dependencies, one can verify that the noise admits the decomposition (3.2). We shall show in Section 5 that the noise model (3.2) arises in practical applications such as linear regression and generalized linear models subject to heavy-tailed data (see also Anantharam and Borkar [2012] for a detailed discussion on this noise model).
In our first result, assuming that the objective function f has a uniformly p-PD Hessian and the noise sequence {ξt}t∈N+ has infinite variance but satisfies Assumption 2, we establish an asymptotic convergence rate in Lp for the SGD algorithm to the unique minimizer x∗. Theorem 3. Suppose Assumptions 1 and 2 hold for some 1 < p 6 2. For step-size satisfying γt t−ρ with ρ ∈ (0, 1), the error of the SGD iterates {xt}t∈N from the minimizer x∗ satisfies
E[|xt − x∗|p] = O ( t−ρ(p−1) ) . (3.4)
Consequently, we have supt∈N+ E[|ξt|p] <∞.
The proof of Theorem 3 is provided in Section B in the supplementary document. We observe that the convergence rate of SGD depends on the highest order finite moment p of the noise sequence, and faster rates are achieved for larger values of p. The fastest convergence rate implied by our result is near O ( t−p+1 ) , which is achieved for ρ ≈ 1. However, SGD converges even for very slowly decaying step-size sequences as ρ gets closer to 0.
If the noise sequence has further integrability properties with a finite p-th moment for all p ∈ [q, α) for some 1 < q < α and if uniform p-PD condition (i.e. Assumption 1) holds, then faster rates are achievable. In particular, the following result is an interesting consequence of Theorem 3, and its proof is provided in Section B in the supplementary document. Corollary 4. For constants q, α satisfying 1 < q < α 6 2, suppose that Assumptions 1 and 2 hold for every p ∈ [q, α). For step-size satisfying γt t−ρ with ρ ∈ (0, 1), the error of the SGD iterates {xt}t∈N from the minimizer x∗ satisfies
E[|xt − x∗|q] = Õ ( t−ρq α−1 α ) . (3.5)
Remark. The additional integrability assumption yields faster rates for any feasible step-size sequence since p(α− 1)/α > p− 1 for p ∈ (1, 2].
Let us briefly compare our results stated above to those in the setting where the noise sequence has a finite variance. A classical convergence result that goes back to Chung [1954, Theorem 5]2 states that
E[|xt − x∗|r] = Θ ( t−ρr/2 ) , (3.6)
where r > 2 is an integer such that the r-th moment exists for the stochastic approximation process, and this is achieved for strongly convex objective functions in one dimension (whose second derivative {f ′′(x) : x ∈ R} satisfies the uniformly 2-PD property) with a step-size choice γt t−ρ for some ρ ∈ (1/2, 1). We point out that our rate (3.5) recovers the rate implied by (3.6) when r = 2, and extends it further to the case 1 6 r < 2.
In the presence of heavy-tailed noise, the folklore is to modify SGD (e.g., clipped gradients) in order to tame the heavy-tails, which considerably simplifies the problem and makes it amenable to classical analysis tools. For instance, to motivate modifying SGD in this regime, in [Zhang et al., 2020, Remark 1] authors prove that E [ |∇f(xt)|2 ] =∞ for vanilla SGD and argue that SGD diverges in this setting. On the contrary, our results show that, without any modifications, SGD can still converge to the optimum in Lp, even when it does not converge in L2 since the second moment is not defined.

heading: 4 Stable Limits for the Polyak-Ruppert Averaging
text: In this section, we establish the limit distribution of the Polyak-Ruppert averaging under infinite noise variance, extending the asymptotic normality result given by Polyak and Juditsky [1992] to α-stable distributions. Let us fix an α ∈ (1, 2] and assume the following throughout this subsection.
2This result, like many other similar studies in the 1950s, concerns only the one-dimensional case. But they generalize easily to higher dimensions.
Assumption 3. The gradient noise sequence {ξt}t∈N+ admits the following decomposition
ξt+1(xt) = mt+1(xt) + ζt+1, (4.1)
where {ζt}t∈N+ is an i.i.d. sequence with E[ζt] = 0, which is in the domain of normal attraction of an n-dimensional symmetric α-stable distribution µ, i.e.,
ζ1 + . . .+ ζt t1/α D−−−→ t→∞ µ. (4.2)
The state dependent component {mt}t∈N+ is a martingale difference sequence with a second-moment satisfying (3.3), and both sequences are adapted to the filtration {Ft}t∈N.
The above assumption is arguably more stringent than Assumption 2, and it implies that E[|ζt|p] <∞ for all p ∈ [1, α); thus, Assumption 2 is satisfied for any p ∈ [1, α). The condition (4.2) is a special case of (2.2) which holds if, for example, ζt’s have power-law tails (i.e. Paretian tail) with index α.
Denoting the Polyak-Ruppert averaging by xt := 1t (x0 + ... + xt−1), we are interested in the asymptotic behavior of
t1−1/α(xt − x∗) = (x0 + . . .+ xt−1)− tx∗
t1/α ,
for α ∈ (1, 2]. In the special case when α = 2, it is known that this limit converges to a multivariate normal distribution (which is a 2-stable distribution), a result proven in the seminal work by Polyak and Juditsky [1992]. Similarly, we begin with a result that considers a quadratic objective where the function ∇f(x) is linear in x, and then building on this result, we establish the limit distribution of Polyak-Ruppert averaging also in the more general non-linear case. Theorem 5 (linear case). Suppose the function ∇f(x) is affine, i.e. ∇f(x) = Ax− b for a real matrix A ∈ Rn×n and a real vector b ∈ Rn and there exist scalars p, ρ satisfying
max
( α+ αρ
1 + αρ , αρ
) 6 p 6 α,
such that A is p-PD and ρ ∈ (0, 1). If the noise sequence satisfies Assumption 3 for the same parameter α, then for the step-size satisfying γt t−ρ, the normalized average t1−1/α(xt − x∗) converges weakly to an n-dimensional α-stable distribution.
The above theorem states that Polyak-Ruppert averaging for any step-size sequence with index ρ ∈ (0, 1] converges weakly to an α-stable limit. Thus, in the linear case, the size of this feasible interval is the same in both heavy- and light-tailed noise settings (see e.g. Polyak and Juditsky [1992] and Ruppert [1988]). Notably, α-stable limit of the averaged iterates does not depend on the index ρ, i.e., limit distribution does not depend on how fast the step-size decays as long as ρ ∈ (0, 1]. The next result generalizes Theorem 5 to the setting where ∇f(x) is non-linear. Theorem 6 (non-linear case). Let 1 < 1/ρ < q < α and suppose Assumption 1 holds for every p ∈ [q, α). Assume further that the gradient ∇f(x) can be approximated using the Hessian matrix ∇2f(x∗) around the minimizer x∗ as∣∣∇f(x)−∇2f(x∗)(x− x∗)∣∣ 6 K|x− x∗|q. (4.3)

heading: If the noise sequence satisfies Assumption 3, for the step-size satisfying γt t−ρ, the normalized
text: average t1−1/α(xt − x∗) converges weakly to an n-dimensional α-stable distribution.
The condition (4.3) is standard (see e.g. Polyak and Juditsky [1992, Assumption 3.2]), which simply imposes a local linearity condition on the gradient of the objective function f , with an order-q polynomial error term. This assumption holds, for example, whenever the Hessian is Lipschitz continuous. We notice that the size of the feasible interval is ρ ∈ (1/α, 1), which is smaller this time compared to the light tailed case; Polyak and Juditsky [1992, Theorem 2] allows ρ ∈ (1/2, 1). The above theorem establishes that, when the noise has diverging variance, the Polyak-Ruppert averaging admits an α-stable limit rather than a standard CLT. This result has potential implications in statistical inference in the presence of heavy-tailed data. Inference procedures that take into account the computational part of the training procedure (instead of drawing conclusions for the minimizer
of the empirical risk) rely typically on variations of Polyak-Ruppert averaging and the CLT they admit [Fang et al., 2018, Su and Zhu, 2018, Chen et al., 2020]. The above theorem simply states this CLT does not hold under heavy-tailed gradient noise. Therefore, many of these procedures require further adaptation, if the gradient has undefined variance. Finally, it is well-known that Polyak-Ruppert averaging achieves the Cramér-Rao lower bound [Polyak and Juditsky, 1992, Gadat and Panloup, 2017], which is a lower bound on the variance of an unbiased estimator. However, it is not clear what this type of optimality means when the variance is not defined. These are important directions that require thorough investigations, and they will be studied elsewhere.

heading: 5 Examples in the Presence of Heavy-tailed Noise
text: In this section, we demonstrate how the stochastic approximation framework discussed in our paper covers several interesting examples. More specifically, we verify the assumptions required for Theorems 3, 5, and 6, for linear regression and generalized linear models (GLMs), where the heavy-tailed noise behavior may naturally arise due to heavy-tailed data.

heading: 5.1 Ordinary Least Squares
text: Let us first consider the following linear model,
y = zTβ0 + ,
where β0 ∈ Rn is the true coefficients, y ∈ R is the response, the random vector z ∈ Rn denotes the covariates with a positive-definite covariance 0 ≺ E[zzT] and a finite fourth moment E[|z|4] <∞, and is the noise with zero conditional mean E[ |z] = 0. In the classical setting, the noise is assumed to be Gaussian whose variance is well-defined. In this case, the population version of the maximum likelihood estimation (MLE) problem corresponds to minimizing
f(x) = 1 2 E [( y − zTx )2] , (5.1)
(where the expectation is over the (y,z) pair), or equivalently solving the following normal equations ∇f(x) := E [ zzT ] x− E[zy] = 0. (5.2)
We easily observe that the true coefficients β0 is the unique zero of the above equation, i.e., x ∗ = β0.
Now, suppose we are given access to a stream of i.i.d. drawn instances of the pair (y,z), denoted by {yt, zt}t∈N+ . In large-scale settings, one generally runs the following stochastic approximation process, which is simply online SGD on the population MLE objective f(x):
xt = xt−1 − γt ( ztz T t xt−1 − ztyt ) . (5.3)
Manifestly, (5.3) is a special case of (2.1), where the gradient noise admitting the decomposition ξt = ζt +mt, for an i.i.d. component ζt and a state-dependent componentmt (cf. (4.1)),{
ζt := E[zy]− ztyt, mt := ( ztz T t − E [ zzT ]) xt−1.
(5.4)
In the presence of heavy-tailed noise, for example when the noise has infinite variance, the population MLE objective f(x) may not be even finite, and one should resort to methods from M-estimation and choose a suitable loss function within robust statistics framework [Huber, 2004, Van der Vaart, 2000]. However, the iterations (5.3) may still be employed to estimate the true coefficients β0 (potentially due to model misspecification), as we demonstrate below. Note that in this case, iterations (5.3) should be seen as solving the root-finding problem (5.2) via stochastic approximation (2.1), rather than a minimization problem stated in (5.1).
First, we verify Assumption 2. We observe from the decomposition given in (5.4) that the i.i.d. component {ζt}t∈N exhibits the heavy-tailed behavior since it contains yt = zTt β0 + t. Assume that this component has the highest order finite moment p satisfying 1 6 p < 2, i.e., E[|ζt|p] <∞. Further, the state dependent component mt defines a martingale difference sequence, and the condition (3.3) is met since the covariates z have finite fourth moment, i.e.,
E [ |mt|2 | xt−1 ] 6 C|xt−1|2.
Hence, Assumption 2 is satisfied. Next, assuming that the second moment of the covariates ∇2f(x) = E[zzT] is p-PD, one can guarantee that Assumption 1 is satisfied. Therefore, the convergence results of our theorems hold. Finally, we note that p-PD assumption is always satisfied if E[zzT] is diagonally dominant, but the condition for p > 1 is weaker.

heading: 5.2 Generalized Linear Models
text: Generalized linear models (GLMs) play a crucial role in numerous problems in statistics, and provide a miscellaneous framework for many regression and classification tasks, with many applications [McCullagh and Nelder, 1989, Nelder and Wedderburn, 1972]. In this section, we consider minimizing the objective function arising from GLMs, for which there are many methods available (see e.g. Erdogdu [2015, 2016] and the references therein). However, we restrict ourselves to the misspecified and online setting. That is, the minimization problem corresponds to a GLM, but the model is misspecified so that the response can be heavy-tailed.
For a response y ∈ R and random covariates z ∈ Rn, the population version of an `2-regularized MLE problem in the canonical GLM framework reads
minimize x
f(x) := E [ ψ ( xTz ) − yxTz ] + λ
2 |x|2 for λ > 0. (5.5)
Here, ψ : R→ R is referred to as the cumulant generating function (CGF), and it is convex. Notable examples include ψ(x) = x2/2 yielding linear regression, ψ(x) = log(1 + ex) yielding logistic regression, and ψ(x) = ex yielding Poisson regression. Gradient of the objective (5.5) is given by
∇f(x) = E [ zψ′ ( zTx )] − E[zy] + λx. (5.6)
We define the unique solution of the population GLM problem as the unique zero of (5.6), which we denote by x∗. To reiterate, we do not assume a model on data, allowing for model misspecification; we simply consider the resulting optimization problem similar to Erdogdu et al. [2016, 2019]. As in the previous section, we assume that the covariates have finite fourth moment and the response yt is contaminated with heavy-tailed noise with infinite variance. In this setting, the objective function is always well-defined, even if the response has infinite variance.
We are given access to a stream of i.i.d. drawn instances of the pair (y,z), denoted by {yt, zt}t∈N+ , and we solve the above non-linear problem using the following stochastic process,
xt = xt−1 − γt ( ztψ ′(zTt xt−1)− ztyt + λxt−1), with gradient noise admitting the decomposition ξt = ζt +mt where{
ζt := E[zy]− ztyt, mt := ztψ ′(zTt xt−1)− E[ztψ′(zTt xt−1)]. In what follows, we verify our assumptions for a CGF satisfying |ψ′(x)| 6 C(1 + |x|), ψ′′(x) > 0, and |ψ′′′(x)| ≤ L for all x ∈ R. These assumptions can be easily verified for any second-order smooth CGF that grows at most quadratically (e.g. if the misspecified model is binomial with k number of trials, we have ψ(x) = k log(1 + ex)). ζt’s are i.i.d. and contain the entire heavy-tailed part of the gradient noise. Assume that this component has the highest defined moment order 1 6 p < 2, i.e., E[|ζt|p] < ∞. Further observe that the state dependent component defines a martingale difference sequence and satisfies the condition (3.3) since the covariates z have finite fourth moment, and |ψ′| grows at most linearly. Therefore, Assumption 2 is satisfied. We note that the Hessian of the objective f is given as
∇2f(x) = E [ zzTψ′′ ( zTx )] + λI.
Since ψ′′(x) > 0, ∇2f(x) is clearly PD for all λ > 0. For sufficiently large λ, this matrix can also be made diagonally dominant, which implies the p-PD condition for any p > 1, further implying Assumption 1. We further note that if ∇2f(x) is Lipshitz (e.g. for the binomail CGF), then (4.3) holds for q = 2 globally; thus it holds for any q < 2 locally. Therefore, for an appropriate step-size sequence, our convergence results on the SGD can be applied to this framework.

heading: 6 Conclusion
text: In this paper, we considered SGD subject to state-dependent and heavy-tailed noise with potentially infinite variance, when the objective belongs to a class of strongly convex functions (termed as p-PD condition). We provided asymptotic Lp convergence rates for vanilla SGD, demonstrating that SGD without any modifications can be still used in the presence of heavy-tailed noise. Furthermore, we provided a generalized central limit theorem for the Polyak-Ruppert averaging, i.e., we proved that the averaged iterates converge to a multivariate α-stable distribution.
We emphasize that p-PD condition is a sufficient condition, and further investigation is needed to see if this condition can be replaced with the standard strong convexity assumption. We also highlight that non-asymptotic Lp rates in the current setting should be achievable, which will be studied elsewhere. Finally, while we leave it for a future study, we emphasized the importance of adapting existing statistical inference techniques that rely on the averaged SGD iterates when the gradient noise is heavy-tailed, which arises naturally in modern statistical learning applications.

heading: Acknowledgements
text: MAE is partially funded by CIFAR AI Chairs program, and CIFAR AI Catalyst grant, NSERC Grant [2019-06167]. MG’s research is supported in part by the grants Office of Naval Research Award Number N00014-21-1-2244, National Science Foundation (NSF) CCF-1814888, NSF DMS2053485, NSF DMS-1723085. LZ is grateful to the partial support from NSF DMS-2053454 and a Simons Foundation Collaboration Grant. UŞ’s research is supported by the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

Review summaries:
This paper concerns the convergence rate of SGD with heavy-tail noise. In particular, the considered gradient noise consists of a finite-variance part (corresponds to the multiplicative noise in SGD for least square) and an infinite-variance part (corresponds to the additive noise in SGD for least square). The main results can be viewed as a law of large number and a central limit theorem in the absence of finite-variance. First the authors show that when the loss is "strongly convex" and "smooth" in the $p$-PD sense and the heavy-tail part noise has finite $p$-moment, SGD (with decaying stepsize) iterates converges to the correct optimum, and the convergence rate is justified in the p-norm sense. Moreover, if the heavy-tail noise is the domain of normal attraction of some stable distribution, then the properly normalized SGD iterates also converges to the stable distribution. Finally, the heavy-tailed noise model is justified from two linear model examples.

The authors analyze SGD for a class of strongly-convex functions and under (possibly) infinite noise variance, i.e., the pth moment of the noise exists for $p \in [1,2)$, but noise variance is not assumed to be bounded. Such noise behavior is defined as heavy-tailed. Additionally, the noise is assumed to be state-dependent such that the noise vector has a component whose variance is bounded by a function of the query point (iterates). For standard analysis of SGD, noise vectors are assumed to be independent of query points and of each other. Specifically, for a class of strongly convex functions (denoted as p-positive definite (p-PD), p=2 being positive definite Hessian) and under heavy-tailed noise, the authors:  -	Quantify the convergence rate with respect to (w.r.t.) distance to unique optimum in $L^p$, which approaches the optimal rate in strongly convex setting as $p \rightarrow 2$. This result is achieved for SGD without any modifications and under state-dependent noise model. -	Prove that normalized Polyak-Ruppert average (uniform average) converges weakly to an $\alpha$-stable distribution, under additional assumptions on the noise and the Hessian. Central Limit Theorem (CLT) does not hold under heavy-tailed noise model.  In order to reinforce the setting they consider and verifiability of their assumptions in practice, authors discuss about the real-world examples that admit heavy-tailed noise and satisfy p-PD assumption.