Estimate Sequences for Stochastic Composite Optimization: Variance Reduction, Acceleration, and Robustness to Noise

Introduction

We consider convex composite optimization problems of the form

More specifically, we focus on stochastic objective functions, which are of utmost importance in machine learning, where $f$ is an expectation or a finite sum of smooth functions

While the finite-sum setting is obviously a particular case of expectation with a discrete probability distribution, the deterministic nature of the resulting cost function drastically changes performance guarantees. In particular, when an algorithm is only allowed to access unbiased measurements of the objective function and gradient—which we assume is the case when $f$ is an expectation—it may be shown that the worst-case convergence rate in expected function value cannot be better than $O(1/k)$ in general, where $k$ is the number of iterations Nemirovski et al. (2009); Agarwal et al. (2012). Such a sublinear rate of convergence is notably achieved by stochastic gradient descent (SGD) algorithms or their variants (see Bottou et al., 2018).

Even though this pessimistic result applies to the general stochastic case, linear convergence rates can be obtained for the finite-sum setting Schmidt et al. (2017). Specifically, a large body of work in machine learning has led to many randomized incremental approaches obtaining linear convergence rates, such as SAG Schmidt et al. (2017), SAGA Defazio et al. (2014a), SVRG Johnson and Zhang (2013); Xiao and Zhang (2014), SDCA Shalev-Shwartz and Zhang (2016), MISO Mairal (2015), Katyusha Allen-Zhu (2017), MiG Zhou et al. (2018), SARAH Nguyen et al. (2017), directly accelerated SAGA Zhou (2019) or the method of Lan and Zhou (2018a). For non-convex objectives, recent approaches have also improved known convergence rates for finding first-order stationary points Fang et al. (2018); Paquette et al. (2018); Lei et al. (2017), which is however beyond the scope of our paper. These algorithms have about the same cost per-iteration as the stochastic gradient descent method, since they access only a single or two gradients $\nabla f_{i}(x)$ at each iteration, and they may achieve lower computational complexity than accelerated gradient descent methods Nesterov (1983, 2004, 2013); Beck and Teboulle (2009) in expectation. A common interpretation is to see these algorithms as performing SGD steps with an estimate of the full gradient that has lower variance Xiao and Zhang (2014).

In this paper, we are interested in providing a unified view of stochastic optimization algorithms, but we also want to investigate their robustness to random perturbations. Specifically, we may consider objective functions with an explicit finite-sum structure such as (2) when only noisy estimates of the gradients $\nabla f_{i}(x)$ are available. Such a setting may occur for various reasons. For instance, perturbations may be injected during training in order to achieve better generalization on new test data Srivastava et al. (2014), perform stable feature selection (Meinshausen and Bühlmann, 2010), improve the model robustness Zheng et al. (2016), or for privacy-aware learning (Wainwright et al., 2012).

Each training point indexed by $i$ is corrupted by a random perturbation $\rho_{i}$ and the resulting function $f$ may be written as

The framework we adopt is that of estimate sequences introduced by Nesterov (2004), which consists of building iteratively a quadratic model of the objective. Typically, estimate sequences may be used to analyze the convergence of existing algorithms, but also to design new ones, in particular with acceleration. Our construction is however slightly different than the original one since it is based on stochastic estimates of the gradients, and some classical properties of estimate sequences are satisfied only approximately. We note that estimate sequences have been used before for stochastic optimization (Lu and Xiao, 2015; Devolder, 2011; Lin et al., 2014), but not for the same generic purpose as ours.

Specifically, our paper makes to the following contributions:

We revisit many stochastic optimization algorithms dealing with composite convex problems; we consider variants of incremental methods such as SVRG, SAGA, SDCA, or MISO. We provide a common convergence proof for these methods and show that they can be modified and become adaptive to the strong convexity constant $\mu$, when only a lower bound is available.

We show that our construction of estimate sequence naturally leads to an accelerated stochastic gradient method for composite optimization as Ghadimi and Lan (2012, 2013); Hu et al. (2009), but simpler as our approach requires to maintain two sequences of iterates instead of three. The resulting complexity in terms of gradient evaluations for $\mu$-strongly convex objectives is

which has also been achieved by Ghadimi and Lan (2013); Aybat et al. (2019); Cohen et al. (2018). When the objective is convex, but non-strongly convex, we also provide a sublinear convergence rate for finite horizon. Given a budget of $K$ iterations, the algorithm returns an iterate $x_{K}$ such that

which is also optimal for stochastic first-order optimization (Ghadimi and Lan, 2012).

We design a new accelerated algorithm for finite sums based on the SVRG gradient estimator, with complexity, for $\mu$-strongly convex functions,

When the problem is convex but not strongly convex, given a budget of $K$ greater than $O(n\log(n))$, the algorithm returns a solution $x_{K}$ such that

This paper is organized as follows. Section 2 introduces the proposed framework based on stochastic estimate sequences; Section 3 is devoted to the convergence analysis and Section 4 introduces accelerated stochastic optimization algorithms; Section 5 presents various experiments to compare the effectiveness of the proposed approaches, and Section 6 concludes the paper.

We note that a short version of this paper was presented by Kulunchakov and Mairal (2019a) at the International Conference on Machine Learning (ICML) in 2019. This paper extends this previous work by (i) providing complexity results for convex but not strongly convex objectives ($\mu=0$), (ii) extending the framework to variants of MISO/Finito/SDCA algorithms, in the context of non-accelerated incremental methods, (iii) providing more experiments with additional baselines and objective functions.

Proposed Framework Based on Stochastic Estimate Sequences

In this section, we present two generic stochastic optimization algorithms to address the composite problem (1). Then, we show their relation to variance-reduction methods.

The iteration (A) is generic and encompasses many existing algorithms, which we review later. Key to our analysis, we are interested in a simple interpretation corresponding to the iterative minimization of strongly convex surrogate functions.

with $\gamma_{0}\geq\mu$ and $d_{0}^{*}$ is a scalar value that is left unspecified at the moment. Then, it is easy to show that $x_{k}$ in (A) minimizes the following quadratic function $d_{k}$ defined for $k\geq 1$ as

where $\delta_{k},\gamma_{k}$ satisfy the system of equations

We note that $\psi^{\prime}(x_{k})$ is a subgradient in $\partial\psi(x_{k})$. By simply using the definition of the proximal operator (7) and considering first-order optimality conditions, we indeed have that $0\in x_{k}-x_{{k\text{--}1}}+\eta_{k}g_{k}+\eta_{k}\partial\psi(x_{k})$ and $x_{k}$ coincides with the minimizer of $d_{k}$. This allows us to write $d_{k}$ in the generic form

The construction (9) is akin to that of estimate sequences introduced by Nesterov (2004), which are typically used for designing accelerated gradient-based optimization algorithms. In this section, we are however not interested in acceleration, but instead in stochastic optimization and variance reduction. One of the main property of estimate sequences that we will use is their ability do behave asymptotically as a lower bound of the objective function near the optimum. Indeed, we have

Relation with existing algorithms.

The iteration (A) encompasses many approaches such as ISTA (proximal gradient descent), which uses the exact gradient $g_{k}=\nabla f(x_{{k\text{--}1}})$ leading to deterministic iterates $(x_{k})_{k\geq 0}$ (Beck and Teboulle, 2009; Nesterov, 2013) or proximal variants of the stochastic gradient descent method to deal with a composite objective (see Lan, 2012, for instance). Of interest for us, the variance-reduced stochastic optimization approaches SVRG (Xiao and Zhang, 2014) and SAGA Defazio et al. (2014a) also follow the iteration (A) but with an unbiased gradient estimator whose variance reduces over time. Specifically, the basic form of these estimators is

2 A Less Classical Iteration with a Different Estimate Sequence

In the previous section, we have interpreted the classical iteration (A) as the iterative minimization of the stochastic surrogate (9). Here, we show that a slightly different construction leads to a new algorithm. To obtain a lower bound, we have indeed used basic properties of the proximal operator to obtain a subgradient $\psi^{\prime}(x_{k})$ and we have exploited the following convexity inequality $\psi(x)\geq\psi(x_{k})+\psi^{\prime}(x_{k})^{\top}(x-x_{k})$. Another natural choice to build a lower bound consists then of using directly $\psi(x)$ instead of $\psi(x_{k})+\psi^{\prime}(x_{k})^{\top}(x-x_{k})$, leading to the construction

where $x_{k\text{--}1}$ is assumed to be the minimizer of the composite function $d_{k\text{--}1}$, $\delta_{k}$ is defined as in Section 2.1, and $x_{k}$ is a minimizer of $d_{k}$. To initialize the recursion, we define then $d_{0}$ as

with $x_{0}=\text{Prox}_{\psi/\gamma_{0}}[\bar{x}_{0}]$ is the minimizer of $d_{0}$ and $d_{0}^{*}=d_{0}(x_{0})=c_{0}+\frac{\gamma_{0}}{2}\|x_{0}-\bar{x}_{0}\|^{2}+\psi(x_{0})$ is the minimum value of $d_{0}$; $c_{0}$ is left unspecified since it does not affect the algorithm. Typically, one may choose $\bar{x}_{0}$ to be a minimizer of $\psi$ such that $x_{0}=\bar{x}_{0}$. Unlike in the previous section, the surrogates $d_{k}$ are not quadratic, but they remain $\gamma_{k}$-strongly convex. It is also easy to check that the relation (11) still holds.

where $c_{k}$ is constant and the inequality on the right is due to the strong convexity of $d_{k}$.

Relation to existing approaches.

The approach (B) is related to several optimization methods. When the objective is a deterministic finite sum, it is possible to relate the update (B) to the MISO Mairal (2015), and Finito Defazio et al. (2014b) algorithms, even though they were derived from a significantly different point of view. This is also the case of a primal variant of SDCA (Shalev-Shwartz, 2016) For instance, SDCA is a dual coordinate ascent approach, whereas MISO and Finito are explicitly derived from the iterative surrogate minimization we adopt in this paper. As the links between (B) and these previous approaches are not obvious at first sight, we detail them in Appendix B.

3 Gradient Estimators and Algorithms

In this paper, we consider the iterations (A) and (B) with the following gradient estimators.

exact gradient with $g_{k}=\nabla f(x_{k\text{--}1})$, when the problem is deterministic and we have an access to the full gradient;

Then, the variables $z_{k}^{i}$ and $\bar{z}_{k}$ also involve noisy estimates of the gradients:

The case with non-uniform sampling is slightly different and is described in Algorithm 3; it requires an additional index $j_{k}$ for updating a variable $z_{k}^{j_{k}}$. The reason for that is to remove a difficulty in the convergence proof, a strategy also adopted by Schmidt et al. (2015) for a variant of SAGA with non-uniform sampling.

SDCA/MISO: To put SAGA, MISO and SDCA under the same umbrella, we introduce a lower bound $\beta$ on the strong convexity constant $\mu$, and a correcting term involving $\beta$ that appears only when the sampling distribution $Q$ is not uniform:

It is then possible to show that when $Q$ is uniform and under the big data condition $L/\mu\leq n$ (used for instanced by Mairal 2015; Defazio et al. 2014b; Schmidt et al. 2017) and with $\beta=\mu$, variant (B) combined with the estimator (15) yields the MISO algorithm, which performs similar updates as a primal variant of SDCA (Shalev-Shwartz, 2016). These links are highlighted in Appendix B.

As we combine different types of iterations and gradient estimators, we recover both known and new algorithms. Specifically, we obtain the following new features:

robustness to noise: we introduce mechanisms to deal with stochastic perturbations and make all these previous approaches robust to noise.

new variants: Whereas SVRG/SAGA were developed with the iterations (A) and MISO in the context of (B), we show that these gradient estimators are both compatible with (A) and (B), leading to new algorithms with similar guarantees.

Convergence Analysis and Robustness

We now present the convergence analysis for iterations (A) or (B). In Section 3.1, we present a generic convergence result. Then, in Section 3.2, we present specific results for the variance-reduction approaches in including strategies to make them robust to stochastic noise. Acceleration is discussed in the next section.

Key to our complexity results, the following proposition gives a first relation between the quantity $F(x_{k})$, the surrogate $d_{k}$, $d_{k\text{--}1}$ and the variance of the gradient estimates.

For either variant (A) or (B), when using the construction of $d_{k}$ from Sections 2.1 or 2.2, respectively, and assuming $\eta_{k}\leq 1/L$, we have for all $k\geq 1$,

Proof We first consider the variant (A) and later show how to modify the convergence proofs to accommodate the variant (B).

where the first inequality comes from Lemma 24—it is in fact an equality when considering Algorithm (A)—and the second inequality simply uses the assumption $\eta_{k}\leq 1/L$, which yields $\delta_{k}=\gamma_{k}\eta_{k}\leq\gamma_{k}/L$. Finally, the last inequality uses a classical upper-bound for $L$-smooth functions presented in Lemma 22. Then, after taking expectations,

where we have defined the following quantity

This relation can now be combined with (11) when $z=x^{*}$, and we obtain (16). It is also easy to see that the proof also works with variant (B). The convergence analysis is identical, except that we take $w_{k\text{--}1}$ to be

and the same result follows. Then, without making further assumption on $\omega_{k}$, we have the following general convergence result, which is a direct consequence of the averaging Lemma 30, inspired by Ghadimi and Lan (2012), and presented in Appendix A.3:

Under the same assumptions as in Proposition 1, we have for all $k\geq 1$, and either variant (A) or (B),

where $\Gamma_{k}=\prod_{t=1}^{k}(1-\delta_{t})$. Then, by using the averaging strategy $\hat{x}_{k}=(1-\delta_{k})\hat{x}_{{k\text{--}1}}+\delta_{k}x_{k}$ of Lemma 30, for any point $\hat{x}_{0}$ (possibly equal to $x_{0}$), we have

Theorem 2 allows us to recover convergence rates for various algorithms. Note that the effect of the averaging strategy is to remove the factor $\delta_{k}$ in front of $F(x_{k})-F^{*}$ on the left part of (17), thus improving the convergence rate by a factor $1/\delta_{k}$. Regarding the quantity $d_{0}(x^{*})-d_{0}^{*}$, we have the following relations

For variant (A), $d_{0}(x^{*})-d_{0}^{*}=\frac{\gamma_{0}}{2}\|x^{*}-x_{0}\|^{2}$;

For variant (B), this quantity may be larger and we may simply say that $d_{0}(x^{*})-d_{0}^{*}=\frac{\gamma_{0}}{2}\|x^{*}-x_{0}\|^{2}+\psi(x^{*})-\psi(x_{0})-\psi^{\prime}(x_{0})^{\top}(x_{0}-x^{*})$ for variant (B), where $\psi^{\prime}(x_{0})=\gamma_{0}(x_{0}-\bar{x}_{0})$ is a subgradient in $\partial\psi(x_{0})$. Note that if $\bar{x}_{0}$ is chosen to be a minimizer of $\psi$, then $d_{0}(x^{*})-d^{*}=\frac{\gamma_{0}}{2}\|x^{*}-x_{0}\|^{2}+\psi(x^{*})-\psi(x_{0})$.

In the next section, we will focus on variance reduction mechanisms, which are able to improve the previous convergence rates by better exploiting the structure of the objective. By controlling the variance $\omega_{k}$ of the corresponding gradient estimators, we will apply Theorem 2 to obtain convergence rates. Before that, we remark that it is relatively straightforward to use this theorem to recover complexity results for proximal SGD, both for the usual variant (A) or the new one (B). Since these results are classical, we present them in Appendix C. As a sanity check, we note that we recover the optimal noise-dependency (see Nemirovski et al., 2009), both for strongly convex cases, or when $\mu=0$.

2 Faster Convergence with Variance Reduction

Stochastic variance-reduced gradient descent algorithms rely on gradient estimates whose variance decreases as fast as the objective function value. Here, we provide a unified proof of convergence for our variants of SVRG, SAGA, and MISO, and we show how to make them robust to stochastic perturbations. Specifically, we consider the minimization of a finite sum of functions as in (3), but, as explained in Section 2, each observation of the gradient $\nabla f_{i}(x)$ is corrupted by a random noise variable. The next proposition extends a proof for SVRG Xiao and Zhang (2014) to stochastic perturbations, and characterizes the variance of $g_{k}$.

where the expectation is with respect to the gradient perturbation, and $Q=\{q_{1},\ldots,q_{n}\}$ is the sampling distribution. As having the variance to be bounded across the domain of $x$ may be a strong assumption, even though classical, we also introduce the quantity

Consider problem (1) when $f$ is a finite sum of functions $f=\frac{1}{n}\sum_{i=1}^{n}f_{i}$ where each $f_{i}$ is convex and $L_{i}$-smooth with $L_{i}\geq\mu$. Then, the gradient estimates $g_{k}$ of the random-SVRG and MISO/SAGA/SDCA strategies defined in Section 2.3 satisfy

where $L_{Q}=\max_{i}L_{i}/(q_{i}n)$, $\rho_{Q}=1/(n\min_{i}q_{i})$, and for all $i$ and $k$, $u_{k}^{i}$ is equal to $z_{k}^{i}$ without noise—that is

and $u^{i}_{*}=\nabla f_{i}(x^{*})-\beta x^{*}$ (with $\beta=0$ for random-SVRG).

Consider the same setting as Proposition 3. For either variant (A) or (B) with the random-SVRG or SAGA/SDCA/MISO gradient estimators defined in Section 2.3, when using the construction of $d_{k}$ from Sections 2.1 or 2.2, respectively, and assuming $\gamma_{0}\geq\mu$ and $(\eta_{k})_{k\geq 0}$ is non-increasing with $\eta_{k}\leq\frac{1}{12L_{Q}}$, we have for all $k\geq 1$,

The proof of the previous proposition is given in Appendix D.2. From the Lyapunov function, we obtain a general convergence result for the variance-reduced stochastic algorithms.

Consider the same setting as Proposition 4, which applies to both variants (A) and (B). Then, by using the averaging strategy of Lemma 30 with any point $\hat{x}_{0}$,

where $\Theta_{k}=\prod_{t=1}^{k}(1-\tau_{t})$. Note that we also have

The proof is given in Appendix D.3. From this generic convergence theorem, we now study particular cases. The first corollary studies the strongly-convex case with constant step size.

Consider the same setting as in Theorem 5, where $f$ is $\mu$-strongly convex, $\gamma_{0}=\mu$, and $\eta_{k}=\frac{1}{12L_{Q}}$. Then, for any point $\hat{x}_{0}$,

with $\tau=\min\left(\frac{\mu}{12L_{Q}},\frac{1}{5n}\right)$, $\Theta_{k}=(1-\tau)^{k}$, and $\alpha=6\min\left(1,\frac{12L_{Q}}{5\mu n}\right)$. Note that $T_{k}\geq\frac{\mu}{2}\|x_{k}-x^{*}\|^{2}$ and for Algorithm (A), we also have $T_{0}\leq(13/12)(F(x_{0})-F^{*})$.

Consider the same setting as Theorem 5, where $f$ is $\mu$-strongly convex, $\gamma_{0}=\mu$, and $\eta_{k}=\eta=\min\left(\frac{1}{12L_{Q}},\frac{1}{5\mu n}\right)$. Then, for all $\hat{x}_{0}$,

The proof follows similar steps as the proof of Corollary 6, after noting that we have $\delta_{k}=\tau_{k}$ for all $k$ for this particular choice of step size. We are now in shape to study a converging algorithm.

Note that $d_{0}(x^{*})-d_{0}^{*}\leq F(x_{0})-F^{*}$ for variant (A).

Accelerated Stochastic Algorithms

We now consider the following iteration, involving an extrapolation sequence $(y_{k})_{k\geq 1}$, which is a classical mechanism from accelerated first-order algorithms Beck and Teboulle (2009); Nesterov (2013). Given a sequence of step-sizes $(\eta_{k})_{k\geq 0}$ with $\eta_{k}\leq 1/L$ for all $k\geq 0$, and some parameter $\gamma_{0}\geq\mu$, we consider the sequences $(\delta_{k})_{k\geq 0}$ and $(\gamma_{k})_{k\geq 0}$ that satisfy

Consider then the stochastic estimate sequence for $k\geq 1$

and $\psi^{\prime}(x_{k})=\frac{1}{\eta_{k}}(y_{k\text{--}1}-x_{k})-g_{k}$ is in $\partial\psi(x_{k})$ by definition of the proximal operator. As in Section 2, $d_{k}(x^{*})$ asymptotically becomes a lower bound on $F^{*}$ since (11) remains satisfied. This time, the iterate $x_{k}$ does not minimize $d_{k}$, and we denote by $v_{k}$ instead its minimizer, allowing us to write $d_{k}$ in the canonical form

The first lemma highlights classical relations between the iterates $(x_{k})_{k\geq 0}$, $(y_{k})_{k\geq 0}$ and the minimizers of the estimate sequences $d_{k}$, which also appears in (Nesterov, 2004, p. 78) for constant step sizes $\eta_{k}$. The proof is given in Appendix D.5.

The sequences $(x_{k})_{k\geq 0}$ and $(y_{k})_{k\geq 0}$ produced by Algorithm (C) satisfy for all $k\geq 0$, with $v_{0}=y_{0}=x_{0}$,

Then, the next lemma is key to prove the convergence of Algorithm (C). Its proof is given in Appendix D.8.

Assuming $(x_{k})_{k\geq 0}$ and $(y_{k})_{k\geq 0}$ are given by Algorithm (C). Then, for all $k\geq 1$,

Finally, we obtain the following convergence result.

Under the assumptions of Lemma 10, we have for all $k\geq 1$,

where, as before, $\Gamma_{t}=\sum_{i=1}^{t}(1-\delta_{i})$.

Proof First, the minimizer $v_{k}$ of the quadratic surrogate $d_{k}$ may be written as

Then, we characterize the quantity $d_{k}^{*}$:

By Lemma 10, we can show that the last term is equal to zero, and we are left with

where we used the fact that $\eta_{k}\leq 1/L$ and $\delta_{k}=\sqrt{\gamma_{k}\eta_{k}}$.

It remains to choose $d_{0}^{*}=F(x_{0})$ and $\xi_{0}=0$ to initialize the induction at $k=0$ and we conclude that

which gives us (30) when noticing that $\xi_{k}=\Gamma_{k}\sum_{t=1}^{k}\frac{\eta_{t}\omega_{t}^{2}}{\Gamma_{t}}$.

Next, we specialize the theorem to various practical cases. For the corollaries below, we assume the variances $\left(\omega_{k}^{2}\right)_{k\geq 1}$ to be upper bounded by $\sigma^{2}$.

Assume that $f$ is $\mu$-strongly convex, and choose $\gamma_{0}=\mu$ and $\eta_{k}=1/L$ with Algorithm (C). Then,

We now show that with decreasing step sizes, we obtain an algorithm with optimal complexity similar to (Ghadimi and Lan, 2013).

The proof is provided in Appendix D.9. We note that despite the “optimal” theoretical complexity, we have observed that Algorithm (C) with the parameters of Corollaries 13 and 14 could be relatively unstable, as shown in Section 5, due to the large radius $\sigma^{2}/\sqrt{\mu L}$ of the noise region. When $\mu$ is small, such a quantity may be indeed arbitrarily larger than $F(x_{0})-F^{*}$. Instead, we have found a minibatch strategy to be more effective in practice. When using a minibatch of size $b=\lceil L/\mu\rceil$, the theoretical complexity becomes the same as SGD, given in Corollary 32, but the algorithm enjoys the benefits of easy parallelization.

Assume that $f$ is convex. Consider a step-size $\eta\leq 1/L$ and run one iteration of Algorithm (A) with a stochastic gradient estimate. Use the resulting point to initialize Algorithm (C) still with constant step size $\eta$, and choose $\gamma_{0}=1/\eta$. Then,

If in addition we choose $\eta=\min\left(\frac{1}{L},\sqrt{\frac{2\|x_{0}-x^{*}\|^{2}}{\sigma^{2}}}\frac{1}{(K+1)^{3/2}}\right)$, then

The proof is given in Appendix D.10. These convergence results are relatively similar to those obtained in Ghadimi and Lan (2013) for a different algorithm and is optimal for convex functions.

2 An Accelerated Algorithm with Variance Reduction

In this section, we show how to combine the previous methodology with variance reduction, and introduce Algorithm 4 based on random-SVRG. Then, we present the convergence analysis, which requires controlling the variance of the estimator in a similar manner to Allen-Zhu (2017), as stated in the next proposition. Note that the estimator does not require storing the seed of the random perturbations, unlike in the previous section.

Consider problem (1) when $f$ is a finite sum of functions $f=\frac{1}{n}\sum_{i=1}^{n}f_{i}$ where each $f_{i}$ is $L_{i}$-smooth with $L_{i}\geq\mu$ and $f$ is $\mu$-strongly convex. Then, the variance of $g_{k}$ defined in Algorithm 4 satisfies

The proof is given in Appendix D.11. Then, we extend Lemma 11 that was used in the previous analysis to the variance-reduction setting.

Consider the iterates provided by Algorithm 4 and call $a_{k}=2L_{Q}\eta_{k}$. Then,

The proof of this lemma is given in Appendix D.12. With this lemma in hand, we may now state our main convergence result.

Consider the iterates provided by Algorithm 4 and assume that the step sizes satisfy $\eta_{k}\leq\min\left(\frac{1}{3L_{Q}},\frac{1}{15\gamma_{k}n}\right)$ for all $k\geq 1$. Then,

Proof Following similar steps as in the proof of Theorem 12, we have

and we obtain (33). We may now derive convergence rates of our accelerated SVRG algorithm under various settings. The proofs of the following corollaries, when not straightforward, are given in the appendix. The first corollary simply uses Lemma 27.

With $\eta_{k}=\min\left(\frac{1}{3L_{Q}},\frac{1}{15\mu n}\right)$ and $\gamma_{0}=\mu$, the iterates produced by Algorithm 4 satisfy

if $\frac{1}{3L_{Q}}\leq\frac{1}{15\mu n}$,

The proof is given in Appendix D.13. Next, we study the case when $\mu=0$.

Experiments

In this section, we evaluate numerically the approaches introduced in the previous sections.

The scalar $\lambda$ is a regularization parameter that acts as a lower bound on the strong convexity constant of the problem. We consider the parameters $\mu=\lambda=1/10n$ in our problems, which is of the order of the smallest values that one would try when doing a parameter search, e.g., by cross-validation. For instance, this is empirically observed for the dataset cifar-ckn described below, where a test set is available, allowing us to check that the “optimal” regularization parameter leading to the lowest generalization error is indeed of this order. We also report an experiment with $\lambda=1/100n$ in order to study the effect of the problem conditioning on the method’s performance.

Following Bietti and Mairal (2017); Zheng and Kwok (2018), we consider DropOut perturbations (Srivastava et al., 2014) to illustrate the robustness to noise of the algorithms. DropOut consists of randomly setting to zero each entry of a data point with probability $\delta$, leading to the optimization problem

where $\rho$ is a binary vector in $\{0,1\}^{p}$ with i.i.d. Bernoulli entries, and $\circ$ denotes the elementwise multiplication between two vectors. We consider two DropOut regimes, with $\delta$ in $\{0.01,0.1\}$, representing small and medium perturbations, respectively.

We consider the following three datasets coming from different scientific fields

alpha is from the Pascal Large Scale Learning Challenge websitehttp://largescale.ml.tu-berlin.de/ and contains $n=250\,000$ points in dimension $p=500$.

gene consists of gene expression data and the binary labels $b_{i}$ characterize two different types of breast cancer. This is a small dataset with $n=295$ and $p=8\,141$.

ckn-cifar is an image classification task where each image from the CIFAR-10 datasethttps://www.cs.toronto.edu/~kriz/cifar.html is represented by using a two-layer unsupervised convolutional neural network (Mairal, 2016). Since CIFAR-10 originally contains 10 different classes, we consider the binary classification task consisting of predicting the class 1 vs. other classes. The dataset contains $n=50\,000$ images and the dimension of the representation is $p=9\,216$.

2 Evaluation of Algorithms without Perturbations

From these experiments, we obtain the following conclusions:

Acceleration for SVRG is effective on the datasets gene and ckn-cifar except on alpha, where all SVRG-like methods perform already well. This may be due to strong convexity hidden in alpha leading to a regime where acceleration does not occur—that is, when the complexity is $O(n\log(1/\varepsilon))$, which is independent of the condition number. Note that this algorithm is now implemented in the open-source Cyanure toolbox Mairal (2019)http://julien.mairal.org/cyanure/.

Acceleration is more effective when the problem is badly conditioned. When $\lambda=1/100n$, acceleration brings several orders of magnitude improvement in complexity.

Accelerated SGD is unstable with the squared hinge loss. During the initial phase with constant step size $1/L$, the expected primal gap is in a region of radius $O(\sigma^{2}/\sqrt{\mu L})\approx\sqrt{n}\sigma^{2}$, which is potentially huge, causing large gradients and instabilities.

Accelerated minibatch SGD performs best among the SGD methods and is competitive with SVRG in the low precision regime. The performance of Adam on these datasets is inconsistent; it performs best among SGD methods on alpha, but is significantly worse on ckn-cifar. Note also that AC-SA performs in general similarly to acc-SGD-d.

3 Evaluation of Algorithms with Perturbations

We now consider the same setting as in the previous section, but we add DropOut perturbations with rate $\delta$ in $\{0.01,0.1\}$. As predicted by theory, all approaches with constant step size do not converge. Therefore, we only report the results for decreasing step sizes in Figures 4, 5, and 6. We evaluate the loss function every $5$ data passes and we estimate the expectation (36) by drawing $5$ random perturbations per data point, resulting in $5n$ samples. The optimal value $F^{*}$ is estimated by letting the methods run for $1000$ epochs and selecting the best point found as a proxy of $F^{*}$.

The conclusions of these experiments are the following:

accelerated minibatch SGD performs the best among SGD approaches in general except on alpha where Adam performs best.

accelerated SVRG performs better than SVRG in general, or they achieve the same performance. As in the deterministic case, the gains are typically more important in ill-conditioned cases.

accelerated SVRG performs better than SGD approaches in the low perturbation regime $\delta=0.01$ and only on the alpha dataset when $\delta=0.1$. Otherwise, the methods perform similarly.

Discussion

Our work is based on stochastic variants of estimate sequences introduced by Nesterov (1983, 2004). The framework leads naturally to many algorithms with relatively generic proofs of convergence, where convergence is proven at the same time as the algorithm’s design. With iterate averaging techniques inspired by Ghadimi and Lan (2013), we show that a large class of variance-reduction stochastic optimization methods can be made robust to stochastic perturbations. Estimate sequences also naturally lead to several accelerated algorithms, some of them we did not present in this paper. For instance, it is possible to show that replacing in (29) the lower bound $\psi(x_{k})+\psi^{\prime}(x_{k})^{\top}(x-x_{k})$ by $\psi(x)$ itself—in a similar way as we proceeded to obtain iteration (B) from iteration (A)—also leads to an accelerated algorithm with similar guarantees as (C).

Possibilities offered by estimate sequences are large, but our framework also admits a few limitations, paving the way for future work. In particular, our results are currently limited to Euclidean metrics—meaning that our convergence rates typically depend on quantities involving the Euclidean norm (e.g., strong convexity or $L$-smooth inequalities), and one may expect extensions of our work to other metrics such as Bregman distances. Estimate sequences admit indeed known extensions to such metrics, and can also deal with higher-order smoothness assumptions than Lipschitz continuity of the gradient (Baes, 2009)—e.g., cubic regularization (Nesterov and Polyak, 2006). We leave such directions for the future.

Another limitation we encountered was the inability to propose robust accelerated variants of SAGA, MISO, or SDCA based on our stochastic estimate sequences framework. To address this problem, after the first version of this manuscript was made publicly available, we investigated in (Kulunchakov and Mairal, 2019b) a significantly different approach based on the Catalyst method (Lin et al., 2018), allowing us to accelerate stochastic first-order methods in a generic fashion, at the price of a logarithmic factor in the optimal complexity—in other words, we were able to obtain for SAGA, MISO, and SDCA a complexity close to (5) up to a logarithmic factor in the condition number $L_{Q}/\mu$. We believe that estimate sequences may be useful to obtain the optimal complexity without this logarithmic term, but the construction would be non-trivial and would rely on a different lower bound than the one we used in Section 4.

This work was supported by the ERC grant SOLARIS (number 714381) and by ANR 3IA MIAI@Grenoble Alpes, (ANR-19-P3IA-0003). The authors would like to thank Anatoli Juditsky and the anonymous reviewers for interesting discussions that greatly improved the quality of this manuscript.

A Useful Mathematical Results

The next three lemmas are classical upper and lower bounds for smooth or strongly convex functions Nesterov (2004).

and again according to Theorem 2.1.5 of Nesterov (2004),

A.2 Useful Results to Select Step Sizes

In this section, we present basic mathematical results regarding the choice of step sizes. The proofs of the first two lemmas are trivial by induction.

$\delta_{k}=\delta$ (constant). Then $\Gamma_{k}=(1-\delta)^{k}$;

$\delta_{k}=1/(k+1)$. Then, $\Gamma_{k}=\delta_{k}=\frac{1}{(k+1)}$;

$\delta_{k}=2/(k+2)$. Then, $\Gamma_{k}=\frac{2}{(k+1)(k+2)}$;

Consider a sequence of weights $(\delta_{k})_{k\geq 0}$ in $(0,1)$. Then,

Consider the same quantities defined in the previous lemma and consider the sequence $\gamma_{k}=(1-\delta_{k})\gamma_{k\text{--}1}+\delta_{k}\mu=\Gamma_{k}\gamma_{0}+(1-\Gamma_{k})\mu$ with $\gamma_{0}\geq\mu$, and assume the relation $\delta_{k}=\gamma_{k}\eta$. Then, for all $k\geq 0$,

when $\gamma_{0}=\mu$, then $\Gamma_{k}=(1-\mu\eta)^{k}$.

when $\mu=0$, $\Gamma_{k}=\frac{1}{1+{\gamma_{0}\eta k}}$.

Proof First, we have for all $k$, $\gamma_{k}\geq\mu$ such that $\delta_{k}\geq\eta\mu$, which leads then to $\Gamma_{k}\leq\left(1-\eta{\mu}\right)^{k}$. Besides, $\gamma_{k}\geq\Gamma_{k}\gamma_{0}$ and thus $\Gamma_{k}=(1-\delta_{k})\Gamma_{k\text{--}1}\leq\left(1-{\Gamma_{k}\gamma_{0}}\eta\right)\Gamma_{k\text{--}1}$. Then, $\frac{1}{\Gamma_{k}}\left(1-\Gamma_{k}{\gamma_{0}}\eta\right)\geq\frac{1}{\Gamma_{k\text{--}1}},$ and

which is sufficient to obtain (38). Then, the fact that $\gamma_{0}=\mu$ leads to $\Gamma_{k}=(1-\mu\eta)^{k}$ is trivial, and the fact that $\mu=0$ yields $\Gamma_{k}=\frac{1}{1+{\gamma_{0}\eta k}}$ can be shown by induction. Indeed, the relation is true for $\Gamma_{0}$ and then, assuming the relation is true for $k-1$, we have for $k\geq 1$,

which leads to $\Gamma_{k}=\frac{1}{1+{\gamma_{0}\eta k}}$.

Consider the same quantities defined in Lemma 27 and consider the sequence $\gamma_{k}=(1-\delta_{k})\gamma_{k\text{--}1}+\delta_{k}\mu=\Gamma_{k}\gamma_{0}+(1-\Gamma_{k})\mu$ with $\gamma_{0}\geq\mu$, and assume the relation $\delta_{k}=\sqrt{\gamma_{k}\eta}$. Then, for all $k\geq 0$,

Besides, when $\gamma_{0}=\mu$, then $\Gamma_{k}=(1-\sqrt{\mu\eta})^{k}$.

Proof see Lemma 2.2.4 of Nesterov (2004).

A.3 Averaging Strategies

Next, we show a generic convergence result and an appropriate averaging strategy given a recursive relation between quantities acting as Lyapunov function.

Assume that there is a sequence $(x_{k})_{k\geq 1}$ generated by an algorithm that minimizes a convex function $F$, and that there exist non-negative sequences $(T_{k})_{k\geq 0}$, $(\delta_{k})_{k\geq 1}$ in $(0,1)$, $(\beta_{k})_{k\geq 1}$ and a scalar $\alpha>0$ such that for all $k\geq 1$,

where the expectation is taken with respect to any random parameter used by the algorithm. Then,

For any point $\hat{x}_{0}$, consider the averaging sequence $(\hat{x}_{k})_{k\geq 0}$,

Uniform averaging strategy.

Assume that $\delta_{k}=\frac{1}{k+1}$ and consider the average sequence $\hat{x}_{k}=\frac{1}{k}\sum_{i=1}^{k}x_{i}$. Then,

Proof Given that $T_{k}\leq(1-\delta_{k})T_{{k\text{--}1}}+\beta_{k}$, we obtain (39) by simply unrolling the recursion. To analyze the effect of the averaging strategies, divide now (39) by $\Gamma_{k}$:

Sum from $t=1$ to $k$ and notice that we have a telescopic sum:

By exploiting the relation (37), we may then use Jensen’s inequality and we obtain (41).

Consider now the specific case $\delta_{k}=\frac{1}{k+1}$, which yields $\Gamma_{k}=\frac{1}{k+1}$. Multiply then Eq. (43) by $\alpha/k$ and use Jensen’s inequality; we obtain Eq. (42).

B Relation Between Iteration (B) and MISO/SDCA

In this section, we derive explicit links between the proximal MISO algorithm Lin et al. (2015), a primal version of SDCA Shalev-Shwartz (2016), and iteration (B) when used with the gradient estimator (15) without stochastic perturbations. Under the big data condition $L/\mu\leq n$, consider indeed $\beta=\mu$, constant step-sizes $\eta_{k}=\eta=\frac{1}{n\mu}$, $\gamma_{k}=\mu$, and a uniform sampling distribution $Q$; then, we obtain the following algorithm

with $\bar{z}_{0}=\bar{x}_{0}=0$. Then, since $\mu\eta=\frac{1}{n}$, it is easy to show that in fact $\bar{z}_{k}=\mu\bar{x}_{k}$ for all $k\geq 0$. This is then exactly the proximal MISO algorithm (see Bietti and Mairal, 2017). For the relation between primal variants of SDCA and MISO, see page 4 and Equation (3) of Bietti and Mairal (2017).

C Recovering Classical Results for Proximal SGD

In this section, we present several corollaries of Theorem 2 to recover classical results for proximal variants of the stochastic gradient descent method. Throughout the section, we assume that the gradient estimates have variance bounded by $\sigma^{2}$:

Convergence results for the deterministic case $\sigma^{2}=0$ can be also recovered naturally from the corollaries. We start by applying Theorem 2 with a constant step-size strategy $\eta_{k}=1/L$, which shows convergence to a noise-dominated region of radius $\sigma^{2}/L$. In all the corollaries below, we use the notation from Theorem 2.

Assume that $f$ is $\mu$-strongly convex, choose $\gamma_{0}=\mu$ and $\eta_{k}=1/L$ with Algorithm (A) or (B). Then, for any point $\hat{x}_{0}$,

when using the averaging strategy from Theorem 2. Note that $d_{k}(x^{*})-d_{k}^{*}\geq\frac{\mu}{2}\|x_{k}-x^{*}\|^{2}$ for all $k\geq 0$ with equality for Algorithm (A).

Next, we show how to obtain converging algorithms by using decreasing step sizes.

Note that $d_{0}(x^{*})-d_{0}^{*}=\frac{\mu}{2}\|x_{0}-x^{*}\|^{2}\leq F(x_{0})-F^{*}$ for Algorithm (A).

where the second inequality uses the fact that $\frac{\mu}{2}\|x_{0}-x^{*}\|^{2}\leq F(x_{0})-F^{*}\leq\frac{2\sigma^{2}}{L}$, and then we use Lemmas 26 and 27. The term on the right is of order $O(\sigma^{2}/\mu k)$ whereas the term on the left becomes of the same order or smaller whenever $k\geq k_{0}=O(L/\mu)$. This leads to the desired iteration complexity. We may now study the case $\mu=0$, first with a constant step size. The next corollary consists of simply applying the uniform averaging strategy of Lemma 30 to Proposition 1, noting that $\delta_{k}=\frac{1}{k+1}$ for all $k\geq 0$ if $\mu=0$ and $\gamma_{0}=1/\eta$.

Assume that $f$ is convex, choose a constant step size $\eta_{k}=\eta\leq\frac{1}{L}$ with Algorithm (A) or (B) with $\gamma_{0}=1/\eta$. Then,

where $\hat{x}_{k}=\frac{1}{k}\sum_{i=1}^{k}x_{i}$. Note that $d_{0}(x^{*})-d_{0}^{*}=\frac{1}{2\eta}\|x_{0}-x^{*}\|^{2}$ for Algorithm (A).

The noise dependency is now illustrated for Algorithm (A) in the next corollary, obtained in a finite horizon setting.

Consider the same setting as in the previous corollary. Assume that we have a budget of $K$ iterations for Algorithm (A). Choose a constant step size

Then, with $\gamma_{0}=1/\eta$ and when using the averaging strategy from Corollary 33,

This corollary is obtained by optimizing the right side of (46) with respect to $\eta$ under the constraint $\eta\leq 1/L$. Considering both cases $\eta=1/L$ and $\eta=\sqrt{{T_{0}}/{K\sigma^{2}}}$, it is easy to check that we have (47) in all cases. Whereas this last result is not a practical one since the step size depends on unknown quantities, it shows that our analysis is nevertheless able to recover the optimal noise-dependency in $O(\sigma\sqrt{{T_{0}}/{K}})$, (see Nemirovski et al., 2009).

D Proofs of the Main Results

We have now two possibilities to control the quantity $A_{k}$ related to $\zeta_{k}$. First, we may simply upper bound it as follows

Then, since $x^{*}$ minimizes $F$, we have $0\in\nabla f(x^{*})+\partial\psi(x^{*})$ and thus $-\nabla f(x^{*})$ is a subgradient in $\partial\psi(x^{*})$. By using as well the convexity inequality $\psi(x)\geq\psi(x^{*})-\nabla f(x^{*})^{\top}(x-x^{*})$, we have

Then, we may combine (50) with (48) and use (49) to obtain (22).

D.2 Proof of Proposition 4

Proof To make the notation more compact, we call

Then, according to Proposition 3, we have

Then, we note that both for the SVRG and SAGA/MISO/SDCA strategies, we have (with $\beta=0$ for SVRG),

By taking a weighted average, this yields

where the second inequality comes from Lemma 25 and the last one uses similar arguments as in the proof of Proposition 3. Then, we add a quantity $\beta_{k}C_{k}$ on both sides of the relation (51) with some $\beta_{k}>0$ that we will specify later:

and then choose $\frac{\beta_{k}}{n}=\frac{5}{2}\eta_{k}\delta_{k}$, which yields

Remember that $\tau_{k}=\min\left(\delta_{k},\frac{1}{5n}\right)$, notice that the sequences $(\beta_{k})_{k\geq 0},(\eta_{k})_{k\geq 0}$ and $(\delta_{k})_{k\geq 0}$ are non-increasing and note that ${4}\leq{5}(1-\frac{1}{5n})$ for all $n\geq 1$. Then,

which immediately yields (23) with the appropriate definition of $T_{k}$, and by noting that $(1-10L_{Q}\eta_{k})\geq\frac{1}{6}$.

D.3 Proof of Theorem 5

where the first inequality uses Lemma 25, and the second one uses the definition of $L_{Q}$, whereas the last one uses (49).

D.4 Proof of Corollary 6

Proof First, notice that $\delta_{k}=\eta_{k}\gamma_{k}=\frac{\mu}{12L_{Q}}$ and that $\alpha=\frac{6\tau_{k}}{\delta_{k}}$. Then, we apply Theorem 5 and obtain

D.5 Proof of Corollary 8

Then, after taking the expectation with respect to the output of the first stage,

Denote now by $k_{0}$ the largest index such that $\frac{2}{\mu(k_{0}+2)}\geq{\eta}$ and thus $k_{0}=\lceil 2/(\mu{\eta})-2\rceil$. Then, according to Lemma 26, for $k\geq k_{0}$,

D.6 Proof of Corollary 9

Then, we consider the main run of the algorithm, and apply Theorem 5, replacing $x_{0}$ by $x_{0}^{\prime}$. With the chosen setup, we have $\delta_{k}=\frac{1}{k+1}$ and since $K\geq 5n$, we have $\delta_{K}=\tau_{K}$, such that (24) becomes

Note that Lemma 26 gives us that $\Theta_{k}=(1-1/5n)^{5n-1}\frac{5n}{k+1}\leq\frac{3n}{k+1}$ for $k\geq 5n$ and since $1+\sum_{t=1}^{K}\frac{\tau_{t}}{\Theta_{t}}=\frac{1}{\Theta_{K}}$ according to Lemma 27,

It remains to optimize it over $\eta$ to get the left side of (28).

D.7 Proof of Lemma 10

Proof Let us assume that the relation $y_{k\text{--}1}=(1-\theta_{k\text{--}1})x_{k\text{--}1}+\theta_{k\text{--}1}v_{k\text{--}1}$ holds and let us show that it also holds for $y_{k}$. Since the estimate sequences $d_{k}$ are quadratic functions, we have

Then note that $\theta_{k\text{--}1}=\frac{\delta_{k}\gamma_{k\text{--}1}}{\gamma_{k\text{--}1}+\delta_{k}\mu}$ and thus, $\frac{\gamma_{k\text{--}1}(1-\theta_{k\text{--}1})}{\gamma_{k}\theta_{k\text{--}1}}=\frac{1}{\delta_{k}}$, and

Then, we note that $x_{k}-x_{k\text{--}1}=\frac{\delta_{k}}{1-\delta_{k}}(v_{k}-x_{k})$ and we are left with

which allows us to conclude that $y_{k}=(1-\theta_{k})x_{k}+\theta_{k}v_{k}$ since the relation holds trivially for $k=0$.

D.8 Proof of Lemma 11

D.9 Proof of Corollary 14

where we use Lemmas 26 and 27. This leads to the desired iteration complexity.

D.10 Proof of Corollary 15

Proof Let us call $x_{0}^{\prime}$ the point obtained by running on step of iteration (A), which according to Theorem 2 satisfies, with $\gamma_{0}=1/\eta$,

Then, we note that according to Lemma 29, we have

and we apply Theorem 12 to obtain the relation

Optimizing with respect to $\eta$ under the constraint $\eta\leq 1/L$ gives (32).

D.11 Proof of Proposition 16

D.12 Proof of Lemma 17

Proof We can show that Lemma 11 still holds and thus,

D.13 Proof of Corollary 20

With the choice of parameters, we have $\gamma_{k}=\mu$ and $\delta_{k}=\sqrt{\frac{5\mu\eta_{k}}{3n}}=\min\left(\sqrt{\frac{5\mu\eta}{3n}},\frac{2}{k+2}\right)$. We may then apply Theorem 18 where the value of $\Gamma_{k}$ is given by Lemma 26. This yields for $k\geq k_{0}=\left\lceil\sqrt{\frac{12n}{5\mu\eta}}-2\right\rceil$,

D.14 Proof of Corollary 21

Then, we note that $\delta_{k}=\min\left(\sqrt{\frac{5\Gamma_{k}}{3n}},\frac{1}{3n}\right)$ such that $\Gamma_{k}=\left(1-\frac{1}{3n}\right)^{k}$ for $k\leq k_{0}$, where $k_{0}$ is the index such that $\left(1-\frac{1}{3n}\right)^{k_{0}+1}\leq\frac{1}{15n}<\left(1-\frac{1}{3n}\right)^{k_{0}}$, which gives us $(3n-1)\log(15n)\leq k_{0}\leq 3n(\log(15n))$. For $k>k_{0}$, we are in a constant step size regime, and we may then use Lemma 29 to obtain

Then, noticing that $K\geq 2k_{0}+1$, we have $K-k_{0}\geq(K+1)/2$, and we conclude that

Then, it remains to optimize with respect to $\eta$, under the constraint $\eta\leq 1/(3L_{Q})$, which provides (35).

References