Optimal mini-batch and step sizes for SAGA

Introduction

Consider the empirical risk minimization (ERM) problem:

where each $f_{i}$ is $L_{i}$-smooth and $f$ is $\mu$-strongly convex. Each $f_{i}$ represents a regularized loss over a sampled data point. Solving the ERM problem is often time consuming for large number of samples $n$, so much so that algorithms scanning through all the data points at each iteration are not competitive. Gradient descent (GD) falls into this category, and in practice its stochastic version is preferred.

Stochastic gradient descent (SGD), on the other hand, allows to solve the ERM incrementally by computing at each iteration an unbiased estimate of the full gradient, $\nabla f_{i}(w^{k})$ for $i$ randomly sampled in $[n]:=\{1,2,\ldots,n\}$ (Robbins & Monro, 1951). On the downside, for SGD to converge one needs to tune a sequence of asymptotically vanishing step sizes, a cumbersome and time-consuming task for the user. Recent works have taken advantage of the sum structure in Eq. 1 to design stochastic variance reduced gradient algorithms (Johnson & Zhang, 2013; Shalev-Shwartz & Zhang, 2013; Defazio et al., 2014; Schmidt et al., 2017). In the strongly convex setting, these methods lead to fast linear convergence instead of the slow ${\cal O}(1/t)$ rate of SGD. Moreover, they only require a constant step size, informed by theory, instead of sequence of decreasing step sizes.

In practice, most variance reduced methods rely on a mini-batching strategy for better performance. Yet most convergence analysis (with the Katyusha algorithm of Allen-Zhu (2017) being an exception) indicates that a mini-batch size of $b=1$ gives the best overall complexity, disagreeing with practical findings, where larger mini-batch often gives better results. Here, we show both theoretically and numerically that $b=1$ is not the optimal mini-batch size for the SAGA algorithm (Defazio et al., 2014).

Our analysis leverages recent results in (Gower et al., 2018), where the authors prove that the iteration complexity and the step size of SAGA, and a larger family of methods called the JacSketch methods, depend on an expected smoothness constant. This constant governs the trade-off between the increased cost of an iteration as the mini-batch size is increased, and the decreased total complexity. Thus if this expected smoothness constant could be calculated a priori, then we could set the optimal mini-batch size and step size. We provide simple formulas for computing the expected smoothness constant when sampling mini-batches without replacement, and use them to calculate optimal mini-batches and significantly larger step sizes for SAGA.

In particular, we provide two bounds on the expected smoothness constant, each resulting in a particular step size formula. We first derive the simple bound and then develop a matrix concentration inquality to obtain the refined Bernstein bound. We also provide substantial theoretical motivation and numerical evidence for practical estimate of the expected smoothness constant. For illustration, we plot in Figure 1 the evolution of each resulting step size as the mini-batch size grows on a classification problem (Section 5 has more details on our experimental settings).

with ${L_{\max}}:=\max_{i\in[n]}L_{i}$, $\bar{L}:=\frac{1}{n}\sum_{i=1}^{n}L_{i}$ and $\epsilon>0$ is the desired precision.

This complexity bound, and others presented in Section 3.3 show that SAGA enjoys a linear speedup as we increase the mini-batch size until an optimal one (as illustrated in Figure 2). After this point, the total complexity increases. We use this observation to develop optimal and practical mini-batch sizes and step sizes.

The rest of the paper is structured as follows. In Section 2 we first introduce variance reduction techniques after presenting our main assumption, the expected smoothnes assumption. We highlight how this assumption is necessary to capture the improvement in iteration complexity, and conclude the section by showing that to calculate the expected smoothness constant we need evaluate an intractable expectation. Which brings us to Section 3 where we directly address this issue and provide several tractable upper-bounds of the expected smoothness constant. We then calculate optimal mini-batch sizes and step sizes by using our new bounds. Finally, we give numerical experiments in Section 5 that verify our theory on artificial and real datasets. We also show how these new settings for the mini-batch size and step size lead to practical performance gains.

Background

We can introduce variance reduced versions of SGD in a principled manner by using a sampling vector.

With a sampling vector we can re-write (1) through the following stochastic reformulation

where $f_{v}(w)$ is called a subsampled function. The stochastic Problem (2) and our original Problem (1) are equivalent :

Consequently the gradient $\nabla f_{v}(w)$ is an unbiased estimate of $\nabla f(w)$ and we could use the SGD method to solve (2). To tackle the variance of these stochastic gradients we can further modify (2) by introducing control variates which leads to the following controlled stochastic reformulation:

That is, starting from a vector $w^{0}$, given a positive step size $\gamma$, we can iterate the steps

where $v^{k}\sim{\cal D}$ are i.i.d. samples at each iteration.

2 The expected smoothness constant

In order to analyze stochastic variance reduced methods, some form of smoothness assumption needs to be made. The most common assumption is

for each $i\in[n]$. That is each $f_{i}$ is uniformly smooth with smoothness constant $L_{\max}$, as is assumed in (Defazio et al., 2014; Hofmann et al., 2015; Raj & Stich, 2018) for variants of SAGAThe same assumption is made in proofs of SVRG (Johnson & Zhang, 2013), S2GD (Konečný & Richtárik, 2017) and the SARAH algorithm (Nguyen et al., 2017).. In the analyses of these papers it was shown that the iteration complexity of SAGA is proportional to $L_{\max},$ and the step size is inversely proportional to $L_{\max}.$

But as was shown in (Gower et al., 2018), we can set a much larger step size by making use of the smoothness of the subsampled functions $f_{v}.$ For this Gower et al. (2018) introduced the notion of expected smoothness, which we extend here to all sampling vectors and control variates.

Note that we refer to any positive constant ${\cal L}$ that satisfies (7) as an expected smoothness constant. Indeed ${\cal L}\rightarrow\infty$ is a valid constant in the extended reals, but as we will see, the smaller ${\cal L}$, the better for our complexity results.

Gower et al. (2018) show that the expected smoothness constant plays the same role that $L_{\max}$ does in the previously existing analysis of SAGA, namely that the step size is inversely proportional to ${\cal L}$ and the iteration complexity is proportional to ${\cal L}$ (see details in Theorem 1). Furthermore, by assuming that $f$ is $L$–smooth, the expected smoothness constant is bounded

as was proven in Theorem 4.17 in (Gower et al., 2018). Also, the bounds $L_{\max}$ and $L$ are attained when using a uniform single element sampling and a full batch, respectively. And as we will show, the constants $L_{\max}$ and $L$ can be orders of magnitude apart on large dimensional problems. Thus we could set much larger step sizes for larger mini-batch sizes if we could calculate ${\cal L}$. Though calculating ${\cal L}$ is not easy, as we see in the next lemma.

If the sampling has a very large combinatorial number of possible realizations — for instance sampling mini-batches without replacement — then this expectation becomes intractable to calculate. This observation motivates the development of functional upper-bounds of the expected smoothness constant that can be efficiently evaluated.

3 Mini-batch without replacement: b𝑏b–nice sampling

Now we will choose a distribution of the sampling vector $v$ based on a mini-batch sampling without replacement. We denote a mini-batch as $B\subseteq[n]$ and its size as $b=|B|$.

$S$ is a $b$-nice sampling if $S$ is a set valued map with a probability distribution given by

where $\mathds{1}_{S}$ denotes the indicator function of the random set $S$. Now taking expectation in (9) gives

using $|\{{{B\subseteq[n]:|B|=b\wedge i\in B}}\}|=\binom{n-1}{b-1}$.

Here we are interested in the mini-batch SAGA algorithm with $b$-nice sampling, which we refer to as the $b$-nice SAGA. In particular, $b$-nice SAGA is the result of using $b$-nice sampling, together with a linear model for the control variate $z_{v}(w)$. Different choices of the control variate $z_{v}(w)$ also recover popular algorithms such as gradient descent, SGD or the standard SAGA method (see Table 1 for some examples).

A naive implementation of $b$-nice SAGA based on the JacSketch algorithm is given in Algorithm 1We also provide a more efficient implementation that we used for our experiments in the appendix in Algorithm 2..

Upper Bounds on the Expected Smoothness

To determine an optimal mini-batch size $b^{*}$ for $b$-nice SAGA, we first state our assumptions and provide bounds of the smoothness of the subsampled function. We then define $b^{*}$ as the mini-batch size that minimizes the total complexity of the considered algorithm, i.e., the total number of stochastic gradients computed. Finally we provide upper-bounds on the expected smoothness constant ${\cal L}$, through which we can deduce optimal mini-batch sizes. Many proofs are deferred to the supplementary material.

We consider that the objective function is a GLM with quadratic regularization controlled by a parameter $\lambda>0$:

We assume that the second derivative of each $\phi_{i}$ is uniformly bounded, which holds for our aforementioned examples.

For a batch $B\subseteq[n]$, we rewrite the subsampled function as

and its second derivative is thus given by

where $I_{d}$ denotes the identity matrix of size $d$.

For a symmetric matrix $M$, we write $\lambda_{\max}(M)$ (resp. $\lambda_{\min}(M)$) for its largest (resp. smallest) eigenvalue. Assumption 1 directly implies the following.

Let $B\subset[n]$, and let $A_{B}=[a_{i}]_{i\in B}$ denote the column concatenation of the vectors $a_{i}$ with $i\in B.$ The smoothness constant of the subsampled loss function $\frac{1}{|B|}\sum_{i\in B}\phi_{i}(a_{i}^{\top}w)$ is given by

Proof. The proof follows from Assumption 1 as

Combined with (11), we get that $f_{B}$ is $(L_{B}+\lambda)$-smooth.

Another key quantity in our analysis is the strong convexity parameter.

The strong convexity parameter is given by

Since we have an explicit regularization term with $\lambda>0$, $f$ is strongly convex and $\mu\geq\lambda>0.$

We additionally define $L_{i}$, resp. $L$, as the smoothness constant of the individual function $\phi_{i}(a_{i}^{\top}w)$, resp. the whole function $\frac{1}{n}\sum_{i=1}^{n}\phi_{i}(a_{i}^{\top}w)$. We also recall the definitions of the maximum of the individual smoothness constants by ${L_{\max}}:=\max_{i\in[n]}L_{i}$ and their average by $\bar{L}:=\frac{1}{n}\sum_{i=1}^{n}L_{i}$. The three constants satisfies

The proof of (13) is given in Lemma 10 in the appendix.

2 Path to the optimal mini-batch size

Our starting point is the following theorem taken from combining Theorem 3.6 and Eq. (103) in (Gower et al., 2018)Note that $\lambda$ has been added to every smoothness constant since the analysis in Gower et al. (2018) depends on the $(L+\lambda)$-smoothness of $f$ and the $(L_{B}+\lambda)$-smoothness of the subsampled functions $f_{B}$..

Consider the iterates $w^{k}$ of Algorithm 1. Let the step size be given by

Through Theorem 1 we can now explicitly see how the expected smoothness constant ${\cal L}$ controls both the step size and the resulting iteration complexity. This is why we need bounds on ${\cal L}$ so that we can set the step size. In particular, we will show that the expected smoothness constant is a function of the mini-batch size $b$. Consequently so is the step size, the iteration complexity and the total complexity. We denote $K_{\text{total}}$ the total complexity defined as the number of stochastic gradients computed, hence with (15),

As we have shown in Lemma 1, computing a precise bound on ${\cal L}$ can be computationally intractable. This is why we focus on finding upper bounds on ${\cal L}$ that can be computed, but also tight enough to be useful. To verify that our bounds are sufficiently tight, we will always have in mind the bounds $L\leq{\cal L}\leq L_{\max}$ given in (8). In particular, after expressing our bounds of ${\cal L}={\cal L}(b)$ as a function of $b$,we would like the bounds (8) to be attained for ${\cal L}(1)=L_{\max}$ and ${\cal L}(n)=L.$

3 Expected smoothness

All bounds we develop on ${\cal L}$ are based on the following lemma, which is a specialization of (1) for $b$-nice sampling.

For the $b$-nice sampling, with $b\in[n]$, the expected smoothness constant is given by

Let $S$ the $b$-nice sampling as defined in Definition 3 and let $v=\frac{n}{b}\sum_{j\in S}e_{j}$ be its corresponding sampling vector. Note that

Taking the maximum over all $i\in[n]$ gives the result. ∎

The first bound we present is technically the simplest to derive, which is why we refer to it as the simple bound.

For a $b$-nice sampling $S$, for $b\in[n]$, we have that

The proof, given in Section A.2, starts by using the that $L_{B}\leq\frac{1}{b}\sum_{j\in B}L_{j}$ for all subsets $B$, which follows from repeatedly applying Lemma 8 in the appendix. The remainder of the proof follows by straightforward counting arguments. ∎

The previous bound interpolates, respectively for $b=1$ and $b=n$, between $L_{\max}$ and $\bar{L}.$ On the one hand, we have that ${\cal L}_{\operatorname{simple}}(b)$ is a good bound for when $b$ is small, since ${\cal L}_{\operatorname{simple}}(1)=L_{\max}$. Though ${\cal L}_{\operatorname{simple}}(b)$ may not be a good bound for large $b$, since ${\cal L}_{\operatorname{simple}}(n)=\overline{L}\quad{\geq}\quad L$, thanks to (13). Thus ${\cal L}_{\operatorname{simple}}(b)$ does not achieve the left-hand side of (8). Indeed $\overline{L}$ can be far from $L$. For instanceWe numerically explore such extreme settings in Section 5., if $f(w)=\frac{1}{n}\sum_{i\in[n]}\frac{1}{2}(a_{i}^{\top}w-b_{i})^{2}$ is a quadratic function, then we have that $\overline{L}=\frac{1}{n}\mbox{Tr}\left(AA^{\top}\right)$ and $L=\frac{1}{n}\lambda_{\max}(AA^{\top})$. Thus if the eigenvalues of $AA^{\top}$ are all equal then $\overline{L}=dL$. Alternatively, if one eigenvalue is significantly larger than the rest then $\overline{L}\approx L$.

Due to this shortcoming of ${\cal L}_{\operatorname{simple}},$ we now derive the Bernstein bound. This bound explicitly depends on $L$ instead of $\bar{L}$, and is developed through a specialized variant of a matrix Bernstein inequality (Tropp, 2012, 2015) for sampling without replacement in Appendix C.

The expected smoothness constant is upper bounded by

Checking again the bounds of ${\cal L}_{\text{Bernstein}}(b)$, we have on the one hand that ${\cal L}_{\text{Bernstein}}(1)=\left(1+\tfrac{4}{3}\log d\right)L_{\max}\geq L_{\max},$ thus there is a little bit of slack for $b$ small. On the other hand, using $\tfrac{1}{n}L_{\max}\leq L$ (see Lemma 10 in appendix), we have that

which depends only logarithmically on $d$. Thus we expect the Bernstein bound to be more useful in the large $d$ domains, as compared to the simple bound. We confirm this numerically in Section 5.1.

The simple bound is relatively tight for $b$ small, while the Bernstein bound is better for large $b$ and large $d$. Fortunately, we can obtain a more refined bound by taking the minimum of the simple and the Bernstein bounds. This is highlighted numerically in Section 5.

Next we propose a practical estimate of ${\cal L}$ that is tight for both small and large mini-batch sizes.

Indeed ${\cal L}_{\text{practical}}(1)=L_{\max}$ and ${\cal L}_{\text{practical}}(n)=L,$ achieving both limits of (8). The downside to ${\cal L}_{\text{practical}}(b)$ is that it is not an upper bound of ${\cal L}$. Rather, we are able to show that ${\cal L}_{\text{practical}}(b)$ is very close to a valid smoothness constant, but it can be slightly smaller. Our theoretical justification for using ${\cal L}_{\text{practical}}(b)$ comes from a mid step in the proof of the Bernstein bound which is captured in the next lemma.

with $N_{i}:=\frac{1}{b}\sum_{j\in S^{i}}a_{j}a_{j}^{\top}-\frac{1}{b}\frac{b-1}{n-1}\sum_{j\in[n]\setminus\{i\}}a_{j}a_{j}^{\top}$.

Lemma 3 shows that the expected smoothness constant is upper-bounded by ${\cal L}_{\text{practical}}(b)$ and an additional term. In this additional term we have the largest eigenvalue of a random matrix. This matrix is zero in expectation, and we also find that its eigenvalues oscillate around zero. Indeed, we provide extensive experiments in Section 5 confirming that ${\cal L}_{\text{practical}}(b)$ is very close to ${\cal L}$ given in (17).

Optimal Mini-Batch Sizes

Now that we have established the simple and the Bernstein bounds, we can minimize the total complexity (16) in the mini-batch size.

For instance for the simple bound, given $\epsilon>0$ and plugging in (18) into (16) gives

where $g_{\text{simple}}(b):=\tfrac{4\left(n\bar{L}-L_{\max}+(n-1)\lambda\right)}{\mu(n-1)}b+\tfrac{4n\left(L_{\max}-\bar{L}\right)}{\mu(n-1)},$ and $h(b):=-\tfrac{4(L_{\max}+\lambda)}{\mu(n-1)}b+n\left(1+\tfrac{1}{n-1}\tfrac{4(L_{\max}+\lambda)}{\mu}\right).$

The right-hand side term $h(b)$ is common to all our bounds since it does not depend on ${\cal L}$. It linearly decreases from $h(1)=n+\frac{4(L_{\max}+\lambda)}{\mu}$ to $h(n)=n$.

We note that $g_{\text{simple}}(b)$ is a linearly increasing function of $b$, because $L_{\max}\leq n\bar{L}$ (as proven in Lemma 10). One can easily verify that $g_{\text{simple}}(b)$ and $h(b)$ cross, as presented in Figure 2, by looking at initial and final values:

At $b\!=\!1$, $g_{\text{simple}}(1)=\tfrac{4}{\mu}(L_{\max}+\lambda)=h(1)-n$. So, $g_{\text{simple}}(1)\leq h(1)$.

At $b\!=\!n$, $g_{\text{simple}}(n)=\tfrac{4(\bar{L}+\lambda)}{\mu}n=\tfrac{4(\bar{L}+\lambda)}{\mu}h(n)$. Since $\bar{L}\geq\mu$, we get $g_{\text{simple}}(n)\geq h(n)$.

Consequently, solving $g_{\text{simple}}(b)=h(b)$ in $b$ gives the optimal mini-batch size

For the Bernstein bound, plugging (19) into (16) leads to

The function $g_{\text{Bernstein}}$ is also linearly increasing in $b$ and its initial and final values are

At $b\!=\!1$, $g_{\text{Bernstein}}(1)=(1+\frac{4}{3}\log d)\tfrac{4L_{\max}}{\mu}+\tfrac{4\lambda}{\mu}.$

At $b\!=\!n$, $g_{\text{Bernstein}}(n)\!=\!n\tfrac{4(2L+\lambda)}{\mu}\!+\!\tfrac{16}{3\mu}{(L_{\max}\!+\!\lambda)}\log(d)$. Since $L\geq\mu$, we get $g_{\text{Bernstein}}(n)\geq h(n)$.

Yet, it is unclear whether $g_{\text{Bernstein}}(1)$ is dominated by $h(1)$. This is why we need to distinguish two cases to minimize the total complexity, which leads to the following solution

Numerical Study

All the experiments were run in Julia and the code is freely available on https://github.com/gowerrobert/StochOpt.jl.

In Figure 4 we see that ${\cal L}_{\text{practical}}$ is arbitrarily close to ${\cal L}$, making it hard to distinguish the two line plots. This was the case in many other experiments, which we defer to Section E.1. For this reason, we use $\gamma_{\text{practical}}$ in our experiments with the SAGA method.

Furthermore, in accordance with our discussion in Section 3.3, we have that ${\cal L}_{\text{simple}}$ and ${\cal L}_{\text{Bernstein}}$ are close to ${\cal L}$ when $b$ is small and large, respectively. In Section E.2 we show, by applying ridge and regularized logistic regression to publicly available datasets from LIBSVMhttps://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ and the UCI repositoryhttps://archive.ics.uci.edu/ml/datasets/, that the simple bound performs better than the Bernstein bound when $n\gg d$, and conversely for $d$ slightly smaller than $n$, larger than $n$ or when scaling the data.

2 Related step size estimation

Different bounds on ${\cal L}$ also give different step sizes (14). Plugging in our estimates ${\cal L}_{\text{simple}}$, ${\cal L}_{\text{Bernstein}}$ and ${\cal L}_{\text{practical}}$ into (14) gives the step sizes $\gamma_{\text{simple}}$, $\gamma_{\text{Bernstein}}$ and $\gamma_{\text{practical}}$, respectively. We compare our resulting step sizes to $\gamma_{{\cal L}}$ where ${\cal L}$ is given by Eq. 17 and to the step size given by Hofmann et al. (2015), which is $\gamma_{\text{Hofmann}}(b)=\frac{K}{2L_{\max}(1+K+\sqrt{1+K^{2}})}$, where $K:=\frac{4bL_{\max}}{n\mu}.$ We can see in Figure 4, that for $b=1$, all the step sizes are approximately the same, with the exceptions of the Bernstein step size. For $b>5$, all of our step sizes are larger than $\gamma_{\text{Hofmann}}(b)$, in particular $\gamma_{\text{practical}}(b)$ is significantly larger. These observations are verified in other artificial and real data examples in Sections E.3 and E.4.

3 Comparison with previous SAGA settings

Here we compare the performance of SAGA when using the mini-batch size and step size $b=1,\gamma_{\text{Defazio}}:=1/{3(n\mu+L_{\max})}$ given in (Defazio et al., 2014), $b=20$ and $\gamma_{\text{Hofmann}}={20}/{n\mu}$ given in Hofmann et al. (2015), to our new practical mini-batch size $b_{\text{practical}}=\left\lfloor 1+\frac{\mu(n-1)}{4(L+\lambda)}\right\rfloor$ and step size $\gamma_{\text{practical}}$. Our goal is to verify how much our parameter setting can improve practical performance. We also compare with a step size $\gamma_{\text{gridsearch}}$ obtained by grid search over odd powers of $2$. These methods are run until they reach a relative error of $10^{-4}$.

We find in Figure 5 that our parameter settings $(\gamma_{\text{practical}},b_{\text{practical}})$ significantly outperforms the previously suggested parameters, and is even comparable to grid search. Finally, In we show in Section E.5 that the settings $(\gamma_{\text{Hofmann}},b=20)$ can lead to very poor performance compared to our settings.

4 Optimality of our mini-batch size

In the last experiment, detailed in Section E.6, we show that our estimation of the optimal mini-batch size $b_{\text{practical}}$ leads to a faster implementation of SAGA. We build a grid of mini-batch sizesOur grid is $\{2^{i},i=0,\dots,14\}$, with $2^{16},2^{18}$ and $n$ being added when needed. and compute the empirical total complexity required to achieve a relative error of $10^{-4}$, as in Section 5.3.

In Figure 6 we can see that the empirical complexity when using $b_{\text{practical}}$ is almost the same as one for the optimal mini-batch size calculated through grid search. Yet, our $b_{\text{practical}}$ being always larger than the mini-batch obtained by grid search, it leads to a faster algorithm for two reasons. Firstly, the corresponding step size is larger and secondly, and secondly, computing the stochastic gradients in parallel improves the running time. What is even more interesting, is that $b_{\text{practical}}$ always predicts a regime change, where using a larger mini-batch size results in a much larger empirical complexity.

Conclusions

We have explained the crucial role of the expected smoothness constant ${\cal L}$ in the convergence of a family of stochastic variance-reduced descent algorithms. We have developed functional upper-bounds of this constant to build larger step sizes and closed-form optimal mini-batch values for the $b$-nice SAGA algorithm. Our experiments on artificial and real datasets showed the validity of our upper-bounds and the improvement in the total complexity using our step and optimal mini-batch sizes. Our results suggest a new parameter setting for mini-batch SAGA, that significantly outperforms previous suggested ones, and is even comparable with a grid search approach, without the computational burden of the later.

Acknowledgements

This work was supported by grants from DIM Math Innov Région Ile-de-France (ED574 - FMJH) and by a public grant as part of the Investissement d’avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH, in a joint call with Gaspard Monge Program for optimization, operations research and their interactions with data sciences.

References

Appendix A Proofs of the Upper Bounds of ℒℒ{\cal L}

Since the $f_{i}$’s are convex, each realization of $f_{v}$ is convex, and it follows from equation 2.1.7 in (Nesterov, 2014) that

Taking expectation over the sampling gives

where in the last equality the full gradient vanishes because it is computed at optimality. The result now follows by comparing the above with the definition of expected smoothness in (7). ∎

A.2 Proof of the simple bound

To derive this bound on ${\cal L}$ we use that

which follows from repeatedly applying Lemma 8. For $b\geq 2$, it follows from Equation 17 and Equation 25 that

Using a double counting argument we can show that

We also verify that this bound is valid for $1$-nice sampling. Indeed, we already have that in this case ${\cal L}=L_{\max}$. ∎

A.3 Proof of the Bernstein bound

To start the proof of Theorem 3, we re-write the expected smoothness constant as the maximum over an expectation. Let $S^{i}$ be a $(b-1)$-nice sampling over $[n]\setminus\{i\}.$ We can write

One can come back to the definition of the subsample smoothness constant Equation 12 and interpret previous expression as an expectation of the largest eigenvalue of a sum of matrices. This insight allows us to apply a matrix Bernstein inequality, see Theorem 7, to bound ${\cal L}$.

For the proof of Theorem 3, we first need the two following results.

This results follows using a double-counting argument at the fourth line of the computation.

We then introduce another two lemmas which give a first intermediate bound.

where in the first inequality we add and remove the mean and then apply Lemma 8. In the second equality we explicit the mean with Lemma 4 and in the last inequality we use again Lemma 8 for the left-hand side term. Finally, we multiply by $U$ on both sides of the inequality. ∎

We recall the following lemma used to introduced the practical estimate given by

The result comes from applying re-writing ${\cal L}$ as an expectation of the largest eigenvalue of a sum of matrices. Then we apply Lemma 5 and then taking the maximum over all $i\in[n]$. Thus, we have

with $N:=\frac{1}{b}\sum_{j\in S^{i}}a_{j}a_{j}^{\top}-\frac{1}{b}\frac{b-1}{n-1}\sum_{j\in[n]\setminus\{i\}}a_{j}a_{j}^{\top}$.

Using this notation, the matrix $N$ which can be further decomposed as

where we have encoded the sampling $S^{i}$ using unit coordinate vectors. The matrices $M_{1},\ldots,M_{b-1}$ are sampled without replacement from the set

Now let $X_{1},\ldots,X_{b}$ be matrices sampled with replacement from (34) and let $X_{k}:=\frac{1}{b}\sum_{j\in[n]\setminus\{i\}}\left(z_{j}^{k}-\frac{1}{n-1}\right)a_{j}a_{j}^{\top}$ and $Y:=\sum_{k=1}^{b-1}X_{k}$ thus the vectors $z^{k}$ are sampled with replacement from $\{e_{1},\ldots,e_{i-1},e_{i+1},\ldots,e_{n}\}.$ Consequently

We are now in a position to apply the Bernstein matrix inequality. To this end we have

Let $k^{*}$ be the unique index such that $z^{k}_{k^{*}}=1.$ We have a uniform bound of the largest eigenvalue of our $X_{k}$

where we applied the Lemma 9 in the first inequality.

Summing in (36), taking the largest eigenvalue and applying Lemma 9 results in

Considering Equations 35 and 37 and applying the matrix Bernstein concentration inequality in Theorem 7 we get

Taking the maximum over $i$ and using $L_{[n]\setminus\{i\}}\leq\frac{n}{n-1}L$ we have that

Combining the above result with (32) leads us to

where in the second inequality we used the inequality $\sqrt{2ab}\leq a+b$. ∎

Appendix B Linear Algebra Tools

This appendix is dedicated to the presentation of useful results to manipulate more easily the smoothness constants.

Let us recall some useful spectral results on Hermitian and positive semi-definite matrices.

whenever $i,j\geq 1$ and $i+j-1\leq n\enspace.$

Moreover, as a direct consequence of the variational characterization of eigenvalues, namely

we have an inequality between the maximum diagonal term of a positive semi-definite matrices and its maximum eigenvalue.

The following lemma is a direct consequence of Weyl’s inequality for $i=j=1.$

Lastly, we present a result arising from previous lemma.

where the first inequality stems from Lemma 8 and the second from $B\succeq 0.$ ∎

B.2 Basic properties of the smoothness constants

The complexity results of Gower et al. (2018) depends on smoothness constants defined in Section 3.1. Here are some inequalities giving an idea of the order of those constants.

Let $\emptyset\neq B\subseteq[n]=\{1,\dots,n\}$ a batch set drawn randomly without replacement. The following inequalities hold

One directly gets that $L_{i}\leq\max_{j=1,\ldots,n}L_{j}={L_{\max}}$.

This inequality states that the smoothness constant $L_{B}$ of the averaged function $f_{B}$ is upper bounded by the average of the corresponding smoothness constants $L_{i}$, over the batch $B$. The proof consists in $|B|$ repetitive calls of Lemma 8.

$(a)$ Direct implication of (ii) for $B=[n]$.

$(c)$ Let us first recall the matrix formulation of our smoothness constants:

Dividing the above by $n$ on both sides gives

Appendix C Matrix Bernstein Inequality: Sampling Without Replacement

In this appendix, we present the matrix Bernstein inequality for independent Hermitian matrices from Tropp (2015). We also provide another version of this theorem for matrices sampled without replacement and prove it as explicitly as possible, taking our inspiration from Tropp (2011). The proof is based the possibility of transferring the results from sampling with to without through the inequality (50) due to Gross & Nesme (2010). The exact same work can be done for the tail bound, which is for instance used in Bach (2013).

We first present Theorem 4 which gives a Bernstein inequality for a sum of random and independent Hermitian matrices whose eigenvalues are upper bounded. If the matrices $X_{k}$ are sampled from a finite set ${\cal X}$, one can interpret this random sampling of independent matrices as a random sampling with replacement.

Consider a finite sequence $\{X_{k}\}_{k=1,\dots,n}$ of $n$ independent, random, Hermitian matrices with dimension $d$. Assume that

Let $v(S_{X})$ be the matrix variance statistic of the sum:

This theorem is the one we extend in Theorem 7 to the case when the random matrices $X_{k}$ are sampled without replacement from a finite set ${\cal X}$. We drew our inspiration from the proof of the matrix Chernoff inequality in Tropp (2011) and the one of the matrix Bernstein tail bound in Bach (2013), both in the case of sampling without replacement.

C.2 Technical random matrices prerequisites

Before proving Theorem 7, which extends the matrix Bernstein inequality to sampling without replacement, we need to introduce the key tools of the matrix Laplace transform technique. This technique is precious to prove tail bounds for sums of random matrices such as Chernoff, Hoeffding or Bernstein bounds, as presented in (Tropp, 2012).

Here, $\left\|\cdot\right\|$ denotes the spectral norm, which is defined for any Hermitian matrix $H$ by

We also introduce the moment generating function (mgf) and the cumulant generating function (cgf) of a random matrix, which are essential in the Laplace transform method approach.

These expectations may not exist for all values of $\theta$.

Let $X$ be a random Hermitian matrix. Then

This proposition is an adaptation of the Laplace transform method to obtain a bound of the expectation of the maximum eigenvalue of a random Hermitian matrix. Contrary to the tail bounds, there is no exact analog of the expectation bounds in the scalar setting.

Fix a positive number $\theta$. Because $\lambda_{\max}(\cdot)$ is a positive-homogeneous map, we have

where in the third line we used the Jensen’s inequality, in the fourth one the spectral mapping theorem and in the last line the domination by the trace of a positive-definite matrix. ∎

Let $H$ be a fixed Hermitian matrix with dimension $d$. The function

is a concave map on the the convex cone of $d\times d$ positive-definite matrices.

Let $H$ be a fixed Hermitian matrix with dimension $d$. Let $X$ be a random Hermitian matrix of same dimension. The following inequality holds

is a concave map on the the convex cone of $d\times d$ positive-definite matrices.

where the inequality comes from the application of Theorem 5 and Jensen’s inequality. ∎

Let $X$ a random Hermitian matrix such that

C.3 Extended results for sampling without replacement

This section is dedicated to the main result, Lemma 13, needed for transferring results from sampling with to without replacement. This lemma is actually the matrix version of a classical result from Hoeffding (1963). We then combine it with previous results of Section C.2 to produce a new master bound in Theorem 6, which is the key inequality of the proof of Theorem 7.

The left-hand side equality directly arises from Definition 6 and the fact that the trace commutes with the expectation because it is a linear operator. For the right-hand side inequality, see the proof in Gross & Nesme (2010). ∎

Consider two finite sequences, of same length $n$, $\{X_{k}\}_{k=1,\dots,n}$ and $\{Y_{k}\}_{k=1,\dots,n}$ of Hermitian random matrices of same size sampled respectively with and without replacement from a finite set ${\cal X}$. Then

This theorem is a modified version of Theorem 3.6.1 in Tropp (2015) for a sum of matrices sampled without replacement.

Consider two finite sequences, of same length, $\{X_{k}\}$ and $\{Y_{k}\}$ of Hermitian random matrices of same size sampled respectively with and without replacement from a finite set ${\cal X}$. Let $\theta$ a positive number.

where we used successively Proposition 2, Lemma 13 and Lemma 11. First, we use the expectation bound for the maximum eigenvalue. We then use the main result of Gross & Nesme (2010) and invoked in Tropp (2011) to extend the matrix Chernoff bound for matrices sampled without replacement. This lemma allows us to transfer our results to sampling with replacement. And finally, we then apply the subadditivity of matrix cgfs to get the desired result. ∎

C.4 Bernstein inequality for sampling without replacement

The following theorem is almost the same than Theorem 4, but in the case of matrices sampled without replacement from a finite set. The proof stems from results established in previous Sections C.3 and C.2.

Let ${\cal X}$ be a finite set of Hermitian matrices with dimension $d$ such that

Sample two finite sequences, of same length $n$, $\{X_{k}\}_{k=1,\dots,n}$ and $\{Y_{k}\}_{k=1,\dots,n}$ uniformly at random from ${\cal X}$ respectively with and without replacement such that

Let $v(S_{X})$ be the matrix variance statistic of the second sum

Consider ${\cal X}$ a finite set of Hermitian matrices of dimension $d$ such that

Sample two finite sequences, of same length, $\{X_{k}\}$ and $\{Y_{k}\}$ uniformly at random from ${\cal X}$ respectively with and without replacement such that

The $\{X_{k}\}$ matrices are thus independent. Introduce the sums $S_{X}=\sum_{k=1}^{n}X_{k}$ and $S_{Y}=\sum_{k=1}^{n}Y_{k}$. Let us bound the expectation of the largest eigenvalue of the latter

Appendix D Miscellaneous

Appendix E Additional Experiments

As described in Section 5, we compute our the simple and Bernstein bounds, our practical estimate and the true ${\cal L}$ for ridge regression applied to small artificial datasets: uniform $(n=24,d=50)$, staircase eigval $(n=d=24)$ and alone eigval $(n=d=24)$. Figure 7 shows first that the practical estimate is a very close approximation of ${\cal L}$. On the one hand, we observe in Figure 7(a) that the Bernstein bound performs poorly since the feature dimension is very small $d=50$. On the other hand, Figure 7(c) shows a regime change for $b\approx 10$, which highlight the usefulness of combining our bounds to approximate the expected smoothness constant. Finally, we observe that for the alone eigval dataset Figure 7(b), which has one very large eigenvalue far from the rest of the spectrum, the simple bound matches ${\cal L}$ because the gap between $\bar{L}$ and $L$ shrinks. Indeed, in this configuration $\bar{L}\approx L\approx\frac{L_{\max}}{n}$. When the spectrum is more concentrated, like for staircase eigval, we get a significant gap between $\bar{L}$ and $L$ as shown in Figure 7(c), where the simple bound is far from ${\cal L}$ when $b=n$.

We also report the influence of changing the value of the regularization parameter $\lambda$. Figure 8 shows that this parameter has little impact on the general shape of the bounds and of ${\cal L}$.

Finally, we study the impact of scaling or standardizing (i.e., removing the mean and dividing by the standard deviation for each feature) our artificial datasets. In order not to benefit from the diagonal shape of the alone eigval and staircase eigval datasets we also give examples of the bounds of ${\cal L}$ after a rotation of the data. The rotation aims at preserving the spectrum while erasing the diagonal structure of the covariance matrix $AA^{\top}$. This rotation procedure consists in transforming $A$ into $Q^{\top}AQ$, where $Q$ is the orthogonal matrix given by the QR decomposition of a random squared matrix (with dimension the same as the one of $A$) with uniformly random coefficients $M$, such that $M=QR$.

We observe in Figure 9 that rotations do not affect our estimates of ${\cal L}$, because they preserve the spectrum. Scaling non-diagonal datasets does not change the general shape neither. As predicted, scaling diagonal matrices leads to a particular case where the spectrum of the covariance matrix is flattened and for all $i\in[n]$, $L_{i}\approx L_{\max}\approx\bar{L}$. This is why we get a flat simple bound in Figures 9(c) and 9(g). Even after those different types of preprocessing (rotation and scaling) and with different values of $\lambda$, we end up with the same strong observation that the practical estimate is a very sharp approximation of the expected smoothness constant.

E.2 Experiment 1: estimates of the expected smoothness constant for real datasets

In what folows, we also used publicly available datasets from LIBSVMhttps://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ provided by Chang & Lin (2011) and from the UCI repositoryhttps://archive.ics.uci.edu/ml/datasets/ provided by Dheeru & Karra Taniskidou (2017). We applied ridge regression to the following datasets: YearPredictionMSD $(n=515,345,d=90)$ from LIBSVM and slice $(n=53,500,d=384)$ from UCI. We also applied regularized logistic regression for binary classification on ijcnn1 $(n=141,691,d=22)$, covtype.binary $(n=581,012,d=54)$, real-sim $(n=72,309,d=20,958)$, rcv1.binary $(n=697,641,d=47,236)$ and news20.binary ($n=19,996,d=1,355,191$) from LIBSVM. When a test set was available, we concatenated it with the train set to have more samples.

One can observe in Figure 10, that for unscaled datasets the Bernstein bound performs better than the simple bound, except for YearPredictionMSD $(n=515,345,d=90)$ and covtype.binary $(n=581,012,d=54)$. From Figure 11, we observe that after feature-scaling, the Bernstein bound is always a tighter upper bound of ${\cal L}$ than the simple bound.

E.3 Experiment 2: step size estimates for artificial datasets

In this section we give the step sizes estimate corresponding to the expected smoothness constant, the simple and Bernstein upper-bounds and the practical estimate for our small artificial datasets. In Figure 12, we show that the practical step size estimate is larger than all others. Moreover, for except for small value sof $b$, our $\gamma_{\text{simple}}$ or $\gamma_{\text{Bernstein}}$ estimates are in most cases larger than the one proposed in (Hofmann et al., 2015).

E.4 Experiment 2: step size estimates for real datasets

Here we show the step sizes estimate corresponding to the simple and Bernstein upper-bounds and the practical estimate for real datasets detailed in Section E.2. On these real data, unscaled in Figure 13 and scaled in Figure 14, we see that the gap between our step size estimates and $\gamma_{\text{Hofmann}}$ is even larger. We observe in Figure 13, accordlingly to previous remarks in Section E.2, that Bernstein bound leads to larger step sizes than the simple one, except for the unscaled YearPredictionMSD and covtype.binary datasets. Yet, as noticed before, Figure 14 seems to show that scaling the data leads to $\gamma_{\text{Bernstein}}$ larger than $\gamma_{\text{simple}}$.

E.5 Experiment 3: comparison with previous SAGA settings

In this section we provide more example of the performance of our practical settings compared to previously known SAGA settings. In Figures 16, 17, 18, 19, 20 and 21 we run our experiments on real datasets introduced in detail in Section E.1. SAGA implementations are run until the suboptimality reaches a relative error $(f(w^{k})-f(w^{*}))/(f(w^{0})-f(w^{*}))$ of $10^{-4}$, except in some cases where the Hofmann’s runs exceeded our maximal number of epochs like in Figure 17. In Figure 19, the curves corresponding to Hofmann’s settings are not displayed because they achieve a total complexity which is too large. Figure 15 shows an example of such a configuration.

These experiments show that our settings ($b_{\text{practical}},\gamma_{\text{practical}}$) most of the time outperforms whether the classical ($b=1,\gamma_{\text{Defazio}}$) or the ($b=20,\gamma_{\text{Hofmann}}$) settings both in terms of epochs and running time.

E.6 Experiment 4: optimality of the mini-batch size

We always observe a change of regime in the empirical complexity. For small values of $b$, the complexity is of the same order of magnitude, then, for values greater than the empirical optimal mini-batch size, the complexity explodes. This experiment shows that our optimal mini-batch size $b_{\text{practical}}$ correctly designates the largest mini-batch achieving the best complexity as large as possible, without reaching the regime where the total complexity explodes.