Mini-Batch Semi-Stochastic Gradient Descent in the Proximal Setting

I Introduction

In this work we are concerned with the problem of minimizing the sum of two convex functions,

where the first component, $F$, is smooth, and the second component, $R$, is possibly nonsmooth (and extended real-valued, which allows for the modeling of constraints).

In the last decade, an intensive amount of research was conducted into algorithms for solving problems of the form (1), largely motivated by the realization that the underlying problem has a considerable modeling power. One of the most popular and practical methods for (1) is the accelerated proximal gradient method of Nesterov , with its most successful variant being FISTA .

In many applications in optimization, signal processing and machine learning, $F$ has an additional structure. In particular, it is often the case that $F$ is the average of a number of convex functions:

Indeed, even one of the most basic optimization problems—least squares regression—lends itself to a natural representation of the form (2).

For problems of the form (1)+(2), and especially when $n$ is large and when a solution of low to medium accuracy is sufficient, deterministic methods do not perform as well as classical stochasticDepending on conventions used in different communities, the terms randomized or sketching are used instead of the word stochastic. In signal processing, numerical linear algebra and theoretical computer science, for instance, the terms sketching and randomized are used more often. In machine learning and optimization, the terms stochastic and randomized are used more often. In this paper, stochasticity does not refer to a data generation process, but to randomization embedded in an algorithm which is applied to a deterministic problem. Having said that, the deterministic problem can and often does arise as a sample average approximation of stochastic problem (average replaces an expectation), which further blurs the lines between the terms. methods. The prototype method in this category is stochastic gradient descent (SGD), dating back to the 1951 seminal work of Robbins and Monro . SGD selects an index $i\in\{1,2,\dots,n\}$ uniformly at random, and then updates the variable ${x}$ using $\nabla f_{i}({x})$ — a stochastic estimate of $\nabla F({x})$. Note that the computation of $\nabla f_{i}$ is $n$ times cheaper than the computation of the full gradient $\nabla F$. For problems where $n$ is very large, the per-iteration savings can be extremely large, spanning several orders of magnitude.

These savings do not come for free, however (modern methods, such as the one we propose, overcome this – more on that below). Indeed, the stochastic estimate of the gradient embodied by $\nabla f_{i}$ has a non-vanishing variance. To see this, notice that even when started from an optimal solution $x_{*}$, there is no reason for $\nabla f_{i}(x_{*})$ to be zero, which means that SGD drives away from the optimal point. Traditionally, there have been two ways of dealing with this issue. The first one consists in choosing a decreasing sequence of stepsizes. However, this means that a much larger number of iterations is needed. A second approach is to use a subset (“minibatch”) of indices $i$, as opposed to a single index, in order to form a better stochastic estimate of the gradient. However, this results in a method which performs more work per iteration. In summary, while traditional approaches manage to decrease the variance in the stochastic estimate, this comes at a cost.

I-B Modern stochastic methods

Very recently, starting with the SAG , SDCA , SVRG and S2GD algorithms from year 2013, it has transpired that neither decreasing stepsizes nor mini-batching are necessary to resolve the non-vanishing variance issue inherent in the vanilla SGD methods. Instead, these modern stochasticThese methods are randomized algorithms. However, the term “stochastic” (somewhat incorrectly) appears in their names for historical reasons, and quite possibly due to their aspiration to improve upon stochastic gradient descent (SGD). method are able to dramatically improve upon SGD in various different ways, but without having to resort to the usual variance-reduction techniques (such as decreasing stepsizes or mini-batching) which carry with them considerable costs drastically reducing their power. Instead, these modern methods were able to improve upon SGD without any unwelcome side effects. This development led to a revolution in the area of first order methods for solving problem (1)+(2). Both the theoretical complexity and practical efficiency of these modern methods vastly outperform prior gradient-type methods.

In order to achieve $\epsilon$-accuracy, that is,

modern stochastic methods such as SAG, SDCA, SVRG and S2GD require only

units of work, where $\kappa$ is a condition number associated with $F$, and one unit of work corresponds to the computation of the gradient of $f_{i}$ for a random index $i$, followed by a call to a prox-mapping involving $R$. More specifically, $\kappa=L/\mu$, where $L$ is a uniform bound on the Lipschitz constants of the gradients of functions $f_{i}$ and $\mu$ is the strong convexity constant of $P$. These quantities will be defined precisely in Section IV.

The complexity bound (4) should be contrasted with that of proximal gradient descent (e.g., ISTA), which requires $O(n\kappa\log(1/\epsilon))$ units of work, or FISTA, which requires $O(n\sqrt{\kappa}\log(1/\epsilon))$ units of workHowever, it should be remarked that the condition number $\kappa$ in these latter methods is slightly different from that appearing in the bound (4).. Note that while all these methods enjoy linear convergence rate, the modern stochastic methods can be many orders of magnitude faster than classical deterministic methods. Indeed, one can have

Based on this, we see that these modern methods always beat (proximal) gradient descent ($n+\kappa\ll n\kappa$), and also outperform FISTA as long as $\kappa\leq{\cal O}(n^{2})$. In machine learning, for instance, one usually has $\kappa\approx n$, in which case the improvement is by a factor of $\sqrt{n}$ when compared to FISTA, and by a factor of $n$ over ISTA. For applications where $n$ is massive, these improvements are indeed dramatic.

For more information about modern dual and primal methods we refer the reader to the literature on randomized coordinate descent methods and stochastic gradient methods , respectively.

I-C Linear systems and sketching.

In the case when $R\equiv 0$, all stationary points (i.e., points satisfying $\nabla F(x)=0$) are optimal for (1)+(2). In the special case when the functions $f_{i}$ are convex quadratics of the form $f_{i}(x)=\tfrac{1}{2}(a_{i}^{T}x-b_{i})$, the equation $\nabla F(x)=0$ reduces to the linear system $A^{T}Ax=A^{T}b$, where $A=[a_{1},\dots,n]$. Recently, there has been considerable interest in designing and analyzing randomized methods for solving linear systems; also known under the name of sketching methods. Much of this work was done independently from the developments in (non-quadratic) optimization, despite the above connection between optimization and linear systems. A randomized version of the classical Kaczmarz method was studied in a seminal paper by Strohmer and Vershynin . Subsequently, the method was extended and improved upon in several ways . The randomized Kaczmarz method is equivalent to SGD with a specific stepsize choice . The first randomized coordinate descent method, for linear systems, was analyzed by Lewis and Leventhal , and subsequently generalized in various ways by numerous authors (we refer the reader to and the references therein). Gower and Richtárik have recently studied randomized iterative methods for linear systems in a general sketch and project framework, which in special cases includes randomized Kaczmarz, randomized coordinate descent, Gaussian descent, randomized Newton, their block variants, variants with importance sampling, and also an infinite array of new specific methods. For approaches of a combinatorial flavour, specific to diagonally dominant systems, we refer to the influential work of Spielman and Teng .

II Contributions

In this paper we equip moderns stochastic methods—methods which already enjoy the fast rate (4)—with the ability to process data in mini-batches. None of the primalBy a primal method we refer to an algorithm which operates directly to solve (1)+(2) without explicitly operating on the dual problem. Dual methods have very recently been analyzed in the mini-batch setting. For a review of such methods we refer the reader to the paper describing the QUARTZ method and the references therein. modern methods have been analyzed in the mini-batch setting. This paper fills this gap in the literature.

While we have argued above that the modern methods, S2GD included, do not have the “non-vanishing variance” issue that SGD does, and hence do not need mini-batching for that purpose, mini-batching is still useful. In particular, we develop and analyze the complexity of mS2GD (Algorithm 1) — a mini-batch proximal variant of semi-stochastic gradient descent (S2GD) . While the S2GD method was analyzed in the $R=0$ case only, we develop and analyze our method in the proximalNote that the Prox-SVRG method can also handle the composite problem (1). setting (1). We show that mS2GD enjoys several benefits when compared to previous modern methods. First, it trivially admits a parallel implementation, and hence enjoys a speedup in clocktime in an HPC environment. This is critical for applications with massive datasets and is the main motivation and advantage of our method. Second, our results show that in order to attain a specified accuracy $\epsilon$, mS2GD can get by with fewer gradient evaluations than S2GD. This is formalized in Theorem 2, which predicts more than linear speedup up to a certain threshold mini-batch size after which the complexity deteriorates. Third, compared to , our method does not need to average the iterates produced in each inner loop; we instead simply continue from the last one. This is the approach employed in S2GD .

III The Algorithm

In this section we first briefly motivate the mathematical setup of deterministic and stochastic proximal gradient methods in Section IIIA, followed by the introduction of semi-stochastic gradient descent in Section IIIB. We will the be ready to describe the mS2GD method in Section IIIC.

The classical deterministic proximal gradient approach to solving (1) is to form a sequence $\{{y}_{t}\}$ via

where $U_{t}(x)\stackrel{{\scriptstyle\text{def}}}{{=}}F({y}_{t})+\nabla F({y}_{t})^{T}({x}-{y}_{t})+\tfrac{1}{2h}\|{x}-{y}_{t}\|^{2}+R(x)$. Note that in view of Assumption 1, which we shall use in our analysis in Section IV, $U_{t}$ is an upper bound on $P$ whenever $h>0$ is a stepsize parameter satisfying $1/h\geq L$. This procedure can be compactly written using the proximal operator as follows:

In a large-scale setting it is more efficient to instead consider the stochastic proximal gradient approach, in which the proximal operator is applied to a stochastic gradient step:

where $G_{t}$ is a stochastic estimate of the gradient $\nabla F(y_{t})$.

III-B Semi-stochastic methods

Of particular relevance to our work are the SVRG , S2GD and Prox-SVRG methods where the stochastic estimate of $\nabla F(y_{t})$ is of the form

where $x$ is an “old” reference point for which the gradient $\nabla F(x)$ was already computed in the past, and $i_{t}\in[n]\stackrel{{\scriptstyle\text{def}}}{{=}}\{1,2,\dots,n\}$ is a random index equal to $i$ with probability $q_{i}>0$. Notice that $G_{t}$ is an unbiased estimate of the gradient of $F$ at $y_{t}$:

Methods such as S2GD, SVRG, and Prox-SVRG update the points $y_{t}$ in an inner loop, and the reference point $x$ in an outer loop (“epoch”) indexed by $k$. With this new outer iteration counter we will have $x_{k}$ instead of $x$, $y_{k,t}$ instead of $y_{t}$ and $G_{k,t}$ instead of $G_{t}$. This is the notation we will use in the description of our algorithm in Section IIIC. The outer loop ensures that the squared norm of $G_{k,t}$ approaches zero as $k,t\to\infty$ (it is easy to see that this is equivalent to saying that the stochastic estimate $G_{k,t}$ has a diminishing variance), which ultimately leads to extremely fast convergence.

III-C Mini-batch S2GD

We are now ready to describe the mS2GD method A more detailed algorithm and the associated analysis (in which we benefit from the knowledge of lower-bound on the strong convexity parameters of the functions $F$ and $R$) can be found in the arXiv preprint . The more general algorithm mainly differs in $t_{k}$ being chosen according to a geometric probability law which depends on the estimates of the convexity constants. (Algorithm 1).

The algorithm includes an outer loop, indexed by epoch counter $k$, and an inner loop, indexed by $t$. Each epoch is started by computing $g_{k}$, which is the (full) gradient of $F$ at $x_{k}$. It then immediately proceeds to the inner loop. The inner loop is run for $t_{k}$ iterations, where $t_{k}$ is chosen uniformly at random from $\{1,\dots,m\}$. Subsequently, we run $t_{k}$ iterations in the inner loop (corresponding to Steps 6–10). Each new iterate is given by the proximal update (5), however with the stochastic estimate of the gradient $G_{k,t}$ in (6), which is formed by using a mini-batch of examples $A_{kt}\subseteq[n]$ of size $|A_{kt}|=b$. Each inner iteration requires $2b$ units of workIt is possible to finish each iteration with only $b$ evaluations for component gradients, namely $\{\nabla f_{i}(y_{k,t})\}_{i\in A_{kt}}$, at the cost of having to store $\{\nabla f_{i}(x_{k})\}_{i\in[n]}$, which is exactly the way that SAG works. This speeds up the algorithm; nevertheless, it is impractical for big $n$..

IV Analysis

In this section, we lay down the assumptions, state our main complexity result, and comment on how to optimally choose the parameters of the method.

Our analysis is performed under the following two assumptions.

Hence, the gradient of $F$ is also Lipschitz continuous with the same constant $L$.

$P$ is strongly convex with parameter $\mu>0$. That is for all ${x},{y}\in\operatorname{dom}(R)$ and $\xi\in\partial{P({x})}$,

where $\partial P({x})$ is the subdifferential of $P$ at ${x}$.

Lastly, by $\mu_{F}\geq 0$ and $\mu_{R}\geq 0$ we denote the strong convexity constants of $F$ and $R$, respectively. We allow both of these quantities to be equal to , which simply means that the functions are convex (which we already assumed above). Hence, this is not necessarily an additional assumption.

IV-B Main result

We are now ready to formulate our complexity result.

Let Assumptions 1 and 2 be satisfied, let ${x}_{*}\stackrel{{\scriptstyle\text{def}}}{{=}}\arg\min_{x}P({x})$ and choose $b\in\{1,2\dots,n\}$. Assume that $0<h\leq 1/L$, $4hL\alpha(b)<1$ and that $m,h$ are further chosen so that

where $\alpha(b)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{n-b}{b(n-1)}$. Then mS2GD has linear convergence in expectation with rate $\rho$:

Notice that for any fixed $b$, by properly adjusting the parameters $h$ and $m$ we can force $\rho$ to be arbitrarily small. Indeed, the second term can be made arbitrarily small by choosing $h$ small enough. Fixing the resulting $h$, the first term can then be made arbitrarily small by choosing $m$ large enough. This may look surprising, since this means that only a single outer loop ($k=1$) is needed in order to obtain a solution of any prescribed accuracy. While this is indeed the case, such a choice of the parameters of the method ($m$, $h$, $k$) would not be optimal – the resulting workload would to be too high as the complexity of the method would depend sublinearly on $\epsilon$. In order to obtain a logarithmic dependence on $1/\epsilon$, i.e., in order to obtain linear convergence, one needs to perform $k=O(\log(1/\epsilon))$ outer loops, and set the parameters $h$ and $m$ to appropriate values (generally, $h={\cal O}(1/L)$ and $m={\cal O}(\kappa)$).

IV-C Special cases: b=1𝑏1b=1 and b=n𝑏𝑛b=n

In the special case with $b=1$ (no mini-batching), we get $\alpha(b)=1$, and the rate given by (8) exactly recovers the rate achieved by Prox-SVRG (in the case when the Lipschitz constants of $\nabla f_{i}$ are all equal). The rate is also identical to the rate of S2GD (in the case of $R=0$, since S2GD was only analyzed in that case). If we set the number of outer iterations to $k=\lceil\log(1/\epsilon)\rceil$, choose the stepsize as $h=\tfrac{1}{(2+4e)L}$, where $e=\exp(1)$, and choose $m=43\kappa$, then the total workload of mS2GD for achieving (3) is $(n+43\kappa)\log(1/\epsilon)$ units of work. Note that this recovers the fast rate (4).

In the batch setting, that is when $b=n$, we have $\alpha(b)=0$ and hence $\rho=1/(mh\mu)$. By choosing $k=\lceil\log(1/\epsilon)\rceil$, $h=1/L$, and $m=2\kappa$, we obtain the rate ${\cal O}\left(n\kappa\log(1/\epsilon)\right)$. This is the standard rate of (proximal) gradient descent.

Hence, by modifying the mini-batch size $b$ in mS2GD, we interpolate between the fast rate of S2GD and the slow rate of GD.

IV-D Mini-batch speedup

In this section we will derive formulas for good choices of the parameter $m,h$ and $k$ of our method as a function of $b$. Hence, throughout this section we shall consider $b$ fixed.

Fixing $0<\rho<1$, it is easy to see that in order for $x_{k}$ to be an $\epsilon$-accurate solution (i.e., in order for (3) to hold), it suffices to choose $k\geq(1-\rho)^{-1}\log(\epsilon^{-1})$. Notice that the total workload mS2GD will do in order to arrive at $x_{k}$ is

units of work. If we now consider $\rho$ fixed (we may try to optimize for it later), then clearly the total workload is proportional to $m$. The free parameters of the method are the stepsize $h$ and the inner loop size $m$. Hence, in order to set the parameters so as to minimize the workload (i.e., optimize the complexity bound), we would like to (approximately) solve the optimization problem

Let $(h_{*}^{b},m_{*}^{b})$ denote the optimal pair (we highlight the dependence on $b$ as it will be useful). Note that if $m_{*}^{b}\leq m_{*}^{1}/b$ for some $b>1$, then mini-batching can help us reach the $\epsilon$-solution with smaller overall workload. The following theorem presents the formulas for $h_{*}^{b}$ and $m_{*}^{b}$.

IV-E Convergence rate

In this section we study the total workload of mS2GD in the regime of small mini-batch sizes.

Fix $\epsilon\in(0,1)$, choose the number of outer iterations equal to

and fix the target decrease in Theorem 2 to satisfy $\rho=\epsilon^{1/k}$. Further, pick a mini-batch size satisfying $1\leq b\leq 29$, let the stepsize $h$ be as in (34) and let $m$ be as in (33). Then in order for mS2GD to find $x_{k}$ satisfying (3), mS2GD needs at most

units of work, where $bm^{b}={\cal O}(\kappa)$, which leads to the overall complexity of

This result shows that as long as the mini-batch size is small enough, the total work performed by mS2GD is the same as in the $b=1$ case. If the $b$ updates can be performed in parallel, then this leads to linear speedup.

IV-F Comparison with Acc-Prox-SVRG

The Acc-Prox-SVRG method of Nitanda, which was not available online before the first version of this paper appeared on arXiv, incorporates both a mini-batch scheme and Nesterov’s acceleration . The author claims that when $b<\lceil b_{0}\rceil$, with the threshold $b_{0}$ defined as $\frac{8\sqrt{\kappa}n}{\sqrt{2}p(n-1)+8\sqrt{\kappa}}$, the overall complexity of the method is

This suggests that acceleration will only be realized when the mini-batch size is large, while for small $b$, Acc-Prox-SVRG achieves the same overall complexity, $\mathcal{O}\left((n+\kappa)\log(1/\epsilon)\right)$, as mS2GD.

We will now take a closer look at the theoretical results given by Acc-Prox-SVRG and mS2GD, for each $\epsilon\in(0,1)$. In particular, we shall numerically minimize the total work of mS2GD, i.e.,

over $\rho\in(0,1)$ and $h$ (compare this with (11)); and compare these results with similar fine-tuned quantities for Acc-Prox-SVRG.$m^{b}$ is the best choice of $m$ for Acc-Prox-SVRG and mS2GD, respectively. Meanwhile, $h$ is within the safe upper bounds for both methods.

Fig. 1 illustrates these theoretical complexity bounds for both ill-conditioned and well-conditioned data. With small-enough mini-batch size $b$, mS2GD is better than Acc-Prox-SVRG. However, for a large mini-batch size $b$, the situation reverses because of the acceleration inherent in Acc-Prox-SVRG.We have experimented with different values for $n,b$ and $\kappa$, and this result always holds. Plots with $b=64$ illustrate the cases where we cannot observe any differences between the methods.

Note however that accelerated methods are very prone to error accumulation. Moreover, it is not clear that an efficient implementation of Acc-Prox-SVRG is possible for sparse data. As shall show in the next section, mS2GD allows for such an implementation.

V Efficient implementation for sparse data

Let us make the following assumption about the structure of functions $f_{i}$ in (2).

Many functions of common practical interest satisfy this assumption including linear and logistic regression. Very often, especially for large scale datasets, the data are extremely sparse, i.e. the vectors $\{a_{i}\}$ contains many zeros. Let us denote the number of non-zero coordinates of $a_{i}$ by $\omega_{i}=\|a_{i}\|_{0}\leq d$ and the set of indexes corresponding to non-zero coordinates by $\operatorname{support}(a_{i})=\{j:a_{i}^{j}\neq 0\}$, where $a_{i}^{j}$ denotes the $j^{th}$ coordinate of vector $a_{i}$.

The regularization function $R$ is separable.

This includes the most commonly used regularization functions as $\frac{\lambda}{2}\|x\|^{2}$ or $\lambda\|x\|_{1}$.

Let us take a brief detour and look at the classical SGD algorithm with $R=0$. The update would be of the form

If evaluation of the univariate function $\phi^{\prime}_{i}$ takes $O(1)$ amount of work, the computation of $\nabla f_{i}$ will account for $O(\omega_{i})$ work. Then the update (12) would cost $O(\omega_{i})$ too, which implies that the classical SGD method can naturally benefit from sparsity of data.

Now, let us get back to the Algorithm 1. Even under the sparsity assumption and structural Assumption 3 the Algorithm 1 suggests that each inner iteration will cost $O(\omega+d)\sim O(d)$ because $g_{k}$ is in general fully dense and hence in Step 9 of Algorithm 1 we have to update all $d$ coordinates.

However, in this Section, we will introduce and describe the implementation trick which is based on “lazy/delayed” updates. The main idea of this trick is not to perform Step 9 of Algorithm 1 for all coordinates, but only for coordinates $j\in\cup_{i\in A_{kt}}\operatorname{support}(a_{i})$. The algorithm is described in Algorithm 2.

To explain the main idea behind the lazy/delayed updates, consider that it happened that during the fist $\tau$ iterations, the value of the fist coordinate in all datapoints which we have used was 0. Then given the values of $y_{k,0}^{1}$ and $g_{k}^{1}$ we can compute the true value of $y_{k,t}^{1}$ easily. We just need to apply the $\operatorname{prox}$ operator $\tau$ times, i.e. $y_{k,\tau}^{1}=\operatorname{prox}^{\tau}_{1}[y_{k,0},g_{k},R,h],$ where the function $\operatorname{prox}^{\tau}_{1}$ is described in Algorithm 3.

The vector $\chi$ in Algorithm 2 is enabling us to keep track of the iteration when corresponding coordinate of $y$ was updated for the last time. E.g. if in iteration $t$ we will be updating the $1^{st}$ coordinate for the first time, $\chi^{1}=0$ and after we compute and update the true value of $y^{1}$, its value will be set to $\chi^{1}=t$. Lines 5-8 in Algorithm 2 make sure that the coordinates of $y_{k,t}$ which will be read and used afterwards are up-to-date. At the end of the inner loop, we will updates all coordinates of $y$ to the most recent value (lines 12-14). Therefore, those lines make sure that the $y_{k,t_{k}}$ of Algorithms 1 and 2 will be the same.

However, one could claim that we are not saving any work, as when needed, we still have to compute the proximal operator many times. Although this can be true for a general function $R$, for particular cases, $R(x)=\frac{\lambda}{2}\|x\|^{2}$ and $R(x)=\lambda\|x\|_{1}^{2}$, we provide following Lemmas which give a closed form expressions for the $\operatorname{prox}^{\tau}_{j}$ operator.

If $R(x)=\frac{\lambda}{2}\|x\|^{2}$ with $\lambda>0$ then

where $\beta\stackrel{{\scriptstyle\text{def}}}{{=}}1/(1+\lambda h)$.

Assume that $R(x)=\lambda\|x\|_{1}$ with $\lambda>0$. Let us define $M$ and $m$ as follows,

and let $[\cdot]_{+}\stackrel{{\scriptstyle\text{def}}}{{=}}\max\{\cdot,0\}$. Then the value of $\operatorname{prox}_{j}^{\tau}[y,g,R,h]$ can be expressed based on one of the 3 situations described below:

If $g^{j}\geq\lambda$, then by letting $p\stackrel{{\scriptstyle\text{def}}}{{=}}\left\lfloor\frac{y^{j}}{M}\right\rfloor$, the operator can be defined as

If $-\lambda<g^{j}<\lambda$, then the operator can be defined as

If $g^{j}\leq-\lambda$, then by letting $q\stackrel{{\scriptstyle\text{def}}}{{=}}\left\lfloor\frac{y^{j}}{m}\right\rfloor$, the operator can be defined as

The proofs of Lemmas 1 and 2 are available in APPENDIX C.

Remark: Upon completion of the paper, we learned that similar ideas of lazy updates were proposed in and for online learning and multinomial logistic regression, respectively. However, our method can be seen as a more general result applied to a stochastic gradient method and its variants under Assumptions 3 and 4.

VI Experiments

In this section we perform numerical experiments to illustrate the properties and performance of our algorithm. In Section VI-A we study the total workload and parallelization speedup of mS2GD as a function of the mini-batch size $b$. In Section VI-B we compare mS2GD with several other algorithms. Finally, in Section VI-C we briefly illustrate that our method can be efficiently applied to a deblurring problem.

In Sections VI-A and VI-B we conduct experiments with $R(x)=\tfrac{\lambda}{2}\|x\|^{2}$ and $F$ of the form (2), where $f_{i}$ is the logistic loss function:

and is used in machine learning for binary classification. In these sections we have performed experiments on four publicly available binary classification datasets, namely rcv1, news20, covtype rcv1, covtype and news20 are available at http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/. and astro-ph Available at http://users.cecs.anu.edu.au/~xzhang/data/..

In the logistic regression problem, the Lipschitz constant of function $\nabla f_{i}$ is equal to $L_{i}=\|a_{i}\|^{2}/4$. Our analysis assumes (Assumption 1) the same constant $L$ for all functions. Hence, we have $L=\max_{i\in[n]}L_{i}$. We set the regularization parameter $\lambda=\frac{1}{n}$ in our experiments, resulting in the problem having the condition number $\kappa=\frac{L}{\mu}=\mathcal{O}(n)$. In Table I we summarize the four datasets, including the sizes $n$, dimensions $d$, their sparsity levels as a proportion of nonzero elements, and the Lipschitz constants $L$.

Mini-batches allow mS2GD to be accelerated on a computer with a parallel processor. In Section IV-D, we have shown in that up to some threshold mini-batch size, the total workload of mS2GD remains unchanged. Figure 2 compares the best performance of mS2GD used with various mini-batch sizes on datasets rcv1 and astro-ph. An effective pass (through the data) corresponds to $n$ units of work. Hence, the evaluation of a gradient of $F$ counts as one effective pass. In both cases, by increasing the mini-batch size to $b=2,4,8$, the performance of mS2GD is the same or better than that of S2GD ($b=1$) without any parallelism.

Although for larger mini-batch sizes mS2GD would be obviously worse, the results are still promising with parallelism. In Figure 3,we show the ideal speedup—one that would be achievable if we could always evaluate the $b$ gradients in parallel in exactly the same amount of time as it would take to evaluate a single gradient.In practice, it is impossible to ensure that the times of evaluating different component gradients are the same..

VI-B mS2GD vs other algorithms

In this part, we implemented the following algorithms to conduct a numerical comparison: 1) SGDcon: Proximal stochastic gradient descent method with a constant step-size which gave the best performance in hindsight. 2) SGD+: Proximal stochastic gradient descent with variable step-size $h=h_{0}/(k+1)$, where $k$ is the number of effective passes, and $h_{0}$ is some initial constant step-size. 3) FISTA: Fast iterative shrinkage-thresholding algorithm proposed in . 4) SAG: Proximal version of the stochastic average gradient algorithm . Instead of using $h=1/16L$, which is analyzed in the reference, we used a constant step size. 5) S2GD: Semi-stochastic gradient descent method proposed in . We applied proximal setting to the algorithm and used a constant stepsize. 6) mS2GD: mS2GD with mini-batch size $b=8$. Although a safe step-size is given in our theoretical analyses in Theorem 1, we ignored the bound, and used a constant step size.

In all cases, unless otherwise stated, we have used the best constant stepsizes in hindsight.

Figure 4 demonstrates the superiority of mS2GD over other algorithms in the test pool on the four datasets described above. For mS2GD, the best choices of parameters with $b=8$ are given in Table II.

VI-C Image deblurring

In this section we utilize the Regularization Toolbox Regularization Toolbox available for Matlab can be obtained from http://www.imm.dtu.dk/~pcha/Regutools/ . We use the blur function available therein to obtain the original image and generate a blurred image (we choose following values of parameters for blur function: $N=256$, band=9, sigma=10). The purpose of the blur function is to generate a test problem with an atmospheric turbulence blur. In addition, an additive Gaussian white noise with stand deviation of $10^{-3}$ is added to the blurred image. This forms our testing image as a vector $b$. The image dimension of the test image is $256\times 256$, which means that $n=d=65,536$. We would not expect our method to work particularly well on this problem since mS2GD works best when $d\ll n$. However, as we shall see, the method’s performance is on a par with the performance of the best methods in our test pool.

Our goal is to reconstruct (de-blur) the original image $x$ by solving a LASSO problem: $\min_{x}\|Ax-b\|_{2}^{2}+\lambda\|x\|_{1}$. We have chosen $\lambda=10^{-4}$. In our implementation, we normalized the objective function by $n$, and hence our objective value being optimized is in fact $\min_{x}\frac{1}{n}\|Ax-b\|_{2}^{2}+\lambda\|x\|_{1}$, where $\lambda=\frac{10^{-4}}{n}$, similarly as was done in .

Figure (5) shows the original test image (left) and a blurred image with added Gaussian noise (right). Figure 6 compares the mS2GD algorithm with SGD+, S2GD and FISTA. We run all algorithms for 100 epochs and plot the error. The plot suggests that SGD+ decreases the objective function very rapidly at beginning, but slows down after 10-20 epochs.

Finally, Fig. 7 shows the reconstructed image after $T=20,60,100$ epochs.

VII Conclusion

We have proposed mS2GD—a mini-batch semi-stochastic gradient method—for minimizing a strongly convex composite function. Such optimization problems arise frequently in inverse problems in signal processing and statistics. Similarly to SAG, SVRG, SDCA and S2GD, our algorithm also outperforms existing deterministic method such as ISTA and FISTA. Moreover, we have shown that the method is by design amenable to a simple parallel implementation. Comparisons to state-of-the-art algorithms suggest that mS2GD, with a small-enough mini-batch size, is competitive in theory and faster in practice than other competing methods even without parallelism. The method can be efficiently implemented for sparse data sets.

Appendix A Technical Results

Note that contractiveness of the proximal operator is a standard result in optimization literature .

Following from the proof of Corollary 3.5 in , by applying Lemma 4 with $\xi_{i}:=\nabla f_{i}(y_{k,t})-\nabla f_{i}({x}_{k})$, we have the bound for variance as follows.

Let $\alpha(b)\stackrel{{\scriptstyle\text{def}}}{{=}}\tfrac{n-b}{b(n-1)}$. Considering the definition of $G_{k,t}$ in Algorithm 1, conditioned on $y_{k,t}$, we have $\operatorname{\mathbf{E}}[G_{k,t}]=\nabla F(y_{k,t})$ and the variance satisfies,

Appendix B Proofs

As in the statement of the lemma, by $\operatorname{\mathbf{E}}[\cdot]$ we denote expectation with respect to the random set $\hat{S}$. First, note that

If we let $C\stackrel{{\scriptstyle\text{def}}}{{=}}\|\bar{\xi}\|^{2}=\frac{1}{n^{2}}\left(\sum_{i,j}\xi_{i}^{T}\xi_{j}\right)$, we can thus write

where in the last step we have used the bound $\tfrac{1}{n}\sum_{i,j}\xi_{i}^{T}\xi_{j}=n\left\|\sum_{i=1}^{n}\tfrac{1}{n}\xi_{i}\right\|^{2}\geq 0.$

B-B Proof of Theorem 1

The proof is following the steps in . For convenience, let us define the stochastic gradient mapping

then the iterate update can be written as $y_{k,t+1}=y_{k,t}-hd_{k,t}.$ Let us estimate the change of $\|y_{k,t+1}-{x}_{*}\|$. It holds that

Applying Lemma 3.7 in (this is why we need to assume that $h\leq 1/L$) with ${x}=y_{k,t}$, $v=G_{k,t}$, ${x}^{+}=y_{k,t+1}$, $g=d_{k,t}$, ${y}={x}_{*}$ and $\Delta=\Delta_{k,t}=G_{k,t}-\nabla F(y_{k,t})$, we get

In order to bound $-\Delta_{k,t}^{T}(y_{k,t+1}-{x}_{*})$, let us define the proximal full gradient update asNote that this quantity is never computed during the algorithm. We can use it in the analysis nevertheless. $\bar{y}_{k,t+1}=\operatorname{prox}_{hR}(y_{k,t}-h\nabla F(y_{k,t})).$ We get

Using Cauchy-Schwarz and Lemma 3, we conclude that

Further, we obtain $\|y_{k,t+1}-{x}_{*}\|^{2}\overset{\eqref{eq:asfaevfwavfa},\eqref{eq:safawvfdwavfda}}{\leq}\left\|y_{k,t}-{x}_{*}\right\|^{2}+2h\left(h\|\Delta_{k,t}\|^{2}-\Delta_{k,t}^{T}(\bar{y}_{k,t+1}-{x}_{*})-[P(y_{k,t+1})-P({x}_{*})]\right).$ By taking expectation, conditioned on $y_{k,t}$For simplicity, we omit the $\operatorname{\mathbf{E}}[\cdot\,|\,y_{k,t}]$ notation in further analysis we obtain

where we have used that $\operatorname{\mathbf{E}}[\Delta_{k,t}]=\operatorname{\mathbf{E}}[G_{k,t}]-\nabla F(y_{k,t})=0$ and hence $\operatorname{\mathbf{E}}[-\Delta_{k,t}^{T}(\bar{y}_{k,t+1}-{x}_{*})]=0$$\bar{y}_{k,t+1}$ is constant, conditioned on $y_{k,t}$. Now, if we substitute (16) into (21) and decrease index $t$ by 1, we obtain

and $\alpha(b)=\frac{n-b}{b(n-1)}$. Note that (22) is equivalent to

Now, by the definition of ${x}_{k}$ in Algorithm 1 we have that

By summing (24) for $1\leq t\leq m$, we get on the left hand side

Combining (26) and (27) and using the fact that $LHS\leq RHS$, we have

Since $\operatorname{\mathbf{E}}\|{y}_{k,m}-{x}_{*}\|^{2}\geq 0$ and ${y}_{k,0}={x}_{k}$, by combining (29) and (28) we get

Notice that in view of our assumption on $h$ and (23), we have $2h>\theta$, and hence

where $\rho=\tfrac{2}{m\mu(2h-\theta)}+\tfrac{\theta(m+1)}{m(2h-\theta)}.$ Applying the above linear convergence relation recursively with chained expectations, we finally obtain $\operatorname{\mathbf{E}}[P({x}_{k})-P(x_{*})]\leq\rho^{k}[P({x}_{0})-P(x_{*})].$

B-C Proof of Theorem 2

Clearly, if we choose some value of $h$ then the value of $m$ will be determined from (8) (i.e. we need to choose $m$ such that we will get desired rate). Therefore, $m$ as a function of $h$ obtained from (8) is

Now, we can observe that the nominator is always positive and the denominator is positive only if $\rho>4\alpha(b)hL(\rho+1)$, which implies $\tfrac{1}{4\alpha(b)L}\cdot\tfrac{\rho}{\rho+1}>h$ (note that $\tfrac{\rho}{\rho+1}\in[0,\frac{1}{2}]$). Observe that this condition is stronger than the one in the assumption of Theorem 1. It is easy to verify that

Also note that $m(h)$ is differentiable (and continuous) at any $h\in(0,\frac{1}{4\alpha(b)L}\cdot\frac{\rho}{\rho+1})=:I_{h}$. The derivative of $m$ is given by $m^{\prime}(h)=\tfrac{-\rho+4\alpha(b)hL(2+(2+h\mu)\rho)}{h^{2}\mu(\rho-4\alpha(b)hL(1+\rho))^{2}}.$ Observe that $m^{\prime}(h)$ is defined and continuous for any $h\in I_{h}$. Therefore there have to be some stationary points (and in case that there is just on $I_{h}$) it will be the global minimum on $I_{h}$. The FOC gives

which is equivalent to $\mu\rho^{2}+4\alpha(b)L(1+\rho)^{2}>2(1+\rho)\sqrt{\alpha(b)L(\mu\rho^{2}+4\alpha(b)L(1+\rho)^{2})}.$ Because both sides are positive, we can square them to obtain the equivalent condition

Claim #2

B-D Proof of Corollary 3

Since $k=\left\lceil\log(1/\epsilon)\right\rceil\geq\log(1/\epsilon)$ if and only if $e^{k}\geq 1/\epsilon$ if and only if $e^{-k}\leq\epsilon,$ we can conclude that $\rho^{k}\overset{\eqref{eqn:convergencerate}}{=}(e^{-1})^{k}=e^{-k}\leq\epsilon$. Therefore, running mS2GD for $k$ outer iterations achieves $\epsilon$-accuracy solution defined in (3). Moreover, since in general $\kappa\gg e,n\gg e$, it can be concluded that

then with the definition $\alpha(b)=\frac{(n-b)}{b(n-1)}$, we derive

so from (11), the total complexity can be translated to $\mathcal{O}\left((n+\kappa)\log(1/\epsilon)\right).$ This result shows that we can reach efficient speedup by mini-batching as long as the mini-batch size is smaller than some threshold $b_{0}\approx 29.75$, which finishes the proof for Corollary 3.

C-B Proof of Lemma 2

For any $s\in\{1,2,\dots,\tau\}$ and $j\in\{1,2,\dots,d\}$,

where $M\stackrel{{\scriptstyle\text{def}}}{{=}}(\lambda+g^{j})h$, $m\stackrel{{\scriptstyle\text{def}}}{{=}}-(\lambda-g^{j})h$ and $M-m=2\lambda h>0$. Now, we will distinguish several cases based on $g^{j}$:

When $-\lambda<g^{j}<\lambda$, then $M=(\lambda+g^{j})h>0,m=-(\lambda-g^{j})h<0$, thus we have that

When $g^{j}\leq-\lambda$, then $m<M=(\lambda+g^{j})h\leq 0$, thus by letting $q=\left\lfloor\frac{y^{j}}{m}\right\rfloor$, we have that: if $y^{j}\leq m$, then

which summarizes the situation as follows:

References