Stability of Stochastic Gradient Descent on Nonsmooth Convex Losses

Raef Bassily, Vitaly Feldman, Cristóbal Guzmán, Kunal Talwar

Introduction

Successful applications of a machine learning algorithm require the algorithm to generalize well to unseen data. Thus understanding and bounding the generalization error of machine learning algorithms is an area of intense theoretical interest and practical importance. The single most popular approach to modern machine learning relies on the use of continuous optimization techniques to optimize the appropriate loss function, most notably the stochastic (sub)gradient descent (SGD) method. Yet the generalization properties of SGD are still not well understood.

where the expectation is taken with respect to the randomness of the sample S\mathbf{S} and internal randomness of A{\mathcal{A}}. A standard way to bound the excess risk is given by its decomposition into optimization error (a.k.a. training error) and generalization error (see eqn. (3) in Sec. 2). The optimization error can be easily measured empirically but assessing the generalization error requires access to fresh samples from the same distribution. Thus bounds on the generalization error lead directly to provable guarantees on the excess population risk.

Classical analysis of SGD allows obtaining bounds on the excess population risk of one pass SGD. In particular, with an appropriately chosen step size, SGD gives a solution with expected excess population risk of O(1/n)O(1/\sqrt{n}) and this rate is optimal (Nemirovsky and Yudin,, 1983). However, this analysis does not apply to multi-pass SGD that is ubiquitous in practice.

In an influential work, Hardt et al., (2016) gave the first bounds on the generalization error of general forms of SGD (such as those that make multiple passes over the data). Their analysis relies on algorithmic stability, a classical tool for proving bounds on the generalization error. Specifically, they gave strong bound on the uniform stability of several variants of SGD on convex and smooth losses (with 2/η2/\eta-smoothness sufficing when all the step sizes are at most η\eta). Uniform stability bounds the worst case change in loss of the model output by the algorithm on the worst case point when a single data point in the dataset is replaced (Bousquet and Elisseeff,, 2002). Formally, for a randomized algorithm A{\cal A}, loss functions f(⋅,z)f(\cdot,z) and S≃S′S\simeq S^{\prime} and z∈Zz\in\mathcal{Z}, let γA(S,S′,z):=f(A(S),z)−f(A(S′),z)\gamma_{\cal A}(S,S^{\prime},z):=f({\mathcal{A}}(S),z)-f({\mathcal{A}}(S^{\prime}),z), where S≃S′S\simeq S^{\prime} denotes that the two datasets differ only in a single data point. We say A{\cal A} is γ\gamma-uniformly stable if

where the expectation is over the internal randomness of A.\cal{A}. Stronger notions of stability can also be considered, e.g., bounding the probability – over the internal randomness of A\cal{A} – that γA(S,S′,z)>γ\gamma_{\cal A}(S,S^{\prime},z)>\gamma. Using stability, Hardt et al., (2016) showed that several variants of SGD simultaneously achieve the optimal tradeoff between the excess empirical risk and stability with both being O(1/n)O(1/\sqrt{n}). Several works have used this approach to derive new generalization properties of SGD (London,, 2017; Chen et al.,, 2018; Feldman and Vondrák,, 2019).

Uniform stability is also closely related to the notion of differential privacy (DP). DP upper bounds the worst case change in the output distribution of an algorithm when a single data point in the dataset is replaced (Dwork et al., 2006b, ). This connection has been exploited in the design of several DP algorithms for SCO. In particular, bounds on the uniform stability of SGD from (Hardt et al.,, 2016) have been crucial in the design and analysis of new DP-SCO algorithms (Wu et al.,, 2017; Dwork and Feldman,, 2018; Feldman et al.,, 2018; Bassily et al.,, 2019; Feldman et al.,, 2020).

We establish tight bounds on the uniform stability of the (stochastic) subgradient descent method on nonsmooth convex losses. These results demonstrate that in the nonsmooth case SGD can be substantially less stable. At the same time we show that SGD has strong stability properties even in the regime when its iterations can be expansive.

This notion is implicit in existing analyses of uniform stability (Bousquet and Elisseeff,, 2002; Shalev-Shwartz et al.,, 2010; Hardt et al.,, 2016) and was explicitly defined by Liu et al., (2017). In this work, we prove stronger – high probability – upper bounds on the random variable δA(S,S′):=∥A(S)−A(S′)∥\delta_{\cal A}(S,S^{\prime}):=\left\|{\mathcal{A}}(S)-{\mathcal{A}}(S^{\prime})\right\|,In fact, for both GD and fixed-permutation SGD we can obtain w.p. 1 upper bounds on δA(S,S′)\delta_{\cal A}(S,S^{\prime}), whereas for sampling-with-replacement SGD, we obtain a high-probability upper bound. and we provide matching lower bounds for the weaker – in expectation – notion of UAS (1). A summary of our bounds is in Table 1. For simplicity, they are provided for constant step size; general step sizes (for upper bounds) are provided in Section 3.

Compared to the smooth case (Hardt et al.,, 2016), the main difference is the presence of the additional ηT\eta\sqrt{T} term. This term has important implications for the generalization bounds derived from UAS. The first one is that the standard step size η=Θ(1/n)\eta=\Theta(1/\sqrt{n}) used in single pass SGD leads to a vacuous stability bound. Unfortunately, as shown by our lower bounds, this is unavoidable (at least in high dimension). However, by decreasing the step size and increasing the number of steps, one obtains a variant of SGD with nearly optimal balance between the UAS and the excess empirical risk.

We highlight two major consequences of our bounds:

Generalization bounds for multi-pass nonsmooth SGD. We prove that the generalization error of multi-pass SGD with KK passes is bounded by O((Kn+K)η)O((\sqrt{Kn}+K)\eta). This result can be easily combined with training error guarantees to provide excess risk bounds for this algorithm. Since training error can be measured directly, our generalization bounds would immediately yield strong guarantees on the excess risk in practical scenarios where we can certify small training error.

Differentially private stochastic convex optimization for non-smooth losses. We show that a variant of standard noisy SGD (Bassily et al.,, 2014) with constant step size and n2n^{2} iterations yields the optimal excess population risk O\big{(}\frac{1}{\sqrt{n}}+\frac{\sqrt{d\log(1/\beta)}}{\alpha n}\big{)} for convex nonsmooth losses under (α,β)(\alpha,\beta)-differential privacy. The best previous algorithm for this problem is substantially more involved: it relies on a multi-phase regularized SGD with decreasing step sizes and variable noise rates and uses O(n2log⁡(1/β))O(n^{2}\sqrt{\log(1/\beta)}) gradient computations (Feldman et al.,, 2020).

2 Overview of Techniques

Upper bounds. When gradient steps are nonexpansive, upper-bounding UAS requires simply summing the differences between the gradients on the neighboring datasets when the replaced data point is used (Hardt et al.,, 2016). This gives the bound of ηT/n\eta T/n in the smooth case.

By contrast, in the nonsmooth case, UAS may increase even when the gradient step is performed on the same function. As a result it may increase in every single iteration. However, we use the fact that the difference in the subgradients has negative inner product with the difference between the iterates themselves (by monotonicity of the subgradient). Thus the increase in distance satisfies a recurrence with a quadratic and a linear term. Solving this recurrence leads to our upper bounds.

Lower bounds. The lower bounds are based on a function with a highly nonsmooth behavior around the origin. More precisely, it is the maximum of linear functions plus a small linear drift that is controlled by a single data point. We show that, when starting the algorithm from the origin, the presence of the linear drift pushes the iterate into a trajectory in which each subgradient step is orthogonal to the current iterate. Thus, if d≥min⁡{T,1/η2}d\geq\min\{T,1/\eta^{2}\}, we get the Tη\sqrt{T}\eta increase in UAS. Our lower bounds are also robust to averaging of the iterates.

3 Other Related Work

Stability is a classical approach to proving generalization bounds pioneered by Rogers and Wagner, (1978); Devroye and Wagner, 1979a ; Devroye and Wagner, 1979b . It is based on analysis of the sensitivity of the learning algorithm to changes in the dataset such as leaving one of the data points out or replacing it with a different one. The choice of how to measure the effect of the change and various ways to average over multiple changes give rise to a variety of stability notions (e.g., (Bousquet and Elisseeff,, 2002; Mukherjee et al.,, 2006; Shalev-Shwartz et al.,, 2010)). ​ Uniform stability was introduced by Bousquet and Elisseeff, (2002) in order to derive general bounds on the generalization error that hold with high probability. These bounds have been significantly improved in a recent sequence of works (Feldman and Vondrák,, 2018, 2019; Bousquet et al.,, 2019). A long line of work focuses on the relationship between various notions of stability and learnability in supervised setting (see (Kearns and Ron,, 1999; Poggio et al.,, 2004; Shalev-Shwartz et al.,, 2010) for an overview). These works employ relatively weak notions of average stability and derive a variety of asymptotic equivalence results. Chen et al., (2018) establish limits of stability in the smooth convex setting, proving that accelerated methods must satisfy strong stability lower bounds. Stability-based data-dependent generalization bounds for continuous losses were studied in (Maurer,, 2017; Kuzborskij and Lampert,, 2018).

First applications of uniform stability in the context of stochastic convex optimization relied on the stability of the empirical minimizer for strongly convex losses (Bousquet and Elisseeff,, 2002). Therefore a natural approach to achieve uniform stability (and also UAS) is to add a strongly convex regularizer and solve the ERM to high accuracy (Shalev-Shwartz et al.,, 2010). Recent applications of this approach can be found for example in (Koren and Levy,, 2015; Charles and Papailiopoulos,, 2018; Feldman et al.,, 2020). In contrast, our approach does not require strong convexity and applies to all iterates of the SGD and not only to a very accurate empirical minimizer.

Classical approach to generalization relies on uniform convergence of empirical risk to population risk. Unfortunately, without additional structural assumptions on convex functions, a lower bound of Ω(d/n)\Omega(\sqrt{d/n}) on the rate of uniform convergence for convex SCO is known (Shalev-Shwartz et al.,, 2010; Feldman,, 2016). The dependence on the dimension dd makes the bound obtained via the uniform-convergence approach vacuous in the high-dimensional settings common in modern applications.

Differentially private convex optimization has been studied extensively for over a decade (see, e.g., (Chaudhuri and Monteleoni,, 2008; Chaudhuri et al.,, 2011; Jain et al.,, 2012; Kifer et al.,, 2012; Smith and Thakurta,, 2013; Bassily et al.,, 2014; Ullman,, 2015; Jain and Thakurta,, 2014; Talwar et al.,, 2015; Bassily et al.,, 2019; Feldman et al.,, 2020)). However, until recently, the research focused on minimization of the empirical risk. Population risk for DP-SCO was first studied by Bassily et al., (2014) who gave an upper bound of max⁡(d14n,dαn)\max\left(\tfrac{d^{\frac{1}{4}}}{\sqrt{n}},\tfrac{\sqrt{d}}{\alpha n}\right) (Bassily et al.,, 2014, Sec. F) on the excess risk. A recent work of Bassily et al., (2019) established that the optimal rate of the excess population risk for (α,β)(\alpha,\beta)-DP SCO algorithms is O\big{(}\frac{1}{\sqrt{n}}+\frac{\sqrt{d\log(1/\beta)}}{\alpha n}\big{)}. Their algorithms are relatively inefficient, especially in the nonsmooth case. Subsequently, Feldman et al., (2020) gave several new algorithms for DP-SCO with the optimal population risk. For sufficiently smooth losses, their algorithms use a linear number of gradient computations. In the nonsmooth case, as mentioned earlier, their algorithm requires O(n2log⁡(1/β))O(n^{2}\sqrt{\log(1/\beta)}) gradient computations and is significantly more involved than the algorithm shown here.

Notation and Preliminaries

Functions with these properties are guaranteed to be subdifferentiable. Moreover, in the convex case, property (2) is “almost” equivalent to having subgradients bounded as ∂f(x)⊆B(0,L)\partial f(x)\subseteq\mathcal{B}(0,L), for all x∈Xx\in\mathcal{X}.For equivalence to hold it is necessary that the function is well-defined and satisfies (2) over an open set containing X\mathcal{X}, see Thm. 3.61 in Beck, (2017). We will assume this is the case, which can be done w.l.o.g.. We denote the class of convex LL-Lipschitz functions as FX0(L)\mathcal{F}_{\mathcal{X}}^{0}(L). With slight abuse of notation, given a function f∈FX0(L)f\in\mathcal{F}_{\mathcal{X}}^{0}(L), we will denote by ∇f(x)\nabla f(x) an arbitrary choice of g∈∂f(x)g\in\partial f(x). In this work, we will focus on the class FX0(L)\mathcal{F}_{\mathcal{X}}^{0}(L) defined over a compact convex set X\mathcal{X}. Since the Euclidean radius of X\mathcal{X} is bounded by RR, we will assume that the range of these functions lies in [−RL,RL][-RL,RL].

Nonsmooth stochastic convex optimization: We study the standard setting of nonsmooth stochastic convex optimization

Here, D\mathcal{D} is an unknown distribution supported on a set Z\mathcal{Z}, and f(⋅,z)∈FX0(L)f(\cdot,z)\in\mathcal{F}_{\mathcal{X}}^{0}(L) for all z∈Zz\in\mathcal{Z}. In the stochastic setting, we assume access to an i.i.d. sample from D\mathcal{D}, denoted as S=(z1,…,zn)∼Dn\mathbf{S}=(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})\sim\mathcal{D}^{n}. Here, we will use the bold symbol S\mathbf{S} to denote a random sample from the unknown distribution. A fixed (not random) dataset from Zn\mathcal{Z}^{n} will be denoted as S=(z1,…,zn)∈ZnS=(z_{1},\ldots,z_{n})\in\mathcal{Z}^{n}.

A stochastic optimization algorithm is a (randomized) mapping A:Zn↦X{\cal A}:\mathcal{Z}^{n}\mapsto\mathcal{X}. When the algorithm is randomized, A(S){\cal A}(\mathbf{S}) is a random variable depending on both the sample S∼Dn\mathbf{S}\sim\mathcal{D}^{n} and its own random coins. The performance of A{\cal A} is quantified by its excess population risk

Note that ε\mboxrisk(A)\varepsilon_{\mbox{\footnotesize{risk}}}({\mathcal{A}}) is a random variable (due to randomness in the sample S\mathbf{S} and any possible internal randomness of the algorithm). Our guarantees on the excess population risk will be expressed in terms of upper bounds on this quantity that hold with high probability over the randomness of both S\mathbf{S} and the random coins of the algorithm.

Empirical risk minimization (ERM) is one of the most standard approaches to stochastic convex optimization. In the ERM problem, we are given a sample S=(z1,…,zn)\mathbf{S}=(\mathbf{z}_{1},\ldots,\mathbf{z}_{n}), and the goal is to find

One way to bound the excess population risk is to solve the ERM problem, and appeal to uniform convergence; however, uniform convergence rates in this case are dimension-dependent, Ω(d/n)\Omega(\sqrt{d/n}) (Feldman,, 2016).

Risk decomposition: Guaranteeing low excess population risk for a general algorithm is a nontrivial task. A common way to bound it is by decomposing it into generalization, optimization and approximation error:

For any θ∈(0,1),\theta\in(0,1), with probability at least 1−θ1-\theta, the approximation error is bounded as

Finally, note that by definition of x∗(S)x^{\ast}(\mathbf{S}), we have FS(x∗(S))−FS(x∗)≤0F_{\mathbf{S}}(x^{\ast}(\mathbf{S}))-F_{\mathbf{S}}(x^{\ast})\leq 0. Combining this with the above bound completes the proof.

We say that two datasets S, S′S,~{}S^{\prime} are neighboring, denoted S≃S′S\simeq S^{\prime}, if they only differ on a single entry; i.e., there exists i∈[n]i\in[n] s.t. for all k≠ik\neq i, zk=zk′z_{k}=z_{k}^{\prime}.

Uniform argument stability (UAS): Given an algorithm A{\cal A} and datasets S≃S′S\simeq S^{\prime}, we define the uniform argument stability (UAS) random variable as

The randomness here is due to any possible internal randomness of A{\cal A}. For any LL-Lipschitz function ff, we have that f(A(S),z)−f(A(S′),z)≤L δA(S,S′).f\left({\mathcal{A}}(S),z\right)-f\left({\mathcal{A}}(S^{\prime}),z\right)\leq L\,\delta_{\cal A}(S,S^{\prime}). Hence, upper bounds on UAS can be easily transformed into upper bounds on uniform stability.

In this work, we will consider two types of bounds on UAS.

In Section 3, we give upper bounds on UAS for three variants of the (stochastic) gradient descent algorithm, namely, (i) full-batch gradient descent, (ii) sampling-with-replacement stochastic gradient descent, and (iii) fixed-permutation stochastic gradient descent. Variant (i) is deterministic (and hence UAS is a deterministic quantity). For variant (ii), for any pair of neighboring datasets S,S′S,S^{\prime}, we give an upper bound on the UAS random variable that holds with high probability over the algorithm’s internal randomness (the sampling with replacement). For variant (iii), we give an upper bound on UAS that holds for an arbitrary choice of permutation; in particular, for any random permutation our upper bound on the UAS random variable that holds with probability 1.

Let A:Zn→X{\mathcal{A}}:\mathcal{Z}^{n}\rightarrow\mathcal{X} be a randomized algorithm. For any pair of neighboring datasets S,S′S,S^{\prime}, suppose that the UAS random variable of A{\mathcal{A}} satisfies:

Then there is a constant cc such that for any distribution D\mathcal{D} over Z\mathcal{Z} and any θ∈(0,1)\theta\in(0,1), we have

2 Expectation guarantees on UAS

In Section 4, we give lower bounds on this quantity for the two variants of the stochastic subgradient method, together with a deterministic lower bound for the full-batch variant.

Upper Bounds on Uniform Argument Stability

We begin by stating a key lemma that encompasses the UAS bound analysis of multiple variants of (S)GD. In particular, all of our UAS upper bounds are obtained by almost a direct application of this lemma. In the lemma we consider two gradient descent trajectories associated to different sequences of objective functions. The degree of concordance of the two sequences, quantified by the distance between the subgradients at the current iterate, controls the deviation between the trajectories. We note that this distance condition is satisfied for all (S)GD variants we study in this work.

Let (xt)t∈[T](x^{t})_{t\in[T]} and (yt)t∈[T](y^{t})_{t\in[T]}, with x1=y1x^{1}=y^{1}, be online gradient descent trajectories for convex LL-Lipschitz objectives (ft)t∈[T−1](f_{t})_{t\in[T-1]} and (ft′)t∈[T−1](f_{t}^{\prime})_{t\in[T-1]}, respectively; i.e.,

for all t∈[T−1]t\in[T-1]. Suppose for every t∈[T−1]t\in[T-1], ∥∇ft(xt)−∇ft′(xt)∥≤at\|\nabla f_{t}(x^{t})-\nabla f_{t}^{\prime}(x^{t})\|\leq a_{t}, for scalars 0≤at≤2L0\leq a_{t}\leq 2L. Then, if t0=inf⁡{t:ft≠ft′},t_{0}=\inf\{t:f_{t}\neq f_{t}^{\prime}\},

Let δt=∥xt−yt∥\delta_{t}=\|x^{t}-y^{t}\|. By definition of t0t_{0} it is clear that δ1=…=δt0=0\delta_{1}=\ldots=\delta_{t_{0}}=0. For t=t0+1t=t_{0}+1, we have that δt0+1=∥ηt0(∇ft0(xt0)−∇ft0′(yt0)∥≤2Lηt0\delta_{t_{0}+1}=\|\eta_{t_{0}}(\nabla f_{t_{0}}(x^{t_{0}})-\nabla f_{t_{0}}^{\prime}(y^{t_{0}})\|\leq 2L\eta_{t_{0}}.

Now, we derive a recurrence for (δt)t∈[T](\delta_{t})_{t\in[T]}:

where at the last step we use the monotonicity of the subgradient. Note that

Now we prove the following bound by induction (notice this claim proves the result):

Indeed, the claim is clearly true for t=t0t=t_{0}. For the inductive step, we assume it holds for some t∈[T−1]t\in[T-1]. To prove the result we consider two cases: first, when δt+1≤max⁡s∈[t]δs\delta_{t+1}\leq\max_{s\in[t]}\delta_{s}, by induction hypothesis we have

In the other case, δt+1>max⁡s∈[t]δs\delta_{t+1}>\max_{s\in[t]}\delta_{s}, we use (4)

Taking square root at this inequality, and using the subadditivity of the square root, we obtain the inductive step, and therefore the result. ∎

2 Upper Bounds for the Full Batch GD

As a direct corollary of Lemma 3.1, we derive the following upper bound on UAS for the batch gradient descent algorithm.

Let X⊆B(0,R)\mathcal{X}\subseteq\mathcal{B}(0,R) and F=FX0(L)\mathcal{F}=\mathcal{F}_{\mathcal{X}}^{0}(L). The full-batch gradient descent (Algorithm 1) has uniform argument stability

The bound of 2R2R is obtained directly from the diameter bound on X\mathcal{X}. Therefore, we focus exclusively on the second term. Let S≃S′S\simeq S^{\prime} be arbitrary neighboring datasets, x1=y1x^{1}=y^{1}, and consider the trajectories (xt)t,(yt)t(x^{t})_{t},(y^{t})_{t} associated with the batch GD method on datasets SS and S′S^{\prime}, respectively. We use Lemma 3.1 with ft=FSf_{t}=F_{S} and ft′=FS′f_{t}^{\prime}=F_{S^{\prime}}, for all t∈[T−1]t\in[T-1]. Notice that

since S≃S′S\simeq S^{\prime}; in particular, ∥∇ft(xt)−∇ft′(xt)∥≤at\|\nabla f_{t}(x^{t})-\nabla f_{t}^{\prime}(x^{t})\|\leq a_{t}, with at=2L/na_{t}=2L/n. We conclude by Lemma 3.1 that for all t∈[T]t\in[T]

Hence, the stability bound holds for all the iterates, and thus for x‾T\overline{x}^{T} by the triangle inequality. ∎

3 Upper Bounds for SGD

Next, we state and prove upper bounds on UAS for two variants of stochastic gradient descent: sampling-with-replacement SGD (Section 3.3.1) and fixed-permutation SGD (Section 3.3.2). Here, we give strong upper bounds that hold with high probability (for sampling-with-replacement SGD) and with probability 1 (for fixed-permutation SGD). In Appendix F, we derive tighter upper bounds for these two variants of SGD in the case where the number of iterations T<T< the number of samples in the data set nn; however, the bounds derived in this case hold only in expectation.

Next, we study the uniform argument stability of the sampling-with-replacement stochastic gradient descent (Algorithm 2). This algorithm has the benefit that each iteration is extremely cheap compared to Algorithm 1. Despite these savings, we will show that same bound on UAS holds with high probability.

We now state and prove our upper bound for sampling-with-replacement SGD.

Let X⊆B(0,R)\mathcal{X}\subseteq\mathcal{B}(0,R) and F=FX0(L)\mathcal{F}=\mathcal{F}_{\mathcal{X}}^{0}(L). The uniform argument stability of the sampling-with-replacement SGD (Algorithm 2) satisfies:

Moreover, if ηt=η>0 ∀t\eta_{t}=\eta>0~{}\forall t then, for any pair (S,S′)(S,S^{\prime}) of neighboring datasets, with probability at least 1−exp⁡(−n/2)1-\exp\left(-n/2\right) (over the algorithm’s internal randomness), the UAS random variable is bounded as

The bound of 2R2R trivially follows from the diameter bound on X\mathcal{X}. We thus focus on the second term of the bound. Let S≃S′S\simeq S^{\prime} be arbitrary neighboring datasets, x0=y0x^{0}=y^{0}, and consider the trajectories (xt)t∈[T],(yt)t∈[T](x^{t})_{t\in[T]},(y^{t})_{t\in[T]} associated with the sampled-with-replacement stochastic subgradient method on datasets SS and S′S^{\prime}, respectively. We use Lemma 3.1 with ft(⋅)=f(⋅,zit)f_{t}(\cdot)=f(\cdot,\mathbf{z}_{\mathbf{i}_{t}}) and ft′(⋅)=f(⋅,zit′)f_{t}^{\prime}(\cdot)=f(\cdot,\mathbf{z}_{{\mathbf{i}_{t}}^{\prime}}). Let us define rt≜1{zit≠zit′}\mathbf{r}_{t}\triangleq\mathbf{1}_{\{\mathbf{z}_{\mathbf{i}_{t}}\neq\mathbf{z}_{\mathbf{i}_{t}}^{\prime}\}}. Note that at every step tt, rt=1\mathbf{r}_{t}=1 with probability 1−1/n1-1/n, and rt=0\mathbf{r}_{t}=0 otherwise. Moreover, note that {rt: t∈[T]}\{\mathbf{r}_{t}:~{}t\in[T]\} is an independent sequence of Bernoulli random variables. Finally, note that ∥∇ft(xt)−∇ft′(xt)∥≤2Lrt\|\nabla f_{t}(x^{t})-\nabla f_{t}^{\prime}(x^{t})\|\leq 2L\mathbf{r}_{t}.

Hence, by Lemma 3.1, for any realization of the trajectories of the SGD method, we have

where ΔT≜2L∑s=1T−1ηs2+4L∑s=1T−1rsηs\Delta_{T}\triangleq 2L\sqrt{\sum_{s=1}^{T-1}\eta_{s}^{2}}+4L\sum_{s=1}^{T-1}\mathbf{r}_{s}\eta_{s}. Taking expectation of (5), we have

This establishes the upper bound on UAS but only in expectation. Now, we proceed to prove the high-probability bound. Here, we assume that the step size is fixed; that is, ηt=η>0\eta_{t}=\eta>0 for all t∈[T−1]t\in[T-1]. Note that each rs,s∈[T],\mathbf{r}_{s},s\in[T], has variance 1n(1−1n)<1n\frac{1}{n}\left(1-\frac{1}{n}\right)<\frac{1}{n}. Hence, by Chernoff’s boundHere, we are applying a bound for (scaled) Bernoulli rvs where the exponent is expressed in terms of the variance., we have

Therefore, with probability at least 1−exp⁡(−n/2)1-\exp\left(-n/2\right), we have

Putting this together with (5), with probability at least 1−exp⁡(−n/2)1-\exp\left(-n/2\right), we have

Finally, by the triangle inequality, we get that with probability at least 1−exp⁡(−n/2)1-\exp\left(-n/2\right), the same stability bound holds for the average of the iterates x‾T\overline{x}^{T}, y‾T\overline{y}^{T}. ∎

3.2 Upper Bounds for the Fixed Permutation SGD

In Algorithm 3, we describe the fixed-permutation stochastic gradient descent. This algorithm works in epochs, where each epoch is a single pass on the data. The order in which data is used is the same across epochs, and is given by a permutation π\pi. The algorithm can be alternatively described without the epoch loop simply by

We show that the same UAS bound of batch gradient descent and sampling-with-replacement SGD holds for the fixed-permutation SGD. We also observe that a slightly tighter bound can be achieved if we consider the expectation guarantee on UAS when π\bm{\pi} is chosen uniformly at random. We leave these details to Theorem F.2 in the Appendix.

In the next result, we assume that the sequence of step sizes (ηt)t∈[T]\left(\eta_{t}\right)_{t\in[T]} is non-increasing, which is indeed the case for almost all known variants of SGD.

Let X⊆B(0,R)\mathcal{X}\subseteq\mathcal{B}(0,R), F=FX0(L)\mathcal{F}=\mathcal{F}_{\mathcal{X}}^{0}(L), and π\pi be any permutation over [n][n]. Suppose the step sizes (ηt)t∈[T]\left(\eta_{t}\right)_{t\in[T]} form a non-increasing sequence. Then the uniform argument stability of the fixed-permutation SGD (Algorithm 3) is bounded as

Again, the bound of 2R2R is trivial. Now, we show the second term of the bound. Let S≃S′S\simeq S^{\prime} be arbitrary neighboring datasets, x1=y1x^{1}=y^{1}, and consider the trajectories (xt)t∈[T],(yt)t∈[T](x^{t})_{t\in[T]},(y^{t})_{t\in[T]} associated with the fixed permutation stochastic subgradient method on datasets SS and S′S^{\prime}, respectively. Since the datasets S≃S′S\simeq S^{\prime} are arbitrary, we may assume without loss of generality that π\pi is the identity, whereas the perturbed coordinate i=i\mathbf{i}=i is arbitrary. We use Lemma 3.1 with ft(⋅)=f(⋅,z(t \mboxmodn))f_{t}(\cdot)=f\left(\cdot,\mathbf{z}_{\left(t~{}\mbox{\footnotesize mod }n\right)}\right) and ft′(⋅)=f(⋅,z(t \mboxmodn)′)f_{t}^{\prime}(\cdot)=f\left(\cdot,\mathbf{z}^{\prime}_{\left(t~{}\mbox{\footnotesize mod }n\right)}\right). It is easy to see then that ∥∇ft(xt)−∇ft′(xt)∥≤at\|\nabla f_{t}(x^{t})-\nabla f_{t}^{\prime}(x^{t})\|\leq a_{t}, with at=2L⋅1{(t \mboxmodn)=i}a_{t}=2L\cdot\mathbf{1}_{\{(t~{}\mbox{\footnotesize mod }n)=i\}}, where 1{condition}\mathbf{1}_{\{\mathsf{condition}\}} is the indicator of condition\mathsf{condition}. Hence, by Lemma 3.1, we have

where at the last step we used the fact that (ηt)t∈[T](\eta_{t})_{t\in[T]} is non-increasing; namely, for any r≥1r\geq 1

Since the bound holds for all the iterates, using triangle inequality, it holds for the output x‾K\overline{x}^{K} averaged over the iterates from the T/nT/n epochs. ∎

4 Discussion of the upper bounds: examples of specific instantiations

The upper bounds on stability from this section all behave very similarly. Let us explore the consequences of the obtained rates in terms of generalization bounds for different choices of the step size sequence. As a case study, we will consider excess risk bounds for the full-batch subgradient method (Algorithm 1), but similar conclusions hold for all the variants that we studied. We emphasize that prior to this work, no dimension-independent bounds on the excess risk were known of this method (specifically, for nonsmooth losses and without explicit regularization).

To bound the excess risk, we will use the risk decomposition, eqn. (3). For simplicity, we will only be studying excess risk bounds in expectation (in Section 6 we consider stronger, high probability, bounds). In this case, the stability implies generalization result (Theorem 2.2) simplifies to Bousquet and Elisseeff, (2002); Hardt et al., (2016)

Finally, the approximation error (Lemma 2.1) simplifies as well: it is upper bounded by 0 in expectation.

Fixed stepsize: Let ηt≡η>0\eta_{t}\equiv\eta>0. By Thm. 3.2, UAS is bounded by 4LTη+4LTηn4L\sqrt{T}\eta+\frac{4LT\eta}{n}. On the other hand, the standard analysis of subgradient descent guarantees that ε\mboxopt(AGD)≤R22ηT+ηL22\varepsilon_{\mbox{\footnotesize opt}}({\mathcal{A}}_{\sf GD})\leq\frac{R^{2}}{2\eta T}+\frac{\eta L^{2}}{2}. Therefore, by the expected risk decomposition (3)

If we consider the standard method choice, η=R/[Ln]\eta=R/[L\sqrt{n}] and T=nT=n, the bound above is at least 4LR4LR (due to the first term). Consequently, upper bounds obtained from this approach are vacuous.

Varying stepsize: For a general sequence of stepsizes the optimization guarantees of Algorithm 1 are the following

In fact, we can show that any choice of step sizes that makes the quantity above O(LR/n)O(LR/\sqrt{n}) must necessarily have T=Ω(n2)T=\Omega(n^{2}). Indeed, notice that in such case

The high iteration complexity required to obtain optimal bounds motivates studying whether it is possible to improve our uniform argument stability bounds. We will show that, unfortunately, they are sharp up to absolute constant factors.

Lower Bounds on Uniform Argument Stability

In this section we provide matching lower bounds for the previously studied first-order methods. These lower bounds show that our analyses are tight, up to absolute constant factors.

We note that it is possible to prove a general purpose lower bound on stability by appealing to sample complexity lower bounds for stochastic convex optimization (Nemirovsky and Yudin,, 1983). This approach in the smooth convex case was first studied in (Chen et al.,, 2018); there, these lower bounds are sharp. However, in the nonsmooth case they are very far from bounds in the previous section. The idea is that for sufficiently small step size, a first-order method must incur Ω(LTη/n)\Omega(LT\eta/n) uniform stability. Details of this lower bound can be found on Appendix C. This reasoning leads to an Ω(LTη/n)\Omega(LT\eta/n) lower bound on uniform argument stability, that can be added to any other lower bound we can prove on specific algorithms that enjoy rates as of gradient descent.

Next we will prove finer lower bounds on the UAS of specific algorithms. For this, note that the objective functions we use are polyhedral, thus the subdifferential is a polytope at any point. Since the algorithm should work for any oracle, we will let the subgradients provided to be extreme points, ∇f(x,z)∈\mboxext(∂f(x,z))\nabla f(x,z)\in\mbox{ext}(\partial f(x,z)). Moreover, we can make adversarial choices of the chosen subgradient.

Let X=B(0,1)\mathcal{X}={\cal B}(0,1), F=FX0(1){\cal F}={\cal F}_{\mathcal{X}}^{0}(1) and d≥min⁡{T,1/η2}d\geq\min\{T,1/\eta^{2}\}. For the full-batch gradient descent (Alg. 1) with constant step size η>0\eta>0, there exist S≃S′S\simeq S^{\prime} such that the UAS is lower bounded as δAGD(S,S′)=Ω(min⁡{1,ηT+ηT/n}).\delta_{{\mathcal{A}}_{\sf GD}(S,S^{\prime})}=\Omega(\min\{1,\eta\sqrt{T}+\eta T/n\}).

The proof of this result is deferred to Appendix D, due to space considerations.

2 Lower Bounds for SGD Sampled with Replacement

We use a similar construction as from the previous result to prove a sharp lower bound on the uniform argument stability for stochastic gradient descent where the sampling is with replacement.

Let D≜min⁡{T,1/η2}≤dD\triangleq\min\{T,1/\eta^{2}\}\leq d, and ν>0\nu>0, K≥DK\geq\sqrt{D}. Consider Z={0,1}{\cal Z}=\{0,1\} and define

where r=(−1,…,−1,0,…,0)r=(-1,\ldots,-1,0,\ldots,0) (i.e., supported on the first DD coordinates). Let the random sequence of indices used by the algorithm: (it)t≥0∼i.i.d.\mboxUnif([n])(\mathbf{i}_{t})_{t\geq 0}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mbox{Unif}([n]). Let S=(1,0,…,0)S=(1,0,\ldots,0) and S′=(0,0,…,0)S^{\prime}=(0,0,\ldots,0) be neighboring datasets, and denote by (xt)t(x^{t})_{t} and (yt)t(y^{t})_{t} the respective stochastic gradient descent trajectories on SS and S′S^{\prime}, initialized at x1=y1=0x^{1}=y^{1}=0. It is easy to see that under S′S^{\prime}, we have yt=0y^{t}=0 for all t∈[T]t\in[T]. Now, suppose that ν<η/K\nu<\eta/K. Then, we only have xt=0x^{t}=0 for all t≤τt\leq\tau, where τ:=inf⁡{t≥1: it=1}\tau:=\inf\{t\geq 1:\,\mathbf{i}_{t}=1\}. After time τ\tau, xτ+1=−ηr/Kx^{\tau+1}=-\eta r/K, and consequently xτ+1+j=−ηk(τ+j)Kr−η∑s=1j−k(τ+j)+1es,x^{\tau+1+j}=-\frac{\eta\mathbf{k}(\tau+j)}{K}r-\eta\sum_{s=1}^{j-\mathbf{k}(\tau+j)+1}e_{s}, for all j∈[D+k(τ+j)−1]j\in[D+\mathbf{k}(\tau+j)-1], where k(t)≜∣{s∈[t]:is=1}∣\mathbf{k}(t)\triangleq|\{s\in[t]:\mathbf{i}_{s}=1\}|. Note that conditioned on any fixed value for τ,\tau, k(τ+j)≤j+1\mathbf{k}(\tau+j)\leq j+1.

where c=(1−exp⁡(−T/4))=Ω(1)c=(1-\exp(-T/4))=\Omega(1). Hence,

We choose KK sufficiently large such that ηDT/K=o(ηmin⁡{T3/2/n,T})\eta\sqrt{D}T/K=o(\eta\min\{T^{3/2}/n,\sqrt{T}\}). Hence, we have

3 Lower Bounds for the Fixed Permutation Stochastic Gradient Descent

Finally, we study fixed permutation SGD. The proof of this result is deferred to Appendix E.

Generalization Guarantees for Multi-pass SGD

One important implication of our stability bounds is that they provide non-trivial generalization error guarantees for multi-pass SGD on nonsmooth losses. Multi-pass SGD is one of the most extensively used settings of SGD in practice, where SGD is run for KK passes (epochs) over the dataset (namely, the number of iterations T=KnT=Kn). To the best of our knowledge, aside from the dimension-dependent bounds based on uniform convergence (Shalev-Shwartz et al.,, 2010), no generalization error guarantees are known for the multi-pass setting on general nonsmooth convex losses. Given our uniform stability upper bounds, we can prove the following generalization error guarantees for the multi-pass setting of sampling-with-replacement SGD. Analogous results can be obtained for fixed-permutation SGD .

Running Algorithm 2 for KK passes (i.e., for T=KnT=Kn iterations) with constant stepsize ηt=η>0\eta_{t}=\eta>0 yields the following generalization error guarantees:

and there exists c>0c>0, such that for any 0<θ<10<\theta<1, with probability ≥1−θ−exp⁡(−n/2)\geq 1-\theta-\exp(-n/2),

First, by the expectation guarantee on UAS given in Theorem 3.3 together with the fact that the losses are LL-Lipschitz, it follows that Algorithm 2 (when run for KK passes with constant stepsize η\eta) is γ\gamma-uniformly stable, where γ=4L2(ηKn+ηK)\gamma=4L^{2}\left(\eta\sqrt{Kn}+\eta K\right). Then, by (Hardt et al.,, 2016, Thm. 2.2), we have

For the high-probability bound, we combine the high-probability guarantee on UAS given in Theorem 3.3 with Theorem 2.2 to get the claimed bound. ∎

These bounds on generalization error can be used to obtain excess risk bounds using the standard risk decomposition (see (3)). In practical scenarios where one can certify small optimization error for multi-pass SGD, Thm. 5.1 can be used to readily estimate the excess risk. In Section 6.2 we provide worst-case analysis showing that multi-pass SGD is guaranteed to attain the optimal excess risk of ≈LR/n\approx LR/\sqrt{n} within nn passes (with appropriately chosen constant stepsize).

Implications of Our Stability Bounds

Using our UAS upper bounds, we show that a simple variant of noisy SGD (Bassily et al.,, 2014), that requires only n2n^{2} gradient computations, yields the optimal excess population risk for DP-SCO. In terms of running time, this is a small improvement over the algorithm of Feldman et al., (2020) for the nonsmooth case, which requires O(n2log⁡1/β)O(n^{2}\sqrt{\log 1/\beta}) gradient computations. More importantly, our algorithm is substantially simpler. For comparison, the algorithm in (Feldman et al.,, 2020) is based on a multi-phase SGD, where in each phase a separate regularized ERM problem is solved. To ensure privacy, the output of each phase is perturbed with an appropriately chosen amount of noise before being used as the initial point for the next phase.

The description of the algorithm is given in Algorithm 4.

Algorithm 4 is (α,β)(\alpha,\beta)-differentially private.

The proof of the theorem follows the same lines of (Bassily et al.,, 2014, Theorem 2.1), but we replace their privacy analysis of the Gaussian mechanism with the tighter Moments Accountant method of Abadi et al., (2016). analysis of Abadi et al., (2016).

In Algorithm 4, let \eta=R/\Big{(}L\cdot n\cdot\max\big{(}\sqrt{n},~{}\frac{\sqrt{d\,\log(1/\beta)}}{\alpha}\big{)}\Big{)}. Then, for any θ∈(6exp⁡(−n/2),1)\theta\in(6\exp(-n/2),1), with probability at least 1−θ1-\theta over the randomness in both the sample and the algorithm, we have

Fix any confidence parameter θ≥6exp⁡(−n/2)\theta\geq 6\exp(-n/2). First, for any data set S∈ZnS\in\mathcal{Z}^{n} and any step size η>0,\eta>0, by Lemma H.1 in Appendix H, we have the following high-probability guarantee on the training error of ANSGD{\mathcal{A}}_{\sf NSGD}:

With probability at least 1−θ/3,1-\theta/3, we have

where the probability is over the sampling in step 4 and the independent Gaussian noise vectors G1,…,Gn2\mathbf{G}_{1},\ldots,\mathbf{G}_{n^{2}}. Given the setting of η\eta in the theorem, we get

Next, it is not hard to show that ANSGD{\mathcal{A}}_{\sf NSGD} attains the same UAS bound as ArSGD{\mathcal{A}}_{\sf rSGD} (Theorem 3.3). Indeed, the only difference is the noise addition in gradient step; however, this does not impact the stability analysis. This is because the sequence of noise vectors {G1,…,Gn2}\{\mathbf{G}_{1},\ldots,\mathbf{G}_{n^{2}}\} is the same for the trajectories corresponding to the pair S, S′S,~{}S^{\prime} of neighboring datasets. Hence, the argument basically follows the same lines of the proof of Theorem 3.3 since the noise terms cancel out. Thus, we conclude that for any pair S≃S′S\simeq S^{\prime} of neighboring datasets, with probability at least 1−exp⁡(n/2)≥1−θ/61-\exp(n/2)\geq 1-\theta/6 (over the randomness of ANSGD{\mathcal{A}}_{\sf NSGD}), the uniform argument stability of ANSGD{\mathcal{A}}_{\sf NSGD} is bounded as: δANSGD≤4Lη(T+Tn),\delta_{{\mathcal{A}}_{\sf NSGD}}\leq 4L\eta\left(\sqrt{T}+\frac{T}{n}\right), where T=n2T=n^{2}. Given the setting of η\eta in the theorem, this bound reduces to 8R/\max\big{(}\sqrt{n},~{}\frac{\sqrt{d\,\log(1/\beta)}}{\alpha}\big{)}.

Hence, by Theorem 2.2, with probability at least 1−θ/31-\theta/3 (over the randomness in both the i.i.d. dataset SS and the algorithm), the generalization error of ANSGD{\mathcal{A}}_{\sf NSGD} is bounded as

where cc in the first bound is a universal constant.

Now, using (7), (8), and Lemma 2.1, we finally conlcude that with probability at least 1−θ1-\theta (over randomness in the sample SS and the internal randomness of ANSGD{\mathcal{A}}_{\sf NSGD}), the excess population risk of ANSGD{\mathcal{A}}_{\sf NSGD} is bounded as

Using the expectation guarantee on UAS given in Theorem 3.3 and following similar steps of the analysis above, we can also show that the expected excess population risk of ANSGD{\mathcal{A}}_{\sf NSGD} is bounded as:

2 Nonsmooth Stochastic Convex Optimization with Multi-pass SGD

Another application of our results concerns obtaining optimal excess risk for stochastic nonsmooth convex optimization via multi-pass SGD. It is known that one-pass SGD is guaranteed to have optimal excess risk, which can be shown via martingale arguments that trace back to the stochastic approximation literature (Robbins and Monro,, 1951; Kiefer and Wolfowitz,, 1952).

Using our UAS bound, we show that Algorithms 2 and 3 can recover nearly-optimal high-probability excess risk bounds by making nn passes over the data. Analogous bounds hold for Algorithm 1, however these are less interesting from a computational efficiency perspective.

In short, we prove the following excess population risk bounds, and their analyses are deferred to Appendix G,

Discussion and Open Problems

In this work we provide sharp upper and lower bounds on uniform argument stability for the (stochastic) subgradient method in stochastic nonsmooth convex optimization. Our lower bounds show inherent limitations of stability bounds compared to the smooth convex case, however we can still derive optimal population risk bounds by reducing the step size and running the algorithms for longer number of iterations. We provide applications of this idea for differentially-private noisy SGD, and for two versions of SGD (sampling-with-replacement and fixed-permutation SGD).

The first open problem regards lower bounds that are robust to general forms of algorithmic randomization. Unfortunately, the methods presented here are not robust in this respect, since random initialization would prevent the trajectories reaching the region of highly nonsmooth behavior of the objective (or doing it in such a way that it does not increase UAS). One may try to strengthen the lower bound by using a random rotation of the objective; however, this leads to an uninformative lower bound. Finding distributional constructions for lower bounds against randomization is a very interesting future direction.

Our privacy application provides optimal risk for an algorithm that runs for n2n^{2} steps, which is impractical for large datasets. Other algorithms, e.g. in (Feldman et al.,, 2020), run into similar limitations. Proving that quadratic running time is necessary for general nonsmooth DP-SCO is a very interesting open problem that can be formalized in terms of the oracle complexity of stochastic convex optimization (Nemirovsky and Yudin,, 1983) under stability and/or privacy constraints.

Part of this work was done while the authors were visiting the Simons Institute for the Theory of Computing during the “Data Privacy: Foundations and Applications” program. RB’s research is supported by NSF Awards AF-1908281, SHF-1907715, Google Faculty Research Award, and OSU faculty start-up support. Work by CG was partially funded by the Millennium Science Initiative of the Ministry of Economy, Development, and Tourism, grant “Millennium Nucleus Center for the Discovery of Structures in Complex Data.” CG would like to thank Nicolas Flammarion and Juan Peypouquet for extremely valuable discussions at early stages of this work.

References

Appendix A Proof of Theorem 3.2

The bound of 2R2R is obtained directly from the diameter bound on X\mathcal{X}. Therefore, we focus exclusively on the second term. Let S≃S′S\simeq S^{\prime} be arbitrary neighboring datasets, x1=y1x^{1}=y^{1}, and consider the trajectories (xt)t,(yt)t(x^{t})_{t},(y^{t})_{t} associated with the batch GD method on datasets SS and S′S^{\prime}, respectively. We use Lemma 3.1 with ft=FSf_{t}=F_{S} and ft′=FS′f_{t}^{\prime}=F_{S^{\prime}}, for all t∈[T−1]t\in[T-1]. Notice that

since S≃S′S\simeq S^{\prime}; in particular, ∥∇ft(xt)−∇ft′(xt)∥≤at\|\nabla f_{t}(x^{t})-\nabla f_{t}^{\prime}(x^{t})\|\leq a_{t}, with at=2L/na_{t}=2L/n. We conclude by Lemma 3.1 that for all t∈[T]t\in[T]

Hence, the stability bound holds for all the iterates, and thus for x‾T\overline{x}^{T} by the triangle inequality. ∎

Appendix B Proof of Theorem 3.4

The stability bound of 2R2R is implied directly by the diameter of the feasible set. Let S≃S′S\simeq S^{\prime}, and let (xt)t∈[T],(yt)t∈[T](x^{t})_{t\in[T]},(y^{t})_{t\in[T]} be the trajectories of Algorithm 3 on SS and S′S^{\prime}, respectively, with x1=y1x^{1}=y^{1}.

Notice that since π\bm{\pi} is a random permutation, we may assume w.l.o.g. that π\bm{\pi} is the identity, whereas the perturbed coordinate between S,S′S,S^{\prime} is i∼\mboxUnif([n])\mathbf{i}\sim\mbox{Unif}([n]). The rest of the proof is a stability analysis conditioned on π\bm{\pi} (which fixes all the randomness of the algorithm), but from the observation above this is the same as conditioning on the random perturbed coordinate i\mathbf{i}.

Let T≤nT\leq n, and δt=∥xt−yt∥\delta_{t}=\|x^{t}-y^{t}\| so that δ1=0\delta_{1}=0. Conditioned on i=i\mathbf{i}=i, we have that for all t≤Tt\leq T,

Indeed, for all t≤it\leq i, δt=0\delta_{t}=0. For t=it=i, we have

where we used xi=yix^{i}=y^{i}, and that both gradients are bounded in norm by LL. Finally, when t>it>i, we have zt=zt′z_{t}=z_{t}^{\prime}, and therefore we can leverage the monotonicity of the subgradients

where the first inequality holds by the Jensen inequality, the second inequality comes from the bound on the conditional expectation, the third inequality from the non-increasing stepsize assumption, and the fourth inequality is from Cauchy-Schwarz. Since averaging can only improve stability, we conclude the result. ∎

Appendix C General Lower Bound on Stability of SGD

Let A{\cal A} be a γ\gamma-uniformly stable stochastic convex optimization algorithm with γ=s(T)/n\gamma=s(T)/n, where s(T)s(T) is increasing and lim⁡T→+∞s(T)=+∞\lim_{T\rightarrow+\infty}s(T)=+\infty. By the lower bound on the optimal risk of nonsmooth convex optimization, ε\mboxrisk≥LRC1n\varepsilon_{\mbox{\footnotesize risk}}\geq\frac{LR}{C_{1}\sqrt{n}}, where C1>0C_{1}>0 is a universal constant (Nemirovsky and Yudin,, 1983). This, combined with the risk decomposition (3), implies that

By our assumption on s(T)s(T), for TT sufficiently large, there always exists nn such that

which leads to ε\mboxopt≥(LR)2C2s(T)\varepsilon_{\mbox{\footnotesize opt}}\geq\frac{(LR)^{2}}{C_{2}s(T)}, where C2>0C_{2}>0 is a universal constant.

If algorithm A{\cal A} is based on TT subgradient iterations with constant step size η>0\eta>0 (these could be either stochastic, batch or minibatch), by standard analysis, the optimization guarantee of such algorithm is ε\mboxopt≤12(R2ηT+ηL2)\varepsilon_{\mbox{\footnotesize opt}}\leq\frac{1}{2}(\frac{R^{2}}{\eta T}+\eta L^{2}). Both bounds in combination give

If we further assume that η≤(R/L)/T\eta\leq(R/L)/\sqrt{T} (notice η=(R/L)/T\eta=(R/L)/\sqrt{T} minimizes the optimization error), then s(T)≥L2ηT/C2s(T)\geq L^{2}\eta T/C_{2}. We also emphasize all the choices of step size that we will make to control generalization error will lie in this range.

Appendix D Proof of Theorem 4.1

Let D≜min⁡{T,1/η2}≤dD\triangleq\min\{T,1/\eta^{2}\}\leq d, and ν,K>0\nu,K>0. We consider Z={0,1}\mathcal{Z}=\{0,1\}, and the objective function

where r=(−1,…,−1,0,…,0)r=(-1,\ldots,-1,0,\ldots,0) (i.e., supported on the first DD coordinates). Notice that for normalization purposes, we need K≥DK\geq\sqrt{D}; furthermore, we will choose KK sufficiently large such that TD/[nK]=o(1)T\sqrt{D}/[nK]=o(1). Consider the data sets S≃S′S\simeq S^{\prime}, leading to the empirical objectives:

Let (xt)t∈[T](x^{t})_{t\in[T]} and (yt)t∈[T](y^{t})_{t\in[T]} be the trajectories of the algorithm over datasets SS and S′S^{\prime}, respectively, initialized from x1=y1=0x^{1}=y^{1}=0. Clearly, yt=0y^{t}=0 for all tt. Now x2=−ηnK rx^{2}=-\frac{\eta}{nK}\,r; choosing ν<η/(nK)\nu<\eta/(nK), we have ∇f(x2,z)=−ηnK r+n−1ne1\nabla f(x^{2},z)=-\frac{\eta}{nK}\,r+\frac{n-1}{n}e_{1}, and hence x3=−2ηnK r−ηn−1ne1.x^{3}=-\frac{2\eta}{nK}\,r-\eta\frac{n-1}{n}e_{1}. Sequentially, the method will perform cumulative subgradient steps on e2,e3…,eDe_{2},e_{3}\ldots,e_{D}. In particular, for any t∈[D+1],t\in[D+1], we have xt+1=−tηnK r−ηn−1n∑s=1t−1es.x^{t+1}=-t\frac{\eta}{nK}\,r-\eta\frac{n-1}{n}\sum_{s=1}^{t-1}e_{s}.

By orthogonality of the subgradients and given our choice of KK, we conclude that

and further subgradient steps t=D+1,…,Tt=D+1,\ldots,T are only given by the linear term, r/[nK]r/[nK], which are negligible perturbations.

We finish by arguing that averaging does not help. First, in the case D=TD=T:

Finally, the additional term Ω(ηT/n)\Omega(\eta T/n) in the lower bound is obtained by Appendix C. ∎

Appendix E Proof of Theorem 4.3

We consider the same function class of Thm. 4.1, and neighbor datasets S′=(0,0,…,0)S^{\prime}=(0,0,\ldots,0), S=(1,0,…,0)S=(1,0,\ldots,0). We will assume in what follows that D=min⁡{T,1/η2}D=\min\{T,1/\eta^{2}\}, KK is sufficiently large and ν<η∥r∥/K\nu<\eta\|r\|/K. Let (xt)t∈[T](x^{t})_{t\in[T]} and (yt)t∈[T](y^{t})_{t\in[T]} be the trajectories of Algorithm 3 over datasets S,S′S,S^{\prime} respectively, both initialized at x1=y1=0x^{1}=y^{1}=0. Let now τ=π−1(1)∼\mboxUnif[n]\tau=\bm{\pi}^{-1}(1)\sim\mbox{Unif}[n]. Arguing as in Thm. 4.1, we have that yt=0y^{t}=0 for all tt, whereas

Later iterations will satisfy ∥xt∥=1−o(1)\|x^{t}\|=1-o(1) if D=1/η2D=1/\eta^{2} (and otherwise the algorithm stops earlier). Therefore, for all t∈[T]t\in[T],

Notice that we used above that KK is such that TD/nK=o(1)T\sqrt{D}/nK=o(1). Analogously as in Thm. 4.1, we can obtain the same conclusion for x‾T\overline{x}^{T}. The lower bound of ηT/n\eta T/n can be added by Appendix C, so the result follows. ∎

Appendix F Upper bounds on UAS of SGD when T≤n𝑇𝑛T\leq n

Let X⊆B(0,R)\mathcal{X}\subseteq\mathcal{B}(0,R) and F=FX0(L)\mathcal{F}=\mathcal{F}_{\mathcal{X}}^{0}(L). Suppose T≤nT\leq n. The UAS of sampling-with-replacement stochastic gradient descent (Algorithm 2) satisfies uniform argument stability

The bound of 2R2R is obtained directly from the diameter bound on X\mathcal{X}. Therefore, we focus exclusively on the second term. Let S≃S′S\simeq S^{\prime}, and let k∈[n]k\in[n] be the entry where both datasets differ. Let (xt)t∈[T],(yt)t∈[T](x^{t})_{t\in[T]},(y^{t})_{t\in[T]} be the trajectories of Algorithm 2 on SS and S′S^{\prime}, respectively, with x1=y1x^{1}=y^{1}.

Let BtB_{t} denote the event that ij=k\mathbf{i}_{j}=k for some j≤tj\leq t; that is, BtB_{t} is the event that the index kk is sampled at least once in the first tt iterations. We note that

Now, conditioned on the past sampled coordinates i1,…,it−1\mathbf{i}_{1},\ldots,\mathbf{i}_{t-1}, we have

where the last inequality is obtained from convexity and Lipschitzness of the objective. Now, squaring we get

From this formula we derive bounds for the two conditional expectations:

where (12) holds by independence of it\mathbf{i}_{t} and Bt−1B_{t-1}, and in (13) we used that δt=0\delta_{t}=0 conditioned on Bt−1‾\overline{B_{t-1}}.

With this last bound we can proceed inductively to show that

proving the inductive step. Finally, putting this together with (10) completes the proof. ∎

Notice that the UAS of Algorithms 2 and 3 are of the same order. The proof is in Appendix B.

Appendix G Generalization in Stochastic Optimization with Data Reuse

Let X⊆B(0,R)\mathcal{X}\subseteq\mathcal{B}(0,R) and F=FX0(L)\mathcal{F}=\mathcal{F}_{\mathcal{X}}^{0}(L). The sampling with replacement stochastic gradient descent (Algorithm 2) with T=n2T=n^{2} iterations and η=R L n3/2\eta=\frac{R}{\,L\,n^{3/2}} satisfies for any 12exp⁡{−n2/32}<θ<1.12\exp\{-n^{2}/32\}<\theta<1.

It should be noted that, similarly to Remark 6.3, if we are only interested in expectation risk bounds, one can shave off the polylogarithmic factor above, which is optimal for the expected excess risk.

Let S∼Dn\mathbf{S}\sim{\cal D}^{n} an i.i.d. random sample for the stochastic convex program, and apply on these data the algorithm ArSGD{\mathcal{A}}_{\sf rSGD} for constant step size η>0\eta>0 and TT iterations.

We consider θ>0\theta>0 such that θ>12exp⁡{−T/32}\theta>12\exp\{-T/32\}. Notice that the sampling-with-replacement stochastic gradient is a bounded first-order stochastic oracle for the empirical objective. It is direct to verify that the assumptions of Lemma H.1 are satisfied with σ=0\sigma=0. Hence, by Lemma H.1, we have that, with probability at least 1−θ/31-\theta/3

On the other hand, Theorem 3.3 together with Theorem 2.2, guarantees that with probability at least 1−θ/31-\theta/3, we have

Finally, Lemma 2.1 ensures that with probability 1−θ/31-\theta/3

By the union bound and the excess risk decomposition (3), we have that, with probability 1−θ1-\theta,

where only at the last step we replaced by the choice of step size and number of iterations from the statement.

G.0.2 Risk Bounds for Fixed-Permutation Stochastic Gradient Descent

As a final application we provide a population risk bound based on the UAS of Algorithm 3. Similarly as in the case of sampling-with-replacement SGD, we need an optimization error analysis, which for completeness is provided in Appendix I, and it is based on the analysis of the incremental subgradient method (Nedic and Bertsekas,, 2001).

Interestingly, the combination of the incremental method analysis for arbitrary permutation (Nedic and Bertsekas,, 2001) and our novel stability bounds that also work for arbitrary permutation, guarantees generalization bounds for fixed permutation SGD without the need of reshuffling, or even any form of randomization. We believe this could be of independent interest.

Algorithm 3 with constant step size ηk≡η=R/[LnK]\eta_{k}\equiv\eta=R/[Ln\sqrt{K}] and K=nK=n epochs is such that for every 0<θ<10<\theta<1,

Similarly to the previous result, we can remove the polylogarithmic factor if we are only interested in expected excess risk guarantees.

by our choice of K,ηK,\eta. On the other hand, Theorem 3.4 guarantees the algorithm is δ\delta-UAS with probability 1, where δ=O(R/n)\delta=O(R/\sqrt{n}). Therefore, by Theorem 2.2, we have that w.p. 1−θ/21-\theta/2

Finally, Lemma 2.1 ensures that with probability 1−θ/21-\theta/2

By the union bound and the excess risk decomposition (3), we have that, with probability at least 1−θ1-\theta,

Appendix H High-probability Bound on Optimization Error of SGD with Noisy Gradient Oracle

It is known that standard online-to-batch conversion technique can be used to provide high-probability bound on the optimization error (i.e., the excess empirical risk) of stochastic gradient descent. For the sake of completeness and to make the paper more self-contained, we re-expose this technique here for stochastic gradient descent with noisy gradient oracle. This is done in the following lemma, which is used in the proofs of our results in Section 6.

Let S=(z1,…,zn)∈ZnS=(z_{1},\ldots,z_{n})\in\mathcal{Z}^{n} be a dataset. Let FS(x)=1n∑i∈[n]f(x,zi)F_{S}(x)=\frac{1}{n}\sum_{i\in[n]}f(x,z_{i}) be the empirical risk associated with SS, where for every z∈Z,z\in\mathcal{Z}, f(⋅,z)f(\cdot,z) is convex, LL-Lipschitz function over X⊆B(0,R)\mathcal{X}\subseteq\mathcal{B}(0,R) for some L,R>0L,R>0. Consider the stochastic (sub)gradient method:

Sub-Gaussian gradient noise: There is σ2≥0\sigma^{2}\geq 0 such that for every x∈X,z∈Zx\in\mathcal{X},z\in\mathcal{Z}, g(x,z)−∇f(x,z)\mathbf{g}(x,z)-\nabla f(x,z) is σ2\sigma^{2}-sub-Gaussian random vector; that is, for every x∈X,z∈Z,x\in\mathcal{X},z\in\mathcal{Z}, ⟨g(x,z)−∇f(x,z), u⟩\langle\mathbf{g}(x,z)-\nabla f(x,z),~{}u\rangle is σ2\sigma^{2}-sub-Gaussian random variable ∀u∈B(0,1)\forall u\in\mathcal{B}(0,1).

Independence of the gradient noise across iterations: conditioned on any fixed realization of (ξt)t∈[T](\xi_{t})_{t\in[T]} the sequence of random maps g(⋅,ξ1),…,g(⋅,ξT)\mathbf{g}(\cdot,\xi_{1}),\ldots,\mathbf{g}(\cdot,\xi_{T}) is independent. (Here, randomness comes only from the gradient oracle.)

Then, for any θ∈(4e−T/32,1)\theta\in(4e^{-T/32},1), with probability at least 1−θ1-\theta, the optimization error (i.e., the excess empirical risk) of this method is bounded as

Let xS∗∈arg⁡min⁡x∈XFS(x)x^{\ast}_{S}\in\arg\min\limits_{x\in\mathcal{X}}F_{S}(x). By convexity of the empirical loss, we have

Since (ξt)t∈[T](\xi_{t})_{t\in[T]} are sampled uniformly without replacement from S,S, we have

for all v∈X,t∈[T]v\in\mathcal{X},t\in[T]. Moreover, since the range of ff lies in [−LR,LR],[-LR,LR], it follows that Yt:=∑j=1tf(xj,ξj), t∈[T]Y_{t}:=\sum_{j=1}^{t}f(x^{j},\xi_{j}),~{}t\in[T] is a martingale with bounded differences (namely, bounded by 2LR2LR). Therefore, by Azuma’s inequality, the first term in (14) satisfies

By Hoeffding’s inequality, the second term in (14) also satisfies the same bound; namely,

Using similar analysis to that of the standard online gradient descent analysis Zinkevich, (2003), the last term in (14) can be bounded as

By the properties of the gradient oracle stated in the lemma, we can see that for any fixed realization of (xt,ξt)t∈[T](x^{t},\xi_{t})_{t\in[T]}, the second term in (17) is (2R+ηL)2σ2T\left(2R+\eta L\right)^{2}\frac{\sigma^{2}}{T}-sub-Gaussian random variable. Hence,

Putting (18) and (19) together, and noticing that T>32log⁡(4/θ)T>32\log(4/\theta), we conclude that with probability at least 1−θ/2,1-\theta/2, the third term of (14) is bounded as

Hence, by the union bound, we conclude that with probability at least 1−θ,1-\theta, the excess empirical risk of the stochastic subgradient method is bounded as

Appendix I Empirical Risk of Fixed-Permutation SGD

Our optimization error analysis is based on (Nedic and Bertsekas,, 2001, Lemma 2.1).

Let us consider the fixed permutation stochastic gradient descent (Algorithm 3), for arbitrary permutation (i.e., not necessarily random) and with constant step size over each epoch (i.e., η(k−1)n+t≡ηk\eta_{(k-1)n+t}\equiv\eta_{k} for all t∈[n]t\in[n], k∈[K]k\in[K]). Then

First, since the permutation is arbitrary, w.l.o.g. π\bm{\pi} is the identity (we make this choice only for notational convenience). Let now y∈Xy\in\mathcal{X}. At each round, the recursion of SGD implies that

Let rt:=∥xt−y∥r_{t}:=\|x^{t}-y\|. Summing up these inequalities from t=1,…,nt=1,\ldots,n

Re-arranging terms we obtain the result. ∎

Using the previous lemma, it is straightforward to derive the optimization accuracy of the method.

The fixed permutation stochastic gradient descent (Algorithm 3), for arbitrary permutation (i.e., not necessarily random) and with constant step size over each epoch (i.e., η(k−1)n+t≡ηk\eta_{(k-1)n+t}\equiv\eta_{k} for all t∈[n]t\in[n], k∈[K]k\in[K]). satisfies