{'fraction_non_alphanumeric': 0.058348117484723044, 'fraction_numerical': 0.016321703134240095, 'mean_word_length': 3.964380076328408, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fine-Grained Analysis of Stability and Generalization for Stochastic Gradient Descent

Recently there are a considerable amount of work devoted to the study of the algorithmic stability and generalization for stochastic gradient descent (SGD). However, the existing stability analysis requires to impose restrictive assumptions on the boundedness of gradients, strong smoothness and convexity of loss functions. In this paper, we provide a fine-grained analysis of stability and generalization for SGD by substantially relaxing these assumptions. Firstly, we establish stability and generalization for SGD by removing the existing bounded gradient assumptions. The key idea is the introduction of a new stability measure called on-average model stability, for which we develop novel bounds controlled by the risks of SGD iterates. This yields generalization bounds depending on the behavior of the best model, and leads to the first-ever-known fast bounds in the low-noise setting using stability approach. Secondly, the smoothness assumption is relaxed by considering loss functions with Holder continuous (sub)gradients for which we show that optimal bounds are still achieved by balancing computation and stability. To our best knowledge, this gives the first-ever-known stability and generalization bounds for SGD with even non-differentiable loss functions. Finally, we study learning problems with (strongly) convex objectives but non-convex loss functions.


Introduction
Stochastic gradient descent (SGD) has become the workhorse behind many machine learning problems. As an iterative algorithm, SGD updates the model sequentially upon receiving a new datum with a cheap per-iteration cost, making it amenable for big data analysis. There is a plethora of theoretical work on its convergence analysis as an optimization algorithm [e.g. 11,21,34,40,46,55].
Concurrently, there are a considerable amount of work with focus on its generalization analysis [9,18,26,43,51]. For instance, using the tool of integral operator the work [9,25,43,50] studied the excess generalization error of SGD with the least squares loss, i.e. the difference between the true risk of SGD iterates and the best possible risk. An advantage of this approach is its ability to capture the regularity of regression functions and the capacity of hypothesis spaces. The results were further extended in Lei & Tang [22], Lin et al. [26] based on tools of empirical processes which are able to deal with general convex functions even without a smoothness assumption. The idea is to bound the complexity of SGD iterates in a controllable manner, and apply concentration inequalities in empirical processes to control the uniform deviation between population risks and empirical risks over a ball to which the SGD iterates belong.
Recently, in the seminal work [18] the authors studied the generalization bounds of SGD via algorithmic stability [5,13] for convex, strongly convex and non-convex problems. This motivates several appealing work on some weaker stability measures of SGD that still suffice for guaranteeing generalization [7,20,57]. An advantage of this stability approach is that it considers only the particular model produced by the algorithm, and can imply generalization bounds independent of the dimensionality.
However, the existing stability analysis of SGD is established under the strong assumptions on the loss function such as the boundedness of the gradient and strong smoothness. Such assumptions are very restrictive which are not satisfied in many standard contexts. For example, the bounded gradient assumption does not hold for the simple least-squares regression, where the model parameter belongs to an unbounded domain. The strong smoothness assumption does not hold for the popular support vector machine. Furthermore, the analysis in the strongly convex case requires strong convexity of each loss function which is not true for many problems such as the important problem of least squares regression.
In this paper, we provide a fine-grained analysis of stability and generalization for SGD. Our new results remove the bounded gradient assumption for differentiable loss functions and remove the strong smoothness assumptions for Lipschitz continuous loss functions, and therefore broaden the impact of the algorithmic stability approach for generalization analysis of SGD. In summary, our main contributions are listed as follows.
• Firstly, we study stability and generalization for SGD by removing the existing bounded gradient assumptions. The key is an introduction of a novel stability measure called on-average model stability, whose connection to generalization is established by using the smoothness of loss functions able to capture the low risks of output models for better generalization. An advantage of on-average model stability is that the corresponding bounds involve a weighted sum of empirical risks instead of the uniform Lipschitz constant. The weighted sum of empirical risks can be bounded via tools in analyzing optimization errors, which implies a key message that optimization is beneficial to generalization. Furthermore, our stability analysis allows us to develop generalization bounds depending on the risk of the best model. In particular, we have established fast generalization bounds O(1/n) for the setting of low noises, where n is the sample size. To our best knowledge, this is the first fast generalization bound of SGD based on stability approach in a low-noise setting.
• Secondly, we consider loss functions with their (sub)gradients satisfying the Hölder continuity which is a much weaker condition than the strong smoothness in the literature. Although stability decreases by weakening the smoothness assumption, optimal generalization bounds can be surprisingly achieved by balancing computation and stability. In particular, we show that optimal generalization bounds can be achieved for the hinge loss by running SGD with O(n 2 ) iterations. Fast learning rates are further derived in the low-noise case.
• Thirdly, we study learning problems with (strongly) convex objectives but non-convex individual loss functions. The nonconvexity of loss functions makes the corresponding gradient update no longer non-expansive, and therefore the arguments in Hardt et al. [18] do not apply. We bypass this obstacle by developing a novel quadratic inequality of the stability using only the convexity of the objective, which shows that this relaxation affects neither generalization nor computation.
The paper is structured as follows. We discuss the related work in Section 2 and formulate the problem in Section 3. The stability and generalization for learning with convex loss functions is presented in Section 4. In Sections 5 and 6, we consider problems with relaxed convexity and relaxed strong convexity, respectively. We conclude the paper in Section 7.

Related Work
In this section, we discuss related work on algorithmic stability, stability of stochastic optimization algorithms and generalization error of SGD.
Algorithmic Stability. The study of stability can be dated back to Rogers & Wagner [42]. A modern framework of quantifying generalization via stability was established in the paper [5], where a concept of uniform stability was introduced and studied for empirical risk minimization (ERM) in the strongly convex setting. This framework was then extended to study randomized learning algorithms [13], transfer learning [20] and privacy-preserving learning [12], etc. The interplay between various notions of stability, learnability and consistency was further studied [39,45]. The power of stability analysis is especially reflected by its ability in deriving optimal generalization bounds in expectation [45]. Very recently, almost optimal high-probability generalization bounds were established via the stability approach [6,14,15]. In addition to the notion of uniform stability mentioned above, various other notions of stability were recently introduced, including uniform argument stability [28] and hypothesis set stability [16].
Stability of Stochastic Optimization Algorithms. In the seminal paper [18], the co-coercivity of gradients was used to study the uniform stability of SGD in convex, strongly convex and non-convex problems. The uniform stability was relaxed to a weaker notion of on-average stability [45], for which the corresponding bounds of SGD can capture the impact of the risk at the initial point [20] and the variance of stochastic gradients [57]. For non-convex learning problems satisfying either a gradient dominance or a quadratic growth condition, pointwise-hypothesis stabilities were studied for a class of learning algorithms that converge to global optima [7], which relaxes and extends the uniform stability of ERM under strongly convex objectives [5]. A fundamental stability and convergence tradeoff of iterative optimization algorithms was recently established, where it was shown that a faster converging algorithm can not be too stable, and vice versa [8]. This together with some uniform stability bounds for several first-order algorithms established there, immediately implies new convergence lower bounds for the corresponding algorithms. Algorithmic stability was also established for stochastic gradient Langevin dynamics with non-convex objectives [24,31] and SGD implemented in a stagewise manner [54].
Generalization Analysis of SGD. A framework to study the generalization performance of largescale stochastic optimization algorithms was established in Bousquet & Bottou [4], where three factors influencing generalization behavior were identified as optimization errors, estimation errors and approximation errors. Uniform stability was used to establish generalization bounds O(1/ √ n) in expectation for SGD for convex and strongly smooth cases [18]. For convex and nonsmooth learning problems, generalization bounds O(n − 1 3 ) were established based on the uniform convergence principle [26]. An interesting observation is that an implicit regularization can be achieved without an explicit regularizer by tuning either the number of passes or the step sizes [26,43]. For the specific least squares loss, optimal excess generalization error bounds (up to a logarithmic factor) were established for SGD based on the integral operator approach [25,38]. The above mentioned generalization results are in the form of expectation. High-probability bounds were established based on either an uniform-convergence approach [22] or an algorithmic stability approach [15]. A novel combination of PAC-Bayes and algorithmic stability was used to study the generalization behavior of SGD, a promising property of which is its applications to all posterior distributions of algorithms' random hyperparameters [29].

Problem Formulation
Let S = {z 1 , . . . , z n } be a set of training examples independently drawn from a probability measure ρ defined over a sample space Z = X × Y, where X ⊆ R d is an input space and Y ⊆ R is an output space. Our aim is to learn a prediction function parameterized by w ∈ Ω ⊆ R d to approximate the relationship between an input variable x and an output variable y. We quantify the loss of a model w on an example z = (x, y) by f (w; z). The corresponding empirical and population risks are respectively given by Here we use E z [·] to denote the expectation with respect to (w.r.t.) z. In this paper, we consider stochastic learning algorithms A, and denote by A(S) the model produced by running A over the training examples S. We are interested in studying the excess generalization error F (A(S)) − F (w * ), where w * ∈ arg min w∈Ω F (w) is the one with the best prediction performance over Ω. It can be decomposed as (3.1) The first term is called the estimation error due to the approximation of the unknown probability measure ρ based on sampling. The second term is called the optimization error induced by running an optimization algorithm to minimize the empirical objective, which can be addressed by tools in optimization theory. A popular approach to control estimation errors is to consider the stability of the algorithm, for which a widely used stability measure is the uniform stability [13,18].
Definition 1 (Uniform Stability). A stochastic algorithm A is ǫ-uniformly stable if for all training datasets S, S ∈ Z n that differ by at most one example, we have The celebrated relationship between generalization and uniform stability is established in the following lemma [18,45].
Lemma 1 (Generalization via uniform stability). Let A be ǫ-uniformly stable. Then Throughout the paper, we restrict our interest to a specific algorithm called projected stochastic gradient descent. It is worth mentioning that our main results in Section 4 hold also when Ω = R d , i.e., no projections.
Definition 2 (Projected Stochastic Gradient Descent). Let Ω ⊆ R d and Π Ω denote the projection on Ω. Let w 1 = 0 ∈ R d be an initial point and {η t } t be a sequence of positive step sizes. Projected SGD updates models by where ∂f (w t , z it ) denotes a subgradient of f w.r.t. the first argument and i t is independently drawn from the uniform distribution over {1, . . . , n}.
Note if f is differentiable, then ∂f denotes the gradient of f w.r.t. the first argument. We say a function g : for all w,w ∈ R d , where ·, · denotes the inner product and w 2 denotes the ℓ 2 norm of w = (3.4) holds with σ = 0, then we say g is convex. We denote B ≍ B if there are absolute constants c 1 and c 2 such that c 1 B ≤ B ≤ c 2 B.

Stability with Convexity
An essential assumption to establish the uniform stability of SGD is the uniform Lipschitz continuity (boundedness of gradients) of loss functions as follows [5,7,18,20,57]. Assumption 1. We assume ∂f (w; z) 2 ≤ G for all w ∈ Ω and z ∈ Z.
Unfortunately, the Lipschitz constant G can be very large or even infinite for some learning problems. Consider the simple least squares loss f (w; z) = 1 2 ( w, x − y) 2 with the gradient ∂f (w; z) = ( w, x − y)x. In this case the G-Lipschitzness of f requires to set G = sup w∈Ω sup z∈Z ( w, x − y)x 2 , which is infinite if Ω is unbounded. As another example, the Lipschitz constant of deep neural networks can be prohibitively large. In this case, existing stability bounds fail to yield meaningful generalization bounds. Furthermore, another critical assumption in the literature is the L-smoothness on f , i.e. for any z and w,w In this section, we will remove the boundedness assumption on the gradients for differentiable loss functions, and establish stability and generalization only under the assumption where loss functions have Hölder continuous (sub)gradients-a condition much weaker than the strong smoothness [23,36,52]. Note that the loss functions can be non-differentiable if α = 0.
If (4.2) holds with α = 1, then f is smooth as defined by (4.1). If (4.2) holds with α = 0, then this amounts to saying that f is Lipschitz continuous as considered in Assumption 1. Examples of loss functions satisfying Definition 3 include the q-norm hinge loss f (w; z) = max(0, 1 − y w, x ) q for classification and the q-th power absolute distance loss f (w; z) = |y − w, x | q for regression [48], whose (sub)gradients are (q−1, C)-Hölder continuous for some C > 0 if q ∈ [1,2]. If q = 1, we get the hinge loss and absolute distance loss with wide applications in machine learning and statistics.

On-average model stability
The key to remove the bounded gradient assumption and the strong smoothness assumption is the introduction of a novel stability measure which we refer to as the on-average model stability. We use the term "on-average model stability" to differentiate it from on-average stability in Kearns & Ron [19], Shalev-Shwartz et al. [45] as we measure stability on model parameters w instead of function values. Intuitively, on-average model stability measures the on-average sensitivity of models by traversing the perturbation of each single coordinate.
Definition 4 (On-average Model Stability). Let S = {z 1 , . . . , z n } and S = {z 1 , . . . ,z n } be drawn independently from ρ. For any i = 1, . . . , n, define S (i) = {z 1 , . . . , z i−1 ,z i , z i+1 , . . . , z n } as the set formed from S by replacing the i-th element withz i . We say a randomized algorithm A is ℓ 1 on- In the following theorem, we build the connection between generalization in expectation and the on-average model stabilities to be proved in Appendix B. Although the generalization by ℓ 1 on-average model stability requires Assumption 1, it is removed for ℓ 2 on-average model stability. We introduce a free parameter γ to tune according to the property of problems. Note we require a convexity assumption in Part (c) by considering non-smooth loss functions. Let c α,1 = (1 + 1/α)    (b) If for any z, the function w → f (w; z) is nonnegative and L-smooth, then (c) If for any z, the function w → f (w; z) is nonnegative, convex and w → ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1], then Remark 1. We explain here the benefit of ℓ 2 on-average model stability. If A is ℓ 2 on-average model ǫ-stable, then we take γ = 2LE F S (A(S)) /ǫ in Part (b) and derive In particular, if the output model has a small empirical risk in the sense of E F S (A(S) That is, our relationship between the generalization and ℓ 2 on-average stability allows us to exploit small risk of output model to get a generalization bound with an improved dependency on the stability measure ǫ. As a comparison, the discussions based on uniform stability (Lemma 1) and the ℓ 1 on-average model stability (Part (a)) only show , which fail to exploit the low-noise condition. We can also take .
If E[F S (A(S))] is small, this also implies an improved dependency of the generalization bound on ǫ.

Strongly smooth case
To justify the effectiveness of the on-average model stability, we first consider its application to learning with smooth loss functions. We first study stability and then generalization.
Stability bounds. The following theorem to be proved in Appendix C.1 establishes on-average model stability bounds in the smooth setting. A key difference from the existing stability bounds is that the uniform Lipschitz constant G is replaced by empirical risks. Since we are minimizing empirical risks by SGD, it is expected that these risks would be significantly smaller than the uniform Lipschitz constant. Actually we will control the weighted sum of empirical risks by tools in analyzing optimization errors. In the optimistic case with F (w * ) = 0, we expect E S, S,A [F S (w t )] = O(1/t), and in this case the discussion based on on-average model stability would imply significantly better generalization bounds. The idea of introducing a parameter p in (4.4) is to make (1 + p/n) t ≤ e by setting p = n/t, where e is the base of the nature logarithm.
Theorem 3 (Stability bounds). Assume for all z ∈ Z, the map w → f (w; z) is nonnegative, convex and L-smooth. Let S, S and S (i) be constructed as Definition 4. Let {w t } and {w (i) t } be produced by (3.3) with η t ≤ 2/L based on S and S (i) , respectively. Then for any p > 0 we have and [20] developed an interesting on-average stability bound O(σ n t j=1 η j ) under the bounded variance assumption E S,z ∂f (w t ; z) − ∂F (w t ; z) 2 2 ≤σ 2 for all t. Although this bound successfully replaces the uniform Lipschitz constant by the milder uniform variance constant σ, the corresponding generalization analysis still requires a bounded gradient assumption. A nice property of the stability bound in Kuzborskij & Lampert [20] is that it depends on the quality of the initialization, i.e., the stability increases if we start with a good model. Our stability bound also enjoys this property. As we can see from Theorem 3, the stability increases if we find good models with small optimization errors in the optimization process. This illustrates a key message that optimization is beneficial to improve the generalization.
Remark 3. The stability bounds in Theorem 3 can be extended to the non-convex case. Specifically, let assumptions of Theorem 3, except the convexity of w → f (w; z), hold. Then for any p > 0 one gets (see Proposition C.3) This result improves the recurrence relationship in Hardt et al. [18] for uniform stability by replacing the uniform Lipschitz constant with empirical risks.
Generalization bounds. We now establish generalization bounds based on ℓ 2 on-average model stability. This approach not only removes a bounded gradient assumption, but also allows us to fully exploit the smoothness of loss functions to derive bounds depending on the behavior of the best model w * . As we will see in Corollary 5, Theorem 4 interpolates between O(1/ √ n) bound in the "pessimistic" case (F (w * ) > 0) and the O(1/n) bound in the "low-noise" case (F (w * ) = 0) [41,47], which is becoming more and more interesting in the deep learning era with possibly more parameters than training examples. To our best knowledge, this is the first optimistic bound for SGD based on a stability approach. Eq. (4.6) still holds if F (w * ) = O(1/n). The proofs are given in Appendix C.2.
Theorem 4 (Generalization bounds). Assume for all z ∈ Z, the function w → f (w; z) is nonnegative, convex and L-smooth. Let {w t } be produced by (3.3) with nonincreasing step sizes satisfying η t ≤ 1/(2L). If γ ≥ 1, then Corollary 5. Assume for all z ∈ Z, the function w → f (w; z) is nonnegative, convex and L-smooth.
from which one can derive the O(1/ √ n) bound at best even if F (w * ) = 0. The improvement of our bounds over (4.7) is due to the consideration of on-average model stability bounds involving empirical risks (we use the same optimization error bounds in these two approaches). Based on the on-average stability bound in Kuzborskij & Lampert [20], one can derive a generalization bound similar to (4.7) with G 2 replaced by Gσ (σ is the uniform variance constant in Remark 2), which also could not yield Remark 5. We compare here our results with some fast bounds for SGD. Some fast convergence rates of SGD were recently derived for SGD under low noise conditions [2,30,47] or growth conditions relating stochastic gradients to full gradients [49]. The discussions there mainly focused on optimization errors, which are measured w.r.t. the iteration number t. As a comparison, our fast rates measured by n are developed for generalization errors of SGD (Part (b) of Corollary 5), for which we need to tradeoff optimization errors and estimation errors by stopping at an appropriate iteration number. Fast generalization bounds are also established for the specific least squares based on an integral operator approach [10,25,32,38]. However, these discussions heavily depend on the structure of the square loss and require capacity assumptions in terms of the decay rate of eigenvalues for the associated integral operator. As a comparison, we consider general loss functions and do not impose a capacity assumption.

Non-smooth case
As a further application, we apply our on-average model stability to learning with non-smooth loss functions, which have not been studied in the literature.

Stability bounds.
The following theorem to be proved in Appendix D.1 establishes stability bounds. As compared to (4.4), the stability bound below involves an additional term O( which is the cost we pay by relaxing the smoothness condition to a Hölder continuity of (sub)gradients. It is worth mentioning that our stability bounds apply to non-differentiable loss functions including the popular hinge loss.
Theorem 6 (Stability bounds). Assume for all z ∈ Z, the map w → f (w; z) is nonnegative, convex and ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1). Let S, S and S (i) be constructed in Definition 4. Let {w t } and {w (i) t } be produced by (3.3) based on S and S (i) , respectively. Then Generalization bounds. We now present generalization bounds for learning by loss functions with Hölder continuous (sub)gradients, which are specific instantiations of a general result (Proposition D.4) stated and proved in Appendix D.2.

Theorem 7 (Generalization bounds). Assume for all
(a) If α ≥ 1/2, we can take θ = 1/2 and T ≍ n to derive E S,A [F (w (1) Remark 6. Although relaxing smoothness affects stability by introducing O( t j=1 η 2 1−α j ) in the stability bound, we achieve a generalization bound similar to the smooth case with a similar computation cost if α ≥ 1/2. For α < 1/2, a minimax optimal generalization bound O(n − 1 2 ) [1] can be also achieved with more computation cost as T ≍ n 2−α 1+α . In particular, if α = 0 we develop the optimal generalization bounds O(n − 1 2 ) for SGD with T ≍ n 2 iterations. To our best knowledge, this gives the first generalization bounds for SGD with non-smooth loss functions (e.g., hinge loss) based on stability analysis. Analogous to the smooth case, we can derive generalization bounds better than O(n − 1 2 ) in the case with low noises. To our best knowledge, this is the first optimistic generalization bound for SGD with non-smooth loss functions.
Remark 7. We can extend our discussion to ERM. If F S is σ-strongly convex and ∂f (w; z) is (α, L)-Hölder continuous, we can apply the on-average model stability to show (see Proposition D.6) where A(S) = arg min w∈R d F S (w). This extends the error bounds developed for ERM with stronglysmooth loss functions [44,47] to the non-smooth case, and removes the G-admissibility assumption in Bousquet & Elisseeff [5]. In a low-noise case with a small E S F (A(S)) , the discussion based on an on-average stability can imply optimistic generalization bounds for ERM.

Stability with Relaxed Convexity
We now turn to stability and generalization of SGD for learning problems where the empirical objective F S is convex but each loss function f (w; z) may be non-convex. For simplicity, we impose Assumption 1 here and use the arguments based on the uniform stability. The proofs of Theorem 8 and Theorem 9 are given in Appendix E.1.
Remark 8. The derivation of uniform stability bounds in Hardt et al. [18] is based on the nonexpansiveness of the operator w → w − ∂f (w; z), which requires the convexity of w → f (w; z) for all z. Theorem 8 relaxes this convexity condition to a milder convexity condition on By the proof, Theorem 8 holds if the convexity condition is replaced by ] for some C and all t ∈ N. As shown below, minimax optimal generalization bounds can be achieved for step sizes η t = η 1 t −θ for all θ ∈ (1/2, 1) as well as the step sizes η t ≍ 1/ √ T with T ≍ n.
Theorem 9. Let Assumption 1 hold. Assume for all z ∈ Z, the function w → f (w; z) is L-smooth.
Let {w t } t be produced by (3.3). Suppose for all S, F S is convex.
. Example: AUC Maximization. We now consider a specific example of AUC (Area under ROC curve) maximization where the objective function is convex but each loss function may be non-convex. As a widely used method in imbalanced classification (Y = {+1, −1}), AUC maximization was often formulated as a pairwise learning problem where the corresponding loss function involves a pair of training examples [17,56]. Recently, AUC maximization algorithms updating models with a single example per iteration were developed [27,33,53]. Specifically, AUC maximization with the square loss can be formulated as the minimization of the following objective function where p = Pr{Y = 1} is the probability of a example being positive. Let An interesting property is that (5.2) involves only a single example z. This observation allows Natole et al. [33] to develop a stochastic algorithm as (3.3) to solve (5.1). However, for each z, the function z → f (w; z) is non-convex since the associated Hessian matrix may not be positively definite. It is clear that its expectation F is convex.
6 Stability with Relaxed Strong Convexity

Stability and generalization errors
Finally, we consider learning problems with strongly convex empirical objectives but possibly nonconvex loss functions. Theorem 10 provides stability bounds, while the minimax optimal generalization bounds O(1/(σn)) are presented in Theorem 11. The proofs are given in Appendix F. Theorem 10. Let Assumptions in Theorem 8 hold. Suppose for all S ⊂ Z, F S is σ S -strongly convex. Then, there exists a constant t 0 such that for SGD with η t = 2/((t + t 0 )σ S ) we have Remark 9. Under the assumption w → f (w, z) is σ-strongly convex and smooth for all z, it was shown that [18]. Indeed, this strong convexity condition is used to show that the operator w → w − ∂f (w; z) is contractive. We relax the strong convexity of f (w; z) to the strong convexity of F S . Our stability bound holds even if w → f (w; z) is non-convex. If t ≍ n, then our stability bound coincides with the one in Hardt et al. [18] up to a constant factor.
Theorem 11. Let Assumption 1 hold. Assume for all z ∈ Z, the function w → f (w; z) is L-smooth. Suppose for all S ⊂ Z, F S is σ S -strongly convex. Then, there exists some t 0 such that for SGD with η t = 2/((t + t 0 )σ S ) and T ≍ n we have

Application: least squares regression
We now consider an application to learning with the least squares loss, where f (w; where Π Ω (w) = min{R/ w 2 , 1}w. Note that each individual loss function f (w t ; z t ) is non-strongly convex. However, as we will show below, the empirical objective satisfies a strong convexity on a subspace containing the iterates {w t }. For any S = {z 1 , . . . , z n } let C S = 1 n n i=1 x i x ⊤ i be the empirical covariance matrix and σ ′ S be the minimal positive eigenvalue of C S . Then it is clear from (6.1) that {w t } t belongs to the range of C S . 1 Let S ⊂ Z n differ from S by a single example. For simplicity, we assume S and S differ by the first example and denote S = {z 1 , z 2 , . . . , z n }. We construct a setS = {0, z 2 , . . . , z n }. Let {w t }, {w t } and {w t } be the sequence by (6.1) based on S, S andS, respectively. Then our previous discussion implies that where we denote by Range(C) the range of a matrix C. It follows that w t −w t andw t −w t are orthogonal to the kernel of C S and C S , respectively. Therefore, . As we will see in the proof, Theorem 10 holds if only the following local strong convexity holds, i.e., Therefore, we can apply Theorem 10 with S =S and A similar inequality also holds for E A w t+1 −w t+1 2 , which together with the subadditivity of · 2 immediately gives the following stability bound on E A [ w t+1 −w t+1 2 ].

Conclusions
In this paper, we study stability and generalization of SGD by removing the bounded gradient assumptions, and relaxing the smoothness assumption and the convexity requirement of each loss function in the existing analysis. We introduce a novel on-average model stability able to capture the risks of SGD iterates, which implies fast generalization bounds in the low-noise case and stability bounds for learning with even non-smooth loss functions. For all considered problems, we show that our stability bounds can imply minimax optimal generalization bounds by balancing optimization and estimation errors. We apply our results to practical learning problems to justify the superiority of our approach over the existing stability analysis. Our results can be extended to stochastic proximal gradient descent, high-probability bounds and SGD without replacement (details are given in Appendix G). In the future, it would be interesting to study stability bounds for other stochastic optimization algorithms, e.g., Nesterov's accelerated variants of SGD [35].

Acknowledgement
The work of Y. Lei is supported by the National Natural Science Foundation of China (Grant Nos. 61806091, 11771012) and the Alexander von Humboldt Foundation. The work of Y. Ying is supported by the National Science Foundation (NSF) under Grant No. #1816227.

A Optimization Error Bounds
For a full picture of generalization errors, we need to address the optimization errors. This is achieved by the following lemma. Parts (a) and (b) consider the convex and strongly convex empirical objectives, respectively (we make no assumptions on the convexity of each f (·, z)). Note in Parts (c) and (d), we do not make a bounded gradient assumption. As an alternative, we require convexity of f (·; z) for all z. An appealing property of Parts (c) and (d) , which is a requirement for developing fast rates in the case with low noises. Our discussion on optimization errors requires to use a self-bounding property for functions with Hölder continuous (sub)gradients, which means that gradients can be controlled by function values. The case α = 1 was established in Srebro et al. [47]. The case α ∈ (0, 1) was established in Ying & Zhou [52]. The case α = 0 follows directly from Definition 3. Define Lemma A.2. (a) Let {w t } t be produced by (3.3) and Assumption 1 hold. If F S is convex, then for all t ∈ N and w ∈ Ω (b) Let F S be σ S -strongly convex and Assumption 1 hold. Let t 0 ≥ 0 and {w t } t be produced by (3.3) with η t = 2/(σ S (t + t 0 )). Then for all t ∈ N and w ∈ Ω (c) Assume for all z ∈ Z, the function w → f (w; z) is nonnegative, convex and L-smooth. Let {w t } t be produced by (3.3) with η t ≤ 1/(2L). If the step size is nonincreasing, then for all t ∈ N and w ∈ Ω independent of the SGD algorithm A (d) Assume for all z ∈ Z, the function w → f (w; z) is nonnegative, convex, and ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1). Let {w t } t be produced by (3.3) with nonincreasing step sizes. Then for all t ∈ N and w ∈ Ω independent of the SGD algorithm A we have Proof. Parts (a) and (b) can be found in the literature [21,34]. We only prove Parts (c) and (d). We first prove Part (c). The projection operator Π Ω is non-expansive, i.e., By the SGD update (3.3), (A.3), convexity and Lemma A.1, we know where the last inequality is due to η t ≤ 1/(2L). It then follows that Multiplying both sides by η t and using the assumption η t+1 ≤ η t , we know Taking a summation of the above inequality gives ( Taking an expectation w.r.t. A gives (note w j is independent of i j ) On the other hand, taking an expectation w.r.t. i t over both sides of (A.4) shows Taking an expectation on both sides followed with a summation, we get where the last step is due to (A.5). The proof is complete since w is independent of A.
We now prove Part (d). Analogous to (A.4), one can show for loss functions with Hölder continuous (sub)gradients (Lemma A.1) By the Young's inequality ab ≤ p −1 |a| p + q −1 |b| q , a, b ∈ R, p, q > 0 with p −1 + q −1 = 1, (A.7) we know (notice the following inequality holds trivially if α = 0) Combining the above two inequalities together, we get Multiplying both sides by η t and using η t+1 ≤ η t , we derive Taking a summation of the above inequality gives According to the Jensen's inequality and the concavity of where in the last step we have used (A.8). Taking an expectation on both sides of (A.6), we know Taking a summation of the above inequality gives where we have used (A.9) and the concavity of x → x 2α 1+α in the last step. The proof is complete by noting the independence between w and A.

B Proofs on Generalization by On-average Model Stability
To prove Theorem 2, we introduce an useful inequality for L-smooth functions w → f (w; z) [35] f (w; z) ≤ f (w; z) + w −w, ∂f (w; z) + L w −w 2 2 2 . (B.1) Proof of Theorem 2. Due to the symmetry, we know where the last identity holds since A(S (i) ) is independent of z i . Under Assumption 1, it is then clear that This proves Part (a). We now prove Part (b). According to (B.1) due to the L-smoothness of f and (B.2), we know According to the Schwartz's inequality we know where the last inequality is due to the self-bounding property of smooth functions (Lemma A.1).
Combining the above two inequalities together, we derive

By the Schwartz's inequality and Lemma A.1 we know
Combining the above two inequalities together, we get 1+α is concave and z i is independent of A(S (i) ), we know A combination of the above two inequalities then gives the stated bound in Part (c). The proof is complete.

C Proof on Learning without Bounded Gradients: Strongly
Smooth Case

C.1 Stability bounds
A key property on establishing the stability of SGD is the non-expansiveness of the gradient-update operator established in the following lemma.
Lemma C.1 (Hardt et al. 18). Assume for all z ∈ Z, the function w → f (w; z) is convex and L-smooth. Then for η ≤ 2/L we know Based on Lemma C.1, we establish stability bounds of models for SGD applied to two sets differing by a single example.
Proof. If i t = i, we know the updates of w t+1 and w (i) t+1 are based on stochastic gradients calculated with the same example z it . By Lemma C.1 we then get where the second inequality follows from the sub-additivity of · 2 and the last inequality is due to Lemma A.1 on the self-bounding property of smooth functions. We first prove Eq. (C.1). Since i t is drawn from the uniform distribution over {1, . . . , n}, we can combine Eqs. (C.3) and (C.5) to derive Since z i andz i follow from the same distribution, we know It then follows that Taking a summation of the above inequality and using w 1 = w (i) 1 then give (C.1). We now turn to (C.2). For the case i t = i, it follows from (C.4) and the standard inequality (a + b) 2 ≤ (1 + p)a 2 + (1 + 1/p)b 2 that where the last inequality is due to Lemma A.1. Combining (C.3), the above inequality together and noticing the distribution of i t , we derive Analogous to (C.6), we have and get Multiplying both sides by (1 + p/n) −(t+1) yields that Taking a summation of the above inequality and using w 1 = w (i) 1 , we get The stated bound then follows.
Proof of Theorem 3. We first prove (4.3). According to Lemma C.2 (Eq. (C.1)), we know It then follows from the concavity of the square-root function and the Jensen's inequality that This proves (4.3). We now turn to (4.4). It follows from Lemma C.2 (Eq. (C.2)) that The proof is complete.

Proposition C.3 (Stability bounds for non-convex learning)
. Let Assumptions of Theorem 3 hold except that we do not require the convexity of w → f (w; z). Then for any p > 0 we have Proof. If i t = i, then by the L-smoothness of f we know If i t = i, then analogous to (C.7), one can get By the uniform distribution of i t ∈ {1, 2, . . . , n} we can combine the above inequality and (C.9) to derive This together with (C.8) implies It then follows that The proof is complete.

C.2 Generalization bounds
We now prove generalization bounds for SGD.
Proof of Theorem 4. According to Part (c) of Lemma A.2 with w = w * , we know the following inequality Let A(S) be the (t + 1)-th iterate of SGD applied to the dataset S. We plug (4.4) into Part (b) of Theorem 2, and derive We can plug (A.5) with w = w * into the above inequality, and derive We choose p = n/T , then (1 + p/n) T −1 = (1 + 1/T ) T −1 ≤ e and therefore the following inequality holds for all t = 1, . . . , Multiplying both sides by η t+1 followed with a summation gives Putting (C.10) into the above inequality then gives The stated inequality then follows from Jensen's inequality. The proof is complete.
Proof of Corollary 5. We first prove Part (a). For the chosen step size, we know T t=1 The stated bound (4.5) then follows from Theorem 4, γ = √ n and (C.11). We now prove Part (b). The stated bound (4.6) then follows from Theorem 4, F (w * ) = 0 and γ = 1. The proof is complete.

D Proof on Learning without Bounded Gradients: Non-smooth Case D.1 Stability bounds
Theorem 6 is a direct application of the following general stability bounds with p = n/t. Therefore, it suffices to prove Theorem D.1.
Theorem D.1. Assume for all z ∈ Z, the function w → f (w; z) is nonnegative, convex and ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1). Let S, S and S (i) be constructed as Definition 4 and c α,3 = respectively. Then for any p > 0 we have We require several lemmas to prove Theorem D.1. The following lemma establishes the co-coercivity of gradients for convex functions with Hölder continuous (sub)gradients. The case α = 1 can be found in Nesterov [35]. The case α ∈ (0, 1) can be found in Ying & Zhou [52]. The case α = 0 follows directly from the convexity of f . Lemma D.2. Assume for all z ∈ Z, the map w → f (w; z) is convex, and w → ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1]. Then for all w,w we have (D. 2) The following lemma controls the expansive behavior of the operator w → w−η∂f (w; z) for convex f with Hölder continuous (sub)gradients.
We now consider the case α > 0. According to Lemma D.2, we know where we have used Young's inequality (A.7). Plugging the above inequality back into (D.3), we derive Combining the above two cases together, we get the stated bound with the definition of c α,3 given in Theorem D.1. The proof is complete.
Proof of Theorem D.1. For the case i t = i, it follows from Lemma D.3 that If i t = i, by (C.4) and the standard inequality (a + b) 2 ≤ (1 + p)a 2 + (1 + 1/p)b 2 , we get Combining the above two inequalities together, using the self-bounding property (Lemma A.1) and noticing the distribution of i t , we derive Analogous to (C.6), we know and therefore Multiplying both sides by (1 + p/n) −(t+1) gives Taking a summation of the above inequality and using w 1 = w (i) 1 , we derive We can take an average over i and get It then follows from the concavity of the function x → x 2α 1+α and the Jensen's inequality that The stated inequality then follows from the definition of F S . The proof is complete.

D.2 Generalization errors
Theorem 7 can be considered as an instantiation of the following proposition on generalization error bounds with specific choices of γ, T and θ. In this subsection, we first give the proof of Theorem 7 based on Proposition D.4, and then turn to the proof of Proposition D.4.
To prove Proposition D.4, we first introduce an useful lemma to address some involved series.
Proof. We first prove (D.6). For the step size sequence η t = cT −θ , we have where we have used the subadditivity of x → x 2α 1+α , the identity in the second step and θ ≥ (1 − α)/2 in the third step (the third term is dominated by the first term). This shows (D.6). We now consider (D.7). Taking an expectation over both sides of (A.8) with w = w * , we get According to the Jensen's inequality and the concavity of where we have used (D.6) in the last step. This shows (D.7). Finally, we show (D.8). Since we consider step sizes η t = cT −θ , it follows from (D.7) that This proves (D.8) and finishes the proof.
Proof of Proposition D. 4. Since E S [F S (w * )] = F (w * ), we can decompose the excess generalization error into an estimation error and an optimization error as follows Our idea is to address separately the above estimation error and optimization error. We first address estimation errors. Plugging (D.1) back into Theorem 2 (Part (c)) with A(S) = w t+1 , we derive By the concavity and sub-additivity of x → x 2α 1+α , we know , 0} for all j ∈ N. It then follows from p = n/T that 1+α .
Solving the above inequality of δ t+1 gives the following inequality for all t ≤ T .
It then follows from the definition of δ t that (note n = O(T )) .
By η t = cT −θ , (D.7) and (D.8), we further get We now consider optimization errors. By Lemma A.2 (Part (d) with w = w * ) and the concavity of where we have used (D.6) in the last step.
Plugging the above optimization error bound and the estimation error bound (D.10) back into the error decomposition (D.9), we finally derive the following generalization error bounds The stated inequality then follows from the convexity of F . The proof is complete.

D.3 Empirical Risk Minimization with Strongly Convex Objectives
In this section, we present an optimistic bound for ERM with strongly convex objectives based on the ℓ 2 on-average model stability. We consider nonnegative and convex loss functions with Hölder continuous (sub)gradients.
Proposition D.6. Assume for any z, the function w → f (w; z) is nonnegative, convex and w → ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1]. Let A be the ERM algorithm, i.e., A(S) = arg min w∈R d F S (w). If for all S, F S is σ-strongly convex, then Proof. Let S and S (i) , i = 1, . . . , n, be constructed as Definition 4. Due to the σ-strong convexity of F S (i) and ∂F S (i) (A(S (i) )) = 0 (necessity condition for the optimality of A(S (i) )), we know . Taking a summation of the above inequality yields According to the definition of S (i) , we know Taking an expectation and dividing both sides by n 2 give (A(S) is independent of S) Furthermore, by symmetry we know Plugging the above identity and (D.12) back into (D.11) gives σ 2n We can now apply Part (c) of Theorem 2 to show the following inequality for all γ > 0 (notice A is a deterministic algorithm) Taking γ = nσ/2, we derive from which we can derive the stated inequality. The proof is complete.

E.1 Stability and generalization errors
For any convex g, we have [35] w −w, ∇g(w) − ∇g(w) ≥ 0, w,w ∈ R d . (E.1) Proof of Theorem 8. Without loss of generality, we can assume that S and S differ by the first example, i.e., z 1 =z 1 and z i =z i , i = 1. According to the update rule (3.3) and (A.3), we know We first study the term ∂f (w t ; z it )−∂f (w t ;z it ) 2 . The event i t = 1 happens with probability 1−1/n, and in this case it follows from the smoothness of f that (z it =z it ) The event i t = 1 happens with probability 1/n, and in this case Therefore, we get Therefore, by (E.1) we derive Plugging (E.3) and the above inequality back into (E.2), we derive and therefore (E.5) By the above recurrence relationship and w 1 =w 1 , we derive Since the above inequality holds for all t ∈ N and the right-hand side is an increasing function of t, we get It then follows that (note Solving the above quadratic function of max 1≤j≤t+1 E A wj −wj 2 The proof is complete.
To prove Theorem 9, we require a basic result on series.
Lemma E.1. We have the following elementary inequalities.
(a) If θ ∈ (0, 1), We denote by ǫ stab (A, n) the infimum over all ǫ for which (3.2) holds, and omit the tuple (A, n) when it is clear from the context. Proof of Theorem 9. For the step sizes considered in both Part (a) and Part (b), one can check that T t=1 η 2 t can be upper bounded by a constant independent of T . Therefore, C t < C for all t = 1, . . . , T and a universal constant C. We can apply Lemma A.2 (Part (a)) on optimization errors to get (E.6) By the convexity of norm, we know . This together with the error decomposition (3.1), Lemma 1 and the optimization error bound (E.6) shows (E.7) For the step sizes η t = η 1 t −θ with θ ∈ (1/2, 1), we can apply Lemma E.1 to show This proves the first part. Part (b) follows by plugging (C.11) into (E.7). The proof is complete.

F Proofs on Stability with Relaxed Strong Convexity
Proof of Theorem 10. Due to the σ S -strong convexity of F S , we can analyze analogously to (E.4) to derive Therefore, analogous to the derivation of (E.5) we can derive We find t 0 ≥ 4L 2 /σ 2 S . Then η t ≤ σ S /(2L 2 ) and it follows that Multiplying both sides by (t + t 0 )(t + t 0 − 1) yields Taking a summation of the above inequality and using w 1 =w 1 then give It then follows The stated bound then follows from the elementary inequality √ a + b ≤ √ a + √ b for a, b ≥ 0. The proof is complete.
Proof of Theorem 11. By the convexity of norm, we know Proof. It suffices to show that the kernel of C S is the orthogonal complement of V = span{x 1 , . . . , x n } (we denote span{x 1 , . . . , x n } the linear span of x 1 , . . . , x n ). Indeed, for any x in the kernel of C S , we know C S x = 0 and therefore x ⊤ C S x = 1 n n i=1 (x ⊤ i x) 2 = 0, from which we know that x must be orthogonal to V . Furthermore, for any x orthogonal to V , it is clear that C S x = 0, i.e., x belongs to the kernel of C S . The proof is complete.

G Extensions
In this section, we present some extensions of our analyses. We consider three extensions: extension to stochastic proximal gradient descent, extension to high probability analysis and extension to SGD without replacement.

G.1 Stochastic proximal gradient descent
Our discussions can be directly extended to study the performance of stochastic proximal gradient descent (SPGD). Let r : R d → R + be a convex regularizer. SPGD updates the models by w t+1 = Prox ηtr w t − η t ∂f (w t ; z it ) , where Prox g (w) = arg minw ∈R d g(w) + 1 2 w −w 2 2 is the proximal operator. SPGD has found wide applications in solving optimization problems with a composite structure [37]. It recovers the projected SGD as a specific case by taking an appropriate r. Our stability bounds for SGD can be trivially extend to SPGD due to the non-expansiveness of proximal operators: Prox g (w) − Prox g (w) 2 ≤ w −w 2 , ∀w,w if g is convex.

G.2 Stability bounds with high probabilities
We can also extend our stability bounds stated in expectation to high-probability bounds, which would be helpful to understand the fluctuation of SGD w.r.t. different realization of random indices.
Proposition G.1. Let Assumption 1 hold. Assume for all z ∈ Z, the function w → f (w; z) is convex and w → ∂f (w; z) is (α, L)-Hölder continuous with α ∈ [0, 1]. Let S = {z 1 , . . . , z n } and S = {z 1 , . . . ,z n } be two sets of training examples that differ by a single example. Let {w t } t and {w t } t be produced by (3.3) based on S and S, respectively, and δ ∈ (0, 1). If we take step size η j = ct −θ for j = 1, . . . , t and c > 0, then with probability at least 1 − δ High-probability generalization bounds can be derived by combining the above stability bounds and the recent result on relating generalization and stability in a high-probability analysis [6,15].
To prove Proposition G.1, we need to introduce a special concentration inequality called Chernoff's fs bound for a summation of independent Bernoulli random variables [3].
Proof of Proposition G.1. Without loss of generality, we can assume that S and S differ by the first example, i.e., z 1 =z 1 and z i =z i for i = 1. If i t = 1, we can apply Lemma D.3 and (A.3) to derive w t+1 −w t+1 2 ≤ w t − η t ∂f (w t ; z it ) −w t + η t ∂f (w t ; z it ) 2 ≤ w t −w t 2 + c α,3 η If i t = 1, we know w t+1 −w t+1 2 ≤ w t − η t ∂f (w t ; z 1 ) −w t + η t ∂f (w t ;z 1 ) 2 ≤ w t −w t 2 + 2η t G.
Combining the above two cases together, we derive w t+1 −w t+1 2 ≤ w t −w t 2 + c α,3 η Since we draw (i k 1 , . . . , i k n ) from the uniform distribution of all permutations, i k takes an equal probability to each 1, . . . , n. Therefore, we can take expectations over A to derive We can take a summation of the above inequality from k = 1 to K to derive the stated bound. The proof is complete.