Stationary Behavior of Constant Stepsize SGD Type Algorithms: An Asymptotic Characterization

Stochastic approximation (SA) and stochastic gradient descent (SGD) algorithms are work-horses for modern machine learning algorithms. Their constant stepsize variants are preferred in practice due to fast convergence behavior. However, constant step stochastic iterative algorithms do not converge asymptotically to the optimal solution, but instead have a stationary distribution, which in general cannot be analytically characterized. In this work, we study the asymptotic behavior of the appropriately scaled stationary distribution, in the limit when the constant stepsize goes to zero. Specifically, we consider the following three settings: (1) SGD algorithms with smooth and strongly convex objective, (2) linear SA algorithms involving a Hurwitz matrix, and (3) nonlinear SA algorithms involving a contractive operator. When the iterate is scaled by $1/\sqrt{\alpha}$, where $\alpha$ is the constant stepsize, we show that the limiting scaled stationary distribution is a solution of an integral equation. Under a uniqueness assumption (which can be removed in certain settings) on this equation, we further characterize the limiting distribution as a Gaussian distribution whose covariance matrix is the unique solution of a suitable Lyapunov equation. For SA algorithms beyond these cases, our numerical experiments suggest that unlike central limit theorem type results: (1) the scaling factor need not be $1/\sqrt{\alpha}$, and (2) the limiting distribution need not be Gaussian. Based on the numerical study, we come up with a formula to determine the right scaling factor, and make insightful connection to the Euler-Maruyama discretization scheme for approximating stochastic differential equations.

1. Introduction. Stochastic approximation (SA) algorithms are the major work-horses for solving large-scale optimization and machine learning problems. Theoretically, to achieve asymptotic convergence, we should use diminishing stepsizes (learning rates) with proper decay rate (Nemirovski et al., 2009). However, constant stepsize SA algorithms are preferred in practice due to their faster convergence (Goodfellow et al., 2016). In that case, instead of converging asymptotically to the desired solution, the iterates of constant stepsize SA have a stationary distribution. Although in many cases such weak convergence to a stationary distribution was established in the literature, it is not possible to fully characterize the limiting distribution. The reason is that, when constant stepsize is used, the distribution of the noise sequence within the SA algorithm plays an important role in the stationary distribution of the iterates. Since the distribution of the noise is in general unknown, the stationary distribution cannot be analytically characterized. In this work, building upon the works on stationary distribution of constant stepsize SA algorithms, we aim at understanding the limiting behavior of the properly scaled stationary distribution as the constant stepsize goes to zero.
More formally, with initialization X (α) 0 ∈ R d , consider the SA algorithm where F : R d → R d is a general nonlinear operator, α is the constant stepsize, and {w k } is the noise sequence. Observe that SA algorithm (1) can be viewed as an iterative algorithm for solving the equation F (x) = 0 in the presence of noise (Robbins and Monro, 1951). A typical example is when F (x) = −c∇f (x) (where c > 0 is a constant) for some objective function f (·), in this case Algorithm (1) becomes the popular stochastic gradient descent (SGD) algorithm for minimizing f (·) (Lan, 2020; Bottou, Curtis and A version of this work was submitted to a conference on 28th May, 2021 Equal contribution between Zaiwei Chen and Shancong Mou Nocedal, 2018). Another example lies in the context of reinforcement learning, where F (x) = T (x) − x, and T (·) is the Bellman operator (Sutton and Barto, 2018). In this case, Algorithm (1) captures popular reinforcement learning algorithms such as TD-learning (Sutton, 1988) and Q-learning (Watkins and Dayan, 1992). Under some mild conditions on the operator F (·), it was shown in the literature that the sequence {X (α) k } converges weakly to some random variable X (α) (Dieuleveut, Durmus and Bach, 2017;Bianchi, Hachem and Schechtman, 2020;Yu et al., 2020;Durmus et al., 2021). However, for a fixed α, it is not possible to fully characterize the distribution of X (α) because it depends on the distribution of the noise sequence {w k }, which is usually unknown. In this work, we further consider letting α go to zero, and study the distribution of a properly centered and scaled iterate. Specifically, let Y (α) k := (X (α) k − x * )/g (α), where x * is the solution of F (x) = 0, and g : R → R is a properly chosen scaling function 1 . When k goes to infinity, we expect that Y (α) k converges weakly to some random variable Y (α) . Then we let α go to zero, and our goal is to further characterize the weak limit of Y (α) . Notice that proper scaling of the iterates is essential for raveling its fine grade behavior because otherwise the limiting distribution of the un-scaled iterates will converge to a singleton as the stepsize α goes to zero, which is analogous to the almost sure convergence result for using diminishing stepsizes in SA literature.
To summarize, we want to find a suitable scaling function g(·) and characterize the following two-step weak convergence of the centered scaled iterate Y where we use the notation ⇒ for weak convergence (or convergence in distribution).

Main Contributions. Our main contributions are twofold.
Characterizing the Distribution of Y . We propose a general framework for characterizing the distribution of Y in the following 3 cases: (1) SGD with a smooth and strongly convex objective, (2) linear SA with a Hurwitz matrix, and (3) SA involving a contractive operator. In particular, we show that in all three cases above the correct scaling function is g(α) = √ α, and the distribution of Y is Gaussian with mean zero and covariance matrix being the unique solution of an appropriate Lyapunov equation. Our proof is to use the characteristic function as a test function to obtain an integral equation of the distribution of Y , and then show that the desired Gaussian distribution solves the integral equation.
Determining the Suitable Scaling Function. For more general SA algorithms, we show empirically that the scaling function need not be g(α) = √ α and the distribution of Y need not be Gaussian. Inspired by this observation, we propose a method to find the the correct scaling function for general SA algorithms. In particular, our results indicate that the scaling function g(α) should be chosen such that (1) lim α→0 1.2. An Illustrative Example. In this section, we provide an example to further illustrate the problem we are going to study. Consider SA algorithm (1). Suppose that F (x) = −x is a scalar valued function, and {w k } is a sequence of i.i.d. standard normal random variables. We make such noise assumption here only for ease of exposition, and it will be relaxed in later sections. In this case, the SA algorithm (1) becomes This algorithm has the following two simple interpretations: (1) it can be viewed as the SGD algorithm for minimizing a quadratic objective function f (x) = x 2 /2, which has a unique minimizer x * = 0, (2) it 1 The scaling function is unique up to a numerical factor can also be viewed as an SA algorithm for solving the fixed-point equation T (x) = x with T (x) being identically equal to zero, therefore x * = 0 is the unique fixed-point.
and V[ · ] denote the mean and variance, respectively, we have converges weakly to a random variable Y (α) , whose distribution is N (0, 1 2−α ). In this case, we are able to analytically characterize the distribution of Y (α) for a fixed α because of the simplicity of the algorithm (3) and the noise {w k } being i.i.d. Gaussian. For the general SA algorithm (1) with limited information on the noise sequence {w k }, it is in general not possible to fully characterize the distribution of Y (α) . Now consider the second weak convergence in Eq.
(2). Since we have already shown that Y (α) ∼ N (0, 1 2−α ). As α goes to zero, we see that Y (α) converges weakly to a random variable Y , whose distribution is N (0, 1 2 ). As opposed to the first weak convergence in Eq.
(2), where the distribution of Y (α) in general cannot be fully characterized, we are able to characterize (in later sections) the distribution of Y for the more general algorithm (1), and for more general noise assumptions. Intuitively, the reason is that as the constant stepsize decreases, the effect of the entire distribution of the noise {w k } on the distribution of Y α is weakened. In the limit, only the mean and the variance of w k play roles in determining the distribution of Y . This is analogous to central limit theorem type of results.
To summarize, we have shown in the special case of Algorithm (3) that the correct scaling function is g(α) = √ α, and the distribution of the limiting random variable Y is a Gaussian distribution with mean zero and variance 1/ √ 2. In Section 2, we extend this result to the more general algorithm (1) with weaker noise assumptions.
1.3. Related Literature. Since proposed in (Robbins and Monro, 1951), SA has been popular for solving large scale optimization in modern machine learning applications. Although require using diminishing stepsizes to achieve convergence (Hu et al., 2017;Xie, Wu and Ward, 2020;Mertikopoulos et al., 2020;Shamir and Zhang, 2013;Li and Orabona, 2019;Fehrman, Gess and Jentzen, 2020;Gower, Sebbouh and Loizou, 2021), constant stepsize is used in practice (Goodfellow et al., 2016). In contrast to the success in machine learning practice, there is little discussion about the stationary distribution of constant step size SGD (Dieuleveut, Durmus and Bach, 2017;Bianchi, Hachem and Schechtman, 2020;Yu et al., 2020;Durmus et al., 2021). Among them, Dieuleveut, Durmus and Bach (2017) bridges Markov chain theory and the constant step size SGD algorithm. They provided an explicit asymptotic expansion of the moments of the averaged SGD iterates. Bianchi, Hachem and Schechtman (2020) studied the asymptotic behavior of constant step size SGD for a nonconvex nonsmooth, locally Lipchitz objective function. It was shown that in a small step size regime, the interpolated trajectory of the algorithm converges in probability towards the solutions of the differential inclusionẋ = ∂F (x) and the invariant distribution of the corresponding Markov chain converges weakly to the set of invariant distribution of the differential inclusion. Yu et al. (2020) established an asymptotic normality result for the constant step size SGD algorithm for a non-convex and non-smooth objective function satisfying a dissipativity property. Durmus et al. (2021) first established non-asymptotic performance bounds under Lyapunov conditions and then proved that for any step size, the corresponding Markov chain admits a unique stationary distribution.
The work mentioned before established the stationary distribution for almost strongly convex and smooth functions for a fixed constant stepsize. Since the SGD iterates will converge to a singleton as the constant step size goes to zero, none of the previously mentioned literates can be applied to study the limiting behavior of SGD in this regime. To understand such behavior, we propose to study the properly centered and scaled iterates. Although not directly related, it shares a similar flavor when studying the limiting behavior of the stationary distribution of the stochastic gradient Langevin dynamics (SGLD) iterates as step size goes to zero.
Another set of related literature is on the diffusion approximation of SGD (Li, Tai and Weinan, 2017;Feng, Li and Liu, 2017;Yang, Hu and Li, 2021;Sirignano and Spiliopoulos, 2020;Latz, 2021). Authors aim to approximate the trajectory of SGD by a diffusion process which solves an SDE. Notice that they also study the scaled version of the diffusion limit of SGD. However, different from our approach, their scale is in temporal domain and cannot be applied in our research.
The Markov chain perspective of studying SGD iterates when step size goes to zero (Dieuleveut, Durmus and Bach, 2017) is related to the heavy traffic analysis in queuing theory (Eryilmaz and Srikant, 2012). It has been studied in the literature using fluid and diffusion limits (Gamarnik and Zeevi, 2006;Harrison, 1988Harrison, , 1998Harrison and López, 1999;Stolyar et al., 2004;Williams, 1998) where the interchange of limit is usually problematic (Eryilmaz and Srikant, 2012). An alternative approach in stochastic networks is based on drift arguments introduced by (Eryilmaz and Srikant, 2012) and further generalized by (Maguluri and Srikant, 2016;Maguluri, Burle and Srikant, 2018;Hurtado-Lange, Varma and Maguluri, 2020;Mou and Maguluri, 2020). We adopt similar techniques in quantifying the limiting distribution of the scaled SGD iterates. Notice that in stochastic networks, people mainly focus on finite state space Markov chains. However, when it comes to SGD iterates, the state space is continuous and thus more challenging.
The rest of this paper is organized as follows. In Section 2, we characterize the distribution of Y (cf. Eq (2)) in the following cases: (1) F (·) is the negative gradient for some smooth and strongly convex function is a contraction operator. In all three cases above, the scaling function is chosen to be g(α) = √ α. Then in Section 3, we first empirically show that for more general SA algorithms beyond these cases, the scaling function need not be g(α) = √ α and the distribution of Y need not be Gaussian. Then we present a method to determine the scaling function for more general SA algorithms and make connection to the Euler-Maruyama discretization scheme for approximating SDE. Finally, we conclude this paper in Section 4.
2. Characterizing the Asymptotic Stationary Distribution. Through out this section, we make the following assumption regarding the noise sequence {w k }.
ASSUMPTION 2.1. The sequence {w k } is independent and identically distributed with mean zero and a positive definite covariance matrix Σ ∈ R d×d .
Note that Assumption 2.1 is much weaker than the assumption used in Section 1.2, where the noise is assumed to obey the Gaussian distribution. Nevertheless, extending our results to the more general noise setting (e.g. martingale difference noise, Markovian noise, etc) is one of our future directions.

SGD for Minimizing a Smooth and Strongly
is an objective function. Then the SA algorithm becomes which is the well-known SGD algorithm for minimizing f (·). To proceed, we require the following definition.
DEFINITION 2.1. A convex differentiable function h : R d → R is said to be L-smooth and σ-strongly convex with respect to the Euclidean norm · 2 if and only if To characterize the asymptotic behavior of Algorithm (4), we make the following assumption.
ASSUMPTION 2.2. The function f : R d → R is twice differentiable, and is L-smooth and σ-strongly convex.
Under Assumption 2.2, the function f (x) has a unique minimizer (or F (x) = 0 has a unique solution), which we have denoted by x * . To proceed, let Y which is obtained by first subtracting both sides of Eq. (4) by x * , and then dividing by √ α. We next characterize the two-step weak convergence (cf. Eq. (2)) in the following theorem. Let H f ∈ R d×d be the Hessian matrix of the objective function f (·) evaluated at x * , which is well-defined because f (·) is twice differentiable (cf. Assumption 2.2). THEOREM 2.1. Consider the iterates {Y (α) k } generated by Algorithm (5). Suppose that Assumptions 2.1 and 2.2 are satisfied, then the following statements hold.
(1) There exists a thresholdᾱ > 0 such that for all α ∈ (0,ᾱ), the sequence of random variables {Y converges weakly to a random variable Y , which satisfies the following integral equation for all t ∈ R d . In addition, suppose that Eq. (6) has a unique solution (in terms of the distribution of Y ), then the distribution of Y is a Multivariate normal distribution with mean zero and covariance matrix Σ Y being the unique solution of the Lyapunov equation REMARK. To establish Theorem 2.1 (2), we require Eq. (6) to have a unique solution in terms of the distribution of Y . Such assumption will be relaxed to some extent in Section 2.4.
Since Σ is positive definite (cf. Assumption 2.1), and H f is also positive definite under Assumption 2.2, it is well established in the literature that the Lyapunov equation H f Σ Y + Σ Y H f = Σ has a unique solution, which is explicitly given by Consider the special case where f (x) = x 2 /2. In this case we have H f = 1, and hence Σ Y = 1 2 Σ by the Lyapunov equation. As a result, the distribution of the limiting random variable Y is Gaussian with mean zero, and covariance matrix being 1 2 Σ. This agrees with the illustrative example (which is for scalar valued iterates) presented in Section 1.2.
From Theorem 2.1, we see that the distribution of Y only depends on the Hessian of f (·) at x * . This makes intuitive sense because we are studying the asymptotic behavior of SA algorithm (4), and only the properties of f (·) around x * should play a role in determining the stationary distribution.

Stochastic Approximation for Solving Linear Systems of Equations
which aims at iteratively solving the linear equation Ax + b = 0. Note that since the matrix A is not necessarily symmetric, F (x) need not be a gradient of any objective function. Such linear SA algorithm arises in many realistic applications. One typical example is the TD-learning algorithm for solving the policy evaluation problem in reinforcement learning, where the goal is to solve a linear Bellman equation. See Bertsekas and Tsitsiklis (1996); Srikant and Ying (2019) for more details.
To study the asymptotic behavior of Algorithm (7), we make the following assumption regarding the matrix A.
ASSUMPTION 2.3. The matrix A is Hurwitz, i.e., all the eigenvalues of the matrix A have strictly negative real part.
REMARK. Since A being Hurwitz implies A being non-singular, Assumption 2.3 implies that the target equation Ax + b = 0 has a unique solution, which we denote by x * . Assumption 2.3 is standard in studying linear SA algorithms. In particular, it was shown in the literature that under Assumption 2.3 and some mild conditions on the noise {w k }, Algorithm (7) converges in the mean square sense to a neighborhood around x * (Bertsekas and Tsitsiklis, 1996).
To study the asymptotic distribution, for a fixed stepsize α, we define the centered scaled iterate Y The full characterization of the two-step weak convergence (cf. Eq. (2) (1) There exists a thresholdᾱ > 0 such that for all α ∈ (0,ᾱ ), the sequence of random variables {Y (2) For any positive sequence {α i } satisfying α i ∈ (0,ᾱ ) for all i and lim i→∞ α i = 0, the sequence of random variables {Y (αi) } converges weakly to a random variable Y , which satisfies the following equation In addition, suppose that Eq. (9) has a unique solution, then the distribution of Y is a Multivariate normal distribution with mean zero and covariance matrix being the unique solution of the Lyapunov equation AΣ Y + Σ Y A + Σ = 0.
Since the matrix A is Hurwitz, and Σ is positive definite, the existence and uniqueness of a positive definition solution to the Lyapunov equation AΣ Y + Σ Y A + Σ = 0 are guaranteed (Haddad and Chellaboina, 2011). Lyapunov equations were used extensively in studying the stability of linear ordinary differential equations (ODE). Interestingly, it also shows up in determining the limit distribution of centered scaled iterates of discrete linear SA algorithms.

Stochastic Approximation under Contraction
which can be interpreted as an SA algorithm for finding the fixed-point of the operator T (·). These type of algorithms arise in the context of reinforcement learning, where T (·) is the Bellman operator. To proceed, we need the following definition.
Next, we state our assumption regarding the operator T (·).
Assumption 2.4 essentially states that the operator T (·) is a contraction mapping with respect to the weighted 2 -norm · µ . By Banach fixed-point theorem (Banach, 1922), the operator T (·) has a unique fixed-point, which we denote by x * .
We next write down the update equation of the centered scaled iterate Y In the next theorem, we characterize the distribution of the limiting random vector Y (cf. Eq. (2)). Let J ∈ R d×d be the Jacobian matrix of T (·) evaluated at x * . THEOREM 2.3. Consider the iterates {Y (α) k } generated by Algorithm (11). Suppose that Assumptions 2.1 and 2.4 are satisfied, then the following statements hold.
(2) For any positive sequence {α i } satisfying α i ∈ (0,ᾱ ) for all i and lim i→∞ α i = 0, the sequence of random variables {Y (αi) } converges weakly to a random variable Y , which satisfies the following equation In addition, suppose that Eq. (12) (12)), to conclude that the distribution of Y is Gaussian, we need to assume that the equation has a unique solution. In this section, we show that such uniqueness assumption can be relaxed to some extend.
where the interchange of integral and differentiation is justified (Flanders, 1973). Now Eq. (13) is an ODE, which has solutions of the form where C is a constant. Since φ Y (t) as a characteristic function satisfies φ Y (0) = 1, we have C = 1 and hence φ Y (t) = exp(− a 4b t 2 ), which is characteristic function for a Gaussian random variable with mean zero and covariance a/(2b). Therefore, the uniqueness assumption can be removed in the uni-dimensional setting.

Multi-Dimensional Setting.
Moving to the multi-dimensional setting, consider Eq. (6) of Theorem 2.1 as a representative example. To show the same result of Theorem 2.1 (2) without imposing the uniqueness assumption, we consider the case where (1) the Hessian matrix H f of the objective function f (·) evaluated at x * is the identity matrix, and (2) the covariance matrix of the noise w k is also an identity matrix. Extending the result to the more general setting where H f and Σ can be any positive definite matrices is a future research direction.
Similarly let φ Y (t) = E[e it Y ] be the characteristic function of the random vector Y . Then in this case Eq. (6) becomes t tφ Y (t) + 2t ∇φ Y (t) = 0, which is equivalent to where ψ Y (t) := log(φ Y (t)). To solve the partial differential equation (PDE) (14), we will first convert the PDE from Cartesian coordinates to spherical coordinates, which then is directly solvable.
2.5. Proof of Theorem 2.1. In this section, we present our proof for Theorem 2.1. The proofs for Theorems 2.2 and 2.3 are similar and hence are omitted.
Before going into details, we first highlight the main ideas for the proof. For Theorem 2.1 (1), to show that {Y (αi) } converges weakly to a some random variable Y , we show that any sub-sequence of {Y (αi) } further contains a weakly convergent sub-sequence, with a common weak limit. We do it by first showing that the sequence of random variables {Y (αi) } is tight. In particular, we show that there exists a constant B independent of α such that E[ Y (α) 2 2 ] ≤ B for all small enough α. This result in conjunction with the Markov inequality implies the tightness of {Y (αi) }. As a result of tightness, {Y (αi) } contains a weakly convergent sub-sequence. Now consider Theorem 2.1 (2). For any positive sequence {α i } such that lim i→∞ α i = 0, since the family of random variables {Y (αi) } is tight, there is a weakly convergent subsequence {Y αi k }. We further show that the weak limit Y of the subsequence {Y (αi k ) } solves Eq. (6). In this case, under the assumption that Eq. (6) has a unique solution, the random variable Y is a Gaussian random variable with mean zero, and covariance matrix Σ Y being the unique solution of the Lyapunov equation H f Σ Y + Σ Y H f = Σ. Since for every sequence {Y (αi) }, there is a weakly convergent subsequence {Y (αi k ) } with a common weak limit, the sequence of random variables {Y (αi) } also converge weakly to the same random variable Y .
2.5.1. Proof of Theorem 2.1 (1). To prove the result, we will apply Proposition 2.1 in Yu et al. (2020). For completeness, we first state this proposition using our notation in the following.
(a) Since the objective function f (·) is assumed to be L-smooth, we have for any (b) Since the objective function is assumed to be σ-strongly convex, we have for any x ∈ R d : (c) This is immediately implied by Assumption 2.1, with L = Trace(Σ) 1/2 . Now apply Proposition 2.1, when the stepsize α satisfies α < σ 2(3L 2 +Trace(Σ)) , the SGD iterates {X (α) k } converge weakly to some random variable X (α) , which is distributed according to the unique stationary distribution π α . In addition, we have k } converges weakly to some random variable Y (α) and E[ Y (α) 2 2 ] < ∞.
2.5.2. Proof of Theorem 2.1 (2). Following the road map described in the beginning of this section, we present and prove a sequence of lemmas in the following. Together they imply the desired result.
LEMMA 2.1. Let α 0 = σ/L 2 . the family of random variables {Y (α) } 0<α≤α0 is tight. PROOF OF LEMMA 2.1. We first show that that there exists an absolute constant C > 0 such that E[ Y (α) 2 ] ≤ C for any α ∈ (0, α 0 ]. Using the update equation (1) The existence and uniqueness of a stationary distribution Y (α) is proved in Part (1) of this theorem. We next show that the family of random variables {Y (α) } 0≤α≤α0 is tight. Using the equation and we have By smoothness, we have By strong convexity, we have . When α ∈ (0, α 0 ], we have from the previous inequality that Hence, for any α > 0, let M = Trace(Σ)/σα, then we have for any α ∈ (0, α 0 ]. It follows that the family of random variables {Y (α) } 0<α≤α0 is tight.
LEMMA 2.2. Let {α i } be a positive sequence of real numbers such that lim i→∞ α i = 0. Suppose that {Y αi } converges weakly to some random variable Y . Then Y satisfies the following equation PROOF OF LEMMA 2.2. For any i ≥ 0, we have which implies for any t ∈ R d : Using the Taylor's theorem and we have Using Theorem 3.3.20 from Durrett (2019) and we have Using the previous two inequalities in Eq. (18) and we have Simplify the above equality and we obtain We next let i go to infinity on both sides of the previous inequality and evaluate the limit of the terms Since {Y (αi) } converges weakly to some random variable Y , we have by continuity theorem (Theorem 3.3.17 in Durrett (2019)) that For the term T 6 , we have by bounded convergence theorem that lim αi→0 T 6 = 0. To evaluate the terms T 2 to T 5 , the following definition and result from Van der Vaart (2000) (2000)). Let f : R d → R be measurable and continuous at every point in a set C. Let X n ⇒ X, where X takes its values in C.
if and only if the sequence of random variables f (X n ) is asymptotically uniformly integrable.

Now consider the term T
which goes to zero as M → ∞, we have by Theorem 2.4 that Using the same line of analysis, we have lim i→∞ T 3 = lim i→∞ T 4 = lim i→∞ T 5 = 0. It follows that Rearranging terms and we obtain the resulting equation.
LEMMA 2.3. Suppose Eq. (17) admits a unique solution. Then the random variable Y given in Lemma 2.2 follows a Gaussian distribution with mean zero, and covariance matrix Σ Y being the unique solution of the Lyapunov equation PROOF OF LEMMA 2.3. Suppose that Eq. (6) has a unique solution, we only need to verify that the multinormal distribution with mean zero and covariance matrix being the solution to the Lyapunov equation (6).
3. Identifying the Suitable Scaling Function for More General SA Algorithms. In the previous section, we have shown that for several particular SA algorithms (e.g. SGD, linear SA, and contractive SA), the scaling function is g(α) = √ α and distribution of the limiting random variable Y is Gaussian. In this section, we consider more general SA algorithms. We first show impirically in the following section that in general the scaling function need not be g(α) = √ α, and the distribution of Y need not be Gaussian. (1) is the SGD algorithm for minimizing the scalar objective f (x) = x 4 /4. That is:

Numerical Experiments. Suppose that Algorithm
Note that f (·) in this case is neither smooth nor strongly convex. It is clear that the unique minimizer of f (·) is zero. Let the centered scaled iterate Y k /g(α). We next use numerical simulation to show that the correct scaling function in this case is g(α) = α 1/4 instead of g(α) = √ α. In Figures 1 and 2, we plot the empirical density function of Y (α) for different α. For the right scaling function, we expect the density function to converge as α decreases, while for the wrong scaling function, we expect the density function to change drastically for order-wise different α. As we see, it is clear that g(α) = √ α is not suitable in this case, and g(α) = α 1/4 seems to be the right scaling. To further verify this result, we plot the logarithmic empirical density function as a function of y 4 in Figure 3. We observe linear growth in Figure 3. This indicates that the density function p Y (y) is proportional to e βy 4 , where β is some numerical constant. Therefore, numerical experiments suggest that the distribution of Y is not Gaussian but super Gaussian in this problem.

3.2.
A Method to Determine the Suitable Scaling Function. Inspired by the numerical simulations provided in the previous section, we here provide a method to determine the correct scaling function for general SA algorithms.
To gain intuition, we consider the centered scaled iterates Y Notably, the factor in terms of the stepsize α in front of the term (Y (α) k ) 3 is α 3/2 , which is equal to the square of the factor α 3/4 in front of the noise term w k . Now for general SA algorithm (1), by rewriting the update equation (1) in terms of the centered scaled iterate Y In view of the previous equation and the empirical observations in the previous section, we see that we need to choose a scaling function g(α) such that the following condition is satisfied.
We next verify the choice of scaling functions in Section 2 using our proposed Condition 3.1. For SGD with a smooth and strong convex objective, since In view of the previous inequality and Condition 3.1, it is clear that the only possible choice of g(α) is g(α) = √ α. For linear SA algorithm studied in Section 2.2, one can also easily show using Condition 3.1 that g(α) = √ α. As for contractive SA studied in Section 2.3, using the contraction property and we have Since all norms are "equivalent" in finite dimensional spaces, the previous inequality implies that we must choose g(α) = √ α. Now to further verify the correctness of the scaling function suggested by Condition 3.1, consider the SGD algorithm with the following two choices of objective functions: (1) f (x) = e x 2 , and (2) f (x) = x 4 4 + sin 2 (x) 2 . Note that in these two cases the function f (·) is not smooth and strongly convex.
Case 1. In the first case where f (x) = e x 2 , since when choosing g(α) = √ α, we haveF (y) = lim α→0 g(α) 2 α 2ye (yg(α)) 2 = 2y. One interesting implication of this example is the following. In the three SA algorithms considered in Section 2, it seems that it is the function F (·) that appears in the algorithm determines the distribution of Y . However, the above example suggests that it is the functionF (·) of Condition 3.1 instead of F (·) that directly impacts the distribution Y . In SGD with a smooth and strongly convex objective, linear SA, and contractive SA, F (·) andF (·) happen to coincides, but this is in general not the case. In fact, since we have x 2i−1 by Taylor series,F (·) in this example is exactly the dominant terms appears in the Taylor series. In addition, this suggests that the distribution of the limiting random variable Y has a density function proportional to e β x 2 , where β is a numerical constant.
We next verify this choice of g(α) and the distribution of Y (α) for small enough α using numerical simulation in the following. We see from Figure 4 that with the scaling function g(α) = √ α, the empirical density function of the random variable Y (α) seems to converge. and Figure 5 further justifies this result by showing that the density function p Y (y) of the distribution of Y in this case is proportional to e β x 2 , where β is a numerical constant.
Since lim x→0 sin(x) x = 1, the only possible choice of the scaling function g(α) to satisfy Condition 3.1 (2) is g(α) = √ α. In this case, we haveF (y) = lim α→0 1 √ α y 3 α 3/2 + sin(y √ α)cos(y √ α) = y by L'Hôpital's rule. This is another example where F (·) =F (·). In fact, since x 4 is dominated by sin 2 (x) as x approaches x * (which is 0), the scaling function and the functionF (·) are determined only by the dominant term. . Similarly, we verify this choice of scaling function via numerical experiments. In Figures 6 and 7, we plot the empirical density function of the random variable Y (α) for different stepsize α, and see if the density function converges as α goes to zero. The results suggest that g(α) = α 1/2 seems to be the correct scaling. To further verify this result, we plot the logarithmic function of the empirical density of Y (α) as a function of y 2 and observe straight lines. Therefore, the distribution of Y (α) is proportional to e β x 2 , where β is a numerical constants.
3.3. Connection to Euler-Maruyama Discretization Scheme for Approximating SDE. The choice of the scaling function suggested by Condition 3.1 has an insightful connection to the Euler-Maruyama discretization scheme for approximate the solution of an SDE, as elaborated below. Let (B t ) t≥0 be a Brownian motion. Consider the following SDE: with initial condition X 0 . The Euler-Maruyama discretization {X k } to the solution (X t ) of SDE (21) is defined as follows. Let ∆t be the discretization accuracy. SetX 0 = X 0 , and recursively defineX k according toX k+1 =X k + ∆tF (X k ) + (B (k+1)∆t − B k∆t ).
Since (B t ) t≥0 is a Brownian motion, we have (B (k+1)∆t − B k∆t ) ∼ N (0, ∆t). Therefore, by letting {Z k } be an i.i.d. sequence of standard normal random variables, we can rewrite the previous equation aŝ The approximation property of the Euler-Maruyama discretization (22) to its corresponding SDE (21) has been studied in the literature, see Wenlong et al. (2019). Specifically, it was shown that under some mild conditions on F (·), the Euler-Maruyama scheme is known to have the first-order accuracy of the SDE (21). As a consequence, intuitively, when (X t ) t≥0 has a stationary distribution µ, the limiting distribution µ ∆t of {X k } as a function of the discretization accuracy ∆t should converge weakly to µ as ∆t tends to zero. If we view the discretization accuracy ∆t as the stepsize in Eq. (22). In order for µ ∆t to converge to some nontrivial distribution µ as ∆t tends to zero, it is important to notice that the scaling factor of the noise Z k in terms of ∆t must be order-wise equal to the square root of the scaling factor of F (X k ). This observation coincides with Eq. (20) in the previous section, which eventually leads to our Condition 3.1.

Conclusion.
In this paper, we characterize the asymptotic stationary distribution of properly centered scaled iterate of SA algorithms. In particular, we show that for (1) SGD with smooth and strongly convex objective, (2) linear SA, and (3) contractive SA, the scaling function is g(α) = √ α and the corresponding stationary distribution is Gaussian. For SA beyond these cases, we empirical show that the stationary distribution need not be Gaussian, and provide a method for determine the suitable scaling function. Since our paper is the first study for this problem, it might open a door for research in this direction.