1

Larger is Better: The Effect of Learning Rates
Enjoyed by Stochastic Optimization with
Progressive Variance Reduction
Fanhua Shang

arXiv:1704.04966v1 [cs.LG] 17 Apr 2017

Department of Computer Science and Engineering, The Chinese University of Hong Kong
fhshang@cse.cuhk.edu.hk
April 18, 2017
Abstract
In this paper, we propose a simple variant of the original stochastic variance reduction gradient (SVRG) [1], where
hereafter we refer to as the variance reduced stochastic gradient descent (VR-SGD). Different from the choices of the
snapshot point and starting point in SVRG and its proximal variant, Prox-SVRG [2], the two vectors of each epoch in VRSGD are set to the average and last iterate of the previous epoch, respectively. This setting allows us to use much larger
learning rates or step sizes than SVRG, e.g., 3/(7L) for VR-SGD vs. 1/(10L) for SVRG, and also makes our convergence
analysis more challenging. In fact, a larger learning rate enjoyed by VR-SGD means that the variance of its stochastic
gradient estimator asymptotically approaches zero more rapidly. Unlike common stochastic methods such as SVRG and
proximal stochastic methods such as Prox-SVRG, we design two different update rules for smooth and non-smooth objective
functions, respectively. In other words, VR-SGD can tackle non-smooth and/or non-strongly convex problems directly without
using any reduction techniques such as quadratic regularizers. Moreover, we analyze the convergence properties of VR-SGD
for strongly convex problems, which show that VR-SGD attains a linear convergence rate. We also provide the convergence
guarantees of VR-SGD for non-strongly convex problems. Experimental results show that the performance of VR-SGD is
significantly better than its counterparts, SVRG and Prox-SVRG, and it is also much better than the best known stochastic
method, Katyusha [3].
Index Terms
Stochastic optimization, stochastic gradient descent (SGD), proximal stochastic gradient, variance reduction, iterate
averaging, snapshot and starting points, strongly convex and non-strongly convex, smooth and non-smooth

F

•

All the codes of VR-SGD and some related variance reduced stochastic methods can be downloaded from the author’s website: https://sites.google.com/site/
fanhua217/publications.

2

1

I NTRODUCTION

In this paper, we focus on the following composite convex optimization problem:
n
X
def 1
fi (x) + g(x)
(1)
min F (x) =
n i=1
x∈Rd
Pn
where f (x) := n1 i=1 fi (x), fi (x) : Rd → R, i = 1, . . . , n are the smooth convex functions, and g(x) is a relatively simple
(but possibly non-differentiable) convex function (referred to as a regularizer). The formulation (1) arises in many places
in machine learning, signal processing, data science, statistics and operations research, such as regularized empirical risk
minimization (ERM). For instance, one popular choice of the component function fi (·) in binary classification problems is
the logistic loss, i.e., fi (x) = log(1 + exp(−bi aTi x)), where {(a1 , b1 ), . . . , (an , bn )} is a collection of training examples, and

bi ∈ {±1}. Some popular choices for the regularizer include the `2 -norm regularizer (i.e., g(x) = (λ1 /2)kxk2 ), the `1 -norm
regularizer (i.e., g(x) = λ2 kxk1 ), and the elastic-net regularizer (i.e., g(x) = (λ1 /2)kxk2 +λ2 kxk1 ), where λ1 ≥ 0 and λ2 ≥ 0
are two regularization parameters. So far examples of some other applications include deep neural networks [1], [4], [5],
[6], group Lasso [7], [8], sparse learning and coding [9], [10], phase retrieval [11], matrix completion [12], [13], conditional
random fields [14], eigenvector computation [15], [16] such as principal component analysis (PCA) and singular value
decomposition (SVD), generalized eigen-decomposition and canonical correlation analysis (CCA) [17].
1.1

Stochastic Gradient Descent

In this paper, we are especially interested in developing efficient algorithms to solve regularized ERM problems involving a
large sum of n component functions. The standard and effective method for solving Problem (1) is the (proximal) gradient
descent (GD) method, including accelerated proximal gradient (APG) [18], [19], [20], [21]. For smooth objective functions,
the update rule of GD is
"

xk+1 = xk − ηk

n

1X
∇fi (xk ) + ∇g(xk )
n i=1

#
(2)

for k = 1, 2, . . ., where ηk > 0 is commonly referred to as the step-size in optimization or the learning rate in machine
learning. When the regularizer g(·) is non-smooth, e.g., the `1 -norm regularizer, we need to introduce the following proximal
operator into (2),

xk+1 = Prox ηk ,g (yk ) := arg min (1/2ηk )·kx − yk k2 + g(x)

(3)

x∈Rd

where yk = xk −(ηk /n)

Pn

i=1 ∇fi (xk ).

The GD methods mentioned above have been proven to achieve linear convergence

for strongly convex problems, and APG attains the optimal convergence rate of O(1/T 2 ) for non-strongly convex problems,
where T denotes the number of iterations. However, the per-iteration cost of all the batch (or deterministic) methods is

O(nd), which is expensive.
Instead of evaluating the full gradient of f (·) at each iteration, an effective alternative is the stochastic (or incremental)
gradient descent (SGD) method [22]. SGD only evaluates the gradient of a single component function at each iteration, thus
it has much lower per-iteration cost, O(d), and has been successfully applied to many large-scale learning problems [4],
[23], [24], [25]. The update rule of SGD is formulated as follows:

xk+1 = xk − ηk [∇fik(xk ) + ∇g(xk )]

(4)

where ηk ∝ 1/k , and the index ik is chosen uniformly at random from {1, . . . , n}. Although the expectation of ∇fik(xk ) is
an unbiased estimation for ∇f (xk ), i.e., E[∇fik(xk )] = ∇f (xk ), the variance of the stochastic gradient estimator ∇fik(xk ) may
be large due to the variance of random sampling [1]. Thus, stochastic gradient estimators are also called “noisy gradients”,
and we need to gradually reduce its step size, leading to slow convergence. In particular, even under the strongly convex
condition, standard SGD attains a slower sub-linear convergence rate of O(1/T ) [26], [27].

3

1.2

Accelerated SGD

Recently, many SGD methods with variance reduction techniques were proposed, such as stochastic average gradient
(SAG) [28], stochastic variance reduced gradient (SVRG) [1], stochastic dual coordinate ascent (SDCA) [29], SAGA [30],
stochastic primal-dual coordinate (SPDC) [31], and their proximal variants, such as Prox-SAG [32], Prox-SVRG [2] and
Prox-SDCA [33]. All these accelerated SGD methods can use a constant step size η instead of diminishing step sizes for
SGD, and fall into the following three categories: primal methods such as SVRG and SAGA, dual methods such as SDCA,
e or
and primal-dual methods such as SPDC. In essence, many of primal methods use the full gradient at the snapshot x
average gradients to progressively reduce the variance of stochastic gradient estimators, as well as dual and primal-dual
methods, which leads to a revolution in the area of first-order methods [34]. Thus, they are also known as the hybrid gradient
descent method [35] or semi-stochastic gradient descent method [36]. In particular, under the strongly convex condition, the
accelerated SGD methods enjoy linear convergence rates and the overall complexity of O((n+L/µ) log(1/)) to obtain an

-suboptimal solution, where each fi (·) is L-smooth and g(·) is µ-strongly convex. The complexity bound shows that they
p
converge significantly faster than deterministic APG methods, whose complexity is O((n L/µ) log(1/)) [36].
SVRG [1] and its proximal variant, Prox-SVRG [2], are particularly attractive because of their low storage requirement
compared with other stochastic methods such as SAG, SAGA and SDCA, which require to store all the gradients of the n
e) is computed
component functions fi (·) or dual variables. At the beginning of each epoch in SVRG, the full gradient ∇f (x
e, which is updated periodically. The update rule for the smooth optimization problem (1) is given by
at the snapshot point x
e i (xk ) = ∇fi (xk ) − ∇fi (x
e) + ∇f (x
e),
∇f
k
k
k

(5a)

e i (xk ) + ∇g(xk )].
xk+1 = xk − η[∇f
k

(5b)

e i (xk ). It is not hard to verify
When g(·) ≡ 0, the update rule in (5b) becomes the original one in [1], i.e., xk+1 = xk −η ∇f
k
2
e
e
that the variance of the SVRG estimator ∇fi (xk ), i.e., Ek∇fi (xk )−∇f (xk )k , can be much smaller than that of the SGD
k

k

2

estimator ∇fik(xk ), i.e., Ek∇fik(xk )−∇f (xk )k . However, for non-strongly convex problems, the accelerated SGD methods
mentioned above converge much slower than batch APG methods such as FISTA [21], namely, O(1/T ) vs. O(1/T 2 ).
More recently, many acceleration techniques were proposed to further speed up those variance-reduced stochastic
methods mentioned above. These techniques mainly include the Nesterov’s acceleration technique used in [24], [37], [38],
[39], [40], reducing the number of gradient calculations in the early iterations [34], [41], [42], the projection-free property of
the conditional gradient method (also known as the Frank-Wolfe algorithm [43]) as in [44], the stochastic sufficient decrease
technique [45], and the momentum acceleration trick in [3], [34], [46]. More specifically, [40] proposed an accelerating
p
Catalyst framework and achieved the complexity of O((n+ nL/µ) log(L/µ) log(1/)) for strongly convex problems. [3]
p
and [46] proved that their accelerated methods can attain the best known complexity of O(n log(1/)+ nL/) for nonstrongly convex problems. The overall complexity matches the theoretical upper bound provided in [47]. Katyusha [3] and
p
point-SAGA [48] achieve the best-known complexity of O((n+ nL/µ) log(1/)) for strongly convex problems, which is
identical to the upper complexity bound in [47]. That is, Katyusha is the best known stochastic optimization method for both
strongly convex and non-strongly convex problems. Its proximal gradient update rules are formulated as follows:
e + (1 − w1 − w2 )zk ,
xk+1 = w1 yk + w2 x


1
2
T e
yk+1 = arg min
ky − yk k + y ∇fik(xk+1 ) + g(y) ,
2η
y∈Rd


3L
2
T e
kz − xk+1 k + z ∇fik(xk+1 ) + g(z)
zk+1 = arg min
2
z∈Rd

(6a)
(6b)
(6c)

where w1 , w2 ∈ [0, 1] are two momentum parameters. To eliminate the need for parameter tuning, η is set to 1/(3w1 L),
and w2 is fixed to 0.5 in [3]. Unfortunately, most of the accelerated methods mentioned above, including Katyusha, require
at least two auxiliary variables and two momentum parameters, which lead to complicated algorithm design and high
per-iteration complexity [34].

4

1.3

Our Contributions

From the above discussion, one can see that most of accelerated stochastic variance reduction methods such as [3], [34],
[37], [42], [44], [45] and applications such as [9], [10], [13], [15], [16], [49] are based on the stochastic variance reduced gradient
(SVRG) method [1]. Thus, any key improvement on SVRG is very important for the research of stochastic optimization.
In this paper, we propose a simple variant of the original SVRG [1], which is referred to as the variance reduced stochastic
gradient descent (VR-SGD). The snapshot point and starting point of each epoch in VR-SGD are set to the average and last
iterate of the previous epoch, respectively. Different from the settings of SVRG and Prox-SVRG [2] (i.e., the last iterate for
the two points of the former, while the average of the previous epoch for those of the latter), the two points in VR-SGD are
different, which makes our convergence analysis more challenging than SVRG and Prox-SVRG. Our empirical results show
that the performance of VR-SGD is significantly better than its counterparts, SVRG and Prox-SVRG. Impressively, VR-SGD
with a sufficiently large learning rate performs much better than the best known stochastic method, Katyusha [3]. The main
contributions of this paper are summarized below.
•

The snapshot point and starting point of VR-SGD are set to two different vectors. That is, for all epochs, except the
s+1
1 Pm
1 Pm−1 s
s
es = m
es = m−1
first one, x
= xsm . In
k=1 xk (denoted by Option I) or x
k=1 xk (denoted by Option II), and x0
particular, we find that the setting of VR-SGD allows us take much larger learning rates or step sizes than SVRG,
e.g., 3/(7L) vs. 1/(10L), and thus significantly speeds up the convergence of SVRG and Prox-SVRG in practice.
Moreover, VR-SGD has an advantage over SVRG in terms of robustness of learning rate selection.

•

Different from proximal stochastic gradient methods, e.g., Prox-SVRG and Katyusha, which have a unified update
rule for the two cases of smooth and non-smooth objectives (see Section 2.2 for details), VR-SGD employs two different
update rules for the two cases, respectively, as in (12) and (13) below. Empirical results show that gradient update
rules as in (12) for smooth optimization problems are better choices than proximal update formulas as in (10).

•

Finally, we theoretically analyze the convergence properties of VR-SGD with Option I or Option II for strongly convex
problems, which show that VR-SGD attains a linear convergence rate. We also give the convergence guarantees of
VR-SGD with Option I or Option II for non-strongly convex objective functions.

2

P RELIMINARY AND R ELATED W ORK

Throughout this paper, we use k · k to denote the `2 -norm (also known as the standard Euclidean norm), and k·k1 is the
P
`1 -norm, i.e., kxk1 = di=1 |xi |. ∇f (·) denotes the full gradient of f (·) if it is differentiable, or ∂f (·) the subgradient if f (·)
is only Lipschitz continuous. For each epoch s ∈ [S] and inner iteration k ∈ {0, 1, . . . , m−1}, isk ∈ [n] is the random chosen
index. We mostly focus on the case of Problem (1) when each component function fi (·) is L-smooth1 , and F (·) is µ-strongly
convex. The two common assumptions are defined as follows.
2.1

Basic Assumptions

Assumption 1 (Smoothness). Each convex function fi (·) is L-smooth, that is, there exists a constant L > 0 such that for all

x, y ∈ Rd ,
k∇fi (x) − ∇fi (y)k ≤ Lkx − yk.

(7)

Assumption 2 (Strong Convexity). The convex function F (x) is µ-strongly convex, i.e., there exists a constant µ > 0 such that for
all x, y ∈ Rd ,

µ
kx − yk2 .
(8)
2
Note that when the regularizer g(·) is non-smooth, the inequality in (8) needs to be revised by simply replacing the
F (y) ≥ F (x) + h∇F (x), y − xi +

gradient ∇F (x) in (8) with an arbitrary sub-gradient of F (·) at x. In contrast, for a non-strongly convex or general convex
function, the inequality in (8) can always be satisfied with µ = 0.
1. Actually, we can extend all the theoretical results in this paper for the case, when the gradients of all component functions have the same
Lipschitz constant L, to the more general case, when some fi (·) have different degrees of smoothness.

5

Algorithm 1 SVRG (Option I) and Prox-SVRG (Option II)
Input: The number of epochs S , the number of iterations m per epoch, and step size η .
e0 .
Initialize: x
1: for

s = 1, 2, . . . , S do

3:

xs0
es

4:

for k = 0, 1, . . . , m − 1 do

2:

5:
6:
7:

8:
9:
10:
11:

es−1 ;
=x
P
es−1 );
µ = n1 ni=1 ∇fi (x

% Initiate the variable xs0
% Compute the full gradient

Pick isk uniformly at random from [n];
e is (xs ) = ∇fis (xs ) − ∇fis (x
es−1 ) + µ
es ;
∇f
k
k
k
k
h k
i
e is (xs ) +∇g(xs ) ,
Option I: xsk+1 = xsk − η ∇f
k
 k k

e is (xs ) ;
or xsk+1 = Prox η, g xsk − η ∇f
k
k
n
e is (xs ) +
Option II: xs = arg miny∈Rd g(y) + y T ∇f
k+1

k

k

% The stochastic gradient estimator
% Smooth case of g(·)
% Non-smooth case of g(·)
1
2η ky

o
− xsk k2 ;

% Proximal update

end for
es = xsm ;
Option I: x
1 Pm
s
es = m
Option II: x
k=1 xk ;

e
% Last iterate for snapshot x
e
% Iterate averaging for snapshot x

12: end for

eS
Output: x

2.2

Related Work

To speed up standard and proximal SGD methods, many variance reduced stochastic methods [28], [29], [30], [35] have
been proposed for some special cases of Problem (1). In the case when each fi (x) is L-smooth, f (x) is µ-strongly convex,
and g(x) ≡ 0, Roux et al. [28] proposed a stochastic average gradient (SAG) method, which attains a linear convergence rate.
However, SAG needs to store all gradients as well as other incremental aggregated gradient methods such as SAGA [30],
so that O(nd) storage is required in general problems [41]. Similarly, SDCA [29] requires storage of all dual variables [1],
which scales as O(n). In contrast, SVRG [1], as well as its proximal variant, Prox-SVRG [2], has the similar convergence
rate to SAG and SDCA but without the memory requirements of all gradients and dual variables. In particular, the SVRG
estimator in (5a) (independently introduced in [1], [35]) may be the most popular choice for stochastic gradient estimators.
Besides, other stochastic gradient estimators include the SAGA estimator in [30] and the stochastic recursive gradient
estimator in [50]. Although the original SVRG in [1] only has convergence guarantees for a special case of Problem (1),
when each fi (x) is L-smooth, f (x) is µ-strongly convex, and g(x) ≡ 0, one can extend SVRG to the proximal setting by
introducing the proximal operator in (3), as shown in Line 7 of Algorithm 1. In other words, when g(·) is non-smooth, the
update rule of SVRG becomes
n
o
e is (xs )]k2 + g(x) .
xsk+1 = arg min (1/2η)·kx − [xsk − η ∇f
k
k

(9)

x∈Rd

Some researchers [8], [51], [52] have borrowed some variance reduction techniques into ADMM for minimizing convex
composite objective functions subject to an equality constraint. [15], [16], [17] applied efficient stochastic solvers to compute
leading eigenvectors of a symmetric matrix or generalized eigenvectors of two symmetric matrices. The first such method
is VR-PCA by Shamir [15], and the convergence properties of the VR-PCA algorithm for such a non-convex problem are
also provided. Garber et al. [16] analyzed the convergence rate of SVRG when f (·) is a convex function that is a sum
of non-convex component functions. Moreover, [6] and [53] proved that SVRG and SAGA with minor modifications can
converge asymptotically to a stationary point of non-convex finite-sum problems. Some distributed variants [54], [55] of
accelerated SGD methods have also been proposed.
An important class of stochastic methods is the proximal stochastic gradient (Prox-SG) method, such as Prox-SVRG [2],
SAGA [30], and Katyusha [3]. Different from standard variance reduction SGD methods such as SVRG, which have a

6

stochastic gradient update as in (5b), the Prox-SG method has a unified update rule for both smooth and non-smooth cases
of g(·). For instance, the update rule of Prox-SVRG [2] is formulated as follows:


e is (xs ) + 1 ky − xs k2 .
xsk+1 = arg min g(y) + y T ∇f
k
k
k
2η
y∈Rd

(10)

For the sake of completeness, the details of Prox-SVRG [2] are shown in Algorithm 1 with Option II. When g(·) is the
widely used `2 -norm regularizer, i.e., g(·) = (λ1 /2)k · k2 , the proximal update formula in (10) becomes
h
i
1
e is (xs ) .
xsk − η ∇f
xsk+1 =
k
k
1 + λ1 η

3

(11)

VARIANCE -R EDUCED S TOCHASTIC G RADIENT D ESCENT

In this section, we propose an efficient variance reduced stochastic gradient descent (VR-SGD) method with iterate averaging.
Different from the choices of the snapshot and starting points in SVRG [1] and Prox-SVRG [2], the two vectors of each
epoch in VR-SGD are set to the average and last iterate of the previous epoch, respectively. Unlike common stochastic
gradient methods such as SVRG and proximal stochastic gradient methods such as Prox-SVRG, we design two different
update rules for smooth and non-smooth objective functions, respectively.
3.1

Iterate Averaging

Like SVRG and Katyusha, VR-SGD is also divided into S epochs, and each epoch consists of m stochastic gradient steps,
es )
where m is usually chosen to be Θ(n) as in [1], [2], [3]. Within each epoch, we need to compute the full gradient ∇f (x
s
s
e i (x ) as in [1]. Unlike SVRG
e and use it to define the variance reduced stochastic gradient estimator ∇f
at the snapshot x
k

k

es of each epoch in VR-SGD is set to
whose snapshot point is set to the last iterate of the previous epoch, the snapshot x
P
m−1 s
1
es = m−1
the average of the previous epoch, e.g., x
k=1 xk in Option I of Algorithm 2, which leads to better robustness to
1 Pm
2
s
es = m
gradient noise , as also suggested in [34], [45], [59]. In fact, the choice of Option II in Algorithm 2, i.e., x
k=1 xk ,
also works well in practice, as shown in Fig. 1. Therefore, we provide the convergence guarantees for both our algorithms
(including Algorithm 2) with Option I and our algorithms with Option II in the next section. In particular, we find that
the one of the effects of the choice in Option I or Option II of Algorithm 2 is to allow taking much larger learning rates
or step sizes than SVRG in practice, e.g., 3/(7L) for VR-SGD vs. 1/(10L) for SVRG (see Fig. 2 and Section 5.3 for details).
This is the main reason why VR-SGD converges significantly faster than SVRG. Actually, a larger learning rate enjoyed by
VR-SGD means that the variance of its stochastic gradient estimator goes asymptotically to zero faster.
Unlike Prox-SVRG [2] whose starting point is initialized to the average of the previous epoch, the starting point xs+1
0
of each epoch in VR-SGD is set to the last iterate xsm of the previous epoch. That is, the last iterate of the previous
epoch becomes the new starting point in VR-SGD, while those of Prox-SVRG are completely different, thereby leading
to relatively slow convergence in general. It is clear that both the starting point and snapshot point of each epoch in the
original SVRG [1] are set to the last iterate of the previous epoch3 , while the two points of Prox-SVRG [2] are set to the
average of the previous epoch (also suggested in [1]). Different from the settings in SVRG and Prox-SVRG, the starting and
snapshot points in VR-SGD are set to the two different vectors mentioned above, which makes the convergence analysis of
VR-SGD more challenging than SVRG and Prox-SVRG, as shown in Section 4.
2. It should be emphasized that the noise introduced by random sampling is inevitable, and generally slows down the convergence speed in
a sense. However, SGD and its variants are probably the most used optimization algorithms for deep learning [56]. In particular, [57] has shown
that by adding noise at each step, noisy gradient descent can escape the saddle points efficiently and converge to a local minimum of non-convex
optimization problems, as well as the application of deep neural networks in [58].
1
3. Note that the theoretical convergence of the original SVRG [1] relies on its Option II, i.e., both x
es and xs+
are set to xsk , where k is randomly
0
chosen from {1, 2, . . . , m}. However, the empirical results in [1] suggest that Option I is a better choice than its Option II, and the convergence
guarantee of SVRG with Option I for strongly convex objective functions is provided in [60].

7

Algorithm 2 VR-SGD for strongly convex objectives
Input: The number of epochs S , the number of iterations m per epoch, and step size η .
e0 .
Initialize: x10 = x
2:

s = 1, 2, . . . , S do
P
es = n1 n
es−1 );
µ
i=1 ∇fi (x

3:

for k = 0, 1, . . . , m − 1 do

1: for

4:
5:
6:

7:
8:
9:
10:

% Compute the full gradient

Pick isk uniformly at random from [n];
e is (xs ) = ∇fis (xs ) − ∇fis (x
es−1 ) + µ
es ;
∇f
k
k
k
i
hk k
e is (xs ) + ∇g(xs ) ,
xsk+1 = xsk − η ∇f
k
k
 k

e is (xs ) ;
or xsk+1 = Prox η, g xsk − η ∇f
k
k

% The stochastic gradient estimator
% Smooth case of g(·)
% Non-smooth case of g(·)

end for
1 Pm−1 s
k=1 xk ;
m−1
1 Pm
s
e = m k=1 xsk ;
Option II: x
xs+1
= xsm ;
0

es =
Option I: x

e
% Iterate averaging for snapshot x
e
% Iterate averaging for snapshot x

% Initiate xs+1
for the next epoch
0

11: end for

PS

es ),
s=1 x

bS = S1
and x

es
s=1 x

VR−SGD−I
VR−SGD−II
VR−SGD−I
VR−SGD−II

−4

10
F (xs ) − F (x* )

PS

−6

10

−8

10

otherwise.

−4

−10

−6

10

−8

10

−10

10

10

−12

10

VR−SGD−I
VR−SGD−II
VR−SGD−I
VR−SGD−II

10
F (xs ) − F (x* )

bS = x
eS if F (x
eS ) ≤ F ( S1
Output: x

−12

0

2

4
6
Gradient evaluations / n

8

10

10

0

(a) Logistic regression

2

4
6
Gradient evaluations / n

8

10

(b) Ridge regression

Fig. 1. Comparison of VR-SGD with Option I (denoted by VR-SGD-I) and Option II (denoted by VR-SGD-II) for solving `2 -norm (i.e., (λ/2)k · k2 )
regularized logistic regression and ridge regression problems on the Covtype data set. In each plot, the vertical axis shows the objective value
minus the minimum, and the horizontal axis is the number of effective passes. Note that the blue lines stand for the results where λ = 10−4 , while
the red lines correspond to the results where λ = 10−5 (best viewed in colors).

3.2

The Algorithm for Strongly Convex Objectives

In this part, we propose an efficient VR-SGD algorithm to solve strongly convex objective functions, as outlined in Algorithm
2. It is well known that the original SVRG [1] only works for the case of smooth and strongly convex objective functions.
However, in many machine learning applications, e.g., elastic net regularized logistic regression, the strongly convex
objective function F (x) is non-smooth. To solve this class of problems, the proximal variant of SVRG, Prox-SVRG [2],
was subsequently proposed. Unlike SVRG and Prox-SVRG, VR-SGD can not only solve smooth objective functions, but
directly tackle non-smooth ones. That is, when the regularizer g(x) is smooth, e.g., the `2 -norm regularizer, the update rule
of VR-SGD is
e is (xs ) + ∇g(xs )].
xsk+1 = xsk − η[∇f
k
k
k

(12)

8

η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L
η=0.5/L
η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L
η=0.5/L
η=6/5L

−4

s

*

F (x ) − F (x )

10

−6

10

−8

10

−10

10

−4

10

−6

10

−8

10

−10

10

−12

10

η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L
η=0.5/L
η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L
η=0.5/L
η=6/5L

10

F (xs ) − F (x* )

−2

10

−2

−12

0

5

10
15
20
Gradient evaluations / n

25

10

30

0

5

10
15
20
25
Gradient evaluations / n

30

35

(a) Logistic regression: λ = 10−4 (left) and λ = 10−5 (right)
−2

10

−8

10

−4

*
s

−6

10

−10

−6

10

−8

10

−10

10

10

−12

10

η=0.05/L
η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L
η=0.05/L
η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L

10
F (x ) − F (x )

−4

10
F (xs ) − F (x* )

−2

10

η=0.05/L
η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L
η=0.05/L
η=0.1/L
η=0.2/L
η=0.3/L
η=0.4/L

−12

0

10

20
30
40
Gradient evaluations / n

50

60

10

0

10

20
30
40
Gradient evaluations / n

50

60

(b) Ridge regression: λ = 10−4 (left) and λ = 10−5 (right)
Fig. 2. Comparison of SVRG [1] and VR-SGD with different learning rates for solving `2 -norm (i.e., (λ/2)k · k2 ) regularized logistic regression and
ridge regression problems on the Covtype data set. In each plot, the vertical axis shows the objective value minus the minimum, and the horizontal
axis is the number of effective passes. Note that the blue lines stand for the results of SVRG with different learning rates, while the red lines
correspond to the results of VR-SGD with different learning rates (best viewed in colors).

When g(x) is non-smooth, e.g., the `1 -norm regularizer, the update rule of VR-SGD becomes


e is (xs ) .
xsk+1 = Prox η, g xsk − η ∇f
k
k

(13)

Different from the proximal stochastic gradient methods such as Prox-SVRG [2], all of which have a unified update rule
as in (10) for both the smooth and non-smooth cases of g(·), VR-SGD has two different update rules for the two cases, as
stated in (12) and (13). This leads to the following advantage over the Prox-SG methods: the stochastic gradient update
rule in (12) usually outperforms the proximal stochastic gradient update rule in (11), as well as the two classes of update
rules for Katyusha [3] (see Section 5.4 for details).
Fig. 2 demonstrates that VR-SGD has a significant advantage over SVRG in terms of robustness of learning rate selection.
That is, VR-SGD yields good performance within the range of the learning rate between 0.1/L and 0.4/L, whereas the
performance of SVRG is very sensitive to the selection of learning rates. Thus, VR-SGD is convenient to apply in various
real-world problems of large-scale machine learning. In fact, VR-SGD can use much larger learning rates than SVRG for
logistic regression problems in practice, e.g., 6/(5L) for VR-SGD vs. 1/(10L) for SVRG, as shown in Fig. 2(a).

9

3.3

The Algorithm for Non-Strongly Convex Objectives

Although many variance reduced stochastic methods have been proposed, most of them, including SVRG and ProxSVRG, only have convergence guarantees for the case of Problem (1) when the objective function F (x) is strongly convex.
However, F (x) may be non-strongly convex in many machine learning applications such as Lasso and `1 -norm regularized
logistic regression. As suggested in [3], [61], this class of problems can be transformed into strongly convex ones by adding
a proximal term (τ /2)kx − xs0 k2 , which can be efficiently solved by Algorithm 2. However, the reduction technique may
degrade the performance of the involved algorithms both in theory and in practice [42]. Thus, we present an efficient
VR-SGD algorithm for directly solving the non-strongly convex problem (1), as outlined in Algorithm 3.
The main difference between Algorithm 2 and Algorithm 3 is the setting for the learning rate. Similar to Algorithm 2,
the learning rate ηs of Algorithm 3 can also be fixed to a constant. Inspired by existing accelerated stochastic algorithms [3],
[34], the learning rate ηs in Algorithm 3 can be gradually increased, which in principle leads to faster convergence (see
Section 5.5 for details). Different from existing stochastic methods such as Katyusha [3], the update rule of the learning rate

ηs for Algorithm 3 is defined as follows: η0 = c, where c > 0 is an initial learning rate, and for any s ≥ 1,
ηs = η0 / max{α, 2/(s + 1)}

(14)

where 0 < α ≤ 1 is a given constant, e.g., α = 0.2.
3.4

Complexity Analysis

From Algorithms 2 and 3, we can see that the per-iteration cost of VR-SGD is dominated by the computation of ∇fisk(xsk ),
es−1 ), and ∇g(xsk ) or the proximal update in (13). Thus, the complexity is O(d), as low as that of SVRG [1] and
∇fisk(x
es−1 ) in the computation of
Prox-SVRG [2]. In fact, for some ERM problems, we can save the intermediate gradients ∇fi (x
es , which generally requires O(n) additional storage. As a result, each epoch only requires (n + m) component gradient
µ
evaluations. In addition, for extremely sparse data, we can introduce the lazy update tricks in [36], [62], [63] to our
algorithms, and perform the update steps in (12) and (13) only for the non-zero dimensions of each example, rather than
all dimensions. In other words, the per-iteration complexity of VR-SGD can be improved from O(d) to O(d0 ), where d0 ≤ d
is the sparsity of feature vectors. Moreover, VR-SGD has a much lower per-iteration complexity than existing accelerated
stochastic variance reduction methods such as Katyusha [3], which have at least two more update rules for additional
variables, as shown in (6a)-(6c).
3.5

Extensions of VR-SGD

It has been shown in [36], [37] that mini-batching can effectively decrease the variance of stochastic gradient estimates. In
this part, we first extend the proposed VR-SGD method to the mini-batch setting, as well as its convergence results below.
Here, we denote by b the mini-batch size and Iks the selected random index set Ik ⊂ [n] for each outer-iteration s ∈ [S] and
inner-iteration k ∈ {0, 1, . . . , m−1}. The variance reduced stochastic gradient estimator in (5a) becomes
X

e I s (xs ) = 1
es−1 ) +∇f (x
es−1 )
∇f
∇fi (xsk )−∇fi (x
k
k
b i∈I s
k

where Iks ⊂ [n] is a mini-batch of size b. If some component functions are non-smooth, we can use the proximal operator
oracle [61] or the Nesterov’s smoothing [64] and homotopy smoothing [65] techniques to smoothen them, and thereby
obtain the smoothed approximations of the functions fi (x). In addition, we can directly extend our algorithms to the
non-smooth setting as in [34], e.g., Algorithm 3 in [34].
Considering that each component function fi (x) maybe have different degrees of smoothness, picking the random index

isk

from a non-uniform distribution is a much better choice than the commonly used uniform random sampling [66], [67],

as well as without-replacement sampling vs. with-replacement sampling [27]. This can be done using the same techniques
Pn
in [2], [3], i.e., the sampling probabilities for all fi (x) are proportional to their Lipschitz constants, i.e., pi = Li / j=1 Lj .

10

Algorithm 3 VR-SGD for non-strongly convex objectives
Input: The number of epochs S , and the number of iterations m per epoch.
e0 , α > 0 and η0 .
Initialize: x10 = x
2:

s = 1, 2, . . . , S do
P
es = n1 n
es−1 );
µ
i=1 ∇fi (x

3:

ηs = η0 / max{α, 2/(s+1)};

4:

for k = 0, 1, . . . , m − 1 do

1: for

% Compute full gradient
% Compute step sizes

Pick isk uniformly at random from [n];
e is (xs ) = ∇fis (xs ) − ∇fis (x
es−1 ) + µ
es ;
∇f
k
k
k
hk k
i
e is (xs ) + ∇g(xs ) ,
xsk+1 = xsk − ηs ∇f
k
k
 k

e is (xs ) ;
or xsk+1 = Prox ηs , g xsk − ηs ∇f
k
k

5:
6:
7:

% The stochastic gradient estimator
% Smooth case of g(·)
% Non-smooth case of g(·)

end for

8:

1 Pm−1 s
k=1 xk ;
m−1
1 Pm
s
e = m k=1 xsk ;
Option II: x
xs+1
= xsm ;
0

es =
Option I: x

9:
10:
11:

e
% Iterate averaging for snapshot x
e
% Iterate averaging for snapshot x

% Initiate xs+1
for the next epoch
0

12: end for

bS = x
eS if F (x
eS ) ≤ F ( S1
Output: x

PS

es ),
s=1 x

bS = S1
and x

PS

es
s=1 x

otherwise.

Moreover, our VR-SGD method can also be combined with other accelerated techniques proposed for SVRG. For instance,
the epoch length of VR-SGD can be automatically determined by the techniques in [42], [68] instead of a fixed epoch
length. We can reduce the number of gradient calculations in the early iterations as in [34], [41], [42], which leads to faster
convergence in general. Moreover, we can also introduce the Nesterov’s acceleration technique as in [24], [37], [38], [39],
[40] and momentum acceleration trick as in [3], [34] to further improve the performance of VR-SGD.

4

C ONVERGENCE A NALYSIS

In this section, we provide the convergence guarantees of VR-SGD for solving both smooth and non-smooth general convex
problems. We also extend these results to the mini-batch setting. Moreover, we analyze the convergence properties of VRSGD for solving both smooth and non-smooth strongly convex objective functions. We first introduce the following lemma
which is useful in our analysis.
Lemma 1 (3-point property, [69]). Let ẑ be the optimal solution of the following problem,

min

z∈Rd

τ
kz − z0 k2 + r(z)
2

where r(z) is a convex function (but possibly non-differentiable), and τ ≥ 0. Then for any z ∈ Rd , the following inequality holds

r(ẑ) +
4.1

τ
τ
τ
kẑ − z0 k2 ≤ r(z) + kz − z0 k2 − kz − ẑk2 .
2
2
2

Convergence Properties for Non-strongly Convex Problems

In this part, we analyze the convergence properties of VR-SGD for solving the more general non-strongly convex problems.
Considering that the two proposed algorithms (i.e., Algorithms 2 and 3) have two different update rules for the smooth
and non-smooth problems, we give the convergence guarantees of VR-SGD for the two cases as follows.
4.1.1

Smooth Objectives

We first give the convergence guarantee of Problem (1) when the objective function F (x) is smooth and non-strongly convex.
In order to simplify analysis, we denote F (x) by f (x), that is, fi (x) := fi (x)+g(x) for all i = 1, 2, . . . , n.

11

Lemma 2 (Variance bound of smooth objectives). Let x∗ be the optimal solution of Problem (1). Suppose Assumption 1 holds, and
e is (xs ) := ∇fis (xs ) − ∇fis (x
es−1 ) + ∇f (x
es−1 ). Then the following inequality holds
let ∇f
k

k

k

k

k

i
e is (xs ) − ∇f (xs )k2 ≤ 4L[f (xs ) − f (x∗ ) + f (x
es−1 ) − f (x∗ )].
E k∇f
k
k
k
k
h

The proof of this lemma is included in APPENDIX A. Lemma 2 provides the upper bound on the expected variance
of the variance reduced gradient estimator in (5a), i.e., the SVRG estimator independently introduced in [1], [35]. For
Algorithm 3 with Option I and a fixed learning rate, we give the following key result for our analysis.
Lemma 3 (Option I and smooth objectives). Let β = 1/(Lη). If each fi (·) is convex and L-smooth, then the following inequality
holds for all s = 1, 2, . . . , S ,


2
1
es ) − f (x∗ )] +
1−
E[f (x
E[f (xsm ) − f (x∗ )]
β −1
m−1




2m
2
Lβ
es−1 )−f (x∗ ) +
≤
E f (x
E[f (xs0 )−f (x∗ )] +
E kx∗ −xs0 k2 − kx∗ −xsm k2 .
(β −1)(m−1)
(β −1)(m−1)
2(m−1)
es−1 )+∇f (x
es−1 ).
Proof. In order to simplify notation, the stochastic gradient estimator is defined as: vks := ∇fisk (xsk )−∇fisk (x
Since each component function fi (x) is L-smooth, which implies that the gradient of the average function f (x) is also Lsmooth, i.e., for all x, y ∈ Rd ,

k∇f (x) − ∇f (y)k ≤ Lkx − yk,
whose equivalent form is

f (y) ≤ f (x) + h∇f (x), y − xi +

L
ky − xk2 .
2

Applying the above smoothness inequality, we have

L s
2
x
− xsk
2 k+1
Lβ s
L(β −1) s
2
= f (xsk ) + ∇f (xsk ), xsk+1 − xsk +
xk+1 − xsk −
xk+1 − xsk
2
2
Lβ s
kxk+1 − xsk k2
= f (xsk ) + vks , xsk+1 − xsk +
2
L(β −1) s
+ ∇f (xsk ) − vks , xsk+1 − xsk −
kxk+1 − xsk k2 ,
2
where β = 1/(Lη) > 3 is a constant. Using Lemma 2, then we get


L(β −1) s
s 2
s
s
s
s
kxk+1 − xk k
E ∇f (xk ) − vk , xk+1 − xk −
2


1
L(β −1) s
L(β −1) s
≤E
k∇f (xsk ) − vks k2 +
kxk+1 −xsk k2 −
kxk+1 −xsk k2
2L(β −1)
2
2

2 
s
∗
s−1
∗
e
≤
f (xk ) − f (x ) + f (x
) − f (x ) ,
β −1
f (xsk+1 ) ≤ f (xsk ) + ∇f (xsk ), xsk+1 − xsk +

2

(15)

(16)

where the first inequality holds due to the Young’s inequality (i.e., y T z ≤ kyk2 /(2γ)+γkzk2 /2 for all γ > 0 and y, z ∈ Rd ),
and the second inequality follows from Lemma 2.
Substituting the inequality in (16) into the inequality in (15), and taking the expectation with respect to the random
choice isk , we have
E[f (xsk+1 )]



Lβ s
2 
s
s
s
s 2
es−1 )−f (x∗ )
+ E vk , xk+1 −xk +
≤
kxk+1 −xk k +
f (xsk )−f (x∗ )+f (x
2
β −1



Lβ
2 
es−1 )−f (x∗ )
≤ E[f (xsk )] + E hvks , x∗ −xsk i +
(kx∗ −xsk k2 − kx∗ −xsk+1 k2 ) +
f (xsk )−f (x∗ )+f (x
2
β −1



Lβ
2 
s
∗
s 2
∗
s
s
∗
s
es−1 )−f (x∗ )
≤ E[f (xk )] + E[h∇f (xk ), x −xk i] + E
(kx −xk k − kx −xk+1 k2 ) +
f (xsk )−f (x∗ )+f (x
2
β −1




Lβ
2
es−1 ) − f (x∗ ) .
≤ f (x∗ ) +
E kx∗ − xsk k2 − kx∗ − xsk+1 k2 +
f (xsk ) − f (x∗ ) + f (x
2
β −1
E[f (xsk )]

12

Here, the first inequality holds due to the inequality in (15) and the inequality in (16); the second inequality follows
from Lemma 1 with ẑ = xsk+1 , z = x∗ , z0 = xsk , τ = Lβ = 1/η , and r(z) := hvks , z − xsk i; the third inequality holds
due to the fact that E[vks ] = ∇f (xsk ); and the last inequality follows from the convexity of the smooth function f (·), i.e.,

f (xsk )+h∇f (xsk ), x∗ −xsk i ≤ f (x∗ ). The above inequality can be rewritten as follows:
E[f (xsk+1 )] − f (x∗ ) ≤

 Lβ  ∗

2 
es−1 ) − f (x∗ ) +
f (xsk ) − f (x∗ ) + f (x
E kx − xsk k2 − kx∗ − xsk+1 k2 .
β −1
2

Summing the above inequality over k = 0, 1, . . . , m−1, then

m
m 
X
X
 Lβ  ∗

2 
es−1 ) − f (x∗ ) +
f (xsk−1 ) − f (x∗ ) + f (x
E kx − xsk−1 k2 − kx∗ − xsk k2 .
E[f (xsk ) − f (x∗ )] ≤
β −1
2
k=1
k=1
1 Pm−1 s
es = m−1
Due to the setting of x
k=1 xk in Option I, and the convexity of f (·), then we have
es ) ≤
f (x

m−1
X
1
f (xsk ).
m − 1 k=1

(17)

(18)

The left and right hand sides of the inequality in (17) can be rewritten as follows:
m
X

E[f (xsk ) − f (x∗ )] =

k=1

m−1
X

E[f (xsk ) − f (x∗ )] + E[f (xsm ) − f (x∗ )] ,

k=1


m 
X
 Lβ  ∗

2 
es−1 ) − f (x∗ ) +
f (xsk−1 ) − f (x∗ ) + f (x
E kx − xsk−1 k2 − kx∗ − xsk k2
β −1
2
k=1
m−1

2 X
Lβ  ∗
2 
es−1 )−f (x∗ )] +
f (xs0 )−f (x∗ ) + m[f (x
E kx −xs0 k2 −kx∗ −xsm k2 .
[f (xsk )−f (x∗ )] +
β −1 k=1
β −1
2
2 Pm−1
s
∗
Subtracting β−1
k=1 [f (xk )−f (x )] from both sides of the inequality in (17), then we obtain

=

 m−1
X
2
1−
E[f (xsk ) − f (x∗ )] + E[f (xsm ) − f (x∗ )]
β −1 k=1


2
Lβ  ∗
es−1 )−f (x∗ )] +
≤
E f (xs0 )−f (x∗ ) + m[f (x
E kx −xs0 k2 − kx∗ −xsm k2 .
β −1
2


Applying the inequality in (18), we have


2
es ) − f (x∗ )] + E[f (xsm ) − f (x∗ )]
1−
(m−1)E[f (x
β −1
 m−1

X
2
E[f (xsk ) − f (x∗ )] + E[f (xsm ) − f (x∗ )]
≤ 1−
β −1 k=1

≤



2
Lβ  ∗
es−1 )−f (x∗ )] +
E f (xs0 )−f (x∗ ) + m[f (x
E kx −xs0 k2 − kx∗ −xsm k2 .
β −1
2

Dividing both sides of the above inequality by (m−1), we arrive at


2
1
es ) − f (x∗ )] +
1−
E[f (x
E[f (xsm ) − f (x∗ )]
β −1
m−1




2
2m
Lβ
es−1 )−f (x∗ ) +
≤
E[f (xs0 )−f (x∗ )]+
E f (x
E kx∗ −xs0 k2 −kx∗ −xsm k2 .
(β −1)(m−1)
(β −1)(m−1)
2(m−1)
This completes the proof.
The first main result is the following theorem, which provides the convergence guarantee of VR-SGD with Option I for
solving smooth and general convex minimization problems.
Theorem 1 (Option I and smooth objectives). Suppose Assumption 1 holds. Then the following inequality holds
h
i
2(m + 1)
β(β − 1)L
bS ) − f (x∗ ) ≤
e0 ) − f (x∗ )] +
e0 − x∗ k2 ,
E f (x
[f (x
kx
[(β −1)(m−1)−4m+2]S
2[(β −1)(m−1)−4m+2]S
P
bS = x
eS if f (x
eS ) ≤ f (xS ), and xS = S1 S
es . Otherwise, x
bS = xS .
where x
s=1 x

13

Proof. Since 2/(β −1) < 1, it is easy to verify that

2
1
{E[f (xsm )] − f (x∗ )} ≤
{E[f (xsm )] − f (x∗ )} .
(β −1)(m−1)
m−1

(19)

Applying the above inequality and Lemma 3, we have


2
2
es ) − f (x∗ )] +
E[f (x
E[f (xsm ) − f (x∗ )]
1−
β −1
(β −1)(m−1)


1
2
es ) − f (x∗ )] +
E[f (x
≤ 1−
E[f (xsm ) − f (x∗ )]
β −1
m−1




2
2m
Lβ
es−1 ) − f (x∗ ) +
≤
E[f (xs0 ) − f (x∗ )] +
E f (x
E kx∗ − xs0 k2 − kx∗ − xsm k2 .
(β −1)(m−1)
(β −1)(m−1)
2(m−1)
Summing the above inequality over s = 1, 2, . . . , S , taking expectation with respect to the history of random variables isk ,
and using the setting of xs+1
= xsm , we obtain
0

S 
X
2
es ) − f (x∗ )]
1−
E[f (x
β
−1
s=1

S 
X
2
2m
es−1 ) − f (x∗ )]
≤
E[f (xs0 ) − f (x∗ ) − (f (xsm ) − f (x∗ ))] +
E[f (x
(β −1)(m−1)
(β −1)(m−1)
s=1
S

Lβ X  ∗
E kx − xs0 k2 − kx∗ − xsm k2 .
2(m−1) s=1
PS−1
2m
es )−f (x∗ )] from both sides of the above inequality, we have
Subtracting (β−1)(m−1)
s=1 [f (x

+


S 
X
4
2
2m
S
∗
e ) − f (x )] +
es ) − f (x∗ )]
E[f (x
1−
−
E[f (x
(β −1)(m−1)
β
−1
(β
−1)(m−1)
s=1
h
i
2m
2
1
∗
∗
e0 ) − f (x∗ )]
E f (x0 ) − f (x ) − (f (xS
E[f (x
≤
m ) − f (x )) +
(β −1)(m−1)
(β −1)(m−1)
h
i
Lβ
2
+
E kx∗ − x10 k2 − kx∗ − xS
.
mk
2(m − 1)
e0 = x10 , we arrive at
Dividing both sides of the above inequality by S , and using the setting of x

X
S
1
2
4
es ) − f (x∗ )]
−
E[f (x
1−
S
β −1 (β −1)(m−1) s=1

X
S
2m
1
2
4
eS ) − f (x∗ )] +
es ) − f (x∗ )]
≤
E[f (x
−
E[f (x
1−
(β −1)(m−1)S
S
β −1 (β −1)(m−1) s=1

2
2m
Lβ
e0 ) − f (x∗ )] +
[f (x10 ) − f (x∗ )] +
[f (x
kx∗ − x10 k2
(β −1)(m−1)S
(β −1)(m−1)S
2(m−1)S
2(m + 1)
Lβ
e0 ) − f (x∗ )] +
e0 − x∗ k2 ,
=
[f (x
kx
(β −1)(m−1)S
2(m−1)S
≤

eS ) − f (x∗ ) ≥ 0; the second inequality holds due to the facts that
where the first inequality holds due to the fact that f (x
e0 = x10 .
f (xSm )−f (x∗ ) ≥ 0 and kx∗ −xSm k2 ≥ 0; and the last equality follows from the setting of x
P
S
es ), and therefore the above inequality
Due to the definition of xS and the convexity of f (·), we have f (xS ) ≤ S1 s=1 f (x
becomes:


 h
i
4
2
−
E f (xS ) − f (x∗ )
β −1 (β −1)(m−1)
2(m + 1)
Lβ
e0 ) − f (x∗ )] +
e0 − x∗ k2 .
≤
[f (x
kx
(β −1)(m−1)S
2(m−1)S

1−

4
2
Dividing both sides of the above inequality by c1 = 1− β−1
− (β−1)(m−1)
> 0, we have

h
i
E f (xS ) − f (x∗ ) ≤

2(m + 1)
Lβ
e0 ) − f (x∗ )] +
e0 − x∗ k2 .
[f (x
kx
c1 (β −1)(m−1)S
2c1 (m−1)S

(20)

14

bS = x
eS if f (x
eS ) ≤ f (xS ). Then
Due to the setting for the output of Algorithm 3, x
h
i
h
i
2(m + 1)
Lβ
bS ) − f (x∗ ) ≤ E f (xS ) − f (x∗ ) ≤
e0 ) − f (x∗ )] +
e0 − x∗ k2 .
E f (x
[f (x
kx
c1 (β −1)(m−1)S
2c1 (m−1)S
bS = xS , and the above inequality still holds. This completes the proof.
eS ) ≥ f (xS ), let x
Alternatively, when f (x
From Lemma 3, Theorem 1, and their proofs, one can see that our convergence analysis is very different from those of
existing stochastic methods, such as SVRG [1] and Prox-SVRG [2]. For Algorithm 3 with Option II and a fixed learning rate,
we give the following key lemma for our convergence analysis.
Lemma 4 (Option II and smooth objectives). If each fi (·) is convex and L-smooth, then the following inequality holds for all

s = 1, 2, . . . , S ,



2
2
es ) − f (x∗ )] +
1−
E[f (x
E[f (xsm ) − f (x∗ )]
β −1
(β −1)m



2
Lβ  ∗
2
es−1 ) − f (x∗ ) +
E f (x
E[f (xs0 ) − f (x∗ )] +
E kx − xs0 k2 − kx∗ − xsm k2 .
≤
(β −1)
(β −1)m
2m
This lemma is a slight generalization of Lemma 3, and we give the proof in APPENDIX B for completeness.
Theorem 2 (Option II and smooth objectives). If each fi (·) is convex and L-smooth, then the following inequality holds
h
i
Lβ(β −1)
2(m+1)
e0 ) − f (x∗ )] +
e0 − x∗ k2 ,
bS ) − f (x∗ ) ≤
[f (x
kx
E f (x
mS(β −5)
2mS(β −5)
P
es . Otherwise, x
bS = xS .
bS = x
eS if f (x
eS ) ≤ f (xS ), and xS = S1 S
where x
s=1 x
The proof of this theorem can be found in APPENDIX C. Clearly, Theorems 1 and 2 show that VR-SGD with Option I
or Option II attains a sub-linear convergence rate for smooth general convex objective functions. This means that VR-SGD
is guaranteed to have a similar convergence rate with other variance reduced stochastic methods such as SAGA [30], and
a slower theoretical rate than accelerated methods such as Katyusha [3]. Nevertheless, VR-SGD usually converges much
faster than the best known stochastic method, Katyusha [3] in practice (see Section 5.3 for details).
4.1.2

Non-Smooth Objectives

Next, we provide the convergence guarantee of Problem (1) when the objective function F (x) is non-smooth (i.e., the
regularizer g(x) is non-smooth) and non-strongly convex. We first give the following lemma.
Lemma 5 (Variance bound of non-smooth objectives). Let x∗ be the optimal solution of Problem (1). If each fi (·) is convex and
e is (xs ) := ∇fis (xs ) − ∇fis (x
es−1 ) + ∇f (x
es−1 ), then the following inequality holds
L-smooth, and ∇f
k

k

k

k

k

h
i
e is (xs ) − ∇f (xs )k2 ≤ 4L[F (xs ) − F (x∗ ) + F (x
es−1 ) − F (x∗ )].
E k∇f
k
k
k
k
This lemma can be viewed as a generalization of Lemma 2, and is essentially identical to Corollary 3.5 in [2] and Lemma
1 in [45], and hence its proof is omitted. For Algorithm 3 with Option II and a fixed learning rate, we give the following
results.
Lemma 6 (Option II and non-smooth objectives). If each fi (·) is convex and L-smooth, then the following inequality holds for all

s = 1, 2, . . . , S ,



2
2
es ) − F (x∗ )] +
E[F (x
E[F (xsm )−F (x∗ )]
β −1
(β −1)m

 Lβ  ∗

2
2
es−1 )−F (x∗ ) +
≤
E[F (xs0 ) − F (x∗ )] +
E F (x
E kx −xs0 k2 −kx∗ −xsm k2 .
(β −1)m
β −1
2m

1−

(21)

The proof of this lemma can be found in APPENDIX D. Using the above lemma, we give the following convergence
result for VR-SGD.

15

Theorem 3 (Option II and non-smooth objectives). Suppose Assumption 1 holds. Then the following inequality holds
h
i
2(m+1)
β(β −1)L
bS ) − F (x∗ ) ≤
e0 ) − F (x∗ )] +
e0 − x∗ k2 .
E F (x
[F (x
kx
(β −5)mS
2(β −5)mS
The proof of this theorem is included in APPENDIX E. Similarly, for Algorithm 3 with Option I and a fixed learning
rate, we give the following results.
Corollary 1 (Option I and non-smooth objectives). If each fi (·) is convex and L-smooth, then the following inequality holds for
all s = 1, 2, . . . , S ,


1
2
es ) − F (x∗ )] +
E[F (x
E[F (xsm ) − F (x∗ )]
1−
β −1
m−1




2m
2
Lβ
es−1 )−F (x∗ ) +
≤
E F (x
E[F (xs0 )−F (x∗ )] +
E kx∗ −xs0 k2 − kx∗ −xsm k2 .
(β −1)(m−1)
(β −1)(m−1)
2(m−1)
Corollary 2 (Option I and non-smooth objectives). Suppose Assumption 1 holds. Then the following inequality holds
h
i
bS ) − F (x∗ )
E F (x

2(m + 1)
β(β − 1)L
e0 ) − F (x∗ )] +
e0 − x∗ k2 ,
[F (x
kx
[(β −1)(m−1)−4m+2]S
2[(β −1)(m−1)−4m+2]S
P
es . Otherwise, x
bS = xS .
bS = x
eS if F (x
eS ) ≤ F (xS ), and xS = S1 S
where x
s=1 x
≤

Corollaries 1 and 2 can be viewed as the generalizations of Lemma 3 and Theorem 1, respectively, and hence their
proofs are omitted. Obviously, Theorems 3 and Corollary 2 show that VR-SGD with Option I or Option II attains a sublinear convergence rate for non-smooth general convex objective functions.
4.1.3

Mini-Batch Settings

e is (xs ) in Lemma 5 is extended to the mini-batch
The upper bound on the variance of the stochastic gradient estimator ∇f
k
k
setting as follows.
e I s (xs ) := 1
Lemma 7 (Variance bound of mini-batch). If each fi (·) is convex and L-smooth, and ∇f
k
b
k

P

i∈Iks



es−1 ) +
∇fi (xsk ) − ∇fi (x

es−1 ), then the following inequality holds
∇f (x
h
i
e I s (xs ) − ∇f (xs )k2 ≤ 4Lδ(b)[F (xs ) − F (x∗ ) + F (x
es−1 ) − F (x∗ )],
E k∇f
k
k
k
k
where δ(b) = (n−b)/[(n−1)b].
It is not hard to verify that 0 ≤ δ(b) ≤ 1. This lemma is essentially identical to Theorem 4 in [36], and hence its proof
is omitted. Based on the variance upper bound in Lemma 7, we further analyze the convergence properties of VR-SGD for
the mini-batch setting. Lemma 6 is first extended to the mini-batch setting as follows.
Lemma 8 (Mini-batch). Using the same notation as in Lemma 7, we have


2δ(b)
2δ(b)
s
∗
es )] − F (x∗ )}
{E[F (xm )] − F (x )} + 1 −
{E[F (x
(β −1)m
β −1


2δ(b) 
2δ(b)
Lβ  ∗
es−1 ) − F (x∗ ) +
≤
F (x
[F (xs0 ) − F (x∗ )] +
E kx − xs0 k2 − kx∗ − xsm k2 .
(β −1)
(β −1)m
2m
Proof. In order to simplify notation, the stochastic gradient estimator of mini-batch is defined as:

1 X
es−1 ) + ∇f (x
es−1 ).
vks :=
∇fi (xsk ) − ∇fi (x
b i∈I s
k

F (xsk+1 ) ≤ g(xsk+1 ) + f (xsk ) + ∇f (xsk ), xsk+1 − xsk +

Lβ s
kxk+1 − xsk k2
2
L(β −1) s
−
kxk+1 − xsk k2 .
2

= g(xsk+1 ) + f (xsk ) + vks , xsk+1 − xsk +
+ ∇f (xsk ) − vks , xsk+1 − xsk

Lβ s
xk+1 − xsk
2

2

−

L(β −1) s
xk+1 − xsk
2

2

(22)

16

Using Lemma 7, then we obtain


L(β −1) s
E ∇f (xsk ) − vks , xsk+1 − xsk −
kxk+1 − xsk k2
2


1
L(β −1) s
L(β −1) s
s
s 2
(23)
≤E
k∇f (xk ) − vk k +
kxk+1 −xsk k2 −
kxk+1 −xsk k2
2L(β −1)
2
2

2δ(b) 
es−1 ) − F (x∗ ) ,
≤
F (xsk ) − F (x∗ ) + F (x
β −1
where the first inequality holds due to the Young’s inequality, and the second inequality follows from Lemma 7. Substituting
the inequality (23) into the inequality (22), and taking the expectation over the random mini-batch set Iks , we have
E[F (xsk+1 )]



Lβ s
2δ(b) 
es−1 )−F (x∗ )
≤ E[g(xsk+1 )] + E[f (xsk )] + E vks , xsk+1 −xsk +
kxk+1 −xsk k2 +
F (xsk )−F (x∗ )+F (x
2
β −1



Lβ
2δ(b) 
∗
s
s
∗
s
∗
s 2
∗
s
2
es−1 )−F (x∗ )
≤ g(x ) + E[f (xk )] + E hvk , x −xk i +
(kx −xk k −kx −xk+1 k ) +
F (xsk )−F (x∗ )+F (x
2
β −1



Lβ
2δ(b) 
∗
s 2
∗
s
s
∗
∗
2
es−1 )−F (x∗ )
≤ g(x ) + f (x ) + E
(kx − xk k − kx − xk+1 k ) +
F (xk )−F (x∗ )+F (x
2
β −1




2δ(b)
Lβ
es−1 )−F (x∗ ) ,
F (xsk )−F (x∗ )+F (x
E kx∗ − xsk k2 − kx∗ − xsk+1 k2 +
= F (x∗ ) +
2
β−1
where the second inequality holds from Lemma 1. Then the above inequality is rewritten as follows:


 Lβ  ∗

2δ(b) 
(24)
es−1 ) − F (x∗ ) +
E F (xsk+1 ) − F (x∗ ) ≤
F (xsk ) − F (x∗ ) + F (x
E kx − xsk k2 − kx∗ − xsk+1 k2 .
β −1
2
Summing the above inequality over k = 0, 1, · · · , (m − 1), then
m−1
X


E[F (xsk+1 )] − F (x∗ )

k=0
m−1
X


 Lβ  ∗

2δ(b) 
s
∗
s−1
∗
s 2
∗
s
2
e
≤
F (xk ) − F (x ) + F (x
) − F (x ) +
E kx − xk k − kx − xk+1 k
.
β −1
2
k=0
1 Pm
1 Pm
s
s
es = m
es ) ≤ m
Since x
k=1 xk , we have F (x
k=1 F (xk ), and


2δ(b)
2δ(b)
s
∗
es )] − F (x∗ )}
{E[F (xm )] − F (x )} + 1 −
{E[F (x
(β −1)m
β −1


2δ(b)
Lβ  ∗
2δ(b) 
es−1 ) − F (x∗ ) +
F (x
[F (xs0 ) − F (x∗ )] +
E kx − xs0 k2 − kx∗ − xsm k2 .
≤
(β −1)
(β −1)m
2m
This completes the proof.
Similar to Lemma 8, we can also extend Lemma 3 and Corollary 1 to the mini-batch setting.
Theorem 4 (Mini-batch). If each fi (·) is convex and L-smooth, then the following inequality holds
h
i




2δ(b)(m + 1)
Lβ(β −1)
bS ) − F (x∗ ) ≤
e0 )−F (x∗ ) +
e0 k2 .
E F (x
E F (x
E kx∗ − x
(β −1−4δ(b))mS
2(β −1−4δ(b))mS

(25)

Proof. Using Lemma 8, we have


2δ(b)
es ) − F (x∗ )]
1−
E[F (x
β −1


2δ(b) 
2δ(b)
Lβ  ∗
es−1 ) − F (x∗ ) +
≤
F (x
[F (xs0 ) − F (x∗ )] +
E kx − xs0 k2 − kx∗ − xsm k2 .
(β −1)
(β −1)m
2m

= xsm , we
Summing the above inequality over s = 1, 2, . . . , S , taking expectation over whole history of Iks , and using xs+1
0
obtain
S 
X

≤

s=1
S
X


2δ(b)
es ) − F (x∗ )]
1−
E [F (x
β −1

S
S
 X

2δ(b) 
2δ(b)
Lβ X  ∗
es−1 ) − F (x∗ ) +
F (x
E[F (xs0 ) − F (xsm )] +
E kx − xs0 k2 − kx∗ − xsm k2 .
(β −1)
(β −1)m
2m s=1
s=1
s=1

17

Subtracting

2δ(b) PS−1
s=1
β−1

es )−F (x∗ )] from both sides of the above inequality, we have
E[F (x


S 
i X
2δ(b) h
4δ(b)
S
∗
es ) − F (x∗ )]
e ) − F (x ) +
E[F (x
E F (x
1−
β −1
β −1
s=1
i
i Lβ h

2δ(b) 
2δ(b) h
∗
1 2
∗
S 2
e0 ) − F (x∗ ) +
≤
F (x
E F (x10 ) − F (xS
)
+
E
kx
−
x
k
−
kx
−
x
k
m
0
m
(β −1)
(β −1)m
2m
P
S
es ), we arrive at
Dividing both sides of the above inequality by S and using E[F (x)] ≤ S1 s=1 F (x


2δ(b)
4δ(b)
eS ) − F (x∗ )] + 1 −
E[F (x) − F (x∗ )]
E[F (x
(β −1)S
β −1
h
i
i

2δ(b)
2δ(b) 
Lβ h ∗
2
e0 ) − F (x∗ ) +
E F (x10 ) − F (xS
E kx − x10 k2 − kx∗ − xS
F (x
.
≤
m) +
mk
(β −1)S
(β −1)mS
2mS
Subtracting

2δ(b)
∗
eS
(β−1)S E[F (x )−F (x )]

from both sides of the above inequality, we have





4δ(b)
E[F (x) − F (x∗ )]
β −1
i
i
h
i
2δ(b) h
2δ(b)
Lβ h ∗
1 2
∗
S 2
e0 ) − F (x
eS ) +
.
≤
E F (x
E F (x10 ) − F (xS
)
+
E
kx
−
x
k
−
kx
−
x
k
m
0
m
(β −1)S
(β −1)mS
2mS
1−

Dividing both sides of the above inequality by (1−

4δ(b)
β−1 ) > 0,

we arrive at

E[F (x)] − F (x∗ )
h
i
h
i


2δ(b)
Lβ(β −1)
2δ(b)
e0 )−F (x
eS ) +
e0 )−F (xS
e0 k2
E F (x
E F (x
)
+
E kx∗ − x
m
(β −1−4δ(b))S
(β −1−4δ(b))mS
2(β −1−4δ(b))mS






2δ(b)
Lβ(β −1)
2δ(b)
e0 )−F (x∗ ) +
e0 )−F (x∗ ) +
e0 k2
≤
E F (x
E F (x
E kx∗ − x
(β −1−4δ(b))S
(β −1−4δ(b))mS
2(β −1−4δ(b))mS




2δ(b)(m + 1)
Lβ(β −1)
e0 )−F (x∗ ) +
e0 k2 .
=
E F (x
E kx∗ − x
(β −1−4δ(b))mS
2(β −1−4δ(b))mS

≤

bS = x
eS , and
eS ) ≤ F (xS ), then x
When F (x
h
i




2δ(b)(m + 1)
Lβ(β −1)
bS ) − F (x∗ ) ≤
e0 )−F (x∗ ) +
e0 k2 .
E F (x
E F (x
E kx∗ − x
(β −1−4δ(b))mS
2(β −1−4δ(b))mS
eS ) ≥ F (xS ), then x
bS = xS , and the above inequality still holds. This completes the proof.
Alternatively, if F (x
From Theorem 4, one can see that when b = n (i.e., the batch setting), we have δ(n) = 0, and the first term on the
right-hand side of (25) diminishes. In other words, our VR-SGD method degenerates to the deterministic method with the
convergence rate of O(1/T ). Furthermore, when b = 1, we have δ(1) = 1, and then Theorem 4 degenerates to Theorem 3.
4.2

Convergence Properties for Strongly Convex Problems

In this part, we analyze the convergence properties of VR-SGD (i.e., Algorithm 2) for solving the strongly convex objective
function (1). According to the above analysis, one can see that if F (·) is convex, and each convex component function fi (·)
is L-smooth, both Algorithms 2 and 3 converge to the optimal solution. In the following, we provide stronger convergence
rate guarantees for VR-SGD under the strongly convex condition. We first give the following assumption.
Assumption 3. For all s = 1, 2, . . . , S , the following inequality holds


es−1 ) − F (x∗ )
E[F (xs0 ) − F (x∗ )] ≤ C E F (x
where C > 0 is a constant4 .
Similar to Algorithm 3, Algorithm 2 also has two different update rules for the two cases of smooth and non-smooth
objective functions. We first give the following convergence result for Algorithm 2 with Option I.
4. This assumption shows the relationship of the gaps between the function values at the starting and snapshot points of each epoch and the
optimal value of the objective function. In fact, both x
es and xsm (i.e., xs+1
) converge to x∗ , and thus C is far less than m, i.e., C  m.
0

18

Theorem 5 (Option I). Suppose Assumptions 1, 2, and 3 hold, and m is sufficiently large so that

2Lη(m+C)
C(1−Lη)
+
< 1.
(m−1)(1−3Lη) µη(m−1)(1−3Lη)
Then Algorithm 2 with Option I has the following geometric convergence in expectation:


es ) − F (x∗ )] ≤ ρsI F (x
e0 ) − F (x∗ ) .
E [F (x
ρI :=

Proof. Since each fi (·) is convex and L-smooth, then Corollary 1 holds, which then implies






2
1
Lβ
∗
es ) − F (x∗ )] +
E[F (x
1−
E F (xs+1
E kxs+1
−x∗ k2
0 ) − F (x ) +
0
β −1
m−1
2(m−1)




2m
2
Lβ
es−1 )−F (x∗ ) +
E F (x
E[F (xs0 )−F (x∗ )]+
E kxs0 −x∗ k2 .
≤
(m−1)(β −1)
(m−1)(β −1)
2(m−1)

(26)

Due to the strong convexity of F (·), we have kxs0 − x∗ k2 ≤ (2/µ)[F (xs0 ) − F (x∗ )]. Then the inequality in (26) can be
rewritten as follows:



2
es ) − F (x∗ )]
E[F (x
β −1




2m
2
Lβ
es−1 )−F (x∗ ) +
E[F (xs0 )−F (x∗ )]
≤
E F (x
+
(m−1)(β −1)
(m−1)(β −1) µ(m−1)






2m
2C
CLβ
es−1 )−F (x∗ ) +
es−1 )−F (x∗ )
≤
E F (x
+
E F (x
(m−1)(β −1)
(m−1)(β −1) µ(m−1)




CLβ
2(m+C)
s−1
e )−F (x∗ )
+
E F (x
=
(m−1)(β −1) µ(m−1)

1−

where the first inequality holds due to the fact that kxs0 − x∗ k2 ≤ (2/µ)[F (xs0 ) − F (x∗ )], and the second inequality follows
from Assumption 3. Dividing both sides of the above inequality by [1−2/(β −1)] > 0 (that is, β is required to be larger than
3) and using the definition of β = 1/(Lη), we arrive at




2(m+C)
CLβ(β −1)
s
∗
e
es−1 )−F (x∗ )
E[F (x ) − F (x )] ≤
+
E F (x
(m−1)(β −3) µ(m−1)(β −3)




C(1−Lη)
2(m+C)Lη
es−1 )−F (x∗ ) .
+
E F (x
=
(m−1)(1−3Lη) µη(m−1)(1−3Lη)
This completes the proof.
Although the learning rate η in Theorem 5 needs to be less than 1/(3L), we can use much larger learning rates in
practice, e.g., η = 3/(7L). Similar to Theorem 5, we give the following convergence result for Algorithm 2 with Option II.
Theorem 6 (Option II). Suppose Assumptions 1, 2, and 3 hold, and m is sufficiently large so that

2Lη(m+C)
C(1−Lη)
+
< 1.
m(1−3Lη)
mµη(1−3Lη)
Then Algorithm 2 with Option II has the following geometric convergence in expectation:


es ) − F (x∗ )] ≤ ρsII F (x
e0 ) − F (x∗ ) .
E [F (x
ρII :=

The proof of Theorem 6 can be found in APPENDIX F. In addition, we can also provide the linear convergence
guarantees of Algorithm 2 with Option I or Option II for solving smooth strongly convex objective functions. From all
the results, one can see that VR-SGD attains a linear convergence rate for both smooth and non-smooth strongly convex
minimization problems.

5

E XPERIMENTAL R ESULTS

In this section, we evaluate the performance of our VR-SGD method for solving various ERM problems, such as logistic
regression, Lasso, and ridge regression, and compare its performance with several related stochastic variance reduced
methods, including SVRG [1], Prox-SVRG [2], and Katyusha [3]. All the codes of VR-SGD and related methods can be
downloaded from the author’s website5 .
5. https://sites.google.com/site/fanhua217/publications

19

Option I
Option II
Option III

−2

10

Option I
Option II
Option III

−4

10

10

−8

10

10

−6

10

−8

10

−6

F (xs ) − F (x* )

F (xs ) − F (x* )

−6

10

−8

10

−10

−10

10

−12

−12

0

10

20
30
Gradient evaluations / n

40

10

50

20

40
60
80
Gradient evaluations / n

10

100

−12

0

(a) Adult: λ = 10−5 (left) and λ = 10−6 (right)
−2

5

10
15
20
Gradient evaluations / n

25

10

30

20

40
60
80
Gradient evaluations / n

100

(b) Protein: λ = 10−5 (left) and λ = 10−6 (right)

10

Option I
Option II
Option III

−4

10

Option I
Option II
Option III

−4

10

Option I
Option II
Option III

−2

10

−8

10

−10

−6

10

−8

10

−10

10

10

−12

5

10
15
20
Gradient evaluations / n

25

30

10

−6

10

−8

10

−10

10

20
30
40
50
Gradient evaluations / n

60

10

−6

10

−8

10
10

−12

0

10

−10

10

−12

0

−4

10

F (xs ) − F (x* )

10

F (xs ) − F (x* )

F (xs ) − F (x* )

−6

Option I
Option II
Option III

−2

10

−4

F (xs ) − F (x* )

0

−2

10

10

−8

10

10

−12

0

−6

10

−10

10

−10

10
10

Option I
Option II
Option III

−4

10

−4

10

F (xs ) − F (x* )

−4

F (xs ) − F (x* )

Option I
Option II
Option III

−2

10

−12

0

(c) Covtype: λ = 10−5 (left) and λ = 10−6 (right)

10

20
30
Gradient evaluations / n

40

10

50

0

50
100
Gradient evaluations / n

150

(d) Sido0: λ = 10−4 (left) and λ = 10−5 (right)

Fig. 3. Comparison of Options I, II, and III for solving ridge regression problems with the regularizer (λ/2)k · k2 . In each plot, the vertical axis shows
the objective value minus the minimum, and the horizontal axis is the number of effective passes.

5.1

Experimental Setup

We used four publicly available data sets in the experiments: Adult (also called a9a), Covtype, Protein, and Sido0, as listed
in Table 1. It should be noted that each example of these date sets was normalized so that they have unit length as in [2], [34], which
leads to the same upper bound on the Lipschitz constants Li , i.e., L = Li for all i = 1, . . . , n. As suggested in [1], [2], [3], the epoch
length is set to m = 2n for the stochastic variance reduced methods, SVRG [1], Prox-SVRG [2], and Katyusha [3], as well as
our VR-SGD method. Then the only parameter we have to tune by hand is the step size (or learning rate), η . Since Katyusha
has a much higher per-iteration complexity than SVRG, Prox-SVRG, and VR-SGD, we compare their performance in terms
of both the number of effective passes and running time (seconds), where computing a single full gradient or evaluating

n component gradients is considered as one effective pass over the data. For fair comparison, we implemented SVRG,
Prox-SVRG, Katyusha, and VR-SGD in C++ with a Matlab interface, and performed all the experiments on a PC with an
Intel i5-2400 CPU and 16GB RAM. In addition, we do not compare with other stochastic algorithms such as SAGA [30] and
Catalyst [40], as they have been shown to be comparable or inferior to Katyusha [3].
TABLE 1
Summary of data sets used for our experiments.
Data sets

Sizes n

Dimensions d

Sparsity

Adult

32,562

123

11.28%

Protein

145,751

74

99.21%

Covtype

581,012

54

22.12%

Sido0

12,678

4,932

9.84%

TABLE 2
The three choices of snapshot and starting points for stochastic optimization.
Option I
1
x
es+1 = xsm and xs+
= xsm
0

Option II
x
es+1 =

1 Pm
s
k=1 xk
m

1
and xs+
=
0

Option III
1 Pm
s
k=1 xk
m

x
es+1 =

1 Pm−1 s
k=1 xk
m−1

1
and xs+
= xsm
0

20

Option I
Option II
Option III

−2

10

−4

Option I
Option II
Option III

−4

10

−8

10

−6

10

−8

10

−6

10

F (xs ) − F (x* )

10

10

F (xs ) − F (x* )

−6

−8

10

−10

−10

10

−12

−12

0

10

20
30
Gradient evaluations / n

10

40

10

20
30
40
Gradient evaluations / n

10

50

−8

10

10

−12

0

−6

10

−10

10

−10

10
10

Option I
Option II
Option III

−4

10

−4

10

F (xs ) − F (x* )

F (xs ) − F (x* )

Option I
Option II
Option III

−2

10

−12

0

5

(a) Adult: λ = 10−4 (left) and λ = 10−5 (right)

10
15
20
25
Gradient evaluations / n

30

10

35

0

20

40
60
Gradient evaluations / n

80

(b) Protein: λ = 10−4 (left) and λ = 10−5 (right)

−2

10

Option I
Option II
Option III

−4

Option I
Option II
Option III

−2

−8

10

−2

10

−2

10
F (xs ) − F (x* )

10

F (xs ) − F (x* )

−6

Option I
Option II
Option III

−1

10

10
F (xs ) − F (x* )

F (xs ) − F (x* )

10

Option I
Option II
Option III

−4

10

−10

−3

10

−4

10

−5

10

10
−6

10
−12

10

−6

0

5

10
15
20
Gradient evaluations / n

25

30

0

50

100
150
Gradient evaluations / n

200

0

(c) Covtype: λ = 10−4 (left) and λ = 10−5 (right)

50

100
150
Gradient evaluations / n

200

10

0

50

100
150
Gradient evaluations / n

200

(d) Sido0: λ = 10−4 (left) and λ = 10−5 (right)

Fig. 4. Comparison of Options I, II, and III for solving Lasso problems with the regularizer λkxk1 . In each plot, the vertical axis shows the objective
value minus the minimum, and the horizontal axis is the number of effective passes.

5.2

Different Choices for Snapshot and Starting Points

es and the starting point xs+1
In the practical implementation of SVRG [1], both the snapshot point x
of each epoch are set
0
to the last iterate xsm of the previous epoch (i.e., Option I in Algorithm 1), while the two vectors in [2] (also suggested
1 Pm
s
es and xs+1
in [1]) are set to the average point of the previous epoch, m
0
k=1 xk (i.e., Option II in Algorithm 1). In contrast, x
P
m−1
1
6
s
s
in this paper are set to m−1 k=1 xk and xm (denoted by Option III, i.e., Option I in Algorithms 2 and 3), respectively.
Note that Johnson and Zhang [1] first presented the choices of Options I and II, while the setting of Option III is suggested
in this paper. In the following, we compare the performance of the three choices (i.e., Options I, II and III in Table 2) for
snapshot and starting points for solving ridge regression and Lasso problems, as shown in Figs. 3 and 4. Except for the
three different settings for snapshot and starting points, we use the update rules in (12) and (13) for ridge regression and
Lasso problems, respectively.
From all the results shown in Figs. 3 and 4, we can see that our algorithms with Option III (i.e., Algorithms 2 and 3 with
their Option I) consistently converge much faster than SVRG with the choices of Options I and II for both strongly convex
and non-strongly convex cases. This indicates that the setting of Option III suggested in this paper is a better choice than
Options I and II for stochastic optimization.
5.3

Logistic Regression

In this part, we focus on the following generalized logistic regression problems for binary classification,
n

min

x∈Rd

1X
λ1
log(1 + exp(−bi aTi x)) + kxk2 + λ2 kxk1 ,
n i=1
2

(27)

where {(ai , bi )} is a set of training examples, and λ1 , λ2 ≥ 0 are the regularization parameters. Note that when λ2 > 0,

fi (x) = log(1+exp(−bi aTi x))+(λ1 /2)kxk2 . The formulation (27) includes the `2 -norm (i.e., λ2 = 0), `1 -norm (i.e., λ1 = 0), and
elastic net (i.e., λ1 6= 0 and λ2 6= 0) regularized logistic regression problems. Figs. 5, 6 and 7 show how the objective gap, i.e.,

F (xs )−F (x∗ ), decreases for the `2 -norm, `1 -norm, and elastic net regularized logistic regression problems, respectively.
From all the results, we make the following observations.
6. As Options I and II in Algorithms 2 and 3 achieve very similar performance, we only report the results of Algorithms 2 and 3 with Option I.

21
•

When the regularization parameters λ1 and λ2 are relatively large, e.g., λ1 = 10−4 or λ2 = 10−4 , Prox-SVRG
usually converges faster than SVRG for both strongly convex (e.g., `2 -norm regularized logistic regression) and
non-strongly convex (e.g., `1 -norm regularized logistic regression) cases, as shown in Figs. 5(a)–5(c) and Figs. 6(a)–
6(b). On the contrary, SVRG often outperforms Prox-SVRG, when the regularization parameters are relatively small,
e.g., λ1 = 10−6 or λ2 = 10−6 , as observed in [45]. The main reason is that they have different initialization settings,
s+1
1 Pm
1 Pm
s
s
s
es = xsm and xs+1
es = m
i.e., x
0 = xm for SVRG vs. x
k=1 xk and x0 = m
k=1 xk for Prox-SVRG.

•

Katyusha converges much faster than SVRG and Prox-SVRG for the cases when the regularization parameters are
relatively small, e.g., λ1 = 10−6 , whereas it often achieves similar or inferior performance when the regularization
parameters are relatively large, e.g., λ1 = 10−4 , as shown in Figs. 5(a)–5(d) and Figs. 6(a)–6(b). Note that we
implemented the original algorithms with Option I in [3] for Katyusha. In other words, Katyusha is an accelerated
proximal stochastic gradient method. Obviously, the above observation matches the convergence properties of
Katyusha provided in [3], that is, only if mµ/L ≤ 3/4, Katyusha attains the best known overall complexities of
p
O((n+ nL/µ) log(1/)) for strongly convex problems.

•

Our VR-SGD method consistently converges much faster than SVRG and Prox-SVRG, especially when the regularization parameters are relatively small, e.g., λ1 = 10−6 or λ2 = 10−6 , as shown in Figs. 5(i)–5(l) and Figs. 6(e)–6(f).
The main reason is that VR-SGD can use much larger learning rates than SVRG (e.g., 6/(5L) for VR-SGD vs. 1/(10L)
for SVRG), which leads to faster convergence. This further verifies that the setting of both snapshot and starting
points in our algorithms (i.e., Algorithms 2 and 3) is a better choice than Options I and II in Algorithm 1.

•

In particular, VR-SGD generally outperforms the best-known stochastic method, Katyusha, in terms of the number
of passes through the data, especially when the regularization parameters are relatively large, e.g., 10−4 and 10−5 ,
as shown in Figs. 5(a)–5(h) and Figs. 6(a)–6(d). Since VR-SGD has a much lower per-iteration complexity than
Katyusha, VR-SGD has more obvious advantage over Katyusha in terms of running time. From the algorithms of
Katyusha proposed in [3], one can see that the learning rate of Katyusha is at least set to 1/(3L). Similarly, the
learning rate used in VR-SGD is comparable to Katyusha, which may be the main reason why the performance of
VR-SGD is much better than that of Katyusha. This also implies that the algorithm that enjoys larger learning rates
can yield better performance.

5.4

Common Stochastic Gradient and Prox-SG Updates

In this part, we compare the original Katyusha algorithm, i.e., Algorithm 1 in [3], with the slightly modified Katyusha
algorithm (denoted by Katyusha-I). In Katyusha-I, only the following two update rules for smooth objective functions are
used to replace the original proximal stochastic gradient update rules in (6b) and (6c).
s
e i (xs ) + ∇g(xs )],
yk+1
= yks − η[∇f
k+1
k+1
k
s
e i (xs ) + ∇g(xs )]/(3L).
zk+1
= xsk+1 − [∇f
k+1
k+1
k

(28)

Similarly, we also implement the proximal versions7 for the original SVRG (also called SVRG-I) and the proposed VR-SGD
(denoted by VR-SGD-I) methods, and denote their proximal variants by SVRG-II and VR-SGD-II, respectively. Here, the
original Katyusha method is denoted by Katyusha-II.
Figs. 8 and 9 show the performance of Katyusha-I and Katyusha-II for solving ridge regression problems on the two
popular data sets: Adult and Covtype. We also report the results of SVRG, VR-SGD, and their proximal variants. It is clear
that Katyusha-I usually performs better than Katyusha-II (i.e., the original proximal stochastic method, Katyusha [3]), and
converges significantly faster for the case when the regularization parameter is 10−4 . This seems to be the main reason
why Katyusha has inferior performance when the regularization parameter is relatively large, as shown in Section 5.3. In
7. Here, the proximal variant of SVRG is different from Prox-SVRG [2], and their main difference is the choices of both the snapshot point and
starting point. That is, the two vectors of the former are set to the last iterate xsm , while those of Prox-SVRG are set to the average point of the
1 Pm
s
previous epoch, i.e., m
k=1 xk .

22

contrast, VR-SGD and its proximal variant have similar performance, and the former slightly outperforms the latter in most
cases, as well as SVRG vs. its proximal variant. All this suggests that stochastic gradient update rules as in (12) and (28)
are better choices than proximal stochastic gradient update rules as in (10), (6b) and (6c) for smooth objective functions. We
also believe that our new insight can help us to design accelerated stochastic optimization methods. Both Katyusha-I and
Katyusha-II usually outperform SVRG and its proximal variant, especially when the regularization parameter is relatively
small, e.g., λ = 10−6 , as shown in Figs. 8(d) and 9(d). Unfortunately, both Katyusha-I and Katyusha-II cannot solve the
convex objectives without any regularization term (i.e., λ = 0, as shown in Figs. 8(f) and 9(f)) or with too large and small
regularization parameters, e.g., 10−3 and 10−7 , as shown in Figs. 9(a) and 8(e).
Moreover, it can be seen that both VR-SGD and its proximal variant achieve much better performance than the other
methods in most cases, and are also comparable to Katyusha-I and Katyusha-II in the remaining cases. This further verifies
that VR-SGD is suitable for various large-scale machine learning.
5.5

Fixed and Varied Learning Rates

Finally, we compare the performance of Algorithm 3 with fixed and varied learning rates for solving `1 -norm regularized
logistic regression and Lasso problems, as shown in Fig. 10. Note that the learning rates in Algorithm 3 are varied according
to the update formula in (14). We can observe that Algorithm 3 with varied step-sizes performs very similar to Algorithm 3
with fixed step-sizes in most cases. In the remaining cases, the former slightly outperforms the latter.

6

C ONCLUSIONS

In this paper, we proposed a simple variant of the original SVRG [1], called variance reduced stochastic gradient descent
(VR-SGD). Unlike the choices of the snapshot point and starting point in SVRG and its proximal variant, Prox-SVRG [2],
the two points of each epoch in VR-SGD are set to the average and last iterate of the previous epoch, respectively. This
setting allows us to use much larger learning rates than SVRG, e.g., 3/(7L) for VR-SGD vs. 1/(10L) for SVRG, and also
makes VR-SGD more robust in terms of learning rate selection. Different from existing proximal stochastic methods such as
Prox-SVRG [2] and Katyusha [3], we designed two different update rules for smooth and non-smooth objective functions,
respectively, which makes our VR-SGD method suitable for non-smooth and/or non-strongly convex problems without
using any reduction techniques as in [3], [61]. In contrast, SVRG and Prox-SVRG cannot directly solve non-strongly convex
objectives [42]. Furthermore, our empirical results showed that for smooth problems stochastic gradient update rules as in
(12) are better choices than proximal stochastic gradient update formulas as in (10).
On the practical side, the choices of the snapshot and starting points make VR-SGD significantly faster than its
counterparts, SVRG and Prox-SVRG. On the theoretical side, the setting also makes our convergence analysis more
challenging. We analyzed the convergence properties of VR-SGD for strongly convex objective functions, which show
that VR-SGD attains a linear convergence rate. Moreover, we also provided the convergence guarantees of VR-SGD for
non-strongly convex problems, which show that VR-SGD achieves a sub-linear convergence rate. All these results imply
that VR-SGD is guaranteed to have a similar convergence rate with other variance reduced stochastic methods such as
SVRG and Prox-SVRG, and a slower theoretical rate than accelerated methods such as Katyusha [3]. Nevertheless, various
experimental results show that VR-SGD significantly outperforms than SVRG and Prox-SVRG, and is still much better than
the best known stochastic method, Katyusha [3].

A PPENDIX A: P ROOF OF L EMMA 2
Before proving Lemma 2, we first give and prove the following lemma.
Lemma 9. Suppose each fi (·) is L-smooth, and let x∗ be the optimal solution of Problem (1) when g(x) ≡ 0, then we have
h
i
2
E k∇fi (x) − ∇fi (x∗ )k ≤ 2L [f (x) − f (x∗ )] .

23

Proof. Following Theorem 2.1.5 in [19] and Lemma 3.4 [2], we have
2

k∇fi (x) − ∇fi (x∗ )k ≤ 2L [fi (x) − fi (x∗ ) − h∇fi (x∗ ), x − x∗ i] .
Summing the above inequality over i = 1, . . . , n, we obtain
n
h
i 1X
2
2
k∇fi (x)−∇fi (x∗ )k ≤ 2L [f (x) − f (x∗ ) − h∇f (x∗ ), x − x∗ i] .
E k∇fi (x)−∇fi (x∗ )k =
n i=1
By the optimality of x∗ , i.e., x∗ = arg minx f (x), we have ∇f (x∗ ) = 0. Then
i
h
2
E k∇fi (x) − ∇fi (x∗ )k ≤ 2L [f (x) − f (x∗ ) − h∇f (x∗ ), x − x∗ i]

= 2L [f (x) − f (x∗ )] .

Proof of Lemma 2:
Proof.


i
h
2
s−1
s−1
s
e is (xs ) − ∇f (xs )k2 = E ∇fis (xs ) − ∇fis (x
e
e
)
+
∇f
(
x
)
−
∇f
(x
)
E k∇f
k
k
k
k
k
k
k


2
2
es−1 )
es−1 )
= E ∇fisk (xsk ) − ∇fisk (x
− ∇f (xsk ) − ∇f (x


2
es−1 )
≤ E ∇fisk (xsk ) − ∇fisk (x



2
s
∗
s
s
es−1 ) − ∇fisk (x∗ )
≤ 2E ∇fik (xk ) − ∇fik (x )
+ 2E ∇fisk (x

2





es−1 ) − f (x∗ ) ,
≤ 4L f (xsk ) − f (x∗ ) + f (x
where the second equality holds due to the fact that E[kx− E[x]k2 ] = E[kxk2 ]−kE[x]k2 ; the second inequality holds due to
the fact that ka − bk2 ≤ 2(kak2 + kbk2 ); and the last inequality follows from Lemma 9.

A PPENDIX B: P ROOF OF L EMMA 5
es−1 )+∇f (x
es−1 ). Since each
Proof. For convenience, the stochastic gradient estimator is defined as: vks := ∇fisk (xsk )−∇fisk (x
component function fi (x) is L-smooth, which implies that the gradient of the average function f (x) is also L-smooth, i.e.,
for all x, y ∈ Rd ,

k∇f (x) − ∇f (y)k ≤ Lkx − yk,
whose equivalent form is

f (y) ≤ f (x) + h∇f (x), y − xi +

L
ky − xk2 .
2

Using the above smoothness inequality, we have

L s
2
x
− xsk
2 k+1
Lβ s
L(β −1) s
2
= f (xsk ) + ∇f (xsk ), xsk+1 − xsk +
xk+1 − xsk −
xk+1 − xsk
2
2
Lβ s
= f (xsk ) + vks , xsk+1 − xsk +
kxk+1 − xsk k2
2
L(β −1) s
kxk+1 − xsk k2 .
+ ∇f (xsk ) − vks , xsk+1 − xsk −
2
Using Lemma 2, then we get


L(β −1) s
E ∇f (xsk ) − vks , xsk+1 − xsk −
kxk+1 − xsk k2
2


1
L(β −1) s
L(β −1) s
s
s 2
≤E
k∇f (xk ) − vk k +
kxk+1 −xsk k2 −
kxk+1 −xsk k2
2L(β −1)
2
2

2 
es−1 ) − f (x∗ ) ,
≤
f (xsk ) − f (x∗ ) + f (x
β −1
f (xsk+1 ) ≤ f (xsk ) + ∇f (xsk ), xsk+1 − xsk +

2

(29)

(30)

24

where the first inequality holds due to the Young’s inequality (i.e., y T z ≤ kyk2 /(2θ)+θkzk2 /2 for all θ > 0 and y, z ∈ Rd ),
and the second inequality follows from Lemma 2.
Taking the expectation over the random choice of isk and substituting the inequality in (30) into the inequality in (29),
we have
E[f (xsk+1 )]



2 
Lβ s
es−1 ) − f (x∗ )
kxk+1 − xsk k2 +
f (xsk ) − f (x∗ ) + f (x
≤ E[f (xsk )] + E vks , xsk+1 − xsk +
2
β −1



Lβ
2 
s
∗
s 2
∗
s
∗
s
es−1 ) − f (x∗ )
≤ E[f (xk )] + E hvk , x − xk i +
(kx − xk k − kx − xsk+1 k2 ) +
f (xsk ) − f (x∗ ) + f (x
2
β −1



Lβ
2 
es−1 ) − f (x∗ )
≤ E[f (xsk )] + h∇f (xsk ), x∗ − xsk i + E
(kx∗ − xsk k2 − kx∗ − xsk+1 k2 ) +
f (xsk ) − f (x∗ ) + f (x
2
β −1



Lβ
2 
s
∗
∗
∗
s 2
∗
s
2
es−1 ) − f (x∗ )
f (xk ) − f (x ) + f (x
(kx − xk k − kx − xk+1 k ) +
≤ f (x ) + E
2
β −1



Lβ
2 
es−1 ) − f (x∗ ) ,
= f (x∗ ) +
E kx∗ − xsk k2 − kx∗ − xsk+1 k2 +
f (xsk ) − f (x∗ ) + f (x
2
β −1
where the first inequality holds due to the inequality in (29) and the inequality in (30); the second inequality follows from
Lemma 1 with ẑ = xsk+1 , z = x∗ , z0 = xsk , τ = Lβ = 1/η , and r(z) := hvks , z −xsk i; the third inequality holds due to the fact
that E[vks ] = ∇f (xsk ); and the last inequality follows from the convexity of f (·), i.e., f (xsk )+h∇f (xsk ), x∗ −xsk i ≤ f (x∗ ). The
above inequality can be rewritten as follows:
E[f (xsk+1 )] − f (x∗ ) ≤

 Lβ  ∗

2 
es−1 ) − f (x∗ ) +
f (xsk ) − f (x∗ ) + f (x
E kx − xsk k2 − kx∗ − xsk+1 k2 .
β −1
2

Summing the above inequality over k = 0, 1, . . . , m−1, we obtain
m−1
X

≤



E[f (xsk+1 )] − f (x∗ )

k=0
m−1
X
k=0


 Lβ  ∗

2 
s
∗
s−1
∗
s 2
∗
s
2
e
f (xk ) − f (x ) + f (x ) − f (x ) +
E kx − xk k − kx − xk+1 k
.
β −1
2

Then
X
m
m
2
2 X
1−
{E[f (xsk )] − f (x∗ )} +
{E[f (xsk )] − f (x∗ )}
β −1 k=1
β −1 k=1

m 
X


Lβ  ∗
2 
es−1 ) − f (x∗ ) +
f (xsk−1 ) − f (x∗ ) + f (x
E kx − xs0 k2 − kx∗ − xsm k2 .
≤
β
−1
2
k=1
1 Pm
s
es = m
Due to the setting of x
k=1 xk in Option II, and the convexity of f (·), then


es ) ≤
f (x

(31)

m
1 X
f (xsk ).
m k=1

Using the above inequality, the inequality in (31) becomes


m
2
2 X
es ) − f (x∗ )] +
m 1−
E[f (x
{E[f (xsk )] − f (x∗ )}
β −1
β −1 k=1

m 
X


2 
Lβ  ∗
es−1 ) − f (x∗ ) +
≤
f (xsk−1 ) − f (x∗ ) + f (x
E kx − xs0 k2 − kx∗ − xsm k2 .
β −1
2
k=1
Pm−1
2
s
∗
Dividing both sides of the above inequality by m and subtracting (β−1)m
k=1 [f (xk )−f (x )] from both sides, we arrive
at



2
2
es ) − f (x∗ )] +
E[f (x
E[f (xsm ) − f (x∗ )]
β −1
(β −1)m



2
2
Lβ  ∗
es−1 ) − f (x∗ ) +
≤
E f (x
E[f (xs0 ) − f (x∗ )] +
E kx − xs0 k2 − kx∗ − xsm k2 .
(β −1)
(β −1)m
2m

1−

This completes the proof.

25

A PPENDIX C: P ROOF OF T HEOREM 2
Proof. Using Lemma 5, we have

2
2
es ) − f (x∗ )] +
)E[f (x
E[f (xsm ) − f (x∗ )]
β −1
(β −1)m

 Lβ
2
2
es−1 ) − f (x∗ ) +
≤
E[f (xs0 ) − f (x∗ )] +
E f (x
E[kx∗ − xs0 k2 − kx∗ − xsm k2 ].
(β −1)m
β −1
2m
(1 −

Summing the above inequality over s = 1, 2, . . . , S , taking expectation with respect to the history of random variables isk ,
and using the setting of xs+1
= xsm in Algorithm 3, we arrive at
0
S
X

≤

(1−

s=1
S 
X
s=1

+

2
es ) − f (x∗ )]
)E[f (x
β −1




2
2
es−1 ) − f (x∗ )
E{f (xs0 ) − f (x∗ ) − [f (xsm ) − f (x∗ )]} +
E f (x
(β −1)m
β −1

S

Lβ X  ∗
E kx − xs0 k2 − kx∗ − xsm k2 .
2m s=1

PS

− f (x∗ )] from both sides of the above inequality, we obtain

S 
X
4
es ) − f (x∗ )]
1−
E[f (x
β
−1
s=1
n
o
2
2
∗
e0 ) − f (x
eS )]
≤
E f (x10 ) − f (x∗ ) − [f (xS
E[f (x
m ) − f (x )] +
(β −1)m
β −1
i
Lβ h ∗
2
+
E kx − x10 k2 − kx∗ − xS
.
mk
2m
e0 )−f (x
eS )] ≤ f (x
e0 )−f (x∗ ). Dividing both sides of the above inequality by S , and using
It is not hard to verify that E[f (x

Subtracting

2
β−1

es )
s=1 [f (x

e0 = x10 , we have
the choice x

X
S
1
4
es ) − f (x∗ )]
E[f (x
1−
S
β −1 s=1
n
o
2
Lβ
2
∗
e0 ) − f (x∗ )] +
E f (x10 ) − f (x∗ ) − [f (xS
[f (x
kx∗ − x10 k2
≤
m ) − f (x )] +
(β −1)mS
(β −1)S
2mS
(32)


2
2
Lβ
0
∗
0
∗
∗
0 2
e ) − f (x ) +
e ) − f (x )] +
e k
≤
f (x
[f (x
kx − x
(β −1)mS
(β −1)S
2mS
2(m+1)
Lβ
e0 ) − f (x∗ )] +
e0 − x∗ k2
=
[f (x
kx
(β −1)mS
2mS


2
e0 )−f (x
eS )] ≤ f (x
e0 )−f (x∗ ) and E kx∗ −xS
where the first inequality holds due to the facts that E[f (x
m k ≥ 0, and the last


∗
e0 = x10 .
inequality uses the facts that E f (xS
m )−f (x ) ≥ 0 and x
P
1 PS
S
es , and using the convexity of f (·), we have f (xS ) ≤ S1 S
es ), and therefore the inequality in
Since x = S s=1 x
s=1 f (x
(32) becomes


 h

X
S
i
4
1
4
S
∗
es ) − f (x∗ )]
1−
E f (x ) − f (x ) ≤
1−
E[f (x
β −1
S
β −1 s=1

≤

Lβ
2(m+1)
e0 ) − f (x∗ )] +
e0 − x∗ k2 .
[f (x
kx
(β −1)mS
2mS

4
Dividing both sides of the above inequality by (1− β−1
) > 0 (i.e., η < 1/(5L)), we arrive at
h
i
Lβ(β −1)
2(m+1)
e0 ) − f (x∗ )] +
e0 − x∗ k2 .
E f (xS ) − f (x∗ ) ≤
[f (x
kx
mS(β −5)
2mS(β −5)

eS ) ≤ f (xS ), then x
bS = x
eS , and
If f (x
h
i
h
i
bS ) − f (x∗ ) ≤ E f (xS ) − f (x∗ )
E f (x

≤

2(m+1)
Lβ(β −1)
e0 ) − f (x∗ )] +
e0 − x∗ k2 .
[f (x
kx
mS(β −5)
2mS(β −5)

26

bS = xS , and the above inequality still holds.
eS ) ≥ f (xS ), then x
Alternatively, if f (x
This completes the proof.

A PPENDIX D: P ROOF OF L EMMA 6
Proof. Since the average function f (x) is L-smooth, then for all x, y ∈ Rd ,

f (y) ≤ f (x) + h∇f (x), y − xi +

L
ky − xk2 ,
2

which then implies

f (xsk+1 ) ≤ f (xsk ) + ∇f (xsk ), xsk+1 − xsk +

L s
x
− xsk
2 k+1

2

.

Using the above inequality, we have

F (xsk+1 ) = f (xsk+1 ) + g(xsk+1 )
≤ f (xsk ) + g(xsk+1 ) + ∇f (xsk ), xsk+1 − xsk +

Lβ s
xk+1 − xsk
2

2

−

L(β −1) s
xk+1 − xsk
2

2

Lβ s
kxk+1 − xsk k2
2
L(β −1) s
−
kxk+1 − xsk k2 .
2

= f (xsk ) + g(xsk+1 ) + vks , xsk+1 − xsk +
+ ∇f (xsk ) − vks , xsk+1 − xsk

According to Lemma 5, then we obtain


L(β −1) s
E ∇f (xsk ) − vks , xsk+1 − xsk −
kxk+1 − xsk k2
2


1
L(β −1) s
L(β −1) s
s
s 2
≤E
k∇f (xk ) − vk k +
kxk+1 −xsk k2 −
kxk+1 −xsk k2
2L(β −1)
2
2

2 
es−1 ) − F (x∗ ) ,
F (xsk ) − F (x∗ ) + F (x
≤
β −1

(33)

(34)

where the first inequality holds due to the Young’s inequality, and the second inequality follows from Lemma 5. Substituting
the inequality (34) into the inequality (33), and taking the expectation over the random choice isk , we arrive at






Lβ s
E F (xsk+1 ) ≤ E[f (xsk )] + E g(xsk+1 ) + E vks , xsk+1 − xsk +
kxk+1 − xsk k2
2

2 
s
∗
s−1
∗
e
F (xk ) − F (x ) + F (x
) − F (x )
+
β −1


Lβ
s
∗
s
∗
s
∗
s 2
∗
s
2
≤ E[f (xk )] + g(x ) + E hvk , x − xk i +
(kx − xk k − kx − xk+1 k )
2

2 
es−1 ) − F (x∗ )
F (xsk ) − F (x∗ ) + F (x
+
β −1


Lβ
≤ f (x∗ ) + g(x∗ ) + E
(kx∗ − xsk k2 − kx∗ − xsk+1 k2 )
2

2 
es−1 ) − F (x∗ )
+
F (xsk ) − F (x∗ ) + F (x
β −1


Lβ  ∗
2 
es−1 ) − F (x∗ ) ,
= F (x∗ ) +
E kx − xsk k2 − kx∗ − xsk+1 k2 +
F (xsk ) − F (x∗ ) + F (x
2
β −1
where the first inequality holds due to the inequality (33) and the inequality (34); the second inequality follows from
Lemma 1 with ẑ = xsk+1 , z = x∗ , z0 = xsk , τ = Lβ = 1/η , and r(z) := hvks , z −xsk i+g(z); and the third inequality holds due to
the fact that E[vks ] = ∇f (xsk ) and the convexity of f (·), i.e., f (xsk )+h∇f (xsk ), x∗ −xsk i ≤ f (x∗ ). Then the above inequality is
rewritten as follows:
E[F (xsk+1 )] − F (x∗ )
 Lβ  ∗

2 
es−1 ) − F (x∗ ) +
≤
F (xsk ) − F (x∗ ) + F (x
E kx − xsk k2 − kx∗ − xsk+1 k2 .
β −1
2

(35)

27

Summing the above inequality over k = 0, 1, . . . , (m−1) and taking expectation over whole history, we have
m−1
X



E[F (xsk+1 )] − F (x∗ )

k=0
m−1
X


 Lβ  ∗

2 
es−1 ) − F (x∗ ) +
F (xsk ) − F (x∗ ) + F (x
E kx − xsk k2 − kx∗ − xsk+1 k2 .
β −1
2
k=0


2 Pm−2
s
∗
Subtracting β−1
k=0 E F (xk+1 )−F (x ) from both sides of the above inequality, we obtain

≤

m−1
X

m−1

2 X 
2
E[F (xsm )−F (x∗ )]
E F (xsk+1 )−F (x∗ ) +
β
−1
β
−1
k=0
k=0

m−2
m−1
X 2 

 Lβ  ∗

2 X 
es−1 )−F (x∗ ) −
F (xsk )−F (x∗ )+F (x
E F (xsk+1 )−F (x∗ ) +
E kx −xs0 k2 −kx∗ −xsm k2 .
≤
β −1
β −1 k=0
2
k=0



E F (xsk+1 ) − F (x∗ ) −

Then


2
β −1

X
m

2
E[F (xsm )−F (x∗ )]
β
−1
k=1
 Lβ  ∗

2
2m 
s
es−1 )−F (x∗ ) +
≤
E[F (x0 ) − F (x∗ )] +
E F (x
E kx −xs0 k2 −kx∗ −xsm k2 .
β −1
β −1
2
P
m
s+1
1
es = m k=1 xsk and x0 = xsm , and the convexity of the objective function F (·), we have F (x
es ) ≤
Due to the settings of x
1 Pm
s
k=1 F (xk ), and
m


2
2
es ) − F (x∗ )] +
E[F (xsm )−F (x∗ )]
E[F (x
m 1−
β −1
β −1
 Lβ  ∗

2
2m 
es−1 )−F (x∗ ) +
≤
E[F (xs0 ) − F (x∗ )] +
E F (x
E kx −xs0 k2 −kx∗ −xsm k2 .
β −1
β −1
2
1−

E[F (xsk ) − F (x∗ )] +

Dividing both sides of the above inequality by m, we arrive at


2
2
es ) − F (x∗ )] +
E[F (x
E[F (xsm )−F (x∗ )]
1−
β −1
(β −1)m

 Lβ  ∗

2
2
es−1 )−F (x∗ ) +
≤
E[F (xs0 ) − F (x∗ )] +
E F (x
E kx −xs0 k2 −kx∗ −xsm k2 .
(β −1)m
β −1
2m
This completes the proof.

A PPENDIX E: P ROOF OF T HEOREM 3
Proof. Summing the inequality in (21) over s = 1, 2, . . . , S , and taking expectation with respect to the history of isk , we have

S 
S
X
X
2
2
es ) − F (x∗ )] +
1−
E[F (x
E[F (xsm ) − F (x∗ )]
β
−1
(β
−1)m
s=1
s=1

S 
X


2
2
es−1 ) − F (x∗ )
≤
E[F (xs0 ) − F (x∗ )] +
E F (x
(β −1)m
β −1
s=1
S

Lβ X  ∗
E kx − xs0 k2 − kx∗ − xsm k2 .
2m s=1
PS
2
2 PS
es )−F (x∗ )] from both sides of the above inequality, and using
Subtracting s=1 (β−1)m
E[F (xsm )−F (x∗ )]+ β−1
s=1 [F (x

+

the setting of xs+1
= xsm , we arrive at
0

S 
X
4
es ) − F (x∗ )]
1−
E[F (x
β
−1
s=1
h
i
h
i Lβ h
i
2
2
2
e0 )−F (x
eS ) +
≤
E F (x10 )−F (xS
E F (x
E kx∗ −x10 k2 − kx∗ −xS
.
m) +
mk
(β −1)m
β −1
2m

28

e0 )−F (x
eS )] ≤ F (x
e0 )−F (x∗ ). Dividing both sides of the above inequality by S , and using the
It is easy to verify that E[F (x
e0 = x10 , we obtain
choice x


S
1X
es ) − F (x∗ )]
E[F (x
S s=1
h
i
2
2
Lβ
e0 ) − F (x
eS )] +
≤
E F (x10 )−F (xS
[F (x
kx∗ − x10 k2
m) +
(β −1)mS
(β −1)S
2mS


2
2
Lβ
e0 ) − F (x∗ ) +
e0 ) − F (x∗ )] +
e0 k2
≤
F (x
[F (x
kx∗ − x
(β −1)mS
(β −1)S
2mS
Lβ
2(m+1)
e0 ) − F (x∗ )] +
e0 − x∗ k2
[F (x
kx
=
(β −1)mS
2mS

1−

4
β −1



(36)



2
1
S
where the first inequality uses the fact that kx∗ −xS
m k ≥ 0; and the last inequality holds due to the facts that E F (x0 )−F (xm )
e0 )−F (x
eS )] ≤ F (x
e0 )−F (x∗ ), and x
e0 = x10 .
≤ F (x10 )−F (x∗ ), E[F (x
P
P
S
es and the convexity of the objective function F (·), we have F (xS ) ≤ S1 S
es ),
Using the definition of xS = S1 s=1 x
s=1 F (x
and therefore we can rewrite the above inequality in (36) as

 h
 X

S
i
4
1
4
S
∗
es ) − F (x∗ )]
E[F (x
1−
E F (x ) − F (x ) ≤ 1 −
β −1
β −1 S s=1

≤

Lβ
2(m+1)
e0 ) − F (x∗ )] +
e0 − x∗ k2 .
[F (x
kx
(β −1)mS
2mS

4
Dividing both sides of the above inequality by (1− β−1
) > 0, we have

i
h
2(m+1)
β(β −1)L
e0 ) − F (x∗ )] +
e0 − x∗ k2 .
E F (xS ) − F (x∗ ) ≤
[F (x
kx
(β −5)mS
2(β −5)mS
eS ) ≤ F (xS ), then x
bS = x
eS , and
When F (x
h
i
bS ) − F (x∗ ) ≤
E F (x

β(β −1)L
2(m+1)
e0 ) − F (x∗ )] +
e0 − x∗ k2 .
[F (x
kx
(β −5)mS
2(β −5)mS

eS ) ≥ F (xS ), then x
bS = xS , and the above inequality still holds.
Alternatively, if F (x
This completes the proof.

A PPENDIX F: P ROOF OF T HEOREM 6
Proof. Since each fi (·) is convex and L-smooth, then we have



 Lβ  s+1

2
2
∗
es ) − F (x∗ )] +
1−
E[F (x
E F (xs+1
E kx0 −x∗ k2
0 )−F (x ) +
β −1
m(β −1)
2m




2
Lβ
2
es−1 )−F (x∗ ) +
E F (x
E[F (xs0 ) − F (x∗ )] +
E kxs0 −x∗ k2
≤
β −1
m(β −1)
2m




2
2
Lβ
s
s−1
∗
e
≤
E F (x
)−F (x ) +
+
E[F (x0 ) − F (x∗ )]
β −1
m(β −1) mµ




2(m+C) CLβ
es−1 ) − F (x∗ )
≤
+
E F (x
m(β −1)
mµ
where the first inequality follows from Lemma 6; the second inequality holds due to the fact that kxs0 −x∗ k2 ≤ (2/µ)[F (xs0 )−

F (x∗ )]; and the last inequality follows from Assumption 3.
Due to the definition of β = 1/(Lη), the above inequality is rewritten as follows:




1−3Lη
2Lη(m+C)
C
s
∗
e ) − F (x )] ≤
es−1 ) − F (x∗ ) .
E[F (x
+
E F (x
1−Lη
m(1−Lη)
mµη
Dividing both sides of the above inequality by (1−3Lη)(1−Lη) > 0, we arrive at




2Lη(m+C)
C(1−Lη)
s
∗
e ) − F (x )] ≤
es−1 ) − F (x∗ ) .
E[F (x
+
E F (x
m(1−3Lη)
mµη(1−3Lη)
This completes the proof.

29

R EFERENCES
[1]

R. Johnson and T. Zhang, “Accelerating stochastic gradient descent using predictive variance reduction,” in Proc. Adv. Neural Inf. Process.
Syst., 2013, pp. 315–323.

[2]

L. Xiao and T. Zhang, “A proximal stochastic gradient method with progressive variance reduction,” SIAM J. Optim., vol. 24, no. 4, pp.
2057–2075, 2014.

[3]
[4]

Z. Allen-Zhu, “Katyusha: The first direct acceleration of stochastic gradient methods,” in Proc. 49th ACM Symp. Theory Comput., 2017.
A. Krizhevsky, I. Sutskever, and G. E. Hinton, “ImageNet classification with deep convolutional neural networks,” in Proc. Adv. Neural Inf.
Process. Syst., 2012, pp. 1097–1105.

[5]

I. Sutskever, J. Martens, G. Dahl, and G. Hinton, “On the importance of initialization and momentum in deep learning,” in Proc. 30th Int.
Conf. Mach. Learn., 2013, pp. 1139–1147.

[6]

Z. Allen-Zhu and E. Hazan, “Variance reduction for faster non-convex optimization,” in Proc. 33rd Int. Conf. Mach. Learn., 2016, pp. 699–707.

[7]

H. Ouyang, N. He, L. Q. Tran, and A. Gray, “Stochastic alternating direction method of multipliers,” in Proc. 30th Int. Conf. Mach. Learn., 2013,
pp. 80–88.

[8]

Y. Liu, F. Shang, and J. Cheng, “Accelerated variance reduced stochastic ADMM,” in Proc. 31st AAAI Conf. Artif. Intell., 2017, pp. 2287–2293.

[9]

C. Qu, Y. Li, and H. Xu, “Linear convergence of SVRG in statistical estimation,” arXiv:1611.01957v2, 2017.

[10] C. Paquette, H. Lin, D. Drusvyatskiy, J. Mairal, and Z. Harchaoui, “Catalyst acceleration for gradient-based non-convex optimization,”
arXiv:1703.10993, 2017.
[11] J. Duchi and F. Ruan, “Stochastic methods for composite optimization problems,” arXiv:1703.08570, 2017.
[12] B. Recht and C. Ré, “Parallel stochastic gradient algorithms for large-scale matrix completion,” Math. Prog. Comp., vol. 5, pp. 201–226, 2013.
[13] X. Zhang, L. Wang, and Q. Gu, “Stochastic variance-reduced gradient descent for low-rank matrix recovery from linear measurements,”
arXiv:1701.00481v2, 2017.
[14] M. Schmidt, R. Babanezhad, M. Ahmed, A. Defazio, A. Clifton, and A. Sarkar, “Non-uniform stochastic average gradient method for training
conditional random fields,” in Proc. 18th Int. Conf. Artif. Intell. Statist., 2015, pp. 819–828.
[15] O. Shamir, “A stochastic PCA and SVD algorithm with an exponential convergence rate,” in Proc. 32nd Int. Conf. Mach. Learn., 2015, pp.
144–152.
[16] D. Garber, E. Hazan, C. Jin, S. M. Kakade, C. Musco, P. Netrapalli, and A. Sidford, “Faster eigenvector computation via shift-and-invert
preconditioning,” in Proc. 33rd Int. Conf. Mach. Learn., 2016, pp. 2626–2634.
[17] Z. Allen-Zhu and Y. Li, “Doubly accelerated methods for faster CCA and generalized eigendecomposition,” arXiv:1607.06017v2, 2017.
[18] Y. Nesterov, “A method of solving a convex programming problem with convergence rate O(1/k2 ),” Soviet Math. Doklady, vol. 27, pp.
372–376, 1983.
[19] ——, Introductory Lectures on Convex Optimization: A Basic Course.

Boston: Kluwer Academic Publ., 2004.

[20] P. Tseng, “On accelerated proximal gradient methods for convex-concave optimization,” Technical report, University of Washington, 2008.
[21] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci., vol. 2, no. 1,
pp. 183–202, 2009.
[22] H. Robbins and S. Monro, “A stochastic approximation method,” Ann. Math. Statist., vol. 22, no. 3, pp. 400–407, 1951.
[23] T. Zhang, “Solving large scale linear prediction problems using stochastic gradient descent algorithms,” in Proc. 21st Int. Conf. Mach. Learn.,
2004, pp. 919–926.
[24] C. Hu, J. T. Kwok, and W. Pan, “Accelerated gradient methods for stochastic optimization and online learning,” in Proc. Adv. Neural Inf.
Process. Syst., 2009, pp. 781–789.
[25] S. Bubeck, “Convex optimization: Algorithms and complexity,” Found. Trends Mach. Learn., vol. 8, pp. 231–358, 2015.
[26] A. Rakhlin, O. Shamir, and K. Sridharan, “Making gradient descent optimal for strongly convex stochastic optimization,” in Proc. 29th Int.
Conf. Mach. Learn., 2012, pp. 449–456.
[27] O. Shamir and T. Zhang, “Stochastic gradient descent for non-smooth optimization: Convergence results and optimal averaging schemes,”
in Proc. 30th Int. Conf. Mach. Learn., 2013, pp. 71–79.
[28] N. L. Roux, M. Schmidt, and F. Bach, “A stochastic gradient method with an exponential convergence rate for finite training sets,” in Proc.
Adv. Neural Inf. Process. Syst., 2012, pp. 2672–2680.
[29] S. Shalev-Shwartz and T. Zhang, “Stochastic dual coordinate ascent methods for regularized loss minimization,” J. Mach. Learn. Res., vol. 14,
pp. 567–599, 2013.
[30] A. Defazio, F. Bach, and S. Lacoste-Julien, “SAGA: A fast incremental gradient method with support for non-strongly convex composite
objectives,” in Proc. Adv. Neural Inf. Process. Syst., 2014, pp. 1646–1654.
[31] Y. Zhang and L. Xiao, “Stochastic primal-dual coordinate method for regularized empirical risk minimization,” in Proc. 32nd Int. Conf. Mach.
Learn., 2015, pp. 353–361.
[32] M. Schmidt, N. L. Roux, and F. Bach, “Minimizing finite sums with the stochastic average gradient,” INRIA, Paris, Tech. Rep., 2013.
[33] S. Shalev-Shwartz and T. Zhang, “Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization,” Math. Program.,
vol. 155, pp. 105–145, 2016.
[34] F. Shang, Y. Liu, J. Cheng, and J. Zhuo, “Fast stochastic variance reduced gradient method with momentum acceleration for machine
learning,” arXiv:1703.07948, 2017.

30

[35] L. Zhang, M. Mahdavi, and R. Jin, “Linear convergence with condition number independent access of full gradients,” in Proc. Adv. Neural
Inf. Process. Syst., 2013, pp. 980–988.
[36] J. Konečný, J. Liu, P. Richtárik, , and M. Takáč, “Mini-batch semi-stochastic gradient descent in the proximal setting,” IEEE J. Sel. Top. Sign.
Proces., vol. 10, no. 2, pp. 242–255, 2016.
[37] A. Nitanda, “Stochastic proximal gradient descent with acceleration techniques,” in Proc. Adv. Neural Inf. Process. Syst., 2014, pp. 1574–1582.
[38] G. Lan and Y. Zhou, “An optimal randomized incremental gradient method,” arXiv:1507.02000v3, 2015.
[39] R. Frostig, R. Ge, S. M. Kakade, and A. Sidford, “Un-regularizing: approximate proximal point and faster stochastic algorithms for empirical
risk minimization,” in Proc. 32nd Int. Conf. Mach. Learn., 2015, pp. 2540–2548.
[40] H. Lin, J. Mairal, and Z. Harchaoui, “A universal catalyst for first-order optimization,” in Proc. Adv. Neural Inf. Process. Syst., 2015, pp.
3366–3374.
[41] R. Babanezhad, M. O. Ahmed, A. Virani, M. Schmidt, J. Konecny, and S. Sallinen, “Stop wasting my gradients: Practical SVRG,” in Proc. Adv.
Neural Inf. Process. Syst., 2015, pp. 2242–2250.
[42] Z. Allen-Zhu and Y. Yuan, “Improved SVRG for non-strongly-convex or sum-of-non-convex objectives,” in Proc. 33rd Int. Conf. Mach. Learn.,
2016, pp. 1080–1089.
[43] M. Frank and P. Wolfe, “An algorithm for quadratic programming,” Naval Res. Logist. Quart., vol. 3, pp. 95–110, 1956.
[44] E. Hazan and H. Luo, “Variance-reduced and projection-free stochastic optimization,” in Proc. 33rd Int. Conf. Mach. Learn., 2016, pp. 1263–1271.
[45] F. Shang, Y. Liu, J. Cheng, K. W. Ng, and Y. Yoshida, “Variance reduced stochastic gradient descent with sufficient decrease,” arXiv:1703.06807,
2017.
[46] L. Hien, C. Lu, H. Xu, and J. Feng, “Accelerated stochastic mirror descent algorithms for composite non-strongly convex optimization,”
arXiv:1605.06892v2, 2016.
[47] B. Woodworth and N. Srebro, “Tight complexity bounds for optimizing composite objectives,” in Proc. Adv. Neural Inf. Process. Syst., 2016,
pp. 3639–3647.
[48] A. Defazio, “A simple practical accelerated method for finite sums,” in Proc. Adv. Neural Inf. Process. Syst., 2016, pp. 676–684.
[49] X. Li, T. Zhao, R. Arora, H. Liu, and J. Haupt, “Nonconvex sparse learning via stochastic optimization with progressive variance reduction,”
in Proc. 33rd Int. Conf. Mach. Learn., 2016, pp. 917–925.
[50] L. Nguyen, J. Liu, K. Scheinberg, and M. Takáč, “SARAH: A novel method for machine learning problems using stochastic recursive
gradient,” arXiv:1703.00102v1, 2017.
[51] L. W. Zhong and J. T. Kwok, “Fast stochastic alternating direction method of multipliers,” in Proc. 31st Int. Conf. Mach. Learn., 2014, pp. 46–54.
[52] S. Zheng and J. T. Kwok, “Fast-and-light stochastic ADMM,” in Proc. 25th Int. Joint Conf. Artif. Intell., 2016, pp. 2407–2613.
[53] S. J. Reddi, S. Sra, B. Poczos, and A. Smola, “Proximal stochastic methods for nonsmooth nonconvex finite-sum optimization,” in Proc. Adv.
Neural Inf. Process. Syst., 2016, pp. 1145–1153.
[54] S. Reddi, A. Hefny, S. Sra, B. Poczos, and A. Smola, “On variance reduction in stochastic gradient descent and its asynchronous variants,” in
Proc. Adv. Neural Inf. Process. Syst., 2015, pp. 2629–2637.
[55] J. D. Lee, Q. Lin, T. Ma, and T. Yang, “Distributed stochastic variance reduced gradient methods and a lower bound for communication
complexity,” arXiv:1507.07595v2, 2016.
[56] Y. Bengio, “Learning deep architectures for AI,” Found. Trends Mach. Learn., vol. 2, no. 1, pp. 1–127, 2009.
[57] R. Ge, F. Huang, C. Jin, and Y. Yuan, “Escaping from saddle points online stochastic gradient for tensor decomposition,” in Proc. 28th Conf.
Learn. Theory, 2015, pp. 797–842.
[58] A. Neelakantan, L. Vilnis, Q. V. Le, I. Sutskever, L. Kaiser, K. Kurach, and J. Martens, “Adding gradient noise improves learning for very
deep networks,” arXiv:1511.06807, 2015.
[59] N. Flammarion and F. Bach, “From averaging to acceleration, there is only a step-size,” in Proc. 28th Conf. Learn. Theory, 2015, pp. 658–695.
[60] C. Tan, S. Ma, Y. Dai, and Y. Qian, “Barzilai-Borwein step size for stochastic gradient descent,” in Proc. Adv. Neural Inf. Process. Syst., 2016,
pp. 685–693.
[61] Z. Allen-Zhu and E. Hazan, “Optimal black-box reductions between optimization objectives,” in Proc. Adv. Neural Inf. Process. Syst., 2016, pp.
1606–1614.
[62] B. Carpenter, “Lazy sparse stochastic gradient descent for regularized multinomial logistic regression,” Tech. Rep., 2008.
[63] J. Langford, L. Li, and T. Zhang, “Sparse online learning via truncated gradient,” J. Mach. Learn. Res., vol. 10, pp. 777–801, 2009.
[64] Y. Nesterov, “Smooth minimization of non-smooth functions,” Math. Program., vol. 103, pp. 127–152, 2005.
[65] Y. Xu, Y. Yan, Q. Lin, and T. Yang, “Homotopy smoothing for non-smooth problems with lower complexity than O(1/),” in Proc. Adv. Neural
Inf. Process. Syst., 2016, pp. 1208–1216.
[66] P. Zhao and T. Zhang, “Stochastic optimization with importance sampling for regularized loss minimization,” in Proc. 32nd Int. Conf. Mach.
Learn., 2015, pp. 1–9.
[67] D. Needell, N. Srebro, and R. Ward, “Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm,” Math.
Program., vol. 155, pp. 549–573, 2016.
[68] J. Konečný and P. Richtárik, “Semi-stochastic gradient descent methods,” ArXiv Preprint: 1312.1666v2, 2015.
[69] G. Lan, “An optimal method for stochastic composite optimization,” Math. Program., vol. 133, pp. 365–397, 2012.

31
−2

−6

10

−8

10

10

−6

10

−8

10

−10

−12

10

5

10
15
20
25
30
Gradient evaluations / n

10

35

−4

−8

10

10

−12

0

−6

10

5

10
15
20
25
Gradient evaluations / n

30

10

35

10

−6

10

−8

10

−10

10

−12

0

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

10

20
30
Gradient evaluations / n

10

40

0

10

20
30
40
Gradient evaluations / n

50

−2

−6

10

−8

10

10

−6

10

−8

10

−10

10

0

0.5

1
1.5
Running time (sec)

10

2

−4

−8

10

10

−12

−12

−6

10

(a) Adult: λ1 = 10−4

0.5

1
1.5
2
Running time (sec)

2.5

10

3

10

−6

10

−8

10

−10

10

−12

0

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

(b) Protein: λ1 = 10−4

5

10
15
Running time (sec)

10

20

0

(c) Covtype: λ1 = 10−4

10

20
30
Running time (sec)

40

(d) Sido0: λ1 = 5∗10−3

−2

−6

10

−8

10

−8

10

−10

10

−12

10

−6

10

20

40
60
Gradient evaluations / n

10

80

−4

−8

10

10

−12

0

−6

10

5

10
15
20
25
Gradient evaluations / n

30

10

35

10

−6

10

−8

10

−10

10

−12

−12

0

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

10
F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

0

5

10
15
20
25
Gradient evaluations / n

30

10

35

0

10

20
30
40
Gradient evaluations / n

50

60

−2

−6

10

−8

10

−8

10

−10

10

−12

10

−6

10

0.5

1
1.5
2
Running time (sec)

2.5

10

3

−4

−8

10

10

−12

0

−6

10

(e) Adult: λ1 = 10−5

1

2

3
4
Running time (sec)

5

10

6

10

−6

10

−8

10

−10

10

−12

0

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

10

−12

0

(f) Protein: λ1 = 10−5

5

10
15
Running time (sec)

10

20

0

(g) Covtype: λ1 = 10−5

5

10
15
Running time (sec)

20

25

(h) Sido0: λ1 = 10−4

−2

−6

10

−8

10

−8

10

−10

10

−12

10

−6

10

50
100
Gradient evaluations / n

10

150

−6

10

−4

−8

10

10

0

50
100
Gradient evaluations / n

10

150

10

−6

10

−8

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−10

10

−12

−12

0

−2

10

−10

10

−10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

10

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

10

20
30
40
50
Gradient evaluations / n

10

60

0

50
100
Gradient evaluations / n

150

−2

−6

10

−8

10

−8

10

−10

10

−12

10

−6

10

1

2

3
4
5
Running time (sec)

(i) Adult: λ1 = 10−6

6

7

10

−8

10

10

0

5
10
Running time (sec)

(j) Protein: λ1 = 10−6

15

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−4

10

−6

10

−8

10

−10

10

−12

−12

0

−6

10

−10

10

−10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

5

10
15
20
Running time (sec)

(k) Covtype: λ1 = 10−6

25

10

0

20

40
60
Running time (sec)

80

100

(l) Sido0: λ1 = 10−5

Fig. 5. Comparison of SVRG [1], Prox-SVRG [2], Katyusha [3], and VR-SGD for solving `2 -norm regularized logistic regression problems (i.e.,
λ2 = 0). In each plot, the vertical axis shows the objective value minus the minimum, and the horizontal axis is the number of effective passes (top)

or running time (bottom).

32
−2

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−4

10

−4

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−8

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−8

10

−10

10

−10

10

10
F (xs ) − F (x* )

−6

10

−6

10

F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

−6

10

−8

10

−4

10

−10

10

−6

10
−12

10

−12

0

10

20
30
40
Gradient evaluations / n

50

10

60

−12

0

20

40
60
Gradient evaluations / n

10

80

0

10

20
30
Gradient evaluations / n

40

50

0

50

100
150
Gradient evaluations / n

200

−2

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−4

10

−4

10

−2

10

−6

10

−8

10

−6

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−8

10

−6

10

−8

10

−10

−12

10

10

−12

0

0.5

1
Running time (sec)

1.5

10

2

SVRG
Prox−SVRG
Katyusha
VR−SGD

−10

10

−10

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

−12

0

2

4
6
Running time (sec)

8

10

10

0

5

10
Running time (sec)

−4

10

−6

10

15

0

20

50
100
Running time (sec)

(a) λ2 = 10−4

150

(b) λ2 = 10−3
−2

−6

10

−8

10

−8

10

−6

10

−8

10

−12

0

50
100
Gradient evaluations / n

10

150

−4

10

−12

10

−6

10

−8

10

−10

10

−12

0

50

100
150
Gradient evaluations / n

10

200

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

10

10

−6

10

−10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

10
F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

20

40
60
Gradient evaluations / n

10

80

0

50

100
150
200
Gradient evaluations / n

250

300

−2

10

−8

10

−6

10

−8

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−10

10

−10

10

−12

−12

0

2

4
Running time (sec)

6

10

8

0

5

10
15
20
Running time (sec)

25

−4

−6

10

−8

10

10

−6

10

−8

10

−10

10

−10

10

−12

10

30

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

F (xs ) − F (x* )

−6

10

−4

10
F (xs ) − F (x* )

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

−4

F (xs ) − F (x* )

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

10

20
30
Running time (sec)

40

10

50

0

50

(c) λ2 = 10−5

−4

−6

10

−8

10

10

−2

−5

10

−6

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

−10

10

−8

50

100
150
200
Gradient evaluations / n

10

300

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−4

10

0

50

100
150
200
250
Gradient evaluations / n

10

−6

10

−8

10

300

350

0

−6

10

−8

10

−10

5

10
Running time (sec)

15

20

10

350

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

0

300

0

0

10

20
30
Running time (sec)

(e) λ2 = 10−6

40

50

100
150
200
250
Gradient evaluations / n

50

0

50
100
Running time (sec)

300

350

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

−10

10

100
150
200
250
Gradient evaluations / n

10

10

−12

50

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

F (xs ) − F (x* )

250

−6

F (xs ) − F (x* )

0

−4

10

10

F (xs ) − F (x* )

−12

10

10

−7

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

10
F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−3

10

200

(d) λ2 = 10−4

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

100
150
Running time (sec)

−4

10

−6

10

150

0

50

100
150
Running time (sec)

200

(f) λ2 = 10−5

Fig. 6. Comparison of SVRG, Prox-SVRG [2], Katyusha [3], and VR-SGD for `1 -norm regularized logistic regression problems (i.e., λ1 = 0) on the
four data sets: Adult (the first column), Protein (the sconced column), Covtype (the third column), and Sido0 (the last column). In each plot, the
vertical axis shows the objective value minus the minimum, and the horizontal axis is the number of effective passes (top) or running time (bottom).

33
−2

−6

s

10

−8

10

−6

10

−8

10

−8

10

10

10

−12

−12

0

10

20
30
Gradient evaluations / n

40

10

50

−4

−6

10

−6

10

−8

10

−10

10

−10

−10

10

−12

0

10

20
30
Gradient evaluations / n

10

40

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

10

−10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

F (xs ) − F (x* )

*

F (x ) − F (x )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

5

10
15
Gradient evaluations / n

10

20

0

10

20
30
40
50
Gradient evaluations / n

60

70

−2

−6

10

−8

10

−6

10

−8

10

−10

−12

10

0

0.5

1

1.5
2
2.5
Running time (sec)

3

10

3.5

−4

−6

10

−8

10

10

−12

1

2
3
Running time (sec)

4

10

5

10

−6

10

−8

10

−10

10

−12

−12

0

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

0

1

2

3
4
5
Running time (sec)

6

10

7

(a) λ1 = 10−5 and λ2 = 10−4

0

10

20
30
40
Running time (sec)

50

60

(b) λ1 = 10−4 and λ2 = 10−5
−2

−6

10

−8

10

−8

10

−10

10

−12

10

10

−6

10

20

40
60
80
Gradient evaluations / n

100

10

120

−4

−6

10

−8

10

10

−12

0

10

−6

10

−8

10

−10

10

−12

−12

0

10

20
30
40
Gradient evaluations / n

10

50

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

10
F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

0

5

10
15
Gradient evaluations / n

10

20

0

50

100
150
Gradient evaluations / n

200

−2

−6

10

−8

10

−8

10

−10

10

−12

10

10

−6

10

1

2
3
Running time (sec)

10

4

−4

−8

10

10

−12

0

−6

10

1

2
3
Running time (sec)

4

10

5

10

−6

10

−8

10

−10

10

−12

0

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

2

4
6
Running time (sec)

8

10

10

(c) λ1 = 10−5 and λ2 = 10−5

0

20

40
60
Running time (sec)

80

100

(d) λ1 = 10−5 and λ2 = 10−4
−2

−6

10

−8

10

−8

10

−10

10

−12

50
100
Gradient evaluations / n

10

150

−4

−8

10

10

10

−6

10

−8

10

−10

10

−12

−12

0

−6

10

0

50
100
Gradient evaluations / n

10

150

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−10

10

−10

10

−6

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

20

40
60
Gradient evaluations / n

10

80

0

50

100
150
Gradient evaluations / n

200

−2

−6

10

−8

10

−8

10

−10

10

−12

2

4
6
Running time (sec)

8

10

10

−8

10

10

0

5

10
15
Running time (sec)

20

(e) λ1 = 10−6 and λ2 = 10−5

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−4

10

−6

10

−8

10

−10

10

−12

−12

0

−6

10

−10

10

−10

10

−6

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10
F (xs ) − F (x* )

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

10

SVRG
Prox−SVRG
Katyusha
VR−SGD

−4

10

F (xs ) − F (x* )

SVRG
Prox−SVRG
Katyusha
VR−SGD

−2

10

−12

0

5

10

15
20
25
Running time (sec)

30

35

10

0

50
100
Running time (sec)

150

(f) λ1 = 10−5 and λ2 = 10−5

Fig. 7. Comparison of SVRG, Prox-SVRG [2], Katyusha [3], and VR-SGD for solving elastic net regularized logistic regression problems on the four
data sets: Adult (the first column), Protein (the sconced column), Covtype (the third column), and Sido0 (the last column). In each plot, the vertical
axis shows the objective value minus the minimum, and the horizontal axis is the number of effective passes (top) or running time (bottom).

34

−4

F (xs ) − F (x* )

−6

10

−8

10

−10

−8

10

10

20
30
Gradient evaluations / n

40

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−4

10

10

−8

10

−10

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

10

0.5
1
Running time (sec)

−8

10

0.5
1
Running time (sec)

10

1.5

−10

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

10

−12

−4

−6

10

−8

10

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

50
100
Gradient evaluations / n

−10

10

−6

10

−8

10

0.5

1
Running time (sec)

1.5

2

(d) λ = 10−6

10

10

−6

10

−8

10

−10

10

−12

0

150

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

−4

10

−12

50
100
Gradient evaluations / n

−2

−10

10

0

10

F (xs ) − F (x* )

−8

10

10

150

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

−4

−6

−8

10

−12

0

−2

10

−6

10

−10

10

F (xs ) − F (x* )

−4

10

10

100

10

10

−12

40
60
80
Gradient evaluations / n

1.5

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

−2

10

−2

0.5
1
Running time (sec)

10

−10

10

10

0

(c) λ = 10−5

F (xs ) − F (x* )

−8

10

20

−8

10

−12

0

−4

−6

0

−6

10

−10

−2

10

10

10

10

10

F (xs ) − F (x* )

−4

10

50

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

(b) λ = 10−4
SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−2

10

20
30
40
Gradient evaluations / n

−4

−6

10

1.5

10

−2

−12

0

0

10

10

(a) λ = 10−3

F (xs ) − F (x* )

10

50

10

−12

F (xs ) − F (x* )

40

−10

10

10

20
30
Gradient evaluations / n

−4

−6

−8

10

−12

0

−2

10

−6

10

−10

10

F (xs ) − F (x* )

−2

10

50

10

10

−12

0

10

F (xs ) − F (x* )

−6

10

10

−12

10

−4

−10

10
10

10

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−2

10

F (xs ) − F (x* )

F (xs ) − F (x* )

−4

10

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−2

10

F (xs ) − F (x* )

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−2

10

−12

0

1

2

3
4
5
Running time (sec)

(e) λ = 10−7

6

7

10

0

1

2

3
4
5
Running time (sec)

6

7

(f) λ = 0

Fig. 8. Comparison of SVRG [1], Katyusha [3], VR-SGD and their proximal versions for solving ridge regression problems with different regularization
parameters on the Adult data set. In each plot, the vertical axis shows the objective value minus the minimum, and the horizontal axis is the number
of effective passes (top) or running time (bottom).

35

−2

−4

10

−6

10

−8

10

−10

5

10
15
20
Gradient evaluations / n

25

10

30

−12

0

10

−6

10

−8

10

−10

4
Running time (sec)

6

−6

10

−8

10

10

8

−8

10

−10

−12

0

2

4
6
8
Running time (sec)

10

10

12

0

2

4
6
8
Running time (sec)

(b) λ = 10−4

10

12

(c) λ = 10−5

−2

−4

10

−6

10

−8

10

10

−4

10

−10

−1

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−6

10

10

F (xs ) − F (x* )

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

F (xs ) − F (x* )

10

−8

10

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

−10

10

10

−12

−2

10

−12

0

10

20
30
40
50
Gradient evaluations / n

10

60

−2

0

50
100
Gradient evaluations / n

150

0

100

200
300
400
Gradient evaluations / n

500

−2

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−4

10

−6

10

−8

10

−10

10

−1

10

−4

10
F (xs ) − F (x* )

10

−6

10

−8

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

10

−10

10

10

−12

−12

0

5
10
Running time (sec)

15

F (xs ) − F (x* )

10

−6

10

10

(a) λ = 10−3
−2

40

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

10

−12

2

20
30
Gradient evaluations / n

−4

−10

0

10

10

10

−12

0

−2

10

10

F (xs ) − F (x* )

10

40

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−4

F (xs ) − F (x* )

10

F (xs ) − F (x* )

20
30
Gradient evaluations / n

−2

−4

−8

10

−10

10

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

−6

10

10

F (xs ) − F (x* )

−2

F (xs ) − F (x* )

−8

10

−12

0

10

10

10

10

−12

10

−6

10

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−4

−10

10

10

10

SVRG−I
SVRG−II
Katyusha−I
Katyusha−II
VR−SGD−I
VR−SGD−II

−4

F (xs ) − F (x* )

F (xs ) − F (x* )

10

−2

10

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

F (xs ) − F (x* )

−2

10

10

0

(d) λ = 10−6

10

20

30
40
50
Running time (sec)

(e) λ = 10−7

60

SVRG−I
SVRG−II
VR−SGD−I
VR−SGD−II

−2

10
70

0

20

40
60
80
Running time (sec)

100

120

(f) λ = 0

Fig. 9. Comparison of SVRG [1], Katyusha [3], VR-SGD and their proximal versions for solving ridge regression problems with different regularization
parameters on the Covtype data set. In each plot, the vertical axis shows the objective value minus the minimum, and the horizontal axis is the
number of effective passes (top) or running time (bottom).

36

Fixed learning rate
Varying learning rate
Fixed learning rate
Varying learning rate

F (xs ) − F (x* )

10

−4

10

*

F (x ) − F (x )

−4

−6

s

10

−8

10

−6

10

−8

10

−10

−10

10

10

−12

10

Fixed learning rate
Varying learning rate
Fixed learning rate
Varying learning rate

−12

0

10
20
30
Gradient evaluations / n

10

40

0

5

10
15
Gradient evaluations / n

20

(a) `1 -norm regularized logistic regression: Adult (left) and Covtype (right)

Fixed learning rate
Varying learning rate
Fixed learning rate
Varying learning rate

F (xs ) − F (x* )

10

−4

10

*

F (x ) − F (x )

−4

−6

s

10

−8

10

−6

10

−8

10

−10

−10

10

10

−12

10

Fixed learning rate
Varying learning rate
Fixed learning rate
Varying learning rate

−12

0

5

10
15
20
Gradient evaluations / n

25

30

10

0

5
10
Gradient evaluations / n

15

(b) Lasso: Adult (left) and Covtype (right)
Fig. 10. Comparison of Algorithm 3 with fixed and varying learning rates for solving `1 -norm (i.e., λkxk1 ) regularized logistic regression and Lasso
problems with λ = 10−4 (blue lines) and λ = 10−5 (red lines, best viewed in colors). Note that the regularization parameter is set to 10−3 and
10−4 for solving Lasso problems on the Covtype data set.

