An Analysis of Constant Step Size SGD in the Non-convex Regime: Asymptotic Normality and Bias
Year: 2022

Abstract:
Structured non-convex learning problems, for which critical points have favorable statistical properties, arise frequently in statistical machine learning. Algorithmic convergence and statistical estimation rates are well-understood for such problems. However, quantifying the uncertainty associated with the underlying training algorithm is not well-studied in the non-convex setting. In order to address this short-coming, in this work, we establish an asymptotic normality result for the constant step size stochastic gradient descent (SGD) algorithm—a widely used algorithm in practice. Specifically, based on the relationship between SGD and Markov Chains [1], we show that the average of SGD iterates is asymptotically normally distributed around the expected value of their unique invariant distribution, as long as the non-convex and non-smooth objective function satisfies a dissipativity property. We also characterize the bias between this expected value and the critical points of the objective function under various local regularity conditions. Together, the above two results could be leveraged to construct confidence intervals for non-convex problems that are trained using the SGD algorithm.

heading: 1 Introduction
text: Non-convex learning problems are prevalent in modern statistical machine learning applications such as matrix and tensor decomposition [2, 3, 4, 5, 6, 7, 8], deep neural networks [9, 10, 11], and robust empirical risk minimization [12, 13, 14]. Developing theoretically principled approaches for tackling such non-convex problems depends critically on the interplay between two aspects. From a computational perspective, variants of stochastic gradient descent (SGD) converge to firstorder critical points [15, 16] or local minimizers [17, 2, 18, 19] of the objective function. From a statistical perspective, oftentimes these critical points or local minimizers have nice statistical properties [20, 3, 12, 21, 22, 5]; see also [23] for a counterexample. For the purpose of uncertainty quantification in such non-convex settings, studying the fluctuations of iterative algorithms used for training becomes extremely important. In this work, we focus on the widely used constant step size SGD, and develop results for quantifying the uncertainty associated with this algorithm for a class of non-convex problems.
We consider minimizing a non-smooth and non-convex objective function f : Rd→ R,
min θ∈Rd f(θ) . (1)
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
The iterations of SGD with a constant step size η > 0, initialized at θ(η)0 ≡ θ0 ∈ Rd, are given by
θ (η) k+1 = θ (η) k − η ( ∇f(θ(η)k ) + ξk+1(θ (η) k ) ) , k ≥ 0 , (2)
where {ξk}k≥1 is a sequence of random functions from Rd to Rd corresponding to the stochasticity in the gradient estimate. Several problems in machine learning and statistics are naturally formulated as the optimization problem in (1), where the function f(θ) is given by
f(θ) := ∫ F (θ, Z) dP (Z) , (3)
for an unknown distribution over the random variable Z ∈ Rp. The function F (θ, Z) is typically the loss function composed with functions from the hypothesis class parametrized by θ ∈ Rd. In online SGD with batch size b, at each iteration k, b independent samples Zj ∼ P (Z) are used to estimate the true gradient with 1b ∑b j=1∇F (θ (η) k , Zj). The above iterates are indeed a special case of the iterates in (2), with the noise sequence {ξk+1(θ(η)k )}k≥0 given by
ξk+1(θ (η) k ) := 1 b ∑b j=1 [ ∇F (θ(η)k , Zj)−∇f(θ (η) k ) ] . (4)
Although proposed in the 1950s by [24], SGD has been the algorithm of choice for training statistical models due to its simplicity, and superior performance in large-scale settings [25, 1, 26, 27]. However, the fluctuations of this algorithm is well-understood only when the objective function f is strongly convex and smooth, and the step size η satisfies a specific decreasing schedule so that the iterates asymptotically converge to the unique minimizer [28, 29, 30]. On the other hand, it is well-known that the SGD iterates in (2) can be viewed as a Markov chain which allows them to converge to a random vector rather than a single critical point [1]. Building on this analogy between SGD and Markov chains, the aforementioned shortcomings can be alleviated by simply relaxing the global smoothness as well as the strong convexity assumptions to the tails of the objective function f , which allows for a flexible non-convex structure around the region of interest. Similar kinds of tail relaxations have been successfully employed in the diffusion theory when the target potential is non-convex [31, 32, 33], but they are not studied in the context of non-convex optimization when the algorithm is SGD. In this work, we study the fluctuations and the bias of the averaged SGD iterates in (2), around the first-order critical points of the minimization problem (1). Our contributions can be summarized as follows.
• For a non-convex and non-smooth objective function f with tails growing at least quadratically, we establish the uniqueness of the stationary distribution of the constant step size SGD iterates in Proposition 2.1, and the asymptotic normality of Polyak-Ruppert averaging in Theorem 2.1. To the best of our knowledge, these are the first uniqueness and normality results for the SGD algorithm when the objective function is non-convex (even not strongly convex) and non-smooth. • We further show in Theorems 3.1 and 3.2 that, with additional local smoothness assumptions on the non-convex objective function f , we can establish a control over the bias in terms of the step size. We further characterize the bias when the objective is (not strongly) convex in Theorem 3.3, providing a thorough bias analysis for the constant step size SGD under various settings that are frequently encountered in statistical learning.
Our results provide algorithm-dependent guarantees for uncertainty quantification, and they could be leveraged to obtain confidence intervals (CIs) for non-convex and non-smooth learning problems. This is contrary to the majority of the existing results in statistics, which only establish normality results for the true stationary point of the non-convex objective function; see for example [12, 34]. While being useful, such results completely ignore the computational hardships associated with non-convex optimization; hence, their practical implications are limited. On the other hand, in the optimization and learning theory literature, a majority of the existing results establish the rate of convergence of an algorithm to a critical point, and do not quantify the fluctuations associated with that algorithm. Our work bridges these separate lines of thought by providing asymptotic normality results directly for the SGD algorithm used for minimizing non-convex and non-smooth functions.
More Related Works. Establishing asymptotic normality results for the SGD algorithm began with the works of [35, 36, 37, 38, 39], with [28] providing a definitive result for strongly convex objectives. In particular, [28] and [38] established that the averaged SGD iterates with an appropriately chosen decreasing step size is asymptotically normal with the variance achieving the Cramer-Rao lower
bound for parameter estimation. Recent works, for example [40, 41, 29, 42, 43, 44], leverage the asymptotic normality analysis of [28], which leads to computing CIs for SGD. The benefits of constant step size SGD for faster convergence under overparametrization has also been demonstrated in the works of [45, 46, 47, 48]. The use of Markov chain theory to study constant step size stochastic approximation algorithms has been considered in several works [49, 50, 51, 25, 52, 53]. Recently, [1, 54] investigated the asymptotic variance of the constant step size SGD. We emphasize here that most of the above works assume strongly convex and smooth objective functions.
The non-linear autoregressive (NLAR) process [55, 56, 57] is a specification of our general framework (2) with the noise sequence {ξk}k≥1 being a collection of i.i.d mean-zero random vectors with continuous density supported on Rd. However, the methodology for establishing the geometric ergodicity of NLAR [55, 56, 57, 58] is by no means straightforward to carry over to the optimization setting, and does not generalize immediately to the state-dependent noise setup considered in our paper (see Assumption 2.3). In contrast, we establish the geometric ergodicity under easily verifiable assumptions on the objective function using tools from Markov chain theory. Moreover, additional steps are needed to go from geometric ergodicity to CLT results (especially if the chain starts with an arbitrary initial distribution), while we directly obtain a CLT result by explicitly leveraging the Markov chain structure.
Furthermore, there exists a vast literature on analyzing Langevin diffusion-based sampling algorithms which relies on the much simpler i.i.d. Gaussian noise sequence. We refer the interested reader to [59, 60, 32, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70] and the references therein, for details.
Finally, recent work [71] establish the trade-off between the stability, statistical accuracy, and computational efficiency for the non-convex optimization algorithms.
Notation. For a, b ∈ R, denote by a ∨ b and a ∧ b the maximum and the minimum of a and b, respectively. We use ‖ · ‖ to denote the Euclidean norm in Rd. We denote the largest eigenvalue of the matrix A as λmax(A), and the smallest one as λmin(A). Let (Ω,F ,P) represent a probability space, and denote by B(Rd), the Borel σ-field of Rd. Let Pk(Rd) := {ν : ∫ Rd ‖θ‖
kν(dθ) < ∞} denote the set of probability measures with finite k-th moments. For a probability distribution π and a function g on X , we define π(g) := ∫ X g(x)dπ(x), and L2(π) := {g : X → R : π(g 2) <∞}.

heading: 2 Central Limit Theorem for The Constant Step Size SGD
text: In this section, we establish an asymptotic central limit theorem (CLT) for the Polyak-Ruppert averaging of the constant step size SGD iterates given in (2) when the objective function is potentially non-convex, non-smooth, and has quadratically growing tails. More specifically, we first prove that there exists a unique stationary distribution πη ∈ P2(Rd) for the Markov chain defined by the SGD algorithm when the objective function is dissipative (see Assumption 2.2) with gradient exhibiting at most linear growth (see Assumption 2.1). Furthermore, under the same conditions, we prove that a CLT holds for the Polyak-Ruppert averaging, and it is independent of the initialization. In what follows, we list and discuss the main assumptions required to establish a CLT for the SGD iterates, and compare them to those existing in the literature. Assumption 2.1 (Linear growth). The gradient of the objective function f has at most linear growth. That is, for some L ≥ 0, we have ‖∇f(θ)‖ ≤ L ( 1 + ‖θ‖ ) for all θ ∈ Rd .
Majority of the results on SGD focus on smooth functions with gradients satisfying ‖∇f(θ) − ∇f(θ′)‖ ≤ ‖θ − θ′‖ for all θ, θ′ ∈ Rd; see e.g. [28, 1]. The above condition allows for non-smooth objectives, and is a significant relaxation of the standard Lipschitz gradient condition. Assumption 2.2 (Dissipativity). The objective function f is (α, β)-dissipative. That is, there exists positive constants α, β such that 〈θ, ∇f(θ)〉 ≥ α ‖θ‖2 − β for all θ ∈ Rd .
The dissipativity assumption has its origins in the analysis of dynamical systems, and is used widely in the analysis of optimization and learning algorithms [72, 31, 33, 73]. It also recently have been used in Bayesian analysis [74]. This assumption could be viewed as a relaxation of strong convexity since it restricts the quadratic growth assumption to the tails of the function f , enforcing no local growth around the first-order critical points. A canonical example for this condition is the sum of a quadratic and any non-convex function with bounded gradient. For example, consider the function x → x2 + 10 sin(x) which is clearly non-convex and (1, 25)-dissipative. It is worth mentioning
that many non-convex problems that arise in statistical learning such as phase retrieval [53] satisfy Assumption 2.2. We provide examples in Section 4. Assumption 2.3 (Noise sequence). Gradient noise sequence {ξk}k≥1 is a collection of i.i.d. random fields satisfying E[ξ1(θ)] = 0 and E1/2[‖ξ1(θ)‖2] ≤ Lξ(1 + ‖θ‖) , for any θ ∈ Rd and a positive constant Lξ. Moreover, for each θ ∈ Rd the distribution of the random variable ξ1(θ) can be decomposed as µ1,θ +µ2,θ where µ1,θ has a density, say pθ, with respect to Lebesgue measure which satisfies infθ∈C pθ(t) > 0 for any bounded set C and any t ∈ Rd.
Assumption 2.3 as formulated above is stronger than what is needed in the proofs. It can easily be seen that the lower bound on the density pθ is only required to hold for a specific set whose form depends on η and various constants from Assumptions 2.1–2.3. The form of this set is complicated, and an exact expression is given in the Appendix – see (12). We also emphasize that Assumption 2.3 does not specify any explicit parametric form for the distribution of the noise sequence contrary to recent works in non-convex settings where dissipitavity condition has been heavily utilized [31, 73, 33].
We now establish the existence and uniqueness of the stationary distribution of the SGD iterates. Proposition 2.1 (Ergodicity of SGD). Let the Assumptions 2.1-2.3 hold, and the step size satisfy
0 < η < α− √
(α2−(3L2+Lξ))∨0 3L2+Lξ .
(a) SGD (2) admits a unique stationary distribution πη ∈ P2(Rd), depending on the step size η.
(b) For a test function φ : Rd → R satisfying |φ(θ)| ≤ Lφ(1 + ‖θ‖), ∀θ ∈ Rd and some Lφ > 0, and for any initialization θ(η)0 = θ0 ∈ Rd of the SGD algorithm, there exists ρ ∈ (0, 1) and κ (both depending on η) such that we have∣∣E[φ(θ(η)k )]− πη(φ)∣∣ ≤ κ ρk(1 + ‖θ0‖2) , where πη(φ) := ∫ φ(x)dπη(x).
The uniqueness of the stationary distribution of the constant step size SGD has been established in [1] for strongly convex and smooth objectives. In Proposition 2.1, we relax both of these assumptions allowing for non-convex and non-smooth objectives. Our proof relies on V -uniform ergodicity [75], which is fundamentally different from the ergodicity analysis in [1]. Under the dissipativity condition (quadratic growth of f ), geometric ergodicity in Proposition 2.1 is not surprising; yet, it is worth highlighting that the function f as well as the noise sequence require significantly less structure than what was assumed in the literature. The above step size assumption is almost standard and it is required to obtain a uniform bound on the moments of SGD iterates. We highlight that similar to the gradient descent algorithm, the step size depends on a quantity that serves as a surrogate condition number in our setting, namely, L/α. Note that ρ depends on η and will typically be converging to one if η → 0. Thus convergence in Proposition 2.1(b) can be expected to be slower when η becomes smaller. However, smaller η leads to a better control of the asymptotic bias (under additional regularity assumptions), see Theorems 3.1-3.3. Both of those statements (slower convergence for smaller η but smaller bias eventually) are confirmed in our numerical experiments, see Figure 1(d,h).
Next, we state our first principal contribution, a central limit theorem for the averaged SGD iterates starting from any initial distribution, for a non-convex objective. For a test function φ : Rd → R, we denote the centered partial sums of φ evaluated at the SGD iterates with Sn(φ), i.e.,
Sn(φ) := ∑n−1 k=0 [ φ ( θ (η) k ) − πη(φ) ] , where πη(φ) := ∫ φ(x)dπη(x) .
Theorem 2.1 (CLT). Let the Assumptions 2.1-2.3 hold. For a step size η and a test function φ satisfying the conditions in Proposition 2.1, we define σ2πη (φ) := limn→∞ 1 nEπη [ S2n(φ) ] . Then,
n−1/2Sn(φ) d−→ N ( 0, σ2πη (φ) ) .
The above result characterizes the fluctuations of a test function φ averaged across SGD iterates, even when the objective function is both non-convex and non-smooth. The asymptotic variance in the above CLT can be equivalently stated in another compact form. If we define the centered test function as h(θ) = φ(θ)− πη(φ), the asymptotic variance can be written as
σ2πη (φ) = 2πη(hĥ)− πη(h 2), where ĥ = ∑∞ k=0 E [ h ( θ (η) k )] .
Indeed, this is the variance we compute at the end of our proof in Section A. However, the expression in Theorem 2.1 is obtained by simply applying [58, Thm 21.2.6]. For the case of strongly convex functions with decreasing step size schedule, it is well-known from the works of [28, 38] that the limiting variance of the averaged SGD iterates achieves the Cramer-Rao lower bound for parameter estimation; see also [76, 30] for non-asymptotic rates in various metrics. The question of providing lower bounds for the limiting variance of the critical points in the non-convex setting is extremely subtle, and is often handled on a case-by-case basis. We refer the interested reader to [77, 78, 12].
There are several important implications of the above CLT for constructing CIs in practice. First note that, following the standard construction in inference, one can write the distribution of the sample mean approximately as n−1Sn(φ) ≈ N ( 0, n−1σ2πη (φ) ) . Here, one needs to estimate the population quantity, the asymptotic variance σ2πη (φ), for the purpose of obtaining CIs. In Section 5, we discuss three strategies for estimating this quantity, which could be eventually used for inference in practice. A theoretical analysis of the proposed approaches in Section 5 is beyond the scope of this work.

heading: 3 Bias of the Constant Step Size SGD
text: Proposition 2.1(b) shows that the expectation of a test function evaluated at the k-th iterate converges exponentially fast to the expected value of the stationary distribution πη. Therefore, a complete characterization of the properties of the SGD requires a control over the asymptotic bias πη(φ)−φ(θ∗) for a critical point θ∗. It turns out that this bias behavior is intimately related to the local properties of the objective around its critical points. Therefore, under the mild assumptions that yield the CLT, one cannot expect a tight control over the bias. This section contains three types of bias analyses under different local growth conditions on the objective function f , characterizing the bias behavior in various non-convex and convex settings. We further note that without local regularity conditions, it is still possible to show that the SGD iterates (2) move towards a compact ball containing all critical points exponentially fast; a formal statement of this result along with a corresponding discussion is provided in Proposition B.1, which is deferred to Appendix B. Throughout this section, we make a slightly stronger assumption on the noise sequence. Assumption 3.1 (Fourth moment of the noise). Gradient noise sequence {ξk}k≥1 satisfies Assumption 2.3, and E [ ‖ξ1(θ)‖4 ] ≤ Lξ(1 + ‖θ‖4) , for any θ ∈ Rd, where Lξ is as in Assumption 2.3.
Localized Dissipativity Condition: We now introduce the generalized dissipativity condition which, in addition to the quadratic tail growth property enforced in Assumption 2.2, imposes a local growth within some compact region, around the unique critical point θ∗. Assumption 3.2 (Localized dissipativity). The objective function f satisfies
〈∇f(θ), θ − θ∗〉 ≥ { α‖θ − θ∗‖2 − β ‖θ − θ∗‖ ≥ R g ( ‖θ − θ∗‖ ) ‖θ − θ∗‖ < R ,
where θ∗ ∈ Rd is the unique minimizer of f , R := δα + √ β α with δ ∈ (0,∞), g : [0,∞)→ [0,∞) is a convex function with g(0) = 0 whose inverse exists.
If g(x) = x2, the objective function is locally strongly convex. However, the above assumption covers a wide range of objectives with different local growth rates depending on the function g. Next, we show that the above assumption along with the assumptions leading to the CLT is sufficient to establish an algorithmic control over the bias with a sufficiently small step size. Theorem 3.1. Let the Assumptions 2.1, 3.1, and 3.2 hold. Then SGD iterates with step size satisfying 0 < η < cL,α for cL,α in (16) admit the stationary distribution θ(η) ∼ πη which satisfies
E [ ‖θ(η) − θ∗‖ ] ≤ Cδ η + g −1(Cη).
Further, for a test function φ : Rd → R that is Lφ-Lipschitz, the bias satisfies∣∣πη(φ)− φ(θ∗)∣∣ ≤ Lφ(Cη/δ + g−1(Cη)), where
C :=3 ( 3L2 + 3L
1/2 ξ (1 + (β/α)
2) )( 1 + ∫ ‖θ‖2πη(dθ) + ‖θ∗‖2 ) . (5)
If the local growth is linear, i.e. g(x) = x, we obtain the bias |πη(φ) − φ(θ∗)| ≤ O(η). If local growth is quadratic, i.e. g(x) = x2, the growth is locally slower than the linear case; thus, we get the bias control |πη(φ) − φ(θ∗)| ≤ O(η1/2), which is worse in step size dependency; it reduces to the bound derived in [1, Lemma 10]. We highlight that [79] proves the following lower bound: lim infk→∞ E[‖θ(η)k − θ∗‖2]1/2 ≥ cη1/2 for some c > 0 under the assumption of Lipschitz gradients. This is in line with our findings since Lipschits gradients imply g(x) ≤ x2 for small x. Generalized Łojasiewicz Condition: In this section we work with a generalization of the commonly used Łojasiewicz condition in optimization. Assumption 3.3 (Generalized Łojasiewicz condition). The objective function f has a critical point θ∗ and it satisfies
‖∇f(θ)‖2 ≥ { γ { f(θ)− f(θ∗) } ‖θ − θ∗‖ ≥ R
g ( f(θ)− f(θ∗) ) ‖θ − θ∗‖ < R,
where γ,R > 0, and g : [0,∞)→ [0,∞) is a convex function with g(0) = 0 whose inverse exists.
In the case g(x) = xκ with κ ∈ [1, 2), for example, the above condition is termed as the Łojasiewicz inequality [80], and for κ = 1, it reduces to the well-known Polyak-Łojasiewicz (PL) inequality [81]. Note that this inequality implies that every critical point is a global minimizer; yet, it does not imply the existence of a unique critical point.
The next result establishes an algorithmically controllable bias bound in terms of the step size. Theorem 3.2. Let the Assumptions 2.1,2.2, 3.1, and 3.3 hold, and the Hessian satisfies ‖∇2f(θ)‖ ≤ L̃(1 + ‖θ‖), ∀θ ∈ Rd and some L̃ > 0. Then, the SGD iterates with a step size satisfying 0 < η < 2
L̃ ∧ cL,α ∧ c†L,α ∧ 1 for cL,α, c † L,α in (16) have the stationary distribution πη ,
πη(f)− f(θ∗) ≤ g−1 ( 2Mη
2−L̃η
) + 2Mη
2−L̃η , where M :=12L̃ ( L+ L
1/2 ξ + L 1/4 ξ
)2( 1 +m+m3/4 + ∫ ‖θ‖2πη(dθ) ) with
m := 87α
[( β + 6L2 + 3L
1/2 ξ + 16 ) ∫ ‖θ‖2πη(dθ) + 16L4 + 2Lξ + 128L6 + 8L3/2ξ ] .
Additionally, if the test function is given as φ = φ̃ ◦ f for a Lφ̃-Lipschitz function φ̃, it holds that∣∣πη(φ)− φ(θ∗)∣∣ ≤ Lφ̃{g−1( 2Mη2−L̃η)+ 2Mη2−L̃η} . For smooth objectives with Lipschitz gradient, [81] provides a linear rate under the PL-inequality (see also [82, Lemma 2]), which yields the asymptotic bias |πη(φ) − φ(θ∗)| ≤ O(η). The above result recovers their findings as a special case, and provides a considerable generalization.
Convexity: To make the analysis of constant step size SGD complete, we digress from the main theme of this paper and consider this algorithm in the (non-strongly) convex regime, for which there is no bias characterization known to authors. We show that, under the convexity assumption, one can achieve the same bias control as in the case of PL-inequality. Theorem 3.3. Let the Assumptions 2.1,2.2, and 3.1 hold for a convex function f . Then, the SGD iterates with a step size 0 < η < cL,α for cL,α in (16) admit the stationary distribution πη , and for a minimizer θ∗ it satisfies
πη(f)− f(θ∗)≤ Cη ,
for C in (5). Further, if the test function is given as φ = φ̃ ◦ f for a Lφ̃-Lipschitz function φ̃ , then,∣∣πη(φ)− φ(θ∗)∣∣ ≤ Lφ̃Cη . Convexity implies that any critical point θ∗ is a global minimizer, which is similar to the PL-inequality; yet, it does not imply a unique minimizer unlike strong convexity. The resulting step size dependency of the bias is the same as in the case of PL-inequality, which is because both of these conditions assert a similar gradient-based domination criterion on the sub-optimality. That is, we have in the convex case 〈∇f(θ), θ−θ∗〉 ≥ f(θ)−f(θ∗), and in the case of PL-inequality γ−1‖f(θ)‖2 ≥ f(θ)−f(θ∗).

heading: 4 Examples and Numerical Studies
text: We now demonstrate the asymptotic normality and bias in non-convex optimization with two examples arising in robust statistics for which our assumptions can be verified. We consider the online SGD setting with the update rule: θ(η)k+1 = θ (η) k − η bk ∑bk j=1∇F (θ (η) k , Zj), for k ≥ 0, with independent samples Zj ∼ P (Z) used to estimate the true gradient in each iteration k ; and also the semistochastic setting, where the noise sequence {ξk(θ)}k≥1 is independent of θ and is simply a sequence of i.i.d. random vectors – such a setting helps verifying our assumptions more explicitly.

heading: 4.1 Regularized MLE for heavy-tailed linear regression
text: While the least-squares loss function is common in the context of linear regression, it is welldocumented that it suffers from robustness issues when the error distribution of the model is heavytailed [83]. Indeed in fields like finance, oftentimes the Student’s t-distribution is used to model the heavy-tailed error [84]. In this case, defining the random vector Z := (X,Y ), the stochastic optimization problem in (3) is given by the expectation of the function F (Z, θ) := log ( 1 + (Y −
〈X, θ〉)2 )
+ λ2 ‖θ‖ 2, which is non-convex (as a function of θ) for small penalty levels λ > 0.
Correspondingly, given n independent and identically distributed samples (xi, yi), the finite-sum version of the optimization problem corresponds to minimizing the following objective function
f(θ) := 12m ∑m i=1 log ( 1 + (yi − 〈xi, θ〉)2 ) + λ2 ‖θ‖ 2. (6)
We consider the finite-sum setup and we verify our assumptions and empirically demonstrate our results on CLT as well as the bias in a clean manner.

heading: 4.1.1 Semi-stochastic Gradient Descent
text: In the experiments, X := (x1, . . . ,xm)> ∈ Rm×d represents a fixed design matrix generated from Xij ∼ Bernoulli(±1)/ √ d, and y := (y1, . . . , ym)> ∈ Rm represents the response vector generated according to the linear model yi = 〈xi, θtrue〉 + ε with (θtrue)i iid∼ Unif(0, 1), and ε is Student-t distributed (df = 10) noise. We choose m = 5000, d = 10, and the Lipschitz test function φ(θ) = ‖θ‖ unless stated otherwise. Asymptotic normality: Fig. 1-(a,b,c,d) demonstrates the normality and the bias of SGD with heavytailed gradient noise distributed as Student-t (df = 5). Each plot has two density curves where red and blue curves in Fig. 1-(a,b) respectively correspond to initializations with θ0 = (1, . . . , 1)> and θ′0 = (1.5, . . . , 1.5)
> with step size η = 0.3; green and orange curves in Fig. 1-c correspond to step sizes η = 0.2 and η′ = 0.3 with initialization θ0. All experiments are based on 4000 Monte Carlo runs. We observe in Fig. 1-a that different initializations have an early impact on the normality when the number of iterations is moderate. However, when SGD is run for a longer time, this effect is removed as in Fig. 1-b. Lastly, Fig.1-c demonstrates the effect of step size on the normality, where the means are different for different step sizes as they depend on the stationary distribution πη. Indeed, the above results are not surprising as one can verify that the assumptions of Theorem 2.1 are satisfied.
Lemma 4.1. The objective function (6) satisfies Assumptions 2.1 and 2.2. Further, Assumption 2.3 is also satisfied with the Student-t distributed (df = 10) noise.
Bias: In order to demonstrate the bias behavior without speculation, one needs the global minimum θ∗ of the non-convex problem. Therefore, we simplify the problem (6) to another non-convex problem
f(θ) := 12 log ( 1 + ‖θ‖2 ) + λ2 ‖θ‖ 2 . (7)
Notice that the general structure is the same, with no data, and θ∗ is known, i.e. θ∗ = 0. We choose the test function φ(θ) = φ̃ ◦ f(θ), where φ̃(x) = 1/(1 + e−x) is Lipschitz. Fig. 1-(d) demonstrates how the bias πη(φ)− φ(θ∗) changes over iterations, where different curves correspond to different step sizes. We notice that larger step size provides fast initial decrease; yet the resulting asymptotic bias is larger which aligns with our theory – indeed, a smaller asymptotic bias for a smaller step size η is predicted by Theorem 3.2 while slower convergence can be expected given the discussion after Proposition 2.1. The following lemma proves that our assumptions are satisfied for this objective.
Lemma 4.2. The objective function (7) is non-convex when λ is sufficiently small, and it satisfies Assumptions 2.1,2.2, and 3.3. Further, Assumption 3.1 is also satisfied for this example.

heading: 4.1.2 Online Stochastic Gradient Descent
text: For our online SGD experiments, we use bk = 2, for all k to obtain the stochastic gradient. We also experimented with mk = 1, 10, 50 and observed similar behavior. The distribution of the random vector Z = (X,Y ) ∈ Rd+1, is as follows: Each coordinate of the vector X ∈ Rd, is generated as Bernoulli(±1)/ √ d and given vector X , the response Y ∈ R is generated according to the linear model Y = 〈X, θtrue〉+ ε with each coordinate of θtrue ∈ Rd generated from Unif(0, 1), and fixed, and ε ∈ R is Student-t (df = 10) noise. We choose d = 10, and set a burn-in period of size 100. Asymptotic normality: Fig. 2-(a,b,c) demonstrates the normality of online SGD. Each plot has two density curves where red and blue curves in Fig. 2-(a,b) respectively correspond to initializations with θ0 = (1, . . . , 1)> and θ ′ 0 = (2.5, . . . , 2.5) > with step size η = 0.3; green and orange curves in Fig. 2-c correspond to step sizes η = 0.2 and η′ = 0.3 with initialization θ0. All experiments are based on 4000 Monte Carlo runs. We observe in Fig. 2-a that different initializations have an early impact on the normality when the number of iterations is moderate. However, when SGD is run for a longer time, this dependence is removed as in Fig. 2-b. Lastly, Fig.2-c demonstrates the effect of step size on the normality, where the means are different for different step sizes as they depend on the stationary distribution πη . Indeed, all these observations are as predicted by our theory.

heading: 4.2 Regularized Blake-Zisserman MLE for corrupted linear regression
text: While the above example was based on linear-regression with heavy-tailed noise, we now consider the case of heavy-tailed regression with corrupted noise. In this setup, the noise model in linear regression is assumed to be Gaussian, but a fraction of the noise vectors are assumed to be corrupted in the sense that they are drawn from a uniform distribution. Such a scenario arises in visual reconstruction problems; see for example [85, 33] for details. In this case, defining the random vector Z := (X,Y ), the stochastic optimization problem in (3) is given by the expectation of the function F (Z, θ) := log ( ν+ e−(Y−〈X, θ〉) 2 )
+ λ2 ‖θ‖ 2, for ν > 0. Similar the previous case, we also
consider the finite-sum version: Given n independent and identically distributed samples (xi, yi), it corresponds to minimizing the following objective function
f(θ) = − 12m ∑m i=1 log ( ν + e−(yi−〈xi, θ〉) 2 ) + λ2 ‖θ‖ 2, ν > 0 . (8)

heading: 4.2.1 Semi-stochastic Gradient Descent
text: Asymptotic normality: In the experiments, we use the same setup and parameters as in Section 4.1.1. Fig 1-(e,f,g) demonstrates the asymptotic normality of the SGD with heavy-tailed gradient noise Student-t(df = 6). The experimental setup is the same as the previous example with the same values for θ0, θ′0, η, η
′. We observe the early impact of initialization in Fig 1-a, the clear normality in Fig. 1-b, and the effect of step size on CLT in Fig.1-c. These observations also align with our theory since this objective also satisfies our assumptions.
Lemma 4.3. The objective function (8) satisfies Assumptions 2.1, 2.2. Further, Assumption 2.3 is also satisfied with the Student-t (df = 10) noise.
Bias: Similar to the previous example, we simplify the problem so that we can compute the bias πη(φ)− φ(θ∗). We consider the function
f(θ) := − 12 log ( ν + e−‖θ‖ 2) + λ2 ‖θ‖ 2, ν > 0 . (9)
We observe in Fig.1-h that smaller step sizes lead to smaller asymptotic bias. One can verify that this can be predicted from Theorem 3.1.
Lemma 4.4. The objective function (9) is non-convex when λ is sufficiently small, and it satisfies Assumptions 2.1 and 3.2. Further, Assumption 3.1 is also satisfied for this example.

heading: 4.2.2 Online Stochastic Gradient Descent
text: Asymptotic normality: In the experiments, we use the same setup as in Section 4.1.2. Fig. 2-(d,e,f) demonstrates the normality of online SGD. Each plot has two density curves where red and blue curves in Fig. 2-(d,e) respectively correspond to initializations with θ0 = (1, . . . , 1) and θ ′′
0 = (1.5, . . . , 1.5) with step size η = 0.3; green and orange curves in Fig. 2-c correspond to step sizes η = 0.2 and η′ = 0.3 with initialization θ0. All experiments are based on 4000 Monte Carlo runs. We observe in Fig. 2-d that different initializations have an early impact on the normality when the number of iterations are moderate. However, when SGD is run for a longer time, this effect is removed as in Fig. 2-e. Lastly, Fig.2-f demonstrates the effect of step size on the normality, where the means are different for different step sizes as they depend on the stationary distribution πη .

heading: 5 Discussions
text: By leveraging the connection between constant step size SGD and Markov chains [1], we provided theoretical results characterizing the fluctuations and bias of SGD for non-convex and non-smooth optimization which arises frequently in statistical learning.
Estimating the Asymptotic Variance: As discussed in Section 2, in order to use the established CLT to compute CIs in practice, the population expectation πη(φ) and asymptotic variance σ2πη (φ) have to be estimated. We suggest the following three ways to do so:
• Estimate them based on sample average of a single trajectory of SGD iterates, i.e., the mean πη(φ) is estimated as n−1 ∑n−1 k=0 φ ( θ (η) k ) , and the asymptotic variance σ2πη (φ) can be estimated by
adopting the online approach of [86] to the constant step size setting. The variance σ2πη (φ) can
also be estimated by the Newey-West long-run variance estimation [87, 88] or empirical variance estimation based on sub-sampling [89, Sections 4.2 and 4.6] for a single trajectory. • First run N parallel SGD trajectories and compute the average of each trajectory, to obtain N independent observations from the stationary distribution πη. Next, use the N observations to compute the sample mean and the sample variance estimators for πη(φ) and σ2πη (φ). • Leverage the online bootstrap and variance estimation approaches proposed in [43, 41, 90] for the constant step size SGD setting in order to obtain estimates for πη(φ) and σ2πη (φ).
As a confirmation of the practicability of constructing CIs, we provide preliminary experimental results for constructing CIs with minibatch SGD. We consider the data generation setup described in Sections 4.1 and 4.2 with step size 0.3, and run online SGD with batch size 2. In each run, the first 200 values are discarded. CIs are constructed for each trajectory based on sub-sampling with empirical CDF (Subsampling quantile) and variance estimation (Subsampling var) [89, Sections 4.2 and 4.6], and Newey-West long-run variance estimation (Long-run var) with data-driven bandwidth selection [87, 88]. Empirical coverage results as a function of iteration numbers (nominal level = 95%, 4000 Monte Carlo replications) for the three methods and different initializations (θ0 = (1, . . . , 1)> and θ′0 = (1.5, . . . , 1.5) >) are reported in Figure 2 (right). A non-asymptotic justification of the relative merits of the above variance estimation approaches are left as future work.

heading: Acknowledgements
text: MAE is partially funded by NSERC Grant [2019-06167], Connaught New Researcher Award, CIFAR AI Chairs program, and CIFAR AI Catalyst grant. KB is partially supported by a seed grant from the Center for Data Science and Artificial Intelligence Research, UC Davis. SV is partially supported by a discovery grant from NSERC of Canada and a Connaught New Researcher Award. The authors thank Yichen Zhang for helpful comments on an earlier version of this manuscript.

Review summaries:
This paper studies the asymptotic normality and the bias for non-convex optimization problem with constant stepsize SGD algorithm. Relaxing previous strongly convex assumption, adopting techniques from Markov chains, they show that under some mild conditions, SGD with constant stepsize will quickly converge to a stationary distribution that is uniquely determined by the stepsize $\eta$. Also, they analyze the bias under various conditions, basically showing that the bias is also (positively) dependent of the stepsize $\eta$. Last, they provide thorough experiments in Regularized MLE and Regularized Blake-Zisserman MLE settings, showing that the numerical results match their theoretical prediction.