Stability and Convergence Trade-off of Iterative Optimization Algorithms

Introduction

For different supervised learning algorithms ranging from classical linear regression, logistic regression, boosting, to modern large-scale deep networks, the overall performance or expected excess risk can always be decomposed into two parts: the empirical error (or the training error) and the generalization error (characterizing the discrepancy between the test error and the training error). A central theme in machine learning is to find an appropriate balance between empirical error and generalization error, because improperly emphasizing one over the other typically results in either overfitting or underfitting. Specifically, in the context of supervised learning models trained by iterative optimization algorithms, the empirical error at each iteration is commonly controlled by convergence rate analysis, and the generalization error can be handled by algorithmic stability analysis (Devroye and Wagner, 1979; Bousquet and Elisseeff, 2002).

Convergence rate of an algorithm portrays how fast the optimization error decreases as the number of iterations grows. Recent years have witnessed a rapid advance on convergence rates analysis of specific optimization methods for a particular class of loss functions that they are optimizing over. In fact, such analysis has been carried out for many gradient methods, including gradient descent (GD), Nesterov accelerated gradient descent (NAG), stochastic gradient descent (SGD), stochastic gradient Langevin dynamics (SGLD) for convex, strongly convex, or even nonconvex functions (see e.g. Boyd and Vandenberghe (2004); Bubeck et al. (2015); Nesterov (2013); Jin et al. (2017); Raginsky et al. (2017)). However, until the optimization error and generalization error of these algorithms are analyzed together, it is not clear whether the fastest converging optimization algorithm is the best for learning.

On the other hand, algorithmic stability (Devroye and Wagner, 1979; Bousquet and Elisseeff, 2002) in learning problems has been introduced as an alternative way to control generalization error instead of uniform convergence results such as classical VC-theory (Vapnik et al., 1994) and Rademacher complexity (Bartlett and Mendelson, 2003). The stability concept has an intuitive appeal: an algorithm is stable if it is robust to small perturbations in the composition of the learning data set. Recently it has been shown that algorithmic stability is well suited for controlling generalization error of stochastic gradient methods (Hardt et al., 2016), as well as stochastic gradient Langevin dynamics algorithm (Mou et al., 2017).

While most previous papers study convergence rate and the algorithmic stability of an optimization algorithm separately, a natural question arises: What is the relationship or trade-off between the convergence rate and the algorithmic stability of an iterative algorithm? Is it possible to design an algorithm that converges the fastest and at the same time most stable? If not, is there any fundamental limit on the trade-off between the two quantities so that a fast algorithm has to be unstable?

This paper shows that there is a fundamental limit on the trade-off. That is, for any iterative algorithms, at any time step, the sum of optimization error and stability is lower bounded by the minimax statistical error over a given loss function class. Therefore, a fast converging algorithm can not be too stable, and a stable algorithm can not converge too fast. This framework therefore provides a new criterion for comparing optimization algorithms by considering jointly convergence rate and algorithm stability. As a consequence, our framework can be immediately applied to provide a new class of convergence lower bounds for algorithms with different stability rates.

In particular, we focus on two settings where the loss functions are either convex smooth or strongly convex smooth. In the first setting, we discuss the stability upper bounds of gradient descent (GD), stochastic gradient descent (SGD) and their variants with decreasing step sizes. New stability upper bounds are provided for Nesterov’s accelerated gradient descent (NAG) and the heavy ball method (HB) under quadratic loss, and we conjecture these upper bounds still hold for the general convex smooth losses. Applying the stability upper bounds for GD and SGD in our trade-off framework, we obtain the convergence lower bounds for them that match the known convergence upper bounds up to constants. Considering jointly convergence rate and algorithm stability for NAG and GD, the trade-off shows that NAG must be less stable than GD even though it converges faster than GD. In the second setting where the loss functions are strongly convex and smooth, we also provide stability upper bound and deduce the convergence lower bound results for GD and NAG via our trade-off framework. Finally, simulations are conducted to show that the stability bounds established have the correct rates as a function of $n$ and iteration $T$. These bounds are demonstrated to be particularly useful in large scale learning settings for understanding the overall performance of an algorithm than the classical uniform convergence bounds because the stability bounds capture better generalization errors at early iterations of these algorithms.

The first quantitative results that focus on generalization error via algorithmic stability date back to (Rogers and Wagner, 1978; Devroye and Wagner, 1979). This line of research was further developed by Bousquet and Elisseeff (2002) to provide guarantees for general supervised learning algorithms and insights for the practice of regularized algorithms. It remains unclear, however, what is the algorithmic stability of general iterative optimization algorithms. Recently, to show the effectiveness of commonly used optimization algorithms in many large-scale learning problems, algorithmic stability has been established for stochastic gradient methods (Hardt et al., 2016), stochastic gradient Langevin dynamics (Mou et al., 2017), as well as for any algorithm in situations where global minima are approximately achieved (Charles and Papailiopoulos, 2017).

Given the importance of efficient optimization methods, many papers have been devoted to understanding the fundamental computational limits of convex optimization. Those lower bounds typically focus on a specific class of algorithms. A classical line of research has been focused on first-order algorithms where only first-order information (i.e. gradients) can be queried through oracle model; see the book (Boyd and Vandenberghe, 2004), the monograph (Bubeck et al., 2015) and references therein for further details. For convex functions, the first lower bound argument given in (Nemirovsky et al., 1982) applies to first-order algorithms whose current iterate lies in the linear span of previous gradients. It has been later extended to any deterministic, then stochastic first-order algorithm (Agarwal and Bottou, 2015; Woodworth and Srebro, 2016).

2 Organization of the paper

The rest of the paper is organized as follows: In Section 2, we set up the necessary backgrounds on the classical excess risk decomposition and introduce the optimization error (or computational bias) and generalization error trade-off. In Section 3, we provide the main theorem on the trade-off between convergence rate (as an upper bound on optimization error) and algorithmic stability (as an upper bound on generalization error). In Section 4, we establish uniform stability bounds for several gradient methods and show that our main theorem applies to these algorithms to obtain their convergence lower bounds. In Section 5, we first provide simulation results validating the correct rates as a function of sample size $n$ and iteration number $T$ of the stability bounds we established, and then illustrate via a simulated logistic regression problem that our stability bounds reflect the generalization errors better than the simple uniform convergence bounds for GD and NAG.

Preliminaries

In this section, we set up the necessary backgrounds on excess risk decomposition and convex optimization. Using classical excess risk decomposition, we introduce the expected optimization error and generalization error trade-off which are crucial to state our main result in the next section.

Given the collection $S$ of $n$ samples and a loss function $l$, the principle of empirical risk minimization is based on the objective function

This empirical risk above serves as a sample-average proxy for the population risk

We denote by ${\hat{\theta}}$ an estimator computed from sample $S$. The statistical question is how to bound the excess risk, measured in terms of the difference between the population risk and the minimal risk over the entire parameter space $\Omega$,

For simplicity, we assume that there exists some ${\theta_{0}}\in\Omega$ such that $R({\theta_{0}})=\inf_{\theta\in\Omega}R(\theta)$.If the infimum is not achieved within $\Omega$ (for example $\Omega$ is an open set), we can choose some ${\theta_{0}}$ where this equality holds up to some arbitrarily small error.

Controlling the excess risk of the estimator ${\hat{\theta}}$ is usually done by decomposing it into three terms as follows:

Term $T_{1}$ is the generalization error of the model ${\hat{\theta}}$. Term $T_{2}$ is the empirical risk difference between the model ${\hat{\theta}}$ and the population risk minimizer ${\theta_{0}}$. Term $T_{3}$ is the generalization error of ${\theta_{0}}$.

Making the optimization error appear in the decomposition is useful for analyzing optimization algorithms in an iterative manner. As noted in Bousquet and Bottou (2008), introducing optimization error allows to analyze algorithms doing approximate optimization. However, our framework is different to that introduced by Bousquet and Bottou (2008). We control the generalization error via iteration-dependent algorithmic stability instead of directly invoking uniform convergence results. As we are going to show, for most iterative optimization algorithms, upper bounding the generalization error by a simple uniform convergence is often loose and algorithmic stability can serve as a tighter bound.

2 Algorithmic Stability

Many forms of algorithmic stability have been introduced to characterize generalization error (Bousquet and Elisseeff, 2002; Kutin and Niyogi, 2002). For the purpose of this paper, we are only interested in the uniform stability notion introduced by Bousquet and Elisseeff (2002).

An algorithm, which outputs a model ${\hat{\theta}_{S}}$ for sample $S$, is $\epsilon$-uniform stable if for all $k\in\{1,...,n\}$, for all data sample pair $S=(z_{1},...,z_{k},...,z_{n})$ and $S^{\prime}=(z_{1},...,z_{k}^{\prime},...,z_{n})$, each $z_{i}$ or $z_{k}^{\prime}$ is i.i.d sampled from $P$, we have

As we did for the generalization error, we use $\mathcal{E}_{\text{stab}}(\hat{\theta},l,P,n)$ to denote the uniform stability of an algorithm $\hat{\theta}$.

A stable algorithm has the property that removing one element in its learning data set does not change much of its outcome. Such a data perturbation scheme is closely related to Jackknife in statistics (Efron, 1982). One can further show that uniform stability implies expected generalization (Bousquet and Elisseeff, 2002) . For completeness, we reformulate this property in the following lemma.

An algorithm, which outputs a model ${\hat{\theta}_{S}}$ for sample $S$, is $\epsilon$-uniformly stable, then its expected generalization error is bounded as follows,

Lemma 2 implies that $\mathcal{E}_{\text{gen}}(\hat{\theta},l,P,n)\leq\mathcal{E}_{\text{stab}}(\hat{\theta},l,P,n)$. The proof provided by Bousquet and Elisseeff (2002) relies on a symmetrization argument and makes use of the i.i.d assumptions of samples in $S$. Combining the expected excess risk decomposition in previous section, we conclude that the sum of uniform stability and expected optimization error (or computational bias) constitutes an upper bound for the expected excess risk,

Note that the result is stated for a fixed loss function $l$ and a fixed data distribution $P$. Equation (2) is a key inequality for our analysis. Not only it provides a way to upper bound the expected excess risk without uniform convergence results, but also it makes the connection between the statistical excess risk and the optimization convergence rate (or computational bias). This can also be seen as reminiscent of the bias-variance trade-off of an algorithm in a computational sense since stability serves as a computational variability term and optimization error as a computational bias term.

3 Convex optimization settings

Throughout the paper, we focus on two types of loss functions: The first type of loss function $l(\cdot,z)$ is $\alpha$-strongly convex and $\beta$-smooth for every $z$; The second type of loss function $l(\cdot,z)$ is convex and $\beta$-smooth for every $z$. We also make use of the $L$-Lipschitz condition. We provide their definitions here. More technical details about convex optimization and relevant results are deferred to Appendix B.

A function $f$ is $L$-Lipschitz if for all $u,v\in\Omega$, we have

A function continuously differentiable $f$ is $\beta$-smooth if for all $u,v\in\Omega$, we have

A function $f$ is convex if for all $u,v\in\Omega$, we have

A function $f$ is $\alpha$-strongly convex if for all $u,v\in\Omega$, we have

Trade-off between stability and convergence rate

In this section, we introduce the trade-off between stability and convergence rate via excess risk decomposition under two settings of loss functions mentioned in the previous section: the convex smooth setting and the strongly convex smooth setting. We show that for any iterative algorithm, at any time step, the sum of optimization error and stability is lower bounded by the minimax statistical error over a given loss function class. Thus algorithms sharing the same stability upper bound can be grouped to obtain convergence rate lower bounds. This provides a new class of convergence lower bounds for algorithms with different stability bounds.

We are interested in distribution independent stability and convergence where we take supremum of these two quantities over distributions and losses. For a fixed iteration algorithm that outputs $\hat{\theta}$ at iteration $T$, we define its uniform stability and optimization error as follows,

Note that in this paper, the supremum is taken over the class of all loss functions $\mathcal{L}$ under either of the two settings considered (convex smooth and strongly convex smooth settings).

Before we state the main theorem, we first define the loss function class of interest in this section. We define the class of all convex smooth loss functions as follows,

In the convex smooth setting, we have the following lower bound on the sum of stability and convergence rate.

Suppose an iterative algorithm outputs $\hat{\theta}_{T}$ at iteration $T$ on an empirical loss built upon a loss $l\in\mathcal{L}_{\text{c}}$ and an i.i.d. sample $S$ of size $n$, and it has uniform stability $\mathcal{E}_{\text{stab}}(T,n,\mathcal{L}_{\text{c}})$ and optimization error $\mathcal{E}_{\text{opt}}(T,n,\mathcal{L}_{\text{c}})$, then there exists a universal constant $C_{1}>0$ such that,

The first inequality of Theorem 7 is a simple outcome of the empirical risk decomposition in Equation (2). This first inequality is not tied to the convex smooth setting and can generalize to a wide class of optimization algorithms. The second inequality is based on an adaptation of the classical Le Cam (1986)’s method for minimax estimation lower bound to the convex smooth loss function class. Further, if we know $\mathcal{E}_{\text{stab}}^{\hat{\theta}}(T,n,\mathcal{L}_{\text{c}})$ precisely, we can obtain an immediate corollary that provide convergence lower bound for stable optimization algorithms.

Under conditions in Theorem 7, if an algorithms has uniform stability

with $s$ a divergent function of $T$, i.e.

then there exists a universal constant $C_{2}>0$, a sample size $n_{0}$ and an iteration number $T_{0}\geq 1$, such that for $T\geq T_{0}$, its convergence rate is lower bounded as follows,

Even though Theorem 7 is valid for any pair of $(T,n)$, Corollary 8 requires to choose a specific sample size $n_{0}$ in construction. However, under the assumption that the optimization algorithm has convergence rate independent of the sample size (i.e. $\mathcal{E}_{\text{opt}}^{\hat{\theta}}(T,n,\mathcal{L}_{\text{c}})$ is not a function of $n$), we can obtain via Corollary 8 a convergence lower bound that is comparable to the lower bounds in the convex optimization literature. We remark that this assumption is satisfied for commonly-used optimization algorithms such as GD and NAG.

Theorem 7 and Corollary 8 provide the trade-off between stability and optimization convergence rate. All iterative optimization methods that are algorithmic uniform stable can not converge too fast. This motivates the idea of grouping optimization methods with their algorithmic stability. Optimization methods that share the same algorithmic stability would have the same optimization lower bound. The proof of Theorem 7 is provided in Appendix A.1 and that of Corollary 8 in Appendix A.2.

2 Trade-off in the strongly convex smooth setting

Similar to the convex smooth setting, we define the class of all strongly convex smooth loss functions as follows,

In the strongly convex smooth setting, we have the following lower bound on the sum of stability and convergence rate.

Suppose an iterative algorithm outputs $\hat{\theta}_{T}$ at iteration $T$ on an empirical loss built upon a loss $l\in\mathcal{L}_{\text{sc}}$ and an i.i.d. sample $S$ of size $n$, and it has uniformly stability $\mathcal{E}_{\text{stab}}^{\hat{\theta}}(T,n,\mathcal{L}_{\text{sc}})$ and has optimization error $\mathcal{E}_{\text{opt}}^{\hat{\theta}}(T,n,\mathcal{L}_{\text{sc}})$, then there exists a universal constant $C_{3}$ such that

The trade-off in the strongly convex smooth setting is similar to that of convex smooth setting, except that the minimax estimation rate is of order $O(\frac{1}{n})$ instead of $O(\frac{1}{\sqrt{n}})$. Theorem 9 provides the trade-off between stability and optimization convergence rate in the strongly convex setting. Note that a similar corollary like Corollary 8. The proof of Theorem 9 is provided in Appendix A.3.

Stability of first order optimization algorithms and implications for convergence lower bounds

This section is devoted to establishing stability bounds of popular first order optimization algorithms and showing that our main theorem can be applied to these algorithms to obtain their convergence lower bounds. In particular, Subsection 4.1 establishes uniform stability for first order iterative methods in the convex smooth setting and Subsection 4.2 discusses the consequence after applying Theorem 7 to various optimization algorithms. Subsection 4.3 provides uniform stability for first order iterative algorithms in the strongly convex smooth setting and Subsection 4.4 discusses the consequence after applying Theorem 9 to GD and NAG.

The goal of proving uniform stability for iteration $T$ is to bound the difference

for the sample $S=(z_{1},\ldots,z_{k},\ldots,z_{n})$ and the perturbed one $S^{\prime}=(z_{1},\ldots,z_{k}^{\prime},\ldots,z_{n})$, uniformly for every $z\in\mathcal{Z}$. $z_{1},\ldots,z_{k},\ldots,z_{n}$ and $z_{k}^{\prime}$ are drawn i.i.d from a distribution $P$. Here ${\hat{\theta}}_{S,T}$ denotes the output model of our optimization algorithm at iteration $T$ based on sample $S$. The optimization algorithm is applied on a pair of data samples $S,S^{\prime}$ to get two sequences of successive models ${\hat{\theta}}_{S,0},{\hat{\theta}}_{S,0},\ldots,{\hat{\theta}}_{S,T}$ and ${\hat{\theta}}_{S^{\prime},0},{\hat{\theta}}_{S^{\prime},1},\ldots,{\hat{\theta}}_{S^{\prime},T}$. For simplicity, we use ${\hat{\theta}}_{t}$ to denote ${\hat{\theta}}_{S,t}$ and ${\hat{\theta}}_{t}^{\prime}$ for ${\hat{\theta}}_{S^{\prime},t}$. We first bound the model estimate difference $\left\|{\hat{\theta}}_{t}-{\hat{\theta}}_{t}^{\prime}\right\|_{2}$, then use the $L$-Lipschitz condition of $l$ to prove stability.

Recall that the empirical loss function for data sample $S=(z_{1},\ldots,z_{n})$ is

where we have replaced $l(\theta;z_{j})$ with $f_{j}(\theta)$ to improve readability. On the other hand, the empirical loss function for the perturbed sample $S^{\prime}=(z_{1},\ldots,z_{k}^{\prime},\ldots,z_{n})$ is

Remark that the two empirical loss functions only differ on one term that is proportional to the inverse of sample size $n$.

We establish uniform stability for gradient descent, stochastic gradient descent, Nesterov accelerated gradient method and heavy ball method with fixed momentum parameter when the loss function is convex smooth.

The gradient descent algorithm is an iterative method for optimization, which uses the full gradient at each iteration (See book by Boyd and Vandenberghe (2004)). Given a convex smooth objective $F$, GD starts at some initial point $\theta_{0}\in\Omega$, and iterates with the following recursion

where $\eta$ is the step-size. Typically, one would choose fixed $\eta\leq\frac{1}{\beta}$ to ensure convergence (Boyd and Vandenberghe, 2004). In the empirical risk minimization setting, the objective $F$ of the optimization is either $R_{S}$ or $R_{S^{\prime}}$.

Given a data distribution $P$, under the assumption that $l(\cdot,z)$ is a convex, $L$-Lipschitz and $\beta$-smooth function for every $z\in\mathcal{Z}$, the gradient method with constant step-size $\eta\leq\frac{1}{\beta}$ on the empirical risk $R_{S}$ with sample size $n$, which outputs $\hat{\theta}_{T}$ at iteration $T$, has the following uniform stability bound for all $T\geq 1$,

We remark that this stability bound does not depend on the exact form of the loss function $l$ and the exact form of the data distribution $P$. The proof of this theorem is provided in Appendix B.1. The key step of our proof is that in such a set-up, the error caused by the difference in empirical loss functions accumulates linearly as the iteration increases. We also show in Appendix B.1 that this stability upper bound can be achieved by a linear loss function.

1.2 Nesterov accelerated gradient methods (NAG)

The Nesterov’s accelerated gradient method attains the optimal convergence rate $O(1/T^{2})$ in the smooth non-strongly convex setting under the deterministic first order oracle (Nesterov, 1983). Given a convex smooth objective $F$, starting at some initial point $\theta_{0}=w_{0}\in\Omega$, NAG uses the following updates,

where $\eta\leq\frac{1}{\beta}$ is the step-size. The parameter $\gamma_{t}$ is defined by the following recursion

satisfying $-1<\gamma_{t}\leq 0$. We only provide a uniform stability bound for NAG when the empirical risk function is quadratic. We conjecture that the same stability bound holds for general convex smooth functions.

Given a data distribution $P$, under the assumption that $l(\cdot,z)$ is a $L$-Lipschitz, $\beta$-smooth convex quadratic loss function defined on a bounded domain for every $z\in\mathcal{Z}$, Nesterov accelerated gradient method with fixed step-size $\eta\leq\frac{1}{\beta}$, which outputs $\hat{\theta}_{T}$ at iteration $T$, has the following uniform stability bound for all $T\geq 1$,

The proof of the theorem is provided in Appendix B.2. We also show in Appendix that this stability upper bound is achieved by a linear loss function. Note that unlike the full gradient method and stochastic gradient descent, the stability bound of Nesterov accelerate gradient method depends quadratically on the iteration $T$. Even though NAG can still have small stability when early stopping is used, its stability grows faster than that of GD at the same iteration.

1.3 The heavy ball method with a fixed momentum

The heavy ball method (HB), like NAG, is also a multi-step extension of the gradient descent method (Polyak, 1964). Fixed step-size and fixed momentum parameter heavy ball method has the following updates. For $t\geq 1$,

with fixed $\gamma\in[0,1),\eta\in\left(0,\frac{2(1-\gamma)}{\beta}\right)$. As for the NAG, we provide only a uniform stability bound for the heavy ball method when the empirical risk function is quadratic. We conjecture that the same stability bound holds for general convex smooth functions.

Given a data distribution $P$, under the assumption that $l(\cdot,z)$ is a $L$-Lipschitz, $\beta$-smooth convex quadratic loss function defined on a bounded domain for every $z$, the heavy ball method with a fixed step-size $\eta\in\left(0,\frac{(1-\gamma)}{\beta}\right)$ and a fixed momentum parameter $\gamma\in[0,1)$, which outputs $\hat{\theta}_{T}$ at iteration $T$, has the following uniform stability bound for all $T\geq 1$,

The proof of this theorem is provided in Appendix B.3. This theorem shows that the Heavy ball method with a fixed step-size and a fixed momentum parameter also uses multi-step gradients, it is more stable than NAG with a stability bound of order $O(T/n)$. This demonstrates that the multi-step setup does not necessarily lead to a similar or worse stability bound than that of NAG.

1.4 Other methods with known stability

In this subsection, we restate the stability bounds of some other gradient methods in this subsection for completeness. The stability bounds stated in this subsection are not new, but they serve as basis of our discussion for their convergence lower bounds implied by Theorem 7 in Subsection 4.2.

The stochastic gradient descent is a randomized iterative algorithm for optimization. Instead of using the full gradient information, it randomly chooses one data sample and updates the parameter estimate according to the gradient on that sample. It starts at some initial point $\theta_{0}\in\Omega$, and iterates with the following recursion with $i$ chosen from the set $\{1,...,n\}$ uniformly at random:

Hardt et al. (2016) adapted the definition of uniform stability to randomized algorithms and showed that the fixed step-size $\eta\leq\frac{1}{\beta}$ stochastic gradient descent has a $\frac{2\eta L^{2}T}{n}$-uniform stability bound in the convex, $L$-Lipschitz and $\beta$-smooth setting. According to Theorem 3.8 in Hardt et al. (2016), we have

for any convex $L$-Lipschitz and $\beta$-smooth loss function $l$. This is a restatement of the result of Hardt et al. (2016) in our notation.

Hardt et al. (2016) further considers stochastic gradient descent with decreasing step-sizes $\eta_{t}=t^{-\alpha}$ and shows that stochastic gradient descent with decreasing step-sizes has $\frac{2\eta L^{2}T^{1-\alpha}}{n}$-uniform stability in the same setting.

SGLD plays an important role in sampling and optimization. It is proposed as a stochastic discrete version of the Langevin Equation $d\theta_{t}=-\nabla f(\theta_{t})dt+\sqrt{\frac{2}{\tau}}dB_{t}$, where $B_{t}$ is the Brownian motion. Recent work by Raginsky et al. (2017) has shown its effective in non-convex learning with optimization and generalization guarantees.

When SGLD is applied to optimization, a decreasing step with $\eta_{t}=O(\eta_{0}/t)$ should be used to ensure convergence to local minima. We study this particular step-size setting of SGLD. It has been shown by Mou et al. (2017) that SGLD has the following uniform stability for $L$-Lipschitz convex loss function,

where $k_{0}=\min\left\{t|\eta_{k}\tau L^{2}<1\right\}$. Plugging in the $O(\eta_{0}/t)$ step-size, we have that SGLD has a uniform stability bound

at iteration $T\geq 1$, for any convex $L$-Lipschitz and $\beta$-smooth loss function $l$. This is an adaptation of the result of Mou et al. (2017) in our notation.

2 Consequences for the convergence lower bound in convex smooth setting

In this section, we apply Theorem 7 and Corollary 8 to obtain convergence lower bounds for a variety of first order optimization algorithms mentioned above. Furthermore, we compare the convergence lower bound we obtain with the known convergence upper bound for each of the optimization methods mentioned in the previous section. The known convergence upper bounds mentioned in this section can be found in the optimization textbooks (See Boyd and Vandenberghe (2004) or Bubeck et al. (2015)). We also discuss how our lower bounds compare to those obtained from classical oracle model of complexity by Nemirovsky et al. (1982).

Note that the assumptions in Theorem 7 are slightly different to what we use when we establish stability bounds in the previous section: the former assume bounded domain $R$ while the latter assume $L$-Lipschitz. To make these two assumptions compatible, in this subsection, we assume that the domain $R=\left|\Omega\right|$ is fixed and for all $z\in\mathcal{Z}$, there exists ${\theta^{*}}\in\Omega$ such that $\nabla l(\cdot,z)=0$. Then we have the loss is $L$-Lipschitz with $L\leq R\beta$. This is because for any $\theta\in\Omega$,

In Table 1, we summarize all the uniform stability results and the corresponding convergence lower bound under convex smooth setting. While exact constants are provided in the main text, we only show the dependency on iteration number $T$ and sample size $n$ in the table.

According to Equation (3) in Theorem 10, the fixed-step-size full gradient method has $\frac{2\eta\left(R\beta\right)^{2}T}{n}$-uniform stability. Applying Corollary 8, knowing that its convergence does not depend on $n$, we obtain that its convergence rate is lower bounded by

The convergence rate lower bound obtained via our stability trade-off thus matches the known upper bound up to constant factors.

2.2 Stochastic gradient descent

According to Hardt et al. (2016), the fixed step-size stochastic gradient descent also has $\frac{2\eta\left(R\beta\right)^{2}T}{n}$-uniform stability. Applying Corollary 8, we obtain a convergence rate lower bound of order $O(1/T)$. However, it is known that fixed-step-size stochastic gradient descent can not converge arbitrarily small error at the rate $O(1/T)$ (Delyon and Juditsky, 1993). The best rate of convergence to minimize a smooth non-strongly convex function with noisy gradients is of order $O(T^{-\frac{1}{2}})$ (Nemirovski et al., 2009). Therefore, in the case of fixed step-size SGD, the convergence lower bound we provide is valid but loose. The fixed step-size SGD is a stable algorithm but is not a convergent algorithm.

On the other hand, it is shown in the same work (Nemirovski et al., 2009) that $O(T^{-\frac{1}{2}})$ convergence rate is achieved by stochastic gradient descent with decreasing step-size of order $O(T^{-\frac{1}{2}})$. Using our stability argument, we provide insights why the stochastic gradient descent with decreasing step-size is not converging too fast. It has also been shown by Hardt et al. (2016) that stochastic gradient descent with decreasing step-size of order $O(T^{-\frac{1}{2}})$ has $O(\sqrt{T}/n)$ uniform stability. Applying Corollary 8, we conclude that when this decreasing step-size is used, gradient descent can not converge as fast as $O(T^{-1})$.

Similar arguments can be used to explain the conjecture by Moulines and Bach (2011) on the optimal convergence rates for stochastic gradient descent of $O(T^{-\alpha})$ step-size. It is shown in Moulines and Bach (2011) that, for $\alpha\in(2/3,1)$, the convergence rate of stochastic gradient descent for the convex $\beta$-smooth case is upper bounded by $O(T^{\alpha-1})$. It is shown by Hardt et al. (2016) that stochastic gradient descent of $O(T^{-\alpha})$ step-size has $O(T^{1-\alpha}/n)$ uniform stability in this set-up. Applying Corollary 8, we provide a proof of this conjecture, confirming the optimality of this convergence rate upper bound.

2.3 Nesterov accelerated gradient descent

According to Theorem 11, the Nesterov accelerated gradient descent with fixed step-size has $\frac{4\eta\left(R\beta\right)^{2}T^{2}}{n}$-uniform stability for quadratic loss functions. Under the conjecture that the same stability holds for convex smooth loss functions, according to Corollary 8, we could obtain that its convergence rate is lower bounded by

This is compatible with its convergence rate upper bound provided in Nesterov (1983). For $f$ convex and $\beta$-smooth function, Nesterov accelerated gradient method with step-size $\eta\leq\frac{1}{\beta}$ satisfies

We can compare our stability based lower bounds to classical ways of getting complexity lower bound using the classical first-order oracle of complexity (Nemirovsky et al., 1982; Nesterov, 2013). The classical oracle model based lower bound provides $O(1/T^{2})$ lower bound for all first order optimization methods that falls into the following black-box framework. It assumes that the optimization methods takes initialization $\theta_{1}=0$ and at iteration t, $\theta_{t}$ is in the linear span of all previous gradients. Whereas our results show that all optimization methods with order $O(T^{2}/n)$ uniform stability in the smooth non-strongly convex setting would have convergence rate lower bounded by $O(1/T^{2})$. The two lower bounds have similar form, but apply under different scenarios. One remarkable property of our result is that it does not depend on how exactly the algorithm is initialized.

2.4 Heavy ball method with fixed step-size

According to Theorem 12, heavy ball method with fixed step-size $\eta\in\left(0,\frac{(1-\gamma)}{\beta}\right)$ and fixed momentum parameter $\gamma\in[0,1)$ has

uniform stability for quadratic loss functions. Under the conjecture that the same stability holds for convex smooth loss functions, applying Corollary 8, we obtain that its convergence rate is lower bounded by $O(1/T)$. First, this lower bound matches the convergence rate upper bound proved in Ghadimi et al. (2015). Second, unlike Nesterov accelerated gradient descent, even though multiple steps of gradients are used, heavy ball method with fixed step-size is not able to achieve the optimal convergence rate $O(1/T^{2})$. Another viewpoint on this result is that the smart choice of weighting coefficients in NAG is necessary to its optimal convergence guarantee.

2.5 Stochastic gradient Langevin dynamics (SGLD)

According to Mou et al. (2017), stochastic gradient Langevin dynamics with temperature $\tau$ and decreasing step-size $O(1/T)$, when used for convex optimization, has

uniform-stability. Applying Corollary 8, we conclude that its convergence rate is lower bounded by $O(1/T^{1/4})$. While the additional noise added in SGLD might be helpful for certain non-convex optimization settings in escaping local minima as stated in Mou et al. (2017), SGLD has a slower worst-case convergence than the GD or SGD based on our stability argument.

3 Stability in the strongly convex smooth setting

In this subsection, we establish uniform stability for gradient descent, Nesterov accelerated gradient method in the strongly convex smooth setting. In the strongly convex smooth setting, the loss function $l(\cdot,z)$ is $\alpha$ strongly-convex, $\beta$-smooth for every $z\in\mathcal{Z}$.

The gradient descent method in the strongly convex setting has exactly the same updates as before, given a strongly convex smooth objective $F$, for $t\geq 0$,

where $\eta\leq 1/\beta$ is the step-size. While the algorithm stays the same, the strongly convex property of the loss function allows the algorithm to have a better stability.

Given a data distribution $P$, under the assumption that $l(\cdot,z)$ is $\alpha$-strongly convex, $\beta$-smooth and $L$-Lipschitz for every $z\in\mathcal{Z}$, the full gradient method with constant step-size $\eta\leq\frac{1}{\beta}$, which outputs $\hat{\theta}_{T}$ at iteration $T\geq 1$, has uniform stability

The proof of this theorem is provided in Appendix C.1.

3.2 Stochastic gradient descent (SGD) with fixed step-size

The stochastic gradient descent in the strongly convex setting has the exactly same updates as before. It starts at some initial point $\theta_{0}\in\Omega$, and iterates with the following recursion with $i$ chosen from the set $\left\{1,...,n\right\}$ uniformly at random,

The stability of SGD under strongly convex setting has been first discussed in Hardt et al. (2016). According to Theorem 3.10 in Hardt et al. (2016), the stability of SGD under strongly convex setting is upper bounded by

at iteration $T\geq 1$, for any $\alpha$-strongly convex, $L$-Lipschitz and $\beta$-smooth loss function $l$.

3.3 Nesterov accelerated gradient descent (NAG)

Unlike in the convex smooth setting, Nesterov’s accelerated gradient descent can take fixed momentum parameter in the strongly convex smooth setting.

where $\eta\leq\frac{1}{\beta}$ is the step-size, $\kappa=\beta/\alpha$.

We prove its uniform stability for $\alpha$ strongly-convex, $\beta$-smooth for quadratic loss function.

Given a data distribution $P$, under the assumption that $l(\cdot,z)$ is $\alpha$-strongly convex, $\beta$-smooth and $L$-Lipschitz for every $z\in\mathcal{Z}$, Nesterov accelerated gradient descent method described above, which outputs $\hat{\theta}_{T}$ at iteration $T\geq 1$, has uniform stability

The proof of this theorem is provided in Appendix C.2.

4 Consequences for the convergence lower bound in the strongly convex setting

In this subsection, we obtain convergence lower bound for GD and NAG in the $\alpha$-strongly convex $\beta$-smooth setting via Theorem 9. In Table 2, we summarize all the uniform stability results and the corresponding convergence lower bounds under strongly convex smooth setting. While exact constants are provided in the main text, we only show the dependency on iteration number $T$ and sample size $n$ in the table.

According to Theorem 13, gradient descent with fixed step-size $\eta$ in the strongly convex smooth setting has

uniform stability. We apply Theorem 9 to obtain a lower bound on the convergence of GD for strongly convex smooth functions.

If the leading constants $\frac{\beta R^{2}}{C_{3}}$ and $\frac{4\left(R\beta\right)^{2}}{\alpha}$ match, we could directly obtain a lower bound on its convergence of order $e^{-O\left(T/\left(1+\kappa\right)\right)}$ as we expect. Unfortunately, due to our proof of the empirical risk minimization lower bound, a couple factors of constants are lost. Thus directly applying the stability bound makes it impossible to match the leading constants. We always have

Therefore, our trade-off result only gives convergence lower bound of GD with an offset of $\frac{\beta R^{2}}{C_{3}n}-\frac{4\left(R\beta\right)^{2}}{\alpha n}$ as stated in Equation (13).

Remark that a similar lower bound can be obtained for stochastic gradient descent using exactly the same argument for GD.

4.2 Nesterov accelerated gradient descent

According to Theorem 14, Nesterov accelerated gradient descent with fixed step-size $\eta$ in the strongly convex smooth setting has

uniform stability for quadratic loss function. Since the construction of the minimax lower bound in Theorem 9 is based on quadratic loss functions, applying Theorem 9 by restricting to quadratic loss functions, we obtain an expected convergence lower bound of order $e^{-O\left(T/\sqrt{\kappa}\right)}$ with an offset,

Simulations Experiments

In this section, we first show via simulation results of a simple logistic regression applied on breast-cancer-wisconsin dataset that the stability bounds established in this paper have the right scaling on the iteration number $T$. Second, we illustrate via a logistic regression problem that the stability bound characterize better the generalization error than simple uniform convergence bound at least for the first iterations of GD and NAG.

We evaluate our stability bounds for all gradient methods mentioned on logistic regression with the binary classification datasets breast-cancer-wisconsin (Wolberg and Mangasarian, 1990). This dataset has sample size $n=699$ and dimension $d=10$. The problem of logistic regression is formulated as follows.

Let $Y=\left(Y_{1},\ldots,Y_{n}\right)^{\top}\in\left\{0,1\right\}^{n}$ and $\mathbf{X}$ be the $n\times d$ matrix with $X_{i}$ as $i^{\text{th}}$-row. The log-likelihood function we optimize over is as follows,

It can be shown that this objective has the Lipschitz constant $L$ equal to $1$ and the smoothness parameter $\beta$ equal to $1/4$ when the covariate matrix $\mathbf{X}$ is normalized to have its maximum eigenvalue equal to $1$. When there is no regularization, each loss function $f_{i}$ is not strongly convex $\mu=0$. In all of our experiments we set constant step-size $\eta=0.1$. To construct samples that differ only on one data point, we first fix a sample $S$ with size $500$ from dataset, then construct a perturbed sample $S^{\prime}$ by changing one data point in $S$ and finally run our optimization algorithm to compute and plot the model difference $\left\|\theta_{t}-\theta_{t}^{\prime}\right\|_{2}$. The norm difference $\left\|\theta_{t}-\theta_{t}^{\prime}\right\|_{2}$ constitute an estimate for the uniform stability up to constants independent of $T$ and $n$. Finally, the perturbation on the sample is repeated $50$ times. Figure 1 shows the estimated uniform stability, averaged over $50$ independent repeats, for all gradient methods methods, Nesterov accelerated gradient, heavy ball method with fixed momentum ($\gamma=0.8$), full gradient method with fixed step-size, full gradient method with decreasing step-size $T^{-\alpha}$ ($\alpha=0.5,0.3$), stochastic gradient method with fixed step-size and stochastic gradient method with decreasing step-size $T^{-\alpha}$ ($\alpha=0.5$). We observe that the estimated uniform stabilities of full gradient method, stochastic gradient method and heavy ball method with fixed step-size all have slope $1$ in log-log plot, while Nesterov accelerated gradient method has slope $2$. As expected, methods with decreasing step-size have a slope smaller than $1$. Even though the stability bounds of NAG and HB are only established for quadratic loss, the estimated stability in the simulation makes us conjecture that the stability bounds of NAG and HB still hold in the general convex smooth setting.

2 Algorithmic stability vs simple uniform convergence bounds

The goal of this simulation is to show that algorithmic stability characterize the generalization error better than the simple uniform convergence bounds, which can not easily take into account of the growth of the function space for iterative algorithms. For $d$-dimensional estimation problem, simple uniform convergence bound would give an generalization error bound of order $O\left(\sqrt{d/n}\right)$. The exact constant in the uniform convergence bound depends on the function space and is hard to characterize for iterative algorithms. We think that more refined uniform convergence bound via Rademacher complexity (Bartlett and Mendelson, 2003) might be possible, but we are not aware of such results for general iterative algorithms. In this section, we show via simulations that the simple uniform convergence bound of order $O\left(\sqrt{d/n}\right)$ is less precise than the stability in characterizing generalization error. More precisely, we can see that when the dimension $d$ and the number of samples $n$ are large and iteration $T$ is small

where $s(T)/n$ is the stability bound for GD or NAG. We show in the next two experiments that this comparison is valid and the stability bound is more relevant in large scale problems.

In the both experiments, we fix the true parameter $\theta^{*}=\left(1,\ldots,1\right)^{\top}$ and we random draw $n$ i.i.d. samples $\left(X_{i},Y_{i}\right)$ according to the following data generation process. Each row of $\mathbf{X}$ is drawn from a standard $d$-dimensional normal distribution, and then $\mathbf{X}$ is renormalized to have row norm $1$. Each label $Y_{i}$, give $X_{i}=x$, is drawn from a Bernoulli distribution with parameter $r(\theta^{*},x)$. We use both the gradient descent and Nesterov accelerated gradient to optimize the empirical log-likelihood objective in Equation $\eqref{eq:logistic_log_likelihood}$. We estimate the stability using its definition in Equation (1) by varying different $z$ from holdout data set. In first experiment, we set $d=20,n=2000$. Figure 2 shows that both the simple uniform convergence bound and estimated stability bound are small compared to optimization error. In this setting, driving optimization error to zero is more important for reducing the test error, as shown in thick red color. We can still observe that the scalings of the estimated stability bound for GD and NAG are different. Our theoretical stability bound follows the estimated stability bound with the same slope, but without the saturation at the end of iterates.

In the second experiment, we set $d=200,n=2000$. Figure 3 shows that the generalization error accounts for a large portion of the test error. Especially, we observe in Figure 3 that the test error of NAG deviates from its training error. Simple uniform convergence bound does not explain the overfitting phenomenon here. The algorithmic stability combined with the training error suggests that early-stopping should be used for NAG in this setting as shown in Figure 4.

Discussions

In this section, we discuss how our stability bound for optimization could served as an early stopping criteria. We also discuss other iterative algorithms such as boosting that could fit into this stability and optimization trade-off framework.

Minimizing empirical risk is often computationally expensive in large scale learning problems. As it has been pointed out in Bousquet and Bottou (2008), optimization algorithms do not need to carry out this minimization with great accuracy since the empirical risk is already an approximation to the expected risk. For example, we can stop an iterative optimization algorithm long before its convergence to reduce computational cost. How early we should stop without deteriorating too much the expected risk becomes the main question we ask in large scale learning problems. The expected excess risk decomposition has been the main theoretical guideline for this kind of early-stopping criteria. Even though in this reasoning we are studying upper bounds of generalization and optimization errors, it is often accepted that these upper bounds give a realistic idea of the actual convergence rates (Vapnik et al., 1994; Bousquet and Elisseeff, 2002; Bartlett et al., 2006; Bousquet and Bottou, 2008).

We would like to stop our optimization algorithm as far as it reaches an optimization error close to its generalization error. However, the uniform convergence bounds are often too pessimistic about the size of the space to search over. Instead, we use our stability based generalization bound as an estimate of the generalization error. Formally, we would choose iteration $T$ such that

As an example, our stability based generalization bound for fixed-step-size full gradient method in the smooth non-strongly convex setting is $\frac{2\eta L^{2}T}{n}$. The first remarkable point is that this generalization error bound is dimension-free. Because it is often hard to access accurate estimates for the uniform convergence bounds based generalization error, it might be advantageous to acquire a theoretical early-stopping criterion via our stability bounds. For the full gradient method trained model, as long as the Lipschitz constant $L$ and smoothness constant $\beta$ can be estimated accurately, we are able to give an early stopping criterion such as $T\approx\sqrt{\frac{n}{\eta^{2}L^{2}R^{2}}}$, given the estimate of $R$ is accurate.

2 Other iterative optimization algorithms such as boosting

Boosting is one of the most successful and practical iterative optimization methods. Unlike gradient method which iterates over parameters, boosting starts with a sensible estimator or classifier, the learner, and seeks its improvements iteratively on the function space. The bias-variance trade-off of L2 boosting discussed in Bühlmann and Yu (2003) shares similar behaviors as the trade-off we discussed in Equation (2). It would be interesting to characterize the stability of boosting algorithms with various kinds of weaker learners and derive precise trade-off results as we did for gradient methods.

This research is supported in part by ONR Grant N00014-16-2664 , NSF Grants DMS-1613002 and IIS 1741340, and the Center for Science of Information (CSoI), a US NSF Science and Technology Center, under grant agreement CCF-0939370. We would like to thank Raaz Dwivedi and Rebecca Barter for fruitful discussions on this topic.

References

A Proof of Main Results

Using Equation (2), Theorem 7 directly follows from the well-known statistical lower bound for empirical risk estimation with adaptation to convex smooth loss functions. For completeness, we restate this lower bound and provide the proof below.

For any fixed sample size $n$, there exists a universal constant $C_{1}>0$ and $\beta$-smooth convex loss function $l$ defined on $\mathcal{Z}\times\Omega$, with $R=\left|\Omega\right|$, such that

The main idea to prove this lemma is to formulate the excess risk minimization problem as binary hypothesis testing problem and then apply Le Cam’s method for lower bound.

For any fixed sample size $n$, define domain $\mathcal{Z}$ be $\left\{-1,1\right\}$ and two probability distributions $P_{1}$ and $P_{2}$ satisfying the following two properties,

We define $P_{1}^{n}$ to be the joint distribution where $Z_{1},\ldots,Z_{n}$ are independent samples from $P_{1}$, and we defin $P_{2}$ accordingly.

Let $\theta_{1}^{*}\in\Omega$ with all other coordinates zero but the first coordinate equals to $-\delta$, and $\theta_{2}^{*}\in\Omega$ with all other coordinates zero but the first coordinate equals to $\delta$, with $0<\delta\leq r$. $\delta$ and $r$ are a constants to be determined later. Let $\theta$ be the first coordinate of $\theta$ and let $\Phi(r)$ be the parameter such that

The exact form of $\Phi(r)$ will be determined after we define the loss function $l$. We have

Le Cam’s method reduce this estimation problem to binary hypothesis testing problem, then we have

where the infimum ranges over all testing functions $\Psi:\mathcal{Z}^{n}\rightarrow\left\{1,2\right\}$.

We have for any $\Psi:\mathcal{Z}^{n}\rightarrow\left\{1,2\right\}$ that the probability of error is

A standard result of Le Cam (1986) gives the exact expression of the minimal possible error in the above hypothesis test. We have

where $\left\|\cdot\right\|_{\text{TV}}$ denotes the total variation distance. Using Pinsker’s inequality, we have

Equality $(i)$ uses the KL divergence formula between two Bernoulli distributions. Inequality $(ii)$ uses the inequality $\delta\log\frac{1+\delta}{1-\delta}\leq 3\delta^{2}$ for $\delta\in\left[0,\frac{1}{2}\right]$. Thus, we show that any test $\Psi$ mistakes one of the probability distribution for the other with probability at least $\frac{1}{4}$.

It remains to design a $\beta$-smooth convex loss function $l$ and determine the exact form of $\Phi$. Without loss of generality, we can assume that $\Omega$ is center around . We define the loss function $l(\theta;z)$ to be

It is easy to verify that the loss function is convex and $\beta$-smooth for each $z$. Then

where $\theta_{1,\text{left}}$ is zero everywhere but $-\frac{3r}{2}$ on the first coordinate. Then

and the same holds for $P_{2}$. Plugging $\Phi(r)=\frac{\beta r^{2}}{\sqrt{96n}}$ into Equation (16), we can conclude that

We remark that we can take $r$ as large as $\frac{R}{2}$. Thus we conclude that

A.2 Proof of Corollary 8

Applying Theorem 7, for any sample size $n$ and $T$, we have

As we only consider optimization method designed for any convex problems, $\mathcal{E}_{\text{optimization}}$ is independent of the sample size $n$. This result is valid for any sample size $n$. We can take $n$ such that the following quadratic function

is maximized to obtain the best lower bound.

$\frac{2C_{1}s(T)}{R^{2}\beta}$ would be the best choice of $\sqrt{n}$, but we have to ensure that $n$ is an integer. Since $s(T)$ is divergent function of $T$, there exists $T_{0}\geq 1$, such that for $T\geq T_{0}$, we can always find integer $n$ satisfying

Plugging $n$, we conclude that there exists universal constant $C_{2}$, and a convex function such that for $T\geq T_{0}$,

A.3 Proof of Theorem 9

We prove the statistical lower bound for empirical risk estimation in the strongly convex case via similar techniques used in the proof of Lemma 15. Le Cam’s argument for reducing an estimation problem to binary hypothesis testing problem is still valid. All we do is to define a $\alpha$-strongly convex $\beta$-smooth loss function $l$ and find the corresponding $\Phi(r)$. We define the loss function $l(\theta;z)$ to be

$l$ is quadratic, so it is $\alpha$-strongly convex and $\beta$ smooth for each $z$. Then

The minimizer $\theta_{1}^{*}$ has the first coordinate equals to $-\frac{r}{\sqrt{6n}}$. And the minimum is $\frac{\beta}{2}\left(r^{2}-\frac{r^{2}}{6n}\right)$.

For $\theta^{\prime}\in\Omega$ such that $\left|\theta^{\prime}-\theta_{1}^{*}\right|\geq r$, we have

The same lower bound holds for $P_{2}$. Plugging $\Phi(r)=\frac{\beta r^{2}}{12n}$ into Equation (16), we can conclude that

We remark that we can take $r$ as large as $\frac{R}{2}$. Thus we conclude that

B Stability Bounds for Convex Smooth Functions

In this section, we prove stability bounds of optimization algorithms (GD, NAG and heavy ball methtod) for convex smooth functions.

Before we proceed to the main proof, we state several well known lemmas about convex optimization which can be found in Boyd and Vandenberghe (2004); Bubeck et al. (2015). The $\beta$-smoothness of a function directly implies the following two lemmas. These two lemmas characterize how well the gradient approximation works for $\beta$-smooth functions in terms of both upper and lower bounds.

Let $f$ be a $\beta$-smooth function on $\Omega$. Then for all $u,v\in\Omega$, we have

Let $f$ be a convex and $\beta$-smooth function on $\Omega$. Then for any $u,v\in\Omega$, we have

An immediate corollary could be obtained by applying from the Lemma 16 to $(u,v)$ and then $(v,u)$. This corollary directly implies the constracting property of the gradient decent method, which is the key component for providing its algorithmic uniform stability.

Let $f$ be a $\beta$-smooth function on $\Omega$. Then for any $u,v\in\Omega$, one has

Recall that in order to prove the uniform stability, we need to bound the loss difference for any fixed sample $z$ at each iteration $t\geq 1$

This quantity is related to the norm difference $\left\|\theta_{t}-\theta^{\prime}_{t}\right\|_{2}$ under the $L$-Lipschitz condition. Using the update rule of full gradient method, we obtain an recursive relation on $\left\|\theta_{t}-\theta^{\prime}_{t}\right\|_{2}$. For $\eta\leq\frac{1}{\beta}$ and $t\geq 1$, we have

The inequality $(i)$ uses triangular inequality. The inequality $(ii)$ follows from the $L$-Lipschitz condition on the perturbed gradient terms. The last inequality $(iii)$ is obtain via the contracting property of gradient descent proved in Lemma 17 and its Corollary 18.

Using the recursive relation, after summing Equation (B.1) from $1$ to $T$, we prove that the fixed-step-size full gradient method at iteration $T$ is $\frac{2\eta L^{2}T}{n}$-uniform stable, for $\eta\leq\frac{1}{\beta}$. That is, for every $z\in\mathcal{Z}$,

We remark that the stability of fixed-step-size full gradient method is linear as a function of iteration $T$. More generally, for gradient descent with varying step-sizes, using the same arguments, we can prove that the stability is upper bounded by the cumulative sum of all previous step-sizes at $T$.

Next, we show that this stability upper bound can be achieved by a linear function. We design the loss function $l(\theta;z)$ such that it is either $L\theta$ or $-L\theta$ depending on $z$. We define the two empirical loss functions on $S$ and $S^{\prime}$,

The two empirical loss functions differ exactly by $\frac{2}{n}L\theta$. We have for iteration $T$,

Then for this linear loss, for any $z\in\mathcal{Z}$,

B.2 Nesterov’s Accelerated Gradient Descent

Recall that the Nesterov’s accelerated gradient method has the following updates for $t\geq 1$:

where $\eta\leq\frac{1}{\beta}$ is the step-size. $\gamma_{t}$ is defined by the following recursion

satisfying $-1<\gamma_{t}\leq 0$. For the updates on the perturbed samples $S^{\prime}$, we have

Denote $\Delta\theta_{t}=\theta_{t}-\theta^{\prime}_{t}$. Taking the difference of Equation (18) and (19), we have

and $\theta_{\text{mid},t}$ is on the path from $\left(1-\gamma_{t-1}\right)\theta_{t}+\gamma_{t-1}\theta_{t-1}$ to $\left(1-\gamma_{t-1}\right)\theta^{\prime}_{t}+\gamma_{t-1}\theta^{\prime}_{t-1}$. Note that we have used the mean value theorem to group two gradient terms.

Because $\nabla R_{S^{\prime}}$ and $\nabla R_{S}$ only differ in one term, using the $L$-Lipschitz gradient property, we obtain an upper bound on the error term

In the case of quadratic objective, we can denote

Using the convex and $\beta$-smooth property, we have

Then we can rewrite Equation 20 as follows,

Writing this equation in matrix form, we have

We have used $\prod_{i=1}^{t}G_{i}$ to denote the matrix product $G_{t}G_{t-1}\ldots G_{1}$. The goal is to bound the norm of $\Delta\theta_{t+1}$. We need the following lemma on the spectral norm of $\prod_{i=1}^{t}G_{i}$ to conclude.

Assuming Lemma 19 as given at the moment, we now complete the proof. According to Equation (22), applying Lemma 19 to $G_{t}$, we can bound the norm of $\Delta\theta_{t+1}$,

We have used the fact that $\left\|\Delta\theta_{0}\right\|_{2}=0$, $\left\|\Delta\theta_{1}\right\|_{2}\leq\frac{2\eta L}{n}$ and $\left\|e_{t}\right\|_{2}\leq\frac{2\eta L}{n}$ in the first inequality. Together with the $L$-Lipschitz condition, we obtain that the Nesterov accelerated gradient method at iteration $T$ is

Since $B$ is symmetric positive-semidefinite, we can diagonalize $B$. There exists a common orthogonal matrix $Q$ and diagonal matrices $D$ such that

with $0\leq h\leq 1$. To bound its spectral norm, we claim the following lemma.

Suppose $H_{i}=\begin{pmatrix}(1-\gamma_{i-1})h&\gamma_{i-1}h\\ 1&0\end{pmatrix}$, where $0\leq h\leq 1$ and $-1<\gamma_{i-1}<1$. Then

Assuming Lemma 20 as given at the moment, the Lemma 19 can be completed.

Note that $\prod_{i=1}^{t}H_{i}$ is a $2\times 2$ matrix. Let $\begin{pmatrix}a_{0}\\ b_{0}\end{pmatrix}$ be a vector with norm $1$. We define

To bound the spectral norm of $\prod_{i=1}^{t}H_{i}$, it is sufficient to bound the norm of $\begin{pmatrix}a_{t}\\ b_{t}\end{pmatrix}$. We going to show by recursion that

For $t=0,t=1$, the statement is easy to verify. Suppose that the statement is true until $t$. We have the following recursion,

We remark that $a_{t+1}$ as a function of $(\gamma_{0},\ldots,\gamma_{t})$ is a multivariate polynomial with degree one. Hence its maximum or minimum value is attained at the extreme values of the variables. Formally,

This is a combinatorial optimization problem. But we observe that there are only four relevant cases.

Applying the assumption of the recursion, we obtain the desired bound for $a_{t+1}$ and $b_{t+1}$.

$\begin{pmatrix}a_{1}\\ b_{1}\end{pmatrix}$ is a vector with norm less than 1. Consider the problem with $\begin{pmatrix}a_{1}\\ b_{1}\end{pmatrix}$ as initialization, we obtain the desired bound for $a_{t+1}$ and $b_{t+1}$.

If there exists $i\in\left\{2,\ldots,t-1\right\}$ such that $\gamma_{i}=1$, then

Since $-1\leq\gamma_{i+1}\gamma_{i-1}\leq 1$, this problem is again reduced to the problem where only $t-2$ matrices are multiplied together: from $H_{t}$ to $H_{i+2}$, then $H_{i+1}H_{i}H_{i-1}$, then from $H_{i-2}$ to $H_{1}$. We apply the assumption of the recursion and obtain the desired bound for $a_{t+1}$.

Otherwise, all $\gamma_{0},...,\gamma_{t}$ should take value $-1$. Then

Let $\begin{pmatrix}H^{11}_{t}&H^{12}_{t}\\ H^{21}_{t}&H^{22}_{t}\end{pmatrix}=\prod_{i=1}^{t}H_{i}$, then we have the following recursion for its entries

We note that $H^{11}_{i}$ satisfies the following second-order recursion

with $H^{11}_{0}=1$ and $H^{11}_{0}=2h$. We observe that $H^{11}_{i}$ is exactly the Chebyshev polynomial Tchebychev (1853); Mason and Handscomb (2002) of second kind with parameter $U_{i}(h)$. It is known that for Chebyshev polynomial of second kind,

Similarly, we show that all entries are less than $t+2$. As a consequence,

This discussion of four relevant cases concludes the recursion part, and thus the proof of Lemma 20.

B.3 Heavy Ball Method with Fixed Momentum

The proof of the fixed momentum heavy ball method proceeds similarly to that of the Nesterov accelerated gradient descent.

Fixed momentum heavy ball method has the following updates.

with fixed momentum $\gamma\in[0,1)$, and fixed step-size $\eta\in\left(0,\frac{(1-\gamma)}{\beta}\right)$. For the updates on the perturbed samples $S^{\prime}$, we have

Denote $\Delta\theta_{t}=\theta_{t}-\theta_{t}^{\prime}$. Taking the difference of Equation (23) and (24), we have

and $\theta_{\text{mid},t}$ is on the path from $\theta_{t}$ to $\theta_{t}^{\prime}$. Here we have used the mean value theorem to group the two gradient terms and to make appear the Hessian terms. Using the $L$-Lipschitz property, we obtain an upper bound on the error term,

In the case of quadratic objective, we can denote

Using the convex and $\beta$-smooth property, we have

We can rewrite Equation (25) in matrix form,

As in the proof of NAG in Appendix B.2, we are going to bound the spectral norm of $\prod_{i=1}^{t}G_{i}$ to conclude. Using diagonalization of the matrices $A$, it is sufficient to consider products of the $2\times 2$ matrices $H=\begin{pmatrix}1+\gamma-a&-\gamma\\ 1&0\end{pmatrix}$, with $0\leq a\leq\eta\beta$. The following lemma characterizes the spectral norm of $\prod_{i=1}^{t}H$.

Suppose $H=\begin{pmatrix}1+\gamma-a&-\gamma\\ 1&0\end{pmatrix}$, where $0<\gamma<1$ and $0\leq a\leq 1-\gamma$. Then

Assuming Lemma 21 as given at the moment, we have

We have used the fact that $\left\|\Delta\theta_{0}\right\|_{2}=0$, $\left\|\Delta\theta_{1}\right\|_{2}\leq\frac{2\eta L}{n}$ and $\left\|e_{t}\right\|_{2}\leq\frac{2\eta L}{n}$ in the first inequality. Together with the $L$-Lipschitz condition, we obtain that the heavy ball method with fixed momentum at iteration $T$ is

Let $\prod_{i=1}^{t}H=\begin{pmatrix}a_{t}&b_{t}\\ c_{t}&d_{t}\end{pmatrix}$. We are going to show by recursion that

For $t=0,1$, the statement is easy to verify. Suppose that the statement is true until $t$. We have by recursion formular

with initialization $a_{1}=1+\gamma-a,c_{1}=1,b_{1}=-\gamma,d_{1}=0$. We remark that $a_{i}$ satisfies the following second-order recursion, for $i\geq 1$,

where $a_{0}=1,a_{1}=1+\gamma-a$. We can also add $a_{-1}=0$.

. We distinguish two cases based on the two roots.

The two roots are distinct. By distinct roots theorem for second order homogeneous system, we have

where $l_{1}$ and $l_{2}$ are constants to be determined by the initial condition. Solving the initial condtion, we have

We have used that $\left|x_{2}\right|\leq\sqrt{\gamma}$. When the two roots have imaginary part, it is clear that $\left|x_{2}\right|=\sqrt{\gamma}$. On the other hand, when the two roots are real, since $\left|x_{1}x_{2}\right|=\gamma$, $\left|x_{2}\right|\leq\left|x_{1}\right|$, we also have $\left|x_{2}\right|\leq\sqrt{\gamma}$.

The two roots are equal. $1+\gamma-a=2\sqrt{\gamma}$.

By single root theorem for second order homogeneous system, We have

Overall, we have proved a bound for $a_{t}$,

We can bound $b_{t},c_{t}$ and $d_{t}$ similarly because they have similar recursion formular.

Using the relationship between spectral norm and Frobenius norm, we have

C Stability Bounds for Strongly Convex Smooth Functions

Recall that in order to prove the uniform stability, we need bound the loss difference for any fixed sample $z$ at each iteration $t\geq 1$

This quantity is related to the norm difference $\left\|\theta_{t}-\theta^{\prime}_{t}\right\|_{2}$ under the $L$-Lipschitz condition. Under $\alpha$-strongly-convex case, we bound $\left\|\theta_{t}-\theta^{\prime}_{t}\right\|_{2}$ slightly different than that in the convex smooth case.

Using the update rule of full gradient method, we obtain an recursive relation on $\left\|\theta_{t}-\theta^{\prime}_{t}\right\|_{2}$. For $\eta\leq\frac{2}{\alpha+\beta}$ and $t\geq 1$, we have

The inequality $(i)$ uses triangular inequality. The inequality $(ii)$ follows from the $L$-Lipschitz condition on the perturbed gradient terms. The inequality $(iii)$ is obtain via the following claim, for $f$ $\alpha$-strongly convex and $\beta$-smooth, we have

This claim can be easily obtain by plugging $f(x)-\frac{\alpha}{2}\left\|x\right\|_{2}^{2}$, which is a convex function into Corollary 18. The inequality $(iv)$ uses the fact $(1-x)^{1/2}\leq 1-x^{1/2}$, for $0\leq x\leq 1$.

Using the recursive relation, after summing Equation (C.1) from $1$ to $T$, we have

Applying the $L$-Lipschitz condition, we have for every $z\in\mathcal{Z}$,

C.2 Nesterov’s Accelerated Gradient Descent

According to the discussion of Equation 21, in the case of quadratic loss, the Nesterov accelerated gradient descent difference term is as follows

As in the proof of NAG in Appendix B.2, we are going to bound the spectral norm of $\prod_{i=1}^{t}G_{i}$ to conclude. Following the proof idea used in Appendix B.2 and Appendix B.3, using diagonalization of the matrices $A$, it is sufficient to consider products of the $2\times 2$ matrices $H=\begin{pmatrix}(1+\gamma)h&-\gamma h\\ 1&0\end{pmatrix}$, with $1-\beta\eta\leq h\leq 1-\alpha\eta$. The following lemma characterizes the spectral norm of $\prod_{i=1}^{t}H$.

Suppose $H=\begin{pmatrix}(1+\gamma)h&-\gamma h\\ 1&0\end{pmatrix}$, where $\gamma=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$ and $1-\beta\eta\leq h\leq 1-\alpha\eta$. Then

Assuming Lemma 22 as given at the moment, we have

We have used the fact that $\left\|\Delta\theta_{0}\right\|_{2}=0$, $\left\|\Delta\theta_{1}\right\|_{2}\leq\frac{2\eta L}{n}$ and $\left\|e_{t}\right\|_{2}\leq\frac{2\eta L}{n}$ in the first inequality. Let $p=\left(\gamma(1-\alpha\eta)\right)^{1/2}$ and

We also have upper and lower bounds on $p$,

Together with the $L$-Lipschitz condition, we obtain that the heavy ball method with fixed momentum at iteration $T$ is

Let $\prod_{i=1}^{t}H=\begin{pmatrix}a_{t}&b_{t}\\ c_{t}&d_{t}\end{pmatrix}$. We are going to show by recursion that

For $t=0,1$, the statement is easy to verify. Suppose that the statement is true until $t$. We have by recursion formular

with initialization $a_{1}=(1+\gamma)h,b_{1}=-\gamma h,c_{1}=1$ and $d_{1}=0$. We remark that $a_{i}$, satisfies the following second-order recursion, for $i\geq 1$,

where $a_{0}=1,a_{1}=(1+\gamma)h$. We can also add $a_{-1}=0$.

because $h\leq 1-\alpha\eta\leq\frac{\kappa-1}{\kappa}$. Hence either we have equal real roots, or we have complex roots with imaginary parts.

The two roots are equal. $(1+\gamma)h=2\sqrt{\gamma h}$. Then

By single root theorem for second order homogeneous system, we have

By distinct roots theorem for second order homogeneous system, we have

where $l_{1}$ and $l_{2}$ are constants to be determined by the initial condition. Solving the initial condtion, we have

We can bound $b_{t},c_{t}$ and $d_{t}$ similarly because they have similar recursion formular.

Using the relationship between spectral norm and Frobenius norm, we have