Parallelizing Stochastic Approximation Through
Mini-Batching and Tail-Averaging
Prateek Jain∗, Sham M. Kakade†, Rahul Kidambi‡, Praneeth Netrapalli§, Aaron Sidford¶

arXiv:1610.03774v2 [stat.ML] 27 Oct 2016

Abstract
This work characterizes the benefits of averaging techniques widely used in conjunction with stochastic
gradient descent (SGD). In particular, this work sharply analyzes: (1) mini-batching, a method of averaging
many samples of the gradient to both reduce the variance of a stochastic gradient estimate and for parallelizing
SGD and (2) tail-averaging, a method involving averaging the final few iterates of SGD in order to decrease
the variance in SGD’s final iterate. This work presents the first tight non-asymptotic generalization error
bounds for these schemes for the stochastic approximation problem of least squares regression.
Furthermore, this work establishes a precise problem-dependent extent to which mini-batching can be
used to yield provable near-linear parallelization speedups over SGD with batch size one. These results are
utilized in providing a highly parallelizable SGD algorithm that obtains the optimal statistical error rate
with nearly the same number of serial updates as batch gradient descent, which improves significantly over
existing SGD-style methods.
Finally, this work sheds light on some fundamental differences in SGD’s behavior when dealing with
agnostic noise in the (non-realizable) least squares regression problem. In particular, the work shows that
the stepsizes that ensure optimal statistical error rates for the agnostic case must be a function of the noise
properties.
The central analysis tools used by this paper are obtained through generalizing the operator view of
averaged SGD, introduced by Défossez and Bach [1] followed by developing a novel analysis in bounding
these operators to characterize the generalization error. These techniques may be of broader interest in
analyzing various computational aspects of stochastic approximation.

1

Introduction and problem setup

With the ever increasing size of modern day datasets, practical algorithms for machine learning are increasingly
constrained to spend less time and use less memory. This makes it particularly desirable for simple streaming
algorithms that generalize well in just a few passes over the dataset. Stochastic gradient descent (SGD) is
perhaps the simplest and most well studied algorithm that meets these constraints. The algorithm simply
repeatedly samples an instance from the stream of data and updates the current parameter estimate using
the gradient of the sampled instance. Nevertheless, SGD has been impressively successful and is the de-facto
method of choice for large scale learning problems.
While a powerful machine learning tool, unfortunately SGD in its simplest forms is inherently serial.
Over the past years as dataset sizes have grown there have been remarkable developments in processing
capabilities, with multi-core/distributed/GPU based computing infrastructure available in relative abundance.
The presence of this computing power has triggered the development of parallel/distributed machine learning
algorithms [2, 3, 4, 5, 6, 7] that possess the capability to utilize multiple cores/machines. However, despite this
exciting work, it is yet unclear how to best parallelize SGD and fully utilize these computing infrastructures.
This paper takes a step towards answering this question, by characterizing the behavior of constant
stepsize SGD for the problem of strongly convex stochastic least square regression (LSR) under two averaging
techniques widely believed to improve the performance of SGD. In particular, this work considers the natural
∗ Microsoft

Research, India, e-mail: prajain@microsoft.com
of Washington, Seattle, e-mail: sham@cs.washington.edu
‡ University of Washington, Seattle, e-mail: rkidambi@uw.edu
§ Microsoft Research, India, e-mail: praneeth@microsoft.com
¶ Stanford University, Palo Alto CA, e-mail: sidford@stanford.edu
† University

1

parallelization technique of mini-batching, where multiple points are processed simultaneously and the current
iterate is updated by the average gradient over these points, and combine it with variance reducing technique
of tail-averaging, where the average of many of the final iterates are returned as SGD’s estimate of the solution.
In this work, the arguments about parallelization are structured through the lens of a work-depth tradeoff:
work refers to the total amount of computation required to reach a certain generalization error, and depth
refers to the amount of serial computation. Depth, defined in this manner, is equal to the runtime of the
algorithm on a large multi-core architecture with shared memory, where there is no communication overhead,
and has strong implications for parallelizability on other architectures as well.
The main contributions of this paper are as follows:
• This work shows that mini-batching yields near-linear parallelization speedups over the standard serial
SGD (i.e. with batch size 1), as long as the mini-batch size is smaller than a problem dependent quantity
(denoted by bthresh ). When batch-sizes increase beyond bthresh , mini-batching is inefficient (owing to the
lack of serial updates), thus obtaining only sub-linear speedups over and above mini-batching with a
batch size bthresh . A by-product of this analysis sheds light on how the step sizes naturally interpolate
from ones used by standard serial SGD (with batch size 1) to ones used by batch gradient descent.
• While the final iterate of SGD decays initial error at a geometric rate but does not obtain optimal
statistical error rates, the average iterate [8, 1] decays the initial error at a sublinear rate while
preserving statistical optimality. This work rigorously shows that tail-averaging obtains the best of both
worlds – decaying initial error at a geometric rate and obtaining near-optimal statistical error rates.
• Combining the above results, this paper provides a mini-batching and tail-averaging version of SGD that
is highly parallelizable – the number of serial steps (which is a proxy for the un-parallelizable time)
of this algorithm nearly matches that of offline gradient descent and is lower than the serial time of all
existing streaming LSR algorithms. See Table 1 for comparison. We note that these results are obtained
by providing a tight finite-sample analysis of the effects of mini-batching and tail-averaging with large
constant learning rate schemes.
• All the results in this paper are established for the general agnostic noise case of the streaming LSR
problem. This reveals a fundamental difference in the behavior of SGD when dealing with the agnostic
case as opposed to the realizable case. In particular, this analysis reveals a surprising insight that the
maximal stepsizes that are sufficient to obtain statistically optimal rates are a function of the properties
of the agnostic noise. The main result of this analysis is that the maximal step sizes for the agnostic case
could be much lower than ones that are used in the realizable case: indeed, an instance that yields such
a separation between the maximal learning rates for the realizable and agnostic case is provided.
The tool employed in obtaining these results is Algorithm 1 and its generalization guarantees as in
Theorem 1, which builds upon and generalizes the analysis of [1] and provides the first tight bound on the
finite-sample generalization error for mini-batching, tail-averaged constant step-size SGD. We note that the
work of [1] does not establish statistically optimal rates while working with large constant step size schemes;
this shortcoming is remedied through coming up with a novel sharp analysis that is used to rigorously establish
statistically optimal rates while working with large constant step sizes (this can be viewed as a special case of
Theorem 2 by setting batch size b = 1 and through averaging iterates from the start by setting s = 0).
Applying Theorem 1 to the general agnostic case of LSR, we obtain Theorems 2 and 3, which shows
an interesting trade-off between the convergence rate or generalization error (work) of SGD, the amount
of parallelization (depth of computation) achieved by increasing the mini-batch size, and the amount of
tail-averaging performed. As with previous SGD-style algorithms, [8, 9, 10, 11, 1, 12], this trade-off stems
from the two terms of the generalization error of SGD: (i) the bias term which represents how quickly the
initial conditions are forgotten, and (ii) the variance term which represents the dependence of the excess
risk on the level of noise (due to streaming data/random sampling/model mismatch) that is present in the
problem. We note that existing results [9, 1] dealing with averaged SGD with batch size 1 applied to the
realizable case of LSR generalize to this paper’s results which establishes (a) the finite-sample effects of
mini-batching, tail averaging in the context of stochastic approximation, and (b) the impact of agnostic noise
on the learning rates employed by SGD to guarantee statistically optimal error rates.
The results of this paper are obtained through building on and generalizing the operator view of averaged
SGD with batch size one, as introduced by [1] to the mini-batched tail-averaged case and through subsequently
2

Algorithm
Batch
gradient
descent
Stochastic
Dual
Co-ordinate
Ascent [13]
Averaged
SGD [1, 9]
Streaming
SVRG (with
initial error
oracle) [12]

Final error
 2 
O σnd

O

O





R4
λ2min n2

σ2 d
n

Runtime/Work
κnd log

R2 t
kHk2 n

κ log

n·init
σ2 d


2

(n+ λRmin d)d·log
2

· init + σnd



t
 κ log(dκ)

σ2 d
n

σ2 d
n

Algorithm 2

O

Algorithm 2
with initial
error oracle




min (H)
· init +
O exp − Rnλ
2 ·log(dκ)

· init +

n·init
σ2 d

2

(n+ λRmin d)·log

nd





(H)
· init +
O exp − nλmin
R2


n·init
σ2 d

Depth

nd

n
2

( λminR(H) ) · log



σ2 d
n

nd


nd

n·init
σ2 d

t log



n·init
σ2 d

nkHk2
tR2

κ log(dκ) log



n·init
σ2 d

Streaming

Agnostic

×

X

×

X

X

×

X

X

X

X

X

X

Table 1: Comparison of Algorithm 2 with existing algorithms given n samples for LSR, with init = R(w0 ) − R(w∗ ). Note that
the bounds have been provided assuming realizable noise; for the algorithms with support for agnostic noise, these can be
appropriately modified. Refer to Section 3 for the definitions of all quantities. We do not consider accelerated variants in this
table. Note that the accelerated variants serve to improve running times of the offline algorithms. Furthermore, offline algorithms
possess the advantage that they do not require the knowledge of R2 . The quantity t in the results for Algorithm 2 can be chosen
arbitrarily and trades off the rate of decay in initial error with the depth of the algorithm. Finally, we note that streaming
SVRG [12] does not conform to the standard first order oracle model [14].
providing a new sharp analysis in order to bound the resulting operators. In particular, straightforward
operator norm bounds of these matrix operators suffices to show convergence of SGD updates, but turn out to
be pretty loose bounds. In order to obtain stronger bounds, this paper employs a much more fine grained
analysis that directly bounds the trace of these operators applied to the relevant matrices, which results in
obtaining much tighter bounds on the performance of mini-batched tail-averaged SGD while simultaneously
shedding light on the role of agnostic noise on the behavior of SGD.
While this paper’s results focus on the case of strongly convex streaming least square regression, we believe
that our techniques extend much more broadly. We hope that this paper will serve as the basis for future work
on analyzing SGD and improving the parallelization of large scale algorithms for machine learning.
Paper organization: A brief overview of related works is presented in Section 2. Then, the problem setup
and notations are introduced in Section 3. The main results of this paper as well as the mini-batched tail
averaging SGD algorithm are presented in Section 5. In Section 6, simulation results are presented using
a synthetic dataset, and these results corroborate well with the main results on the mini-batch thresholds
established by the paper as well as the advantages of mini-batch tail averaged SGD algorithm (i.e. algorithm 1).
The proofs of all the claims and theorems are provided in the appendix.

2

Related work

We split the work related to this paper’s contributions into several algorithm classes:
Stochastic Approximation: Seminal works of Ruppert [15], and Polyak and Juditsky [8] showed that
averaging iterates of SGD applied to stochastic optimization problems like stochastic least squares regression
achieves asymptotically optimal statistical rates. A more recent line of work [16, 9, 10, 11, 17, 18, 1] has
provided rigorous non-asymptotic analysis (with explicit constants) of SGD for stochastic approximation
problems including linear and logistic regression. Equipped with averaging, these algorithms achieve (up to
constants) minimax optimal statistical rates. Note that all these results assume that the noise is realizable when
stating the generalization error bounds. We are unaware of any results that provide sharp non-asymptotic
analysis of SGD and the related step size issues in the general agnostic noise setting. Several specialized
3

stochastic algorithms such as SVRG [19], SDCA [13] and SAGA [20] have been proposed for the agnostic
setting, among which, only the streaming version of SVRG [12] is known to provide minimax optimal
statistical rates and is highly parallelizable. Note that offline SVRG [19] is parallelizable as well, but is not
streaming and incurs additional log factors without carefully employing doubling tricks. Note that streaming
SVRG [12] does not take advantage of minibatching and its depth depends on a stronger fourth moment
condition on the input whereas our mini-batching bounds depend only on the condition number of the
covariance matrix. Moreover, our constant-step size SGD algorithm is arguably simpler than streaming SVRG
and can be run in a more general setting requiring less information from a single sample. In particular, SGD
based algorithms conform to the standard first order oracle model [14], whereas the streaming SVRG [12]
does not, since it requires gradient information at two different points from a single sample.
Another key aspect of this paper’s main result is a precise characterization of the effects of mini-batching,
which is also studied extensively in this literature [21, 22, 23, 6, 24]. Previous results on mini-batching within
the context of stochastic approximation yielded results that suggest increasing minibatch sizes tend to hurt
the convergence rate [22, 21]. [6, 24, 21, 23] study mini-batching for various algorithms and show faster
decay in the bias term but do not show optimal statistical error rates. Moreover, all of these results claim
speed-up from mini-batching by comparing the upper bounds on the training/generalization error. In contrast,
this paper’s result for mini-batching compares the exact generalization error for different batch sizes and
demonstrates an improvement with larger sizes.
The other scheme for parallelizing stochastic approximation is through model averaging [2, 3, 25, 26].
These schemes are communication efficient since they average the final output of independent runs of the
algorithm, while being minimax optimal. However, as noted in [26], these schemes do not yield improvements
to the bias and it is unclear how these results relate to this work’s depth bounds. Other efforts [25, 2, 3]
provide guarantees on upperbounds of the error, with only [25] guaranteeing minimax optimality.
Non-Stochastic Approximation: In order to utilize the recent developments in computing frameworks,
there has been a number of efforts in developing distributed/parallel machine learning algorithms that utilize
multi-core [4, 5, 21, 22, 27, 24] and multi-machine [2, 3, 6, 23, 7] computing frameworks. As in the stochastic
case, parallelization speedups are argued only in terms of upper bounds of training/generalization error. On
the other hand, a distinctive feature of this paper’s mini-batching result is that the prescribed parallelization
speedups are on the exact generalization error on a per problem basis, and do not rely on providing linear
speedups on weaker upper bounds. Moreover, none of the existing parallelization results characterize regimes
where there are unconditional linear speedups on the overall generalization error while maintaining minimax
optimal statistical error rates. Indeed, model averaging results are known to be minimax optimal [25, 26],
however, the bias term is unimprovable using these schemes, as noted by [26].
Asynchronous Stochastic Approximation: There exists a recent line of work [5, 28, 29] that provides
linear speedups of SGD in the context of asynchronous stochastic optimization. However, in contrast to
our work, these efforts deal with the much harder case of asynchronous optimization, and moreover, these
results provide speedups on upper bounds on the training/generalization error while relying on specific input
structures such as hard sparsity. Moreover, while [29] proves statistical optimality in an asymptotic sense,
none of the other efforts claim statistical optimality.

3

Preliminaries

Problem formulation: Let D be a distribution on tuples (x, y) ∈ Rd × R. Let R : Rd → R be defined as the
expected square loss over the pairs (x, y):
R(w) =

1
2
E(x,y)∼D (y − hw, xi) ∀ w ∈ Rd .
2

(1)

Our goal is to compute the predictor w∗ = arg minw R(w). In this paper, we provide the maiden sharpest
analysis of SGD with batch sizes ≥ 1, applied to the general agnostic case of LSR (also referred to as the
def
unrealizable case) where y = hw∗ , xi +  where  ∈ R is the noise.
Notation: We use boldface small letters (x, w etc.) for vectors and boldface capital letters (A, H etc.)
for matrices. We use ⊗ to denote the outer product of two vectors or matrices. S(d) denotes the space of
real symmetric d × d matrices and M(d) denotes the space of linear transformations from Rd×d to Rd×d .
Every element in M(d) is a d × d × d × d tensor, which is represented using normal font scripted letters
4

Algorithm 1 Minibatch-TailAveraging-SGD
Input: Initial point w0 , stepsize γ, minibatch size b, initial iterations s, total number of samples m.
1: for t = 1, 2, .., s do
2:
Sample “b” tuples {(xti , yti )}bi=1 ∼ Db
Pb
3:
wt ← wt−1 − γb i=1 (hwt−1 , xti i − yti )xti
4: w̄ = 0
5: for t = s + 1, s + 2, .., m
b do
6:
Sample “b” tuples {(xti , yti )}bi=1 ∼ Db
Pb
7:
wt ← wt−1 − γb i=1 (hwt−1 , xti i − yti )xti
1
8:
w̄ = t−s
wt + t−s−1
t−s w̄
Output: w̄
(M, T etc.).PA tensor M applied to a matrix A gives a matrix MA such that its ij th element is given by
(MA)ij = k,`∈[d] Mik`j Ak` . Any d × d matrix A can be lifted to AL , AR ∈ M(d), which represents the
left/right multiplication operator of the matrix A respectively. More concretely, AL , AR ∈ M(d) are defined
so that AL W = AW and AR W = WA ∀ W ∈ Rd×d . We let I ∈ M(d) represent the identity map.
We let k·k or k·k2 denote `2 norm for vectors and spectral or operator norm for matrices and tensors. More
concretely, we define
def

kAk = max
x∈Rd

kT AkF
kAxk
def
∀ T ∈ M(d),
∀ A ∈ Rd×d , and kT k = max
kxk
A∈S(d) kAkF

where k·kF denotes the Frobenius norm of a matrix. When we refer to the largest/smallest eigenvalues
λmax (S) or λmin (S) of a tensor S, this is obtained by viewing the operator S as a member of M(d).
Relevant quantities: We denote by H and M the population second moment matrix and fourth moment
tensor of the covariates x, i.e.:
def

def

H = E(x,y)∼D [x ⊗ x], M = E(x,y)∼D [x ⊗ x ⊗ x ⊗ x]
We assume
that both H and M are finite and are related by MI  R2 H. This in particular implies that

2
E kxk ≤ R2 . We further assume that H is non-singular, with λmin (H) > 0 denoting the smallest eigenvalue
of H. This makes R(w) strongly convex and implies the existence of a unique minimizer w∗ of the objective
R(w). The condition number κ of the population covariance matrix H is defined to be the ratio of the largest
def λmax (H)
λmin (H) .

to smallest eigenvalues of H, i.e., κ =

Next, we define the noise covariance matrix Σ as:



def
Σ = E(x,y)∼D (y − hw∗ , xi)2 · x ⊗ x
In the realizable case, the noise  is independent of the covariate x, thus implying Σ = σ 2 H.
 Note that M

can also be thought of as an operator in M(d) which acts on matrix A to produce MA = E (x> Ax)x ⊗ x .
For a given γ > 0 and b ∈ Z+ , a key quantity in our results turns out to be the operator Tbγ defined as
def

Tbγ = HL + HR − γb M −
Tbγ .

γ(b−1)
HL HR .
b

In most cases we drop γ for notational convenience and use Tb for

def

We let T = HL + HR − γHL HR denote the limit of Tb as b → ∞. For a given γ > 0, we also use G to
def

denote the operator G = HL −1 + HR −1 − γI.
Assumptions: For the realizable case, the generalization error estimate in Theorem 1 is exact on a
per-problem basis if E [yi |xi ] = hw∗ , xi i (which is true if we have a constant component
inx, as is typically

the case in practice). For the agnostic case, we note that this result is exact if E x(i) x(j) x(k)  = 0 ∀ (i, j, k) ∈
⊗3 {1, 2, .., d} with x(i) referring to the ith dimension of x. If this is not the case, the generalization error in
Theorem 1 are off the exact value by a factor of at most 2 (as noted by [9, 1], owing to the application of
Minkowski’s inequality).

4

Viewing Stochastic Gradient Descent as a Stochastic Process

In this section, we consider a stochastic process view of Stochastic Gradient Descent updates and shed some
light on why this view allows us to provide a sharp characterization of SGD as a tool to solve the general LSR
5

Algorithm 2 MinibatchDoublingPartialAveragingSGD
Input: Initial point w0 , stepsize γ, initial minibatch size b, number of iterations in each epoch s, number of
samples n.
1: /*Run logarithmic number of epochs where each epoch runs t iterations of minibatch SGD (with out
averaging). Double minibatch size after each epoch.*/
n
2: for ` = 1, 2, · · · , log bt
− 1 do
`−1
3:
b` ← 2 b
4:
w` ← Minibatch-TailAveraging-SGD(w`−1 , γ, b` , t − 1, t · b` )
5: /*For the last epoch, run tail averaged minibatch SGD with initial point wt , stepsize γ, minibatch size
n
2log bt −1 · b = n/2t, number of initial iterations t/2 and number of samples n/2.*/
6: w ← Minibatch-TailAveraging-SGD(ws , γ, n/2t, t/2, n/2)
Output: w
problem. In particular, let us consider the mini-batch SGD update that moves iterate wt−1 to wt :
wt = wt−1 −

b
γX
(hwt−1 , xti i − yti )xti
b i=1

def

By defining ηt = wt − w∗ , we write the recursion relating ηt to ηt−1 as:

γX
γX
xti ⊗ xti ηt−1 +
ti xti
ηt = I −
b i
b i
Now, since we deal with the square loss case, the generalization error is written as:
E [R(wt )] − R(w∗ ) =

1
Tr H E [ηt ⊗ ηt ]
2

We thus consider E [ηt ⊗ ηt ], which can be written in terms of E [ηt−1 ⊗ ηt−1 ] as:
E [ηt ⊗ ηt ] = (I − γTb ) E [ηt−1 ⊗ ηt−1 ] + γ 2 Σ
Divergent Learning Rates: This operator view of mini-batch SGD immediately sheds light on the maximal
div
stepsizes γb,max
that prevents divergence of the resulting mini-batch SGD updates: in particular, a necessary
condition that the set of non-divergent learning rates γ ≥ 0 should satisfy is Tbγ  0. This condition yields the
div
following generalized eigenvector problem for computing γb,max
:
2
div
γb,max

=

hW, MWi + (b − 1) · Tr WHWH
b · Tr WHW
W∈S(d)
sup

This characterization generalizes the divergent stepsize characterization of [1] for batch sizes > 1. We note
that such a characterization automatically sheds light on how the maximal divergent learning rates interpolate
batch sizes of 1 (which is ≤ Tr2H [1]) to the batch gradient descent learning rate (setting b to ∞), which
div
turns out to be λmax2(H) . A property of γb,max
that is worth noting is that it does not depend on properties of
the noise covariance Σ, and depends solely on the statistics of the covariate x.
div
We note that in this paper, our interest does not lie in the non-divergent stepsizes 0 ≤ γ ≤ γb,max
, but in
div
the set of (maximal) step sizes 0 ≤ γ ≤ γb,max (< γb,max ) that are sufficient to guarantee statistically optimal
error rates of O( n1 Tr H−1 Σ). We show that our characterization of γb,max is indeed tight in the sense that
div
there are problem instances, where γb,max and γb,max
are off by a factor of d and γb,max is still the largest step
size which provides statistically optimal error rates. A critical property of γb,max is that it is dependent on the
properties of the noise covariance Σ. We return to the characterization of γb,max in section 5.4.

5

Main Results

This section presents the main Algorithms analyzed by this paper and their associated guarantees. In
particular, the standard mini-batch tail-averaged SGD Algorithm for the streaming LSR problem is presented
6

as Algorithm 1. Theorem 1 presents an exact characterization of the error achieved by Algorithm 1. Theorem 2
establishes a scalar estimate of the generalization error while dealing with the general agnostic case of LSR,
in order to better illustrate the performance of Algorithm 1. Theorem 2 shows that increasing minibatch sizes
(up to some threshold) and tail averaging achieves the same work (number of first order oracle calls to get
gradient at a point) and statistical error performance as the standard averaged SGD but decreases the depth
significantly, thereby increasing its parallelizability.
Theorem 3 presents the guarantees of Algorithm 2, which is a doubling based Algorithm that decays
the initial error at a superpolynomial rate (and a geometric rate in the presence of an initial error oracle)
while performing the same work and achieving same statistical accuracy (up to constant factors) compared to
standard sequential SGD (with batch size 1). Algorithm 2’s depth nearly matches the depth of standard offline
gradient descent, while providing the same generalization error (upto constants) as offline gradient descent
on every problem instance, while being a single pass streaming Algorithm, thus performing substantially
lower work compared to the standard offline gradient descent. Refer to table 1 for the associated bounds.
The paper’s final result is lemma 4, which presents an instance that shows learning rates necessary to retain
statistical optimality in the agnostic case could be much lower than ones required by realizable case.

5.1

Characterization of SGD with minibatching and tail averaging

This section presents Theorem 1, which provides an exact characterization of the error incurred by minibatched tail-averaged SGD, as described in Algorithm 1, working with a mini-batch size b and step size
γ. Algorithm 1 consists of two phases. The first phase (lines 1-3) is the burn-in phase that consists of “s”
iterations where mini-batch SGD updates are performed without averaging the iterates. The second phase
(lines 5-8) is the tail-averaging phase, where, iterates of mini-batch SGD are averaged. Intuitively, the burn-in
phase is used to decrease the bias which is dependent on initial error kw0 − w∗ k and the averaging phase is
to decrease the variance. Note that the bias exhibits a geometric decay during the unaveraged phase.
−1

def

Σ
2b
R)
, 2b Tr (HL +H
Theorem 1. For all b ≥ 1, let γb,max be defined as γb,max = min{ (b−1)kHk
−1 Σk R2 d }. Suppose we
2 k(HL +HR )
2
run Algorithm 1 with initial point w0 , stepsize 0 < γ < γb,max , minibatch size b, initial iterations s, total number
of samples n, and ∃ universal constants C1 , C2 such that the final error w − w∗ satisfies:

i
h
>
E (w − w∗ )(w − w∗ ) =

n/b

ρb
1
 kE0 kF · M1
GTb −1 (I − γTb )s E0 + C1 ·
− s)2 γ 2
γλmin (H) nb − s
{z
} |
{z
}

( nb
|

"

+

1
I−
n − sb
γ
|

T1

T2

#

1
 (HL −1 + HR −1 ) GTb −1 Σ
−s
{z
}

n
b

T3

n

2

+
γb
|
def

n
b

−1

2 HL HR
−s
{z

−1

Tb

−1

−s+1

ρbb
kΣkF · M2 ,
Σ + C2 ·
(n − bs)λmin (H)λmin (Tb )
{z
}
} |

(2)

T5

T4

>

def

where E0 = (w0 − w∗ )(w0 − w∗ ) and ρb = max(kI − γHk , kI − γTb k) < 1 and M1 , M2 are matrices
whose frobenius norm is less than or equal to 1.
Under the fourth moment assumption on the covariates MI  R2 H, the stepsize γb,max is the maximal
step size we can use with SGD for the general agnostic LSR problem while being able to guarantee minimax
optimal statistical error rates. We return to this step size characterization after presenting this paper’s main
results. Lemma 7 (in the appendix) provides details regarding why this learning rate is a maximal; briefly put,
this stepsize limit characterizes the regime where the operators that govern the SGD update rules possess
favorable contractive properties while guaranteeing minimax optimal statistical rates.
Though (2) looks quite complicated, each of the terms in the equation has a very intuitive conceptual
meaning. The terms dependent on E0 (i.e., T1 and T2 ) are the bias terms while those depending on Σ (i.e.,
T3 , T4 and T5 ) are the variance terms. We observe the geometric decay of the bias terms T1 , T2 during the “s”
initial iterations when the iterates are not averaged. The term T3 (which is the leading order variance term)
1
contains a coefficient n−sb
, which is the number of samples that have been used in the averaging phase.
7

5.2

Benefits of minibatching and tail-averaging

This section presents a means to combine minibatching and tail averaging in order to obtain better performance
than standard averaged SGD for the general agnostic case of LSR. We note that these results naturally
interpolate to results for the realizable case by setting Σ = σ 2 H. Tail averaging allows a burn-in phase
when iterates are not averaged, thereby decaying the initial error geometrically during this phase, while
minibatching decreases the depth (i.e. enables parallelizing SGD) while performing the same amount of work
(up to constants). By the same amount of work, we mean that after using the same number of samples n, the
generalization error is the same as SGD with a batch-size 1 (up to constant factors): i.e., we obtain near linear
parallelization speedups via mini-batching.
Theorem 2. Consider the general agnostic case of LSR in the setting of Theorem 1. If a minibatch size
2
−1
def R dk(HL +HR ) Σk
γ
1
≤ b,max
of b ≥ bthresh = kHk Tr(HL +HR )−1 Σ2 and stepsize γ = 3kHk
is chosen, we obtain the following
6
2
2
generalization error bound on the performance of Algorithm 1:

h
i
>
R(w) − R(w∗ ) = Tr HE (w − w∗ )(w − w∗ )
√
 s 
 n  3κ2 √d
36 dκ2
∗
{R(w0 ) − R(w )} + exp −
≤
{R(w0 ) − R(w∗ )}
2 exp −
n
3κ
3bκ nb − s
−
s
{z
}
{z
} |
|b
T1

+

T2

√

3 Tr H−1 Σ 9κ Tr H−1 Σ 2 dκ2 Tr H−1 Σ
+
exp −
2 +
b( nb − s)
b( nb − s)
b nb − s
{z
} |
{z
|
{z
} |
T3



( nb

− s + 1)
3κ

T5

T4


.

(3)

}

Tail averaging enables the bias terms T1 and T2 to decay exponentially with number of initial iterations
s, which is burn-in phase when iterates of SGD are not averaged. Compared to standard SGD which reduces
the initial error
R(w0 ) as O( s12 ) [1], Algorithm 1 reduces the dependence on R(w0 ) − R(w∗ ) geometrically

s
as exp − κ . We note that this contraction factor is reminiscent of the contraction factor that is observed
with deterministic batch gradient descent. This paper’s analysis of minibatching thus characterizes how the
behavior of mini-batch SGD smoothly
from vanilla
SGD (with a batch size 1) to batch gradient

 interpolates
R2 dk(HL +HR )−1 Σk
2
descent with mini-batch sizes of O kHk Tr((HL +HR )−1 Σ) .
2

Minibatching with a batchsize b = bthresh allows us reduce the depth (serial work) while achieving (up to
−1
constants) the statistically optimal rate of Tr Hn Σ . Note that while standard averaged SGD has a depth of n,
Algorithm 1 has a depth of n/b. We note that this reduction in depth does not incur an extra overhead in terms
of work (the extra work overhead implies that the Algorithm is wasteful of samples, or in other words, we no
longer achieve linear parallelization speedups as a function of increasing the mini-batch size). We refer the
reader to Table 2, which compares the work, depth and accuracy achieved by minibatching and tail-averaging
R2 dk(HL +HR )−1 Σk
as compared to averaged SGD with n samples where minibatch size b = bthresh = kHk Tr(HL +HR )−1 Σ2 and
2
n
number of initial steps s = 2b
.

5.3

(Nearly) Matching the depth of Batch Gradient Descent

The result of section 5.2 establishes a scalar generalization error bound of Algorithm 1 for the general agnostic
case of LSR and showed that the depth (number of sequential updates in our algorithm) is decreased to n/b.
This section builds upon this result to present a simple and intuitive doubling based streaming algorithm that
works in epochs and processes a total of n/2 points. In each epoch, the minibatch size is increased by a factor
of 2 while thus applying Algorithm 1 with twice as many samples as the previous epoch. After running over
n/2 samples using this epoch based Algorithm, we run Algorithm 1 with the remaining n/2 points. Intuitively
each of these epochs decreases the bias linearly and then resets the bias of the next epoch as the total risk of
previous epoch, and then again decays the bias. The final application of the tail-averaging algorithm is to
ensure that the variance is small.

8

Final error
Sequential
averaged SGD
(b = 1)
MinibatchTailAveragingSGD (Algorithm 1)

O

O

√

√

dR4

λ2min n2

dR4
{R(w0 )
λ2min n2

− R(w∗ )} +

σ2 d
n

Runtime

Depth

O (nd)

n

O (nd)

n
bthresh




min
exp − nλ
{R(w0 ) − R(w∗ )} +
R2

σ2 d
n



Table 2: Comparison of minibatching and tail averaging in Algorithm 1 with with averaged SGD, both using
n samples (assuming realizable noise; the agnostic noise bound for Algorithm 1 can be derived similarly).
R2 dk(HL +HR )−1 Σk
1
n
R2
Algorithm 1 uses minibatch size bthresh = kHk Tr(HL +HR )−1 Σ2 = kHk
, γ = 3kHk
and s = 2b
while averaged
2

SGD uses the n samples sequentially. We assume that n >
make any progress [12].

R2
λmin

2

2

which is the minimum number of samples to

The next theorem formalizes 
this intuition and shows that Algorithm 2 indeed improves the depth
∗
0 )−R(w )}
) in the presence of an error oracle that
exponentially from n/bthresh to O κ log(dκ) log( n{R(w
Tr H−1 Σ
provides us with the initial excess risk R(w0 ) − R(w∗ ) and the noise level Tr H−1 Σ.
Theorem 3. Consider the general Agnostic case of LSR. Suppose in MinibatchDoublingPartialAveragingSGD (AlR2 dk(HL +HR )−1 Σk
γ
1
gorithm 2), we use initial minibatchsize of b = kHk Tr(HL +HR )−1 Σ2 , stepsize γ = 3kHk
≤ b,max
and number of
6
2
2
iterations in each epoch being t ≥ 16κ log(dκ). We obtain the following bound on the excess risk of w:
E [R(w)] − R(w∗ ) ≤



2bt
n

t
 3κ log(dκ)

{R(w0 ) − R(w∗ )} + 40

Tr H−1 Σ
.
n

Remarks: The final error again has two parts – the bias term that depends on the initial error R(w0 ) and
the variance term that depends on the statistical noise Tr H−1 Σ.
Algorithm 2 decays the bias at a superpolynomial rate by choosing t large enough. If Algorithm 2 has
access to an initial error oracle that provides R(w0 ) − R(w∗ ), we can run Algorithm 2 with a batch size bthresh
until the excess risk drops to the noise level Tr H−1 Σ and subsequently begin doubling the batch size. Such
an algorithm indeed gives geometric convergence with a generalization error bound as:


Tr(HL + HR )−1 Σ
Tr H−1 Σ
nλmin
E [R(w)] − R(w∗ ) ≤ exp −( 2
)·
{R(w0 ) − R(w∗ )} + 40
,
−1
R · log(dκ) d k(HL + HR ) Σk2
n


∗
0 )−R(w )}
with a depth of O κ log(dκ) log n{R(w
. The proof of this claim follows relatively straightforwardly
Tr H−1 Σ
from the proof of Theorem 3. We note that this depth nearly matches (up to log factors), the depth of standard
offline gradient descent. Our final observation is that the statistical error decays at a rate of O Tr H−1 Σ/n
which is statistically optimal up to constant factors.

5.4

Tradeoffs between mini-batch sizes, learning rates and the role of agnostic noise

While existing sharp non-asymptotic analyses of SGD with large learning rates [9, 10, 17, 1] focus on the
realizable case (i.e. the case with independent noise, so that Σ = σ 2 H), this paper provides results in the
agnostic noise case of LSR that point towards fundamental differences in the behavior of SGD between the
realizable and agnostic noise setting. Moreover, this paper’s results in general agnostic noise case naturally
specialize to existing results in the realizable case with batch size 1. Before discussing these differences, we
first understand mini-batching and related learning rate issues better. We start by writing out the maximal
learning rate that permits statistical optimality:
γb,max = min

n

2b
2b Tr (HL + HR )−1 Σ o
,
(b − 1) kHk2 k(HL + HR )−1 Σk2 R2 d

9

We note that the maximal learning rates γb,max tends to increase linearly as a function of the batch size b, until
R2 dk(HL +HR )−1 Σk
the point that the batch sizes hit a problem dependent limit bthresh = kHk Tr(HL +HR )−1 Σ2 ; this is precisely the
2
range of mini-batch sizes that yield a provably near-linear parallelization speedup over standard SGD with
batch size 1, while retaining statistical optimality. Moreover, when the batch sizes hit bthresh , the first term in
the definition of γb,max becomes active, and we begin using learning rates that mimic batch gradient descent,
i.e., with γb,max = O( λmax1(H) ). Beyond this mini-batching threshold bthresh , mini-batching is inefficient in that
while the variance term still improves as the batch size increases, the rate of decay of the bias does not show
a linear improvement as a function of the batch size (owing to the lack of serial updates), thus leading to
work inefficiency.
2b
Next, we note that in the realizable case, γb,max = min{ (b−1)kHk
, R2b2 }. This fact, coupled with the case of
2

SGD with batch size 1 (as considered in [9, 10, 17]) provides γ1,max = R22 , which matches the best learning
rates known to provide statistical optimality. The corresponding mini-batching threshold for the realizable
R2
case is bthresh = kHk
.
2
This paper’s final result deals with shedding light on some fundamental differences in the behavior of
SGD when dealing with the realizable case as opposed to the agnostic noise case. In particular, the following
lemma presents an instance working with agnostic noise case, wherein, SGD cannot employ the large learning
rates that typically are used in the realizable case while retaining statistical optimality; the proof of this lemma
can be found in the appendix. Note that this behavior is in stark contrast to algorithms such as streaming
SVRG [12], which work with the same large learning rates in the agnostic case as in the realizable case, while
guaranteeing statistically optimal rates.
Lemma 4. Consider a Streaming LSR example with Gaussian covariates with a diagonal second moment matrix
defined by:
(
1
if i = 1
Hii =
1/d
if i > 1
Further, let the noise covariance matrix Σ be diagonal as well, with the following entries:
(
1
if i = 1
Σii =
1/(d − 1)d
if i > 1
For this problem instance, γ1,max ≤

4
1
(d+2)(1+ d
)

is necessary for retaining statistical optimality, while the realizable

version of this problem permits a maximal learning rate ≤
learning rates between the realizable and agnostic case.

d
1 ,
(d+2)(1+ d
)

thus implying an O(d) separation in

The proof of lemma 4 can be found in the appendix in section A.4. As a final remark, note that in the
realizable case, we are unaware of whether there is a separation between O( R12 ) learning rate and the
divergent step size schema as introduced by [1]. Furthermore, it is unclear whether it is possible if we can
prove statistical optimality when working with step sizes that are larger than O( R12 ). However, we believe this
question is relatively of less consequence given the separation that exists between the realizable and agnostic
case, where a much smaller learning rate is necessary for obtaining statistically optimal rates.

6

Experimental simulations

We conduct experiments using a synthetic example to illustrate the implications of our theoretical results. The
data is sampled from a 50− dimensional Gaussian with eigenvalues decaying as { k1 }50
k=1 (condition number
κ = 50), and the variance σ 2 of the noise is 0.01. In this case, our estimated batch size according to Theorem 2
is bthresh = 11. In all our plots, the x-axis is the log of number of iterations (at each iteration, we use “b”
samples in the SGD update rule), and y-axis is the log of the generalization error. We assume that each
iteration takes constant time for all batch sizes; this is just in order to present evidence regarding the tightness
of our mini-batching characterization of the regime of mini-batch sizes where we obtain problem dependent
linear speedups while achieving the same generalization error as SGD with mini-batch size of 1.
We first consider the effect of mini-batching with batch sizes of 1, 4, bthresh = 11, 2 ∗ bthresh = 22 and d = 50.
In figure 1, we observe the generalization error curves obtained by averaging 100 runs of the Algorithm, each
10

(a) Bias term

(b) Variance term

(c) Total Risk

Figure 1: Effect of increasing batch sizes on the Algorithm’s generalization error. As expected, variance
decreases monotonically with increasing batch size. But the bias term shows that the rate of decay increases
till the optimal bthresh . Total risk shows that with optimal batch size, our method indeed obtains the optimal
statistical rate but using smaller number of iterations compared to sub-optimal batch sizes.
run utilizing the same set of samples for all mini-batch sizes. Averaging begins after observing a fixed number
of samples (set as 5κ). We see that the rate of bias decay (figure 1a) increases until reaching a mini-batch size
of bthresh , saturating thereafter; this implies we are wasteful in terms of sample size. As expected, the rate of
decay of variance (figure 1b) is monotonic as a function of mini-batch size. Finally, the overall generalization
error (figure 1c) shows the tightness of our mini-batching characterization. This is because, with a batch size
of bthresh , we obtain a generalization error that is the same as using batch size of 1 with the number of (serial)
iterations that is an order of magnitude smaller. Subsequently, we note that larger mini-batch sizes tend to
worsen generalization error thus depicting the tightness of our characterization of bthresh .
In the next experiment, we fix the batch size (as bthresh ) as well as the number of iterations that we
run mini-batch SGD, and thus hold the total number of samples as 10κ. We then consider the effect of
when averaging begins; in particular, we consider averaging iterates right from the start, as prescribed by
Dèfossez and Bach [1]. We also consider averaging iterates after a quarter/half/three-fourth of total number
of iterations, and consider unaveraged SGD as well. In figure 2, we see that the bias (figure 2a) exhibits
a geometric decay in the unaveraged phase while switching to an inferior O( t12 ) rate with averaging. The
σ2
variance (figure 2b) tends to increase and stabilize at O( bthresh
) in the absence of averaging, while switching to
1
a O( N ) decay rate when averaging begins. The overall generalization error (figure 2c) shows the superiority
of the scheme where averaging after a burn-in period allows us decay the bias towards the noise level at a
geometric rate, following which tail-averaging allows us to decay the variance term, thus providing credence
to our theoretical results that averaging after a burn-in time allows us to obtain better generalization error as
a function of sample size.
Finally, we consider the effect of step size choice; in particular, we fix the number of samples to be 10κ
and begin averaging after observing 2κ samples. We employ a mini-batch size of bthresh = 11 and compute
the corresponding γb,max . We then run mini-batch SGD with a step size γ = c · γbthresh ,max with c varied
as 0.25/0.5/0.75/1/1.1. The results, as observed in figure 3 presents the effect of increasing stepsizes. In
particular, for the bias term (figure 3a), the bias is decays progressively faster during the unaveraged phase as
the step size is increased, until we hit γb,max , where the behavior tends to be composed of a combination of
divergent and non-divergent steps. The trend is very similar with regard to the behavior of the variance term
(figure 3b) and the total risk (figure 3c).

7

Concluding remarks

The framework established as a part of Theorem 1 can be extended in a very straightforward manner to provide
sharp finite sample generalization error bounds for parameter mixing/model averaging schemes [2, 3, 25, 26].
Note that parameter mixing/model averaging is a parallelization scheme that can be used in conjunction with
mini-batching and tail-averaging. Furthermore, as noted by [26], these schemes do not offer any improvement
on the bias part of the error, while improving the variance term owing to the effect of averaging: this suggests

11

(a) Bias term

(b) Variance term

(c) Total Risk

Figure 2: Effect of tail-averaging with mini-batch size of bthresh = 11.

(a) Bias term

(b) Variance term

(c) Total Risk

Figure 3: Effect of increasing step sizes, i.e. with γ = c · γbthresh ,max with c = 0.25/0.5/0.75/1/1.1.
the effectiveness of these schemes when the bias error is roughly at the noise level.
While this paper provided the first rigorous analysis of the effect of mini-batching with SGD for the general
agnostic case of the streaming LSR problem, we believe the effect of mini-batching holds in much more
generality, with corresponding algorithm/problem dependent thresholds on the mini-batch size. In particular,
possible future directions could include understanding mini-batching thresholds and effects for stochastic
approximation with the Logistic Loss using the framework of [16], more general convex functions [9] and for
PCA using the framework of [30]. Moreover, similar mini-batching thresholds could be derived for algorithms
other than SGD, such as that of SVRG, by building on the framework established by [12].

8

Acknowledgments

Sham Kakade acknowledges funding from the Washington Research Foundation Fund for Innovation in
Data-Intensive Discovery. Rahul Kidambi thanks James Saunderson for useful discussions on matrix operator
theory.

12

References
[1] Alexandre Défossez and Francis R. Bach. Averaged Least-Mean-Squares: Bias-Variance Trade-offs and Optimal
Sampling Distributions. In AISTATS, volume 38, 2015.
[2] Gideon Mann, Ryan T. McDonald, Mehryar Mohri, Nathan Silberman, and Dan Walker. Efficient Large-Scale
Distributed Training of Conditional Maximum Entropy Models. In NIPS 22, 2009.
[3] Martin A. Zinkevich, Alex Smola, Markus Weimer, and Lihong Li. Parallelized stochastic gradient descent. In NIPS
24, 2011.
[4] Joseph K. Bradley, Aapo Kyrola, Danny Bickson, and Carlos Guestrin. Parallel Coordinate Descent for L1-Regularized
Loss Minimization. In ICML, 2011.
[5] Feng Niu, Benjamin Recht, Christopher Re, and Stephen J. Wright. Hogwild: A Lock-Free Approach to Parallelizing
Stochastic Gradient Descent. In NIPS 24, 2011.
[6] Mu Li, Tong Zhang, Yuqiang Chen, and Alexander J. Smola. Efficient mini-batch training for stochastic optimization.
In KDD, 2014.
[7] Yuchen Zhang and Lin Xiao. DiSCO: Distributed Optimization for Self-Concordant Empirical Loss. In ICML, 2015.
[8] B. T. Polyak and A. B. Juditsky. Acceleration of Stochastic Approximation by Averaging. SIAM J Control Optim,
volume 30, 1992.
[9] Francis Bach and Eric Moulines. Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Machine
Learning. In NIPS 24, 2011.
[10] Francis R. Bach and Eric Moulines. Non-strongly-convex smooth stochastic approximation with convergence rate
O(1/n). In NIPS 26, 2013.
[11] Francis R. Bach. Adaptivity of averaged stochastic gradient descent to local strong convexity for logistic regression.
JMLR, volume 15, 2014.
[12] Roy Frostig, Rong Ge, Sham M. Kakade, and Aaron Sidford. Competing with the Empirical Risk Minimizer in a
Single Pass. In COLT, 2015.
[13] Shai Shalev-Shwartz and Tong Zhang. Stochastic Dual Coordinate Ascent Methods for Regularized Loss Minimization.
CoRR, abs/1209.1873, 2012.
[14] Alekh Agarwal, Peter L. Bartlett, Pradeep Ravikumar, and Martin J. Wainwright. Information-theoretic lower bounds
on the oracle complexity of stochastic convex optimization. IEEE Transactions on Information Theory, 2012.
[15] D. Ruppert. Efficient Estimations from a Slowly Convergent Robbins-Monro Process. Tech. Report, ORIE, Cornell
University, 1988.
[16] Francis R. Bach. Self-concordant analysis for logistic regression. Elec. J. of Statistics, volume 4, 2010.
[17] Aymeric Dieuleveut and Francis Bach. Non-parametric Stochastic Approximation with Large Step sizes. The Annals
of Statistics, 2015.
[18] Nicolas Flammarion and Francis R. Bach. From Averaging to Acceleration, There is Only a Step-size. In COLT,
volume 40, 2015.
[19] Rie Johnson and Tong Zhang. Accelerating Stochastic Gradient Descent using Predictive Variance Reduction. In
NIPS 26, 2013.
[20] Aaron Defazio, Francis R. Bach, and Simon Lacoste-Julien. SAGA: A Fast Incremental Gradient Method With Support
for Non-Strongly Convex Composite Objectives. In NIPS 27, 2014.
[21] Andrew Cotter, Ohad Shamir, Nati Srebro, and Karthik Sridharan. Better Mini-Batch Algorithms via Accelerated
Gradient Methods. In NIPS 24, 2011.
[22] Ofer Dekel, Ran Gilad-Bachrach, Ohad Shamir, and Lin Xiao. Optimal Distributed Online Prediction Using MiniBatches. JMLR, volume 13, 2012.
[23] Shai Shalev-Shwartz and Tong Zhang. Accelerated Mini-Batch Stochastic Dual Coordinate Ascent. In NIPS 26, 2013.
[24] Martin Takác, Peter Richtárik, and Nati Srebro. Distributed Mini-Batch SDCA. CoRR, abs/1507.08322, 2015.
[25] Yuchen Zhang, John C. Duchi, and Martin Wainwright. Divide and Conquer Ridge Regression: A Distributed
Algorithm with Minimax Optimal Rates. JMLR, volume 16, 2015.
[26] Jonathan Rosenblatt and Boaz Nadler. On the Optimality of Averaging in Distributed Statistical Learning. CoRR,
abs/1407.2724, 2014.
[27] Martin Takác, Avleen Singh Bijral, Peter Richtárik, and Nati Srebro. Mini-Batch Primal and Dual Methods for SVMs.
In ICML, volume 28, 2013.
[28] Horia Mania, Xinghao Pan, Dimitris Papailiopoulos, Benjamin Recht, Kannan Ramchandran, and Michael I. Jordan.
Perturbed Iterate Analysis for Asynchronous Stochastic Optimization. CoRR, abs/1507.06970, 2015.
[29] John C. Duchi, Sorathan Chaturapruek, and Christopher Ré. Asynchronous stochastic convex optimization. CoRR,
abs/1508.00882, 2015.
[30] Prateek Jain, Chi Jin, Sham M. Kakade, Praneeth Netrapalli, and Aaron Sidford. Streaming PCA: Matching Matrix
Bernstein and Near-Optimal Finite Sample Guarantees for Oja’s Algorithm. In COLT, 2016.

13

A

Appendix

A brief note on the organization of the appendix:
Theorem 1 establishes the generalization error of Algorithm 1. The proof of Theorem 1 can be found
in A.1 and consists of three parts, namely:
• Lemma 7, which rigorously establishes convergence (in section A.1.6) and provides an estimate of the
(scalar) contraction factor (in section A.1.7).
• Theorem 5 provides the asymptotic expansion of the finite-sample generalization error of Algorithm 1
on the noiseless problem (also referred to as the bias term), in section A.1.4.
• Theorem 6 provides the asymptotic expansion of the finite-sample generalization error of running
Algorithm 1 by initializing it at the solution w∗ (also referred to as the variance term), in section A.1.5.
The proofs of Theorem 2 and Theorem 3 are established in section A.3.
• Theorem 2 decomposes the tensor bound of Theorem 1 and provides a scalar estimate of the generalization performance of the mini-batch tail averaged SGD based Algorithm (Algorithm 1). The proof of
Theorem 2 can be found in section A.3.1.
• Theorem 3 provides the generalization error guarantee of the doubling based Algorithm (Algorithm 2).
The proof of Theorem 3 can be found in section A.3.2.
The proofs of both Theorem 2 and Theorem 3 utilize the claims established by lemmas in section A.2
along with the tensor bound established in Theorem 1.
Finally, the proof of lemma 4, which describes an instance that yields separation between learning rates
required for the realizable and agnostic case is described in section A.4.
Definition A.1. We introduce some notations relating the current iterate wt and the averaged iterate
Pt+N −1
def
def
wt,N to the optimum w∗ , namely: ηt = wt − w∗ , and η̄t,N = wt,N − w∗ = N1 s=t
ηs . Let Lt =
E [(wt − w∗ ) ⊗ (wt − w∗ )], Lt,N = E [(wt,N − w∗ ) ⊗ (wt,N − w∗ )] denote the covariance of the respective
def

def

parameter errors. Furthermore, we define the following contraction factors, ρT b = kI − γTb (γ)k2 , ρb =
max{ρT b , kI − γHk}.

Theorem 5. (Asymptotic Expansion of the error covariance for the bias term) Let E0 = η0 ⊗ η0 be the covariance
of the initial parameter error. If 0 < γ < γb,max and if t,i = 0 ∀ i = {1, ..., b}, t > 0 (i.e. we deal with the
noiseless problem), then, by averaging the iterates of “b” sample mini-batch SGD for N iterations starting from
the sth iteration, the covariance term of the error has the following asymptotic expansion:
Lbias
s,N = E[η̄s,N ⊗ η̄s,N ] =

 −1
ρs+N
1  −1
−1
s
b
H
+
H
kE0 kF · M1 (4)
−
γI
T
(I
−
γT
)
E
+
C
·
b
b
0
1
L
R
N 2γ2
γλmin (H)N

Where C1 is a universal constant and M1 is a matrix with Frobenius norm ≤ 1.
Theorem 6. (Asymptotic Expansion of the error covariance for the variance term) Let Σ = E[2 x ⊗ x] and let
η0 = 0 (i.e. we start at the solution). If 0 < γ < γb,max , then, by averaging the iterates of “b” sample mini-batch
SGD for N iterations starting from the sth iteration, the covariance term of the error has the following asymptotic
expansion:




1
1
−1
−1
variance
Lt,N
= E[η̄s,N ⊗ η̄s,N ] =
I−
(HL + HR ) (HL )−1 + (HR )−1 − γI Tb −1 Σ
bN
γN
+

+1
ρN
2
b
HL −1 HR −1 Tb −1 Σ + C2 ·
kΣkF · M2
2
γbN
bN λmin (H)λmin (Tb )

Where C2 is a universal constant and M2 is a matrix with Frobenius norm ≤ 1.

14

(5)

Lemma 7. If 0 < γ < γb,max , then, Tb  0 and ρb < 1. Furthermore, for dimension d ≥ 2, the contraction factor
ρb can be upper bounded as:
ρb ≤ max(1 − 2γ(1 −

γ
γ
γ 2
)λmin (H), −1 + 2γ(1 −
)λmax (H) + 2(
) , 1 − γλmin (H))
γb,max
γb,max
γb,max

while in dimension d = 1, a very similar statement holds, i.e.:
ρb ≤ max(1 − 2γ(1 −

γ
γb,max

)λmin (H), −1 + 2γ(1 −

γ
γb,max

)λmax (H) + 2(

γ
γb,max

)2 , |1 − γλmin (H)|)

Remark 1. We note that γ ≤ γb,max is sufficient to ensure Algorithm 1 with constant step sizes achieves
statistically optimal rates. Note, our characterization of γb,max ensures that Tb  0 (refer to the proof of
lemma 7), thus ensuring contraction of all operators that characterize both the burn-in phase as well as
averaging phases of the Algorithm.

A.1

Proofs of Theorems 5, 6 and Lemma 7

We begin by detailing the notation and setup required for proofs of Theorems 5, 6 and lemma 7.
A.1.1

Deriving and setting up the basic recursion

At each iteration t of Algorithm 1, we are provided with b fresh samples {(xti , yti )}bi=1 drawn i.i.d. from the
distribution D. We start by recounting the mini-batch gradient descent update rule that allows us to move
from iterate wt−1 to wt :
b
γX
(hwt−1 , xti i − yti )xti
wt = wt−1 −
b i=1
where, 0 < γ < γb,max is the constant step size that is set to a value less than the maximum allowed learning
rate γb,max . We also recount the definition of wt,N which is the iterate obtained by averaging for N iterations
starting from the tth iteration, i.e.,
t+N −1
1 X
wt,N =
ws
N s=t
Let us first denote the residual error term by i = yi − hw∗ , xi i. By the first order optimality conditions of
w∗ , we observe that  and x are uncorrelated, i.e, E(x,y)∼D [x] = 0. However, independence of  and x is
guaranteed only iff E(x,y)∼D [y|x] = hw∗ , xi, i.e., in other words, if the noise is well specified. Note that the
def

risk of the iterate wt,N is written as R̄t,N = R(w̄t,N ). We can now write out the the generalization error as:
E[R̄t,N ] − R∗ =

1
1
Tr H E η̄t,N ⊗ η̄t,N = Ehη̄t,N , H η̄t,N i
2
2

We now write out the main recursion governing the mini-batch SGD updates in terms of η. :
b
b


γX
γX
ηt = I −
xti ⊗ xti ηt−1 +
ti xti
b i=1
b i=1
b
b


γX
γX
= I−
xti ⊗ xti ηt−1 +
ξti
b i=1
b i=1

def

Where, Ptb =



I−

γ
b

= Ptb ηn−1 + γζtb

Pb
Pb
def
x
⊗
x
and ζtb = 1b i=1 ξti =
ti
ti
i=1

(6)
1
b

Pb

i=1 ti xti .

Equation 6 automatically brings

out the “operator” view of analyzing the (expected) covariance of the centered estimate E [ηt ⊗ ηt ]. We now
note the following about the covariance of ζtb :
E ζtb ⊗ ζt0 b =

1 X
E ξti ⊗ ξt0 j
b2 i,j
15

=

b
h1 X
i
1
E
ξ
⊗
ξ
1[t = t0 ] = Σ
ti
ti
2
b i=1
b

1[t = t0 ]

(7)

Where, 1[.] is the indicator function, and equals 1 if the argument inside [.] is true and 0 otherwise. We
note that the expectation of the cross terms in equation 7 is zero owing to independence of the samples
{xti , yti }bi=1 as well as between {xti , yti }bi=1 , {xt0 i , yt0 i }bi=1 ∀ t 6= t0 and owing to the first order optimality
conditions. Owing to the invariance of ζtb on the iteration t, context permitting, we sometimes drop the
iteration index t from ζtb and simply refer to it as ζb . Q
t
Next we expand out the recurrence (6). Let Qj,t = ( k=j Pkb )T with the convention that Qt0 ,t = I ∀ t0 > t.
With this notation we have:
ηt = Ptb ηn−1 + γζtb
= Ptb Pt−1,b ...P1,b η0 + γ

t−1
X

{Ptb ....Pt−j+1,b }ζt−j,b

j=0

= Q1,t η0 + γ

t−1
X

Qt−j+1,t ζt−j,b

j=0

= Q1,t η0 + γ

t
X

Qj+1,t ζj,b .

(8)

j=1

A.1.2

The terms in the averaged iterate

Let us consider beginning averaging after a certain number of iterations “s”, i.e for t > s. In particular,
after averaging the updates for the next “N ” steps, we examine the quantity η̄s+1,N = ws+1,N − w∗ , where
ws+1,N is the iterate we obtain by averaging for “N ” iterations starting from iteration s + 1. We write out the
expression for η̄s+1,N starting out from equation 8:
η̄s+1,N =

A.1.3

s+N
1 X
ηt
N t=s+1

=

s+N
s+N
t
o
γ X X
1n X
Q1,t η0 +
Qj+1,t ζj,b
N t=s+1
N t=s+1 j=1

=

s+N
s s+N
o
  s+N
o
X s+N
X
1n X
γ n X X
Q1,t η0 +
Qj+1,t ζj,b +
Qj+1,t ζj,b
N t=s+1
N
j=1 t=s+1
j=s+1 t=j

=

s+N
s
s+N
s+N
s+N
o
o
o
1n X
γ Xn X
γ X nX
Q1,t η0 +
Qj+1,t ζj,b +
Qj+1,t ζj,b
N t=s+1
N j=1 t=s+1
N j=s+1 t=j

(9)

Computing the outerproduct

Firstly, the outer product of η̄s+1,N with itself contains three terms, each involving the outerproduct of each
term in equation 9 with itself. The 6 other terms are the cross terms all of whose expectation are zero. We
note that E[ξt,i ] = E[t,i xt,i ] = 0 ∀ i, t, and using this, it is clear to see the outer product of the cross term
between terms 2 and 3 is zero, since they deal with disjoint samples j = 1, ..., s and j = s + 1, ..., s + N and
the q.,. doesn’t deal with the same sample index j. The outer product of the 1st term with the 2nd and 3rd
term, E[x(i) x(j) x(k) ] = 0 (as noted by [9, 1], and if this is not satisfied we lose at most a constant factor two
in our generalization error estimate), and as noted in section 4. This leaves us to deal with the 3 of the terms,
involving outer product of the each of the terms with itself. Finally, we note that the bound in the theorem 5
refers to the outerproduct computation of the first term, and that of theorem 6 refers to the outerproduct
computation of the second and third term.

16

A.1.4

Proof of Theorem 5

We begin by noting that this proof generalizes the proof of Défossez and Bach [1] to account for mini-batching
and tail-averaged SGD. In this proof, we assume that the problem is noiseless, i.e., we assume j,k = 0 a.s.
∀ j, k.
" s+N
#
n X
o
n s+N
oT
X
1
E[η̄s+1,N ⊗ η̄s+1,N ] = 2 E
Q1,t η0 ⊗ η0
Q1,t0
N
t=s+1
t0 =s+1


s+N
X
1
Q1,t η0 ⊗ η0 QT1,t0 
= 2E
N
t,t0 =s+1
"
#
s+N
s+N
o
X n
1 X
T
T
T
T
E Q1,t η0 ⊗ η0 Q1,t +
= 2
Q1,t η0 ⊗ η0 Q1,t Qt+1,t0 + Qt+1,t0 Q1,t η0 ⊗ η0 Q1,t
N t=s+1
t0 =t+1
"
#
s+N
s+N
X
1 X
(t0 −t)
(t0 −t)
I+
(I − γH)R
+ (I − γH)L
E[Q1,t η0 ⊗ η0 QT1,t ] (∗)
= 2
N t=s+1
0
t =t+1
"
#
s+N
s+N
X
X
1
(t0 −t)
(t0 −t)
I+
= 2
(I − γH)R
+ (I − γH)L
(I − γTb )t E0 (∗∗)
N t=s+1
0
t =t+1

=

1
N2

s+N
X

h
I + (γHR )−1 {(I − γH) − (I − γH)s+N −t+1 }R

t=s+1

i
+ (γHL )−1 {(I − γH) − (I − γH)s+N −t+1 }L (I − γTb )t E0

(∗ ∗ ∗)

0

where, (∗) follows since E [Qt+1,t0 ] = (I − γH)t −t , and (∗ ∗ ∗) follows by summing up the geometric series
1
(which converges when γ ≤ kHk
, which is the case as discussed in section A.1.6), and finally, (∗∗) follows
2
through the following argument:
E[Q1,t η0 ⊗ η0 QT1,t ] = E[Ptb Pt−1,b ...P1b η0 ⊗ η0 P1b ...Pt−1,b Ptb ]
= E[E[Ptb (Pt−1,b ...P1b η0 ⊗ η0 P1b ...Pt−1,b )Ptb |Ft−1 ]]
where, Ft−1 refers to the sigma-field defined by samples {{(xt0 ,i , yt0 ,i )}bi=1 }t−1
t0 =1 . We note that (Pt−1,b ...P1b η0 ⊗
η0 P1b ...Pt−1,b ) is a constant w.r.t. the inner expectation because of the conditioning Ft−1 . Now, the observation can be seen clearly using the argument below (with C being a matrix that is independent w.r.t. the
expectation):
b
b
h
i
γX
γX
E[Ptb CPttb ] = E (I −
xti ⊗ xti )C(I −
xti ⊗ xti )
b i=1
b i=1

γ2
γ 2 (b − 1)
E[xT Cx x ⊗ x] +
E[xT1 Cx2 x1 ⊗ x2 ]
b
b
γ(b − 1)
γ
HL HR ))C
= (I − γ(HL + HR − M −
b
b
= (I − γTb )C
= C − γCH − γHC +

(10)

where, we apply the fact that the samples {(xti , yti )}bi=1 are independent and identically distributed draws
from the distribution D. The final expression is obtained by applying the argument recursively for all
t0 = t, t − 1, .., 1.
Next we bound the exponential term, in particular:
AN =

s+N
i
−1 X h
−t+1
s+N −t+1
(γHR )−1 (I − γH)s+N
+ (γHL )−1 (I − γH)L
(I − γTb )t E0
R
2
N t=s+1

=⇒ kAN kF ≤

2
ρs+N +1 kE0 kF
γλmin (H)N b
17

The remaining terms are bounded as the following:
E[η̄s+1,N

s+N
1 X
⊗ η̄s+1,N ] =
[HL −1 + HR −1 − γI](I − γTb )t E0 + AN
γN 2 t=s+1

=

1
[HL −1 + HR −1 − γI]Tb−1 [(I − γTb )s+1 − (I − γTb )s+N +1 ]E0 + AN
γ2N 2

Considering the exponential terms separately, we write out the covariance operator as:
E[η̄s+1,N ⊗ η̄s+1,N ] =

1
[HL −1 + HR −1 − γI]Tb−1 (I − γTb )s+1 E0 + AN + BN
γ2N 2

with BN defined as:
−1
[HL −1 + HR −1 − γI]Tb−1 (I − γTb )s+N +1 E0
γ2N 2


2
1
+1
=⇒ kBN kF ≤ 2 2
− γ ρs+N
kE0 kF
Tb
γ N λmin (Tb ) λmin (H)
BN =

The last equation above assumes that γ < λmin2(H) . Finally, using triangle inequality that kAN + BN kF ≤
kAN kF + kBN kF and hence, when we begin averaging starting at iteration s + 1 and average for N iterations
subsequently, we get a geometric decay of the bias term. The final covariance operator of the bias term thus
has the following asymptotic expansion:
E[η̄s+1,N ⊗ η̄s+1,N ] =

ρbs+N +1
1
−1
−1
−1
s+1
[H
+
H
−
γI]T
(I
−
γT
)
E
+
O(
kE0 kF )
L
R
b
0
b
γ2N 2
γλmin (H)N

The above bound requires the following conditions to hold: (i) Tb (γ)  0, (ii) γ < λmax1(H) . This wraps up
the proof of Theorem 5 and indicates the geometric decay of the bias term during the unaveraged phase and
the O( N12 ) decay during the averaging phase.
A.1.5

Proof of Theorem 6

We begin by noting that this proof generalizes the proof of Défossez and Bach [1] to account for mini-batching
and tail-averaged SGD. For bounding the variance part of the generalization error, we go over a sharp analysis
of the stochastic process defined by the mini-batch tail averaged SGD initialized at the solution w∗ , i.e.
w0 = w∗ , and thus the process is driven solely by the noise. We note that we require bounding the behavior
of two terms, the first of which is defined by the noise introduced during the unaveraged phase, and the
second being one that deals with the noise introduced during the tail averaging phase.


s
s+N
s+N
s+N
2
X
X
X
X
γ
E[η̄s+1,N ⊗ η̄s+1,N ] = 2 
E[Qj+1,t ζjb ⊗ ζjb QTj+1,t0 ] +
E[Qj+1,t ζjb ⊗ ζjb QTj+1,t0 ]
N
j=1 0
j=s+1 0
t,t =j

t,t =s+1

2

=

γ
(K1 + K2 )
N2

(11)

We consider the first term, and specifically the inner two terms of the summation:
K1j =

s+N
X

s+N
X

E[Qj+1,t ζjb ⊗ ζjb QTj+1,t0 ]

t=s+1 t0 =s+1

=

s+N
X
t=s+1

s+N
t−1
h
i
X
X
E Qj+1,t ζjb ⊗ ζjb QTj+1,t +
Qj+1,t ζjb ⊗ ζjb QTj+1,t QTt+1,t0 +
Qt0 +1,t Qj+1,t0 ζjb ⊗ ζjb QTj+1,t0
t0 =t+1

t0 =s+1

18

=

s+N
X

s+N
t−1
i
h
X
X
0
0
E Qj+1,t ζjb ⊗ ζjb QTj+1,t +
Qj+1,t ζjb ⊗ ζjb QTj+1,t (I − γH)t −t +
(I − γH)t−t Qj+1,t0 ζjb ⊗ ζjb QTj+1,t0
t0 =t+1

t=s+1

=

s+N
X

t0 =s+1

h
i
I + (γHR )−1 {(I − γH) − (I − γH)s+N −t+1 }R + (γHL )−1 {(I − γH) − (I − γH)s+N −t+1 }L (I − γTb )t−j E[ζjb ⊗ ζjb ]

t=s+1

=

s+N
i
1 X h
I + (γHR )−1 {(I − γH) − (I − γH)s+N −t+1 }R + (γHL )−1 {(I − γH) − (I − γH)s+N −t+1 }L (I − γTb )t−j Σ
b t=s+1

Note, the second to last line above is obtained by seeing that we can swap the summations for the third term
and then replace the roles of t and t0 and applying equation 10, following which we note that E[ζjb ⊗ζjb ] = 1b Σ
from equation 7 for the equation in the last line above. And hence, the first term K1 is:
K1 =

s
X

K1j

j=1

=

s+N
s
i
1X X h
I + (γHR )−1 {(I − γH) − (I − γH)s+N −t+1 }R + (γHL )−1 {(I − γH) − (I − γH)s+N −t+1 }L (I − γTb )t−j Σ
b j=1 t=s+1

(12)

In a very similar manner, it is easy to observe the second term is:
K2 =

s+N s+N
i
1 X Xh
I + (γHR )−1 {(I − γH) − (I − γH)s+N −t+1 }R + (γHL )−1 {(I − γH) − (I − γH)s+N −t+1 }L (I − γTb )t−j Σ
b j=s+1 t=j

(13)
2

γ
If K1 and K2 can be bounded, then it is clear that the variance decays as N
2 (K1 + K2 ) as is seen from
equation 11. Let us begin with K1 (from equation 12). We first exchange the two summations and sum up
the resulting series w.r.t. j.
s+N
1 X h
I + (γHR )−1 {(I − γH) − (I − γH)s+N −t+1 }R + (γHL )−1 {(I − γH)−
K1 =
b t=s+1
i
h
i
(I − γH)s+N −t+1 }L (γTb )−1 (I − γTb )t−s − (I − γTb )t Σ
N

1 Xh
I + (γHR )−1 [(I − γH) − (I − γH)N −t+1 ]R + (γHL )−1 [(I − γH)−
b t=1
i
h
i
(I − γH)N −t+1 ]L (γTb )−1 (I − γTb )t − (I − γTb )t+s Σ
=

Next we bound the exponential terms, i.e.
N
i
h
i
1 Xh
(γHL )−1 + (γHR )−1 − I (γTb )−1 (I − γTb )t − (I − γTb )t+s Σ + A0N
b t=1
i
h
ih
i
1h
= (γHL )−1 + (γHR )−1 − I (γTb )−2 (I − γTb ) − (I − γTb )s+1 I − (I − γTb )N Σ + A0N
b

K1 =

where A0N is:
A0N = −

N
i
h
i
1 Xh
−t+1
−t+1
(γHR )−1 (I − γH)N
+ (γHL )−1 (I − γH)N
(γTb )−1 (I − γTb )t − (I − γTb )t+s Σ
R
L
b t=1

whose norm is bounded as:
kA0N kF ≤

2N
γ 2 bλmin (H)λmin (Tb )

19

+1
+s+1
(ρN
+ ρN
)kΣkF
b
b

0
Again, we accumulate the exponential terms existing in K1 in BN
, leading to:

i
h
i
1h
(γHL )−1 + (γHR )−1 − I (γTb )−2 (I − γTb ) − (I − γTb )s+1 Σ + B0N + A0N
b
h
i
i
1h
= (γHL )−1 + (γHR )−1 − I (γTb )−2 (I − γTb ) Σ + C0N + B0N + A0N
b

K1 =

(14)

where, B0N is:
B0N = −

i
h
i
1h
(γHL )−1 + (γHR )−1 − I (γTb )−2 (I − γTb ) − (I − γTb )s+1 (I − γTb )N Σ
b

whose norm is bounded using triangle inequality as:
kB0N kF ≤

1
bγ 3 (λmin (Tb ))2

(

2
N
− γ)(ρT b + ρs+1
T b )ρT b kΣkF
λmin (H)

finally, C0N is:
C0N = −

i
1h
(γHL )−1 + (γHR )−1 − I (γTb )−2 (I − γTb )s+1 Σ
b

whose norm is bounded as:
kC0N k ≤

2
1
(
− γ)ρs+1
T b kΣkF
bγ 3 (λmin (Tb ))2 λmin (H)

which rounds up bounding all terms of equation 14, except for the leading order term. This term ends up
canceling with a term involved in bounding K2 in equation 17.
We now begin bounding the terms in K2 (from equation 13). Towards this effect, we begin by interchanging
the summations in K2 :
t
s+N
1 X X h
I + (γHR )−1 [(I − γH) − (I − γH)s+N −t+1 ]R +
b t=s+1 j=s+1
i
(γHL )−1 [(I − γH) − (I − γH)s+N −t+1 ]L (I − γTb )t−j Σ

K2 =

s+N
1 X h
I + (γHR )−1 [(I − γH) − (I − γH)s+N −t+1 ]R +
b t=s+1
i
(γHL )−1 [(I − γH) − (I − γH)s+N −t+1 ]L (γTb )−1 [I − (I − γTb )t−s ]Σ

=

N

1 Xh
I + (γHR )−1 [(I − γH) − (I − γH)N −t+1 ]R +
b t=1
i
h
i
(γHL )−1 [(I − γH) − (I − γH)N −t+1 ]L (γTb )−1 I − (I − γTb )t Σ
=

We begin considering the exponential terms, i.e.,:
N

K2 =

i
h
i
1 Xh
(γHR )−1 + (γHL )−1 − I (γTb )−1 I − (I − γTb )t Σ + A00N
b t=1

where A00N is:
N

A00N = −

i
h
i
1 Xh
N −t+1
−t+1
−1
t
(γHR )−1 (I − γH)R
+ (γHL )−1 (I − γH)N
(γT
)
I
−
(I
−
γT
)
Σ
b
b
L
b t=1

= B00N + C00N
20

where,
N
i
1 Xh
N −t+1
−t+1
(γHR )−1 (I − γH)R
+ (γHL )−1 (I − γH)N
(γTb )−1 Σ
L
b t=1
h
i
h
ii
1h
+1
+1
+ (γHL )−2 (I − γH)L − (I − γH)N
(γTb )−1 Σ
= − (γHR )−2 (I − γH)R − (I − γH)N
R
L
b
= F00N + G00N

B00N = −

Where F00N is:

1
(γHR )−2 (I − γH)R + (γHL )−2 (I − γH)L (γTb )−1 Σ
b

1 
= − 3 (HR −2 − γHR −1 ) + (HL −2 − γHL −1 ) Tb −1 Σ
bγ

1 
= − 3 (HL −1 + HR −1 − γI)(HL −1 + HR −1 ) − HL −1 HR −1 − HR −1 HL −1 Tb −1 Σ
bγ

1 
= − 3 (HL −1 + HR −1 − γI)(HL −1 + HR −1 ) − 2HL −1 HR −1 Tb −1 Σ
bγ

F00N = −

(15)

The norm of F00N is bounded along with the leading order variance term. G00N is then written out as:
G00N =


1
+1
+1
(γHR )−2 (I − γH)N
+ (γHL )−2 (I − γH)N
(γTb )−1 Σ
R
L
b

whose norm is bounded as:
kG00N kF ≤

+1
ρN
Hb
kΣkF
bγ 3 (λmin (H))2 λmin (Tb )

N

C00N =

i
1 Xh
−t+1
N −t+1
−1
(γHR )−1 (I − γH)N
+
(γH
)
(I
−
γH)
(γTb )−1 (I − γTb )t Σ
L
R
L
b t=1

whose norm is bounded as:
kC00N kF ≤

2
+1
(N ρN
)kΣkF
b
γ 2 bλmin (H)λmin (Tb )

Now, K2 is written as the following:
i
Nh
K2 =
(γHR )−1 + (γHL )−1 − I (γTb )−1 Σ + D00N + A00N
b

(16)

where the N upfront is the time for which the iterates have been averaged, and this is the term that leads to
γ2
1
00
N rate of decay for the variance term (which has the factor of N 2 K2 ), and DN is obtained through summing
the resulting geometric series:
i
h
i
1h
D00N = − (γHR )−1 + (γHL )−1 − I (γTb )−2 (I − γTb ) − (I − γTb )N +1 Σ
b
i
1h
= − (γHR )−1 + (γHL )−1 − I (γTb )−2 (I − γTb )Σ + E00N
(17)
b
where, E00N is:
E00N =

i
1h
(γHR )−1 + (γHL )−1 − I (γTb )−2 (I − γTb )N +1 Σ
b

whose norm is bounded as:
kE00N kF ≤

1
γ 3 b(λ

min (Tb

))2

(

2
+1
− γ)ρN
T b kΣkF
λmin (H)

21

We note that the leading term in D00N in equation 17 is exactly the negative of the leading order term in
equation 14, thus canceling out without requiring us to bound their behavior. This allows us to provide the
γ2
final theorem statement, combining the leading order terms from equations 15 and 16 along with the N
2
outside the summation,

1 
(γHR )−1 + (γHL )−1 − I (γTb )−1 Σ
bN

1 
−
(HL −1 + HR −1 − γI)(HL −1 + HR −1 ) − 2HL −1 HR −1 Tb −1 Σ
bγN 2

E[η̄s,N ⊗ η̄s,N ] =

+1
ρN
b
kΣkF )
bN λmin (H)λmin (Tb )




1
1
=
(HL −1 + HR −1 ) (HR )−1 + (HL )−1 − γI Tb −1 Σ
I−
bN
γN

+ O(

+

+1
ρN
2
−1
−1
−1
b
H
H
T
Σ
+
O(
kΣkF )
L
R
b
γbN 2
bN λmin (H)λmin (Tb )

with which the proof of Theorem 6 concludes.
A.1.6

Some intermediate Lemmas required to prove Lemma 7

We establish a sequence of lemmas in this subsection in order to prove that the operators I − γH and I − γTb
contract with the choice of our stepsizes as described in lemma 7. We then provide a scalar bound for the
contraction factor, and this is done by generalizing the lemma appearing in Défossez and Bach [1] to the
mini-batching case.
Lemma 8. Without loss of generality, let us assume that the eigenvalues {λi }di=1 > 0 of H be arranged in
decreasing order of their values. If γ > 0 is chosen such that Tb (γ)  0, i.e.:
hA, Tb Ai = 2 Tr AT HA −

γ
b−1
Tr AT MA − γ
Tr AT HL HR A > 0 ∀ A ∈ S(d)
b
b

then the following upper bound on the stepsize γ holds: γ <

λi +λj
b
b+1 ( λi λj )

holds ∀ {i, j} ∈ {1, 2, ..., d}×{1, 2, ..., d}

Proof. This proof of this lemma generalizes the proof appearing in Défossez and Bach [1] to establish bounds
on the stepsize for the mini-batching case.
Expanding out the inner product, we have:
γ
b−1
M − γ(
)HL HR )A)
b
b
b−1
γ
= 2 Tr(AT HA) − Tr(AT MA) − γ(
) Tr(AT HL HR A)
b
b
γ
b−1
= 2 Tr(AT HA) − E(xT Ax)2 − γ(
)E(xT1 Ax2 )2
b
b

hA, Tb Ai = Tr(AT (HL + HR −

Where, (x1 , y1 ), (x2 , y2 ) and (x, y) are independent draws from the distribution D. Next, we consider the
following:
2 Tr (AT HA) −

γ
b−1
γ
b−1
[Tr (AH)]2 − γ(
) Tr(AT HL HR A) = 2 Tr (AT HA) − [Tr AE [x ⊗ x]]2 − γ(
) Tr(AHAH)
b
b
b
b
γ
b−1
≥ 2 Tr (AT HA) − E(xT Ax)2 − γ(
) Tr(AT HAH)
b
b
= hA, Tb Ai > 0

Where, we note that the second line follows through the application of Jensen’s inequality. Further, Let
H = UDUt to be the eigen decomposition of H so that D contains the eigenvalues {λi }di=1 in its diagonal. For
any (i, j) ∈ {1, ..., d} × {1, ..., d}, we set the matrix A = ci ui ⊗ ui + cj uj ⊗ uj so as to observe the following:
2 Tr (AT HA) −

By setting ci =

γ
b−1
γ
b−1 2 2
[Tr AH]2 − γ(
) Tr(AT HL HR A) = 2(c2i λi + c2j λj ) − (ci λi + cj λj )2 − γ(
)(ci λi + c2j λ2j )
b
b
b
b
>0

q

λj
λi

and cj =

q

λi
λj ,

the result follows.
22

b
Corollary 8.1. By setting i = j = 1 in lemma 8, we observe that the condition Tb  0 allows γ < ( b+1
)·

2
λmax (H) .

Remark 2. We note that an immediate consequence of corollary 8.1 is that by setting γ such that Tb  0,
−I ≺ I − γH ≺ I.
Lemma 9. With Tb = HL + HR − γb M − γ( b−1
b )HL HR defined as before, if γ > 0 is chosen such that
hZ, Tb Zi > 0 ∀ Z ∈ S(d), then, ∀ d ≥ 2, we have:
hZ, (I − γTb ) Zi ≥ −hZ, Zi ∀ Z ∈ S(d).
while in dimension d = 1 we have the following condition:
I − γTb ≥ 0
Before we prove the above lemma, we prove a simpler version of this lemma where a stronger conclusion
holds, but only ∀ Z ∈ S + (Rd ) (as opposed to the statement of lemma 9 where Z ∈ S(d)).
Lemma 10. If γ <

1
kHk2 ,

the following statement holds:
hZ, (I − γ(HL + HR − γHL HR )) Zi ≥ 0 ∀ Z ∈ S(d).

Proof. Assuming the eigenvalues of H are {λi }di=1 arranged in decreasing order, the proof requires the
observation that the eigenvalues of (HL + HR − γbHL HR ) are respectively {λi + λj − γλi λj }di,j=1 . This
is because if {ui }di=1 are the eigenvectors of H, then, {ui ⊗ uj }di,j=1 are the eigenvectors for HL , HR and
HL HR with corresponding eigenvalues being λi , λj and λi λj respectively. This is because:
HL ui ⊗ uj = Hui ⊗ uj = λi ui ⊗ uj
HR ui ⊗ uj = ui ⊗ uj H = λj ui ⊗ uj
HL HR ui ⊗ uj = Hui ⊗ uj H = λi λj ui ⊗ uj
If we prove that the eigenvalues of (HL + HR − γHL HR ) are all upper bounded by γ1 , the lemma is proven
since this implies HL + HR − γHL HR ≺ γ1 I.
We consider two cases, (i) when i = j, in which case we deal with eigenvalues of the form 2λi − γλ2i . If we
require this to be less than γ1 , then, (γλi − 1)2 > 0 which holds automatically ∀ γ. (ii) When i 6= j, we require
λi + λj − γλi λj < γ1 , or, γ 2 λi λj − γ(λi + λj ) + 1 > 0. If we assume without loss of generality, λi > λj , then,
the solution of interest is γ < λ1i . Taking the intersection of all the solutions obtained from solving equations
∀ i, j, we obtain, γ < λmax1(H) . We consider lemma 8 with i = 1 and j 6= 1. In this case, we know that
b(λ

(H)+λ )

b(λ

(H)+λ )

λ +λ

(H)

λ

max
j
max
j
j
max
γ < (b+1)λ
. The proof completes iff λmax1(H) < (b+1)λ
or b+1
implying b > λmax j(H) ,
b <
λj
max (H)λj
max (H)λj
which is true ∀ b ≥ 1. We also note that the result holds even in the case that there exists eigenvalues with
multiplicity ≥ 1 since that is covered through case(i).

Lemma 11. For any PSD matrix A ∈ S + (Rd ), the following inequality holds:
hA, HL HR Ai ≤ hA, MAi
Proof. The claim is proven through the following argument:
hA, HL HR Ai = hA, HAHi


= Ex1 ,x2 (x1 > Ax2 )2

2
≤ Ex x> Ax


≤ Ex (x> Ax)2
= hA, MAi.

23

Lemma 12. Under the hypothesis of Lemma 9, we have:
hZ, (I − γTb ) Zi ≥ 0 ∀ Z ∈ S + (Rd ).
Proof of Lemma 12. The first observation is that for any PSD matrix Z ∈ S + (Rd ), using lemma 11, hZ, HL HR Zi ≤
hZ, MZi.
This implies that hZ, Tb Zi = hZ, (HL + HR − γb M − γ( b−1
b )HL HR )Zi ≤ hZ, (HL + HR − γHL HR )Zi.
The proof rounds up by by employing Lemma 10.
Proof of Lemma 9. The proof goes via a contradiction. Let us assume that there exists Z ∈ S(d) such that the
conclusion of Lemma 9 is false. So we have hZ, (I − γTb ) Zi < −hZ, Zi, which implies that
2

hZ, Tb Zi > 2 kZkF /γ

(18)

Using the spectral decomposition theorem, let Z = Z1 −Z2 where Z1 , Z2 ∈ S + (Rd ) and Z1 Z2 = Z2 Z1 = 0. Let
2
2
2
2
b def
b Tb Zi
b > 0,
Z
= Z1 + Z2 . This means that kZk = Zb
= kZ1 k + kZ2 k . Since by hypothesis we have hZ,
F

adding this to (18), we obtain

F

F

F

b Tb Zi
b > 2 kZk2 /γ
hZ, Tb Zi + hZ,
F
2

2

⇒ 2hZ1 , Tb Z1 i + 2hZ2 , Tb Z2 i > 2(kZ1 kF + kZ2 kF )/γ
2

⇒ either

hZ1 , Tb Z1 i > kZ1 kF /γ

⇒ either

hZ1 , (I − γTb )Z1 i < 0

or
or

2

hZ2 , Tb Z2 i > kZ2 kF /γ
hZ2 , (I − γTb )Z2 i < 0.

This contradicts with the conclusion of Lemma (12) which shows that for every Z ∈ S + (Rd ), we have
hZ, (I − γTb )Zi ≥ 0. This proves that there does not exist Z ∈ S(d) for which (18) holds, thus concluding the
proof of Lemma 9.
Note that in dimension d = 1, we have:
  γ  
b − 1  2 2
I − γTb = 1 − γ(2E x2 − E x4 − γ(
)E x )
b
b
  γ2  
b − 1  2 2
= 1 − 2γE x2 + E x4 + γ 2 (
)E x
b
b
 
 2
≥ 1 − 2γE x2 + γ 2 E x2
 
= (1 − γE x2 )2 > 0
 2
 
with the inequality relating the second to the fourth moment namely, E x2 ≤ E x4 being through the
application of Jensen’s inequality.
T
Lemma 13. With γ > 0, Tb = HL + HR − γb M − γ( b−1
b )HL HR , we have ∀A ∈ S(d), 2 Tr(A HA) −
γ
b−1
T
2
T
b E(x Ax) − γ( b ) Tr A HL HR A > 0, then,

• kI − γTb k < 1
• kI − γHk < 1
Proof. This lemma generalizes the lemma governing contraction of all associated operators appearing in DéfosT
sez and Bach [1]. We note that the condition ∀ A ∈ S(d), 2 Tr(AT HA)− γb E(xT Ax)2 −γ( b−1
b ) Tr A HL HR A >
b
2
2
0 implies Tb  0, with the step size set as 0 < γ < γb,max < b+1 ( λmax (H) ) < λmax (H) .
Next, to prove kI − γTb k < 1, we require −I ≺ I − γTb ≺ I, the right hand side of which follows from
the fact that Tb  0. The left hand side follows using lemma 9.
We know that H  0. Next, to ensure kI − γHk < 1, we require −I ≺ I − γH ≺ I. The right side of the
inequality holds owing to the fact that H  0. The left side is equivalent to having H ≺ γ2 I, or equivalently,
γ < λmax2(H) which follows from corollary 8.1 derived from lemma 8.

24

A.1.7

Proof of Lemma 7

We first note that with our choice of γb,max , we have Tb ≥ 0. This is because
 γ(b − 1)
γ 
E hx, Axi2 −
hH, AHAi
b
b
h
i
γ
γ(b − 1)
2
2
≥ 2 Tr (AHA) − E kxk kAxk −
kHk Tr (AHA)
b
b
h
i
γ
γ(b − 1)
2
≥ 2 Tr (AHA) − R2 E kAxk −
kHk Tr (AHA)
b
b



γ
R2 + (b − 1) kHk Tr (AHA) ,
≥ 2−
b

hA, Tb Ai = 2 Tr (AHA) −

and using the definition of γb,max finishes the claim.
The rest of the proof is a generalization of the bound on the contraction factor provided by Défossez and
Bach [1]. As we operate with a step size 0 < γ < γb,max , all the operators contract, i.e., ρb < 1 and Tbγ  0 as
per lemma 13. Furthermore, given the largest and smallest eigenvalues of H as λmax (H) and λmin (H), then,
we denote the largest and smallest eigenvalues of HL + HR as 2λmax (H) and 2λmin (H).
We first provide bounds on the operator norm of Tbγ . Let λmax (Tbγ ) and λmin (Tbγ ) represent the largest and
smallest eigenvalues of the operator Tbγ . Then,
ρTbγ = max{1 − γλmin (Tbγ )), γλmax (Tbγ ) − 1}
For a particular step size γ, we define α =
Tbγ = HL + HR −

γ
γb,max

implying the following:

b−1
γ
M − γ(
)HL HR
b
b

= (1 − α)(HL + HR ) + α(HL + HR −

γb,max
b−1
M − γb,max (
)HL HR )
b
b

γb,max

= (1 − α)(HL + HR ) + αTb
γ

Since Tb b,max  0, Tbγ  2(1 − α)λmin (H)I, implying λmin (Tbγ ) ≥ 2(1 − α)λmin (H).
γ
γ
2
I. This implies that,
Next, −I ≺ I − γb,max Tb b,max implying Tb b,max ≺ γb,max
γ

Tbγ = (1 − α)(HL + HR ) + αTb b,max
2α
≺ 2(1 − α)λmax (H)I +
I
γb,max
implying λmax (Tbγ ) ≤ 2(1 − α)λmax (H) +

2α
γb,max .

This when plugged into the operator norm of Tbγ yields,

α
)}
γb,max
γ
γ
γ
= max{1 − 2γ(1 −
)λmin (H), −1 + 2γ(1 −
)λmax (H) + 2(
)2 }
γb,max
γb,max
γb,max

ρTbγ = max{1 − 2α(1 − α)γb,max λmin (H), −1 + 2αγb,max ((1 − α)λmax (H) +

Similarly, the operator norm of I − γH is bounded as: ρHb = 1 − γλmin (H), since in dimension d ≥ 2,
1 − γλmin (H) ≥ γλmax (H) − 1. In d=1, ρHb = |1 − γλmin (H)|. This rounds up the claim in lemma 7.
A.1.8

Proof of Theorem 1

The proof of Theorem 1 builds on the results of lemmas 7, 13 and Theorems 5, 6.
By instantiating Equation 4 of Theorem 5 and Equation 5 of Theorem 6 with N =
the output of the Algorithm 1 is w, the result follows, i.e.:
h

∗

∗ >

E (w − w )(w − w )

i

1
= n
GTb −1 (I − γTb )s E0 + O
( b − s)2 γ 2

25

n
b

− s, and noting that
!

n/b

ρb
γλmin (H)

n
b

 kE0 kF
−s

"
1
+
I−
n − sb
γ

#
1
−1
−1
 (HL + HR ) GTb −1 Σ
n
−
s
b
n

2

+
γb

A.2

n
b

2 HL
−s

−1

HR

−1

Tb

−1

Σ+O

−s+1

ρbb
kΣkF
(n − bs)λmin (H)λmin (Tb )

!

Proofs of some useful Lemmas

Lemma 14.
i
h
2
E kxk xx> 

2
· H.
γ1,max

Proof. The proof of the claim above is straightforward given that Ekxk2 x ⊗ x  R2 H and owing to the
definition of γ1,max = R22 .
Lemma 15. Denoting the assumption (A) γ ≤ γb,max /2,
1. With (A) in place, Tb −1 W  0 for every W ∈ S(d), W  0

2
2. With (A) in place, Tr Tb −1 W ≤ λmin
· Tr (W) ∀ W ∈ S(d), W  0

3
· kWk∗ for every W ∈ Rd×d
3. With (A) in place, HL −1 + HR −1 − γI W ∗ ≤ λmin


4. With (A) in place, Tr HL −1 + HR −1 − γI Tb −1 W ≤ λ26 · Tr (W) ∀ W ∈ S(d), W  0
min

5. (I − γTb ) W  0 ∀ W ∈ S(d), W  0
6. With (A) in place, Tr ((I − γTb ) W) ≤ (1 − γλmin /2) Tr (W) ∀ W ∈ S(d), W  0



−1
7. Tr (HR + HL ) A = 21 Tr H−1 A ∀ A ∈ S + (Rd )
8. With (A) in place,
 3

Tr Tb −1 Σ ≤ Tr H−1 Σ
2
Proof. Proof of claim 1 in Lemma 15: We require to prove Tb −1 operating on a PSD matrix produces a PSD
matrix, or in other words, Tb −1 is a PSD map.
Tb −1 = [HL + HR −

γ
(M + (b − 1)HL HR )]−1
b
1

1

= (HL + HR )− 2 (HL + HR ) 2 [HL + HR −
1

= (HL + HR )− 2 [I −

1
1
γ
(M + (b − 1)HL HR )]−1 (HL + HR ) 2 (HL + HR )− 2
b

1
1
1
γ
(HL + HR )− 2 (M + (b − 1)HL HR )(HL + HR )− 2 ]−1 (HL + HR )− 2
b
1

1

(19)

Now, we prove that k γb (HL + HR )− 2 (M + (b − 1)HL HR )(HL + HR )− 2 k < 1. Given γ < γb,max /2, we employ
lemma 8 to note that Tb  0.
Tb  0
⇒ HL + HR −

γ
(M + (b − 1)HL HR )  0
b

γ
(M + (b − 1)HL HR ) ≺ HL + HR
b
1
1
γ
⇒ (HL + HR )− 2 (M + (b − 1)HL HR )(HL + HR )− 2 ≺ I
b
1
1
γ
⇒ k (HL + HR )− 2 (M + (b − 1)HL HR )(HL + HR )− 2 k < 1
b
⇒

26

With this fact in place, we employ Taylor series to expand Tb −1 in equation 19, i.e.:
1

Tb −1 = (HL + HR )− 2
=

∞
X
1
1
1
γ
( (HL + HR )− 2 (M + (b − 1)HL HR )(HL + HR )− 2 )i (HL + HR )− 2
b
i=0

∞
X
γ
( (HL + HR )−1 (M + (b − 1)HL HR ))i (HL + HR )−1
b
i=0

The proof completes by employing the following facts: Using Lyapunov’s theorem, we know (HL + HR )−1
is a PSD map, i.e. if (HL + HR )−1 (A) = B, then, if A is PSD =⇒ B is PSD. Furthermore, M is also a PSD
map, i.e. if A1 is PSD, M(A1 ) = E[(xT A1 x)x ⊗ x] is PSD as well. Finally, HL HR is also a PSD map, since, if
A2 is PSD, then, HL HR (A2 ) = HA2 H which is PSD as well. With all these facts in place, we note that each
term in the Taylor’s expansion above is a PSD map implying the overall map is PSD as well, thus rounding up
the proof to claim 1 in Lemma 15.
Proof of claim 2 in Lemma 15: Since we know that Tb −1 is a PSD map, it suffices to prove the following
stronger statement where we replace the condition W  0 with Tb −1 W  0.

Tr Tb −1 W ≤

1
λmin

· Tr (W) ∀ W ∈ S(d), Tb −1 W  0.

This is equivalent to showing that
Tr (U) ≤
Since Tb U = HU + UH − γb MU −

1
· Tr (Tb U) ∀ U ∈ S(d), U  0.
λmin

γ(b−1)
HUH,
b

(20)

we bound the following quantities:

h
i
h
i
2
2
Tr (MU) = E kxk x> Ux = hE kxk xx> , Ui ≤

2
γ1,max

hH, Ui,

where the last step follows from Lemma 14 and the fact that U  0. We also have


Tr (HUH) = hH, H1/2 UH1/2 i ≤ kHk2 Tr H1/2 UH1/2 = kHk2 Tr (HU) .
Using (21) and (22), we have
γ(b − 1)
γ
Tr (MU) −
Tr (HUH)
b
b
γ kHk2 (b − 1)
γ
Tr (HU) −
Tr (HU)
≥ 2 Tr (HU) −
bγ1,max
b


γ kHk2 (b − 1)
γ
= 2−
−
Tr (HU) .
bγ1,max
b

Tr (Tb U) = 2 Tr (HU) −

Since

γ
bγ1,max

<

γb,max
2bγ1,max

≤ 12 , and γ kHk2 ≤ γb,max · kHk2 /2 ≤
Tr (Tb U) ≥

2
kHk2

· kHk2 /2 = 1, we have:

1
λmin
Tr (HU) ≥
· Tr (U) .
2
2

This proves the claim.
Proof of claim 3 in Lemma 15: Again using Lemma 16, we know that

HL −1 + HR −1 − γI W ∗ = H−1 W + WH−1 − γW ∗


3
2
−1
−1
≤ H W ∗ + WH
+ γ kWk∗ ≤
kWk∗ ,
+ γ kWk∗ ≤
∗
λmin
λmin
1
where we used the fact that γ < γb,max /2 ≤ kHk
.
2
Proof of claim 4 in Lemma 15: This follows easily from combining Claims 1, 2 and 3.



Tr HL −1 + HR −1 − γI Tb −1 W ≤ HL −1 + HR −1 − γI Tb −1 W ∗ ( from Lemma 16)

27

(21)

(22)

3
· Tb −1 W ∗ (from Claim 3)
λmin

3
=
· Tr Tb −1 W (from Claim 1)
λmin
6
≤ 2 · Tr (W) (from Claim 2).
λmin

≤

Proof of claim 5 in Lemma 15: Since (I − γTb ) W = (I − γH) W (I − γH) + γ 2 /b (MW − HWH), it suffices
to show that (MW − HWH)  0. So we need to show
(MW − HWH)  0
⇔MW  HWH
⇔u> (MW)u ≥ u> HWHu ∀ u ∈ Rd




⇔Ex (u> x)2 (x> Wx) ≥ Ex1 ,x2 (x1 > Wx2 )(u> x1 )(u> x2 ) ∀ u ∈ Rd .
This last statement is proved using Cauchy-Schwarz inequality as follows:


E (x1 > Wx2 )(u> x1 )(u> x2 )
h
i
≤ E W1/2 x1 W1/2 x2 u> x1 u> x2

1/2


≤ E (x1 > Wx1 )(x2 > Wx2 )(u> x1 )2 (u> x2 )2
= E (u> x)2 (x> Wx) .
This finishes the proof.
Proof of claim 6 in Lemma 15: Since in the proof of Claim 2 we proved (20) which says that
Tr (Tb U) ≥ λmin · Tr (U) ∀ U ∈ S(d), U  0,
Claim 6 follows.
Proof of claim 7 in Lemma 15:
−1
We know that the operator (HR + HL ) is a PSD map, i.e, it maps PSD matrices to PSD matrices. Since
−1
A  0, we replace this condition with U = (HR + HL ) A  0 implying, we need to show the following:
Tr (U) =


1
Tr H−1 U ∀ U  0
2

Examining the right hand side, we see the following:
 1

1
Tr H−1 U = Tr H−1 (HL + HR )U
2
2

1
= Tr H−1 HU + H−1 UH
2
= Tr (U)
thus wrapping up the proof of claim 8.
Proof of claim 8 in Lemma 15:
We first denote U = HL + HR − γ b−1
b HL HR . Then, we note that the taylor’s expansion holds since
γ < γb,max and hence Tb  0.

γ −1
Σ
Tb −1 Σ = U − M
b
γ
= U−1/2 (I − U−1/2 MU−1/2 )−1 U−1/2 Σ
b
∞
X
γ
= U−1 Σ +
( )i (U−1 M)i (U−1 Σ)
(23)
b
i=1
Before we proceed further, we provide some upper bounds on the spectral properties of the map U−1 as
applied to any psd matrix B  0:

−1
b−1
−1
U B = HL + HR − γ
HL HR
B
b
28


−1
b−1
−1/2
−1/2
(HL + HR )
(HL + HR )−1/2 B
= (HL + HR )
I −γ
HL HR (HL + HR )
b
i
∞ 
X
b−1
=
(HL + HR )−1 HL HR B̂
γ
b
i=0
−1/2

where, B̂ = (HL + HR )−1 B. Next,
(HL + HR )−1 HL HR B̂ = (HL + HR )−1 HB̂H
1
= B̂H
2
kHk2

B̂
2
implying,
i
∞ 
X
γ b−1
( )(
) kHk2 B̂
2
b
i=0
!
2b
1
B̂ , γ <
kHk2
=
b−1
b−1
1 − γ 2b kHk2

U−1 B 

Next, we establish some lower bounds on the spectral properties of the map U−1 when applied to any psd
matrix B  0. In particular, our main claim is:
(HL + HR )−1 B  U−1 B
In order to prove this claim, let us consider C1 = (HL + HR )−1 B, and C2 = U−1 B. Given that B  0, we
know that C1  0 (by Lyapunov’s Theorem), and C2  0 (using an argument similar to claim 1, with the step
2b
size γ < (b−1)kHk
chosen in order to have U  0 for the arguments similar to claim 1 to hold). Then, we
2
establish the claim using the following series of equations:
b−1
HL HR )C2
b
b−1
HC2 H
HC1 + C1 H = HC2 + C2 H − γ
b
b−1
γ
HC2 H = H(C2 − C1 ) + (C2 − C1 )H
b
b−1
γ
HC2 H = (HL + HR )(C2 − C1 )
b

B = (HL + HR )C1 = (HL + HR − γ

γ

b−1
(HL + HR )−1 HC2 H = C2 − C1
b

We know that the left hand side is psd (using Lyapunov’s theorem), and hence, C2  C1 which establishes
the claim. Now, we combine the resulting lower bounds with the upper bounds to obtain:
Tr(HL + HR )−1 B ≤ Tr U−1 B ≤ {
(HL + HR )−1 B
with the condition that γ <

2

≤ U−1 B

2b
(b−1)kHk2 .

2

1
1−

γ(b−1)kHk2
2b

1

≤{
1−

γ(b−1)kHk2
2b

} Tr(HL + HR )−1 B

} (HL + HR )−1 B

2

(24)
(25)

Next, denoting A = U−1 Σ, we have by equation 23:
Tb −1 Σ = A +

∞
X
(U−1 M)i A
i=1

29

(26)

With this in place, we resort to bounding the behavior of equation 26 in a psd sense, with the following
arguments:
U−1 MA = U−1 E(xt Ax)x ⊗ x
 kAk2 R2 U−1 H
= kAk2 R2 A0
where, A0 = U−1 H. We note that U−1 is a psd map (following arguments similar to claim 1 in lemma 15
2b
with the step size γ < (b−1)kHk
), which is necessary for the above sequence of inequalities to hold. This
2
implies,
(U−1 M)2 A = U−1 M(U−1 MA)
 kAk2 R2 U−1 MA0
 kAk2 R2 (kA0 k2 R2 )A0
implying, the ith term in equation 26 is bounded as:
(U −1 M)i A  kAk2 R2 (kA0 k2 R2 )i−1 A0
Using these bounds in equation 26 and summing up the resulting geometric series and applying trace on both
sides, we get the following bound:
Tr Tb −1 Σ ≤ Tr U−1 Σ +

γ
b

U−1 Σ
1−

γ
b

R2 Tr U−1 H

kU−1 Hk2 R2
2

,γ<

b
kU−1 Hk2 R2

(27)

To bound the behavior of the first term in equation 27, we employ the bound in equation 24 (which ensures
2b
statistically optimal rates), and is satisfied if γ < (b−1)kHk
. Next, the denominator of the second term
2

in equation 27 holds iff the geometric series converges, which is if γ <
equation 25 is upperbounded as γ <

b
R2 kU−1 Hk2 ,

which by employing

= R2b2 . Finally, for the numerator of the second term
b Tr (HL +HR )−1 Σ
γ < kU
−1 Σk R2 Tr U−1 H , which by employing the bounds
2
2b Tr(HL +HR )−1 Σ
k(HL +HR )−1 Σk2 R2 d . Combining these three bounds, the

b
R2 k(HL +HR )−1 Hk2

in equation 27 to be statistically optimal, we require

in equations 25 and 24, we obtain a bound as γ <
step size that suffices for the statistically optimal rates are:


2b 2b Tr (HL + HR )−1 Σ
2b
,
,
γ < min
(b − 1) kHk2 R2 k(HL + HR )−1 Σk2 R2 d
We also note that the third term is always ≤ the second term, implying,


2b
2b Tr (HL + HR )−1 Σ
γ < min
,
(b − 1) kHk2 k(HL + HR )−1 Σk2 R2 d

In particular, with the choice of step size that is half of this maximum, by employing claim 7, we have
Tr Tb −1 Σ ≤

3
Tr H−1 Σ
2

Lemma 16. Recall that the nuclear norm kUk∗ of a matrix U is defined to be the sum of its singular values. The
following properties are satisfied by the nuclear norm:
• kUAk∗ ≤ kAk2 kUk∗ for all U, A ∈ Rd×d .
• kUk∗ = Tr (U) for all PSD matrices U.
• kUk∗ = max A∈Rd×d hU, Ai for all U ∈ Rd×d . Moreover, the optimizer A∗ satisfies A∗ U = UA∗  0.
kAk2 ≤1

30

A.3

Proofs of Theorems 2 and 3

A.3.1

Proof of Theorem 2

Proof. Since γ ≤

γb,max
3 ,

we first invoke lemma 7 to get the contraction factor ρb :

ρb ≤ max(1 − 2γ(1 −

γ
γb,max

)λmin (H), −1 + 2γ(1 −

γ
γb,max

)λmax (H) + 2(

γ
γb,max

)2 , 1 − γλmin (H))

We now wish to get a finer bound on the contraction factor, for which, we first define α =
equation 28 we subtract the second term from the first, to yield:

γ
γb,max

(28)
in the

1 − 2γ(1 − α)λmin (H) − (−1 + 2γ(1 − α)λmax (H) + 2α2 )
= 2 − 2αγb,max (1 − α)(λmin (H) + λmax (H)) − 2α2
≥ 2 − 4αγb,max (1 − α)λmax (H) − 2α2
= 2(1 − 2α(1 − α)γb,max λmax (H) − α2 )
≥ 2(1 − 4α(1 − α) − α2 )
= 2(1 − 4α + 3α2 )
= 2(1 − α)(1 − 3α)
≥0
implies that α ≤ 1 and α ≤ 13 , which is indeed true from the statement of the theorem. Once this is indeed
true, we bound the contraction factor from equation 28 as:
ρb ≤ max(1 − 2γ(1 − α)λmin (H), 1 − γλmin (H))
By setting α = 61 , we get,

5
1
1
,1 −
)=1−
9κ
3κ
3κ
With a bound on the contraction factor established explicitly, we now turn towards providing a scalar bound
by taking the innerproduct of each term in equation 2 with H, which helps in establishing the statement of
theorem 2.
ρb ≤ max(1 −

hH, T1 i =
=
≤
=
≤
≤
≤
≤
≤

1
hH, GTb −1 (I − γTb )s E0 i
− s)2 γ 2
1
hTb −1 (I − γTb )s E0 , (HL −1 + HR −1 − γI)Hi
n
( b − s)2 γ 2
2
hTb −1 (I − γTb )s E0 , Ii
( nb − s)2 γ 2

2
Tr Tb −1 (I − γTb )s E0
( nb − s)2 γ 2
√
2 d
Tb −1 (I − γTb )s E0 F
( nb − s)2 γ 2
√
2 d
s
Tb −1 kI − γTb k kE0 kF
( nb − s)2 γ 2
√
s
2 d
e(− 3κ ) kTb −1 kkE0 kF
n
2
2
( b − s) γ
√
s
2 d
e(− 3κ ) kTb −1 k Tr (E0 )
n
2
2
( b − s) γ
√


s
2 d
1
(− 3κ
) kT −1 k
hH,
E
i
e
b
0
γ 2 ( nb − s)2
λmin (H)
( nb

31

√


s
1
4 d
(− 3κ
) kH−1 k
e
hH,
E
i
0
γ 2 ( nb − s)2
λmin (H)
√
s
4 d
= 2 n
e(− 3κ ) hH, E0 i
2
2
γ ( b − s) (λmin (H))
≤

which wraps up bounding hH, T1 i by plugging in the value of γ. We switch our attention to the second term.
hH, T2 i = kHk Tr (T2 )
√
≤ d kHk kT2 kF
√

n/b

ρb
 kE0 kF
γλmin (H) nb − s
√
3 exp(−n/3bκ)κ2 d
≤
hH, E0 i,
( nb − s)

≤

d kHk ·

which bounds the second term. Next, we proceed to bounding the third term:
"
#
1
1
−1
−1
 (HL + HR ) GTb −1 Σi
hH, I −
hH, T3 i =
n − sb
γ nb − s
"
#
1
1
−1
−1
−1
 (HL + HR ) Hi
=
hΣ, Tb G I −
n − sb
γ nb − s
1
hΣ, Tb −1 (HL −1 + HR −1 − γI)Hi
n − sb
2
hΣ, Tb −1 Ii
≤
n − sb
2
hI, Tb −1 Σi
=
n − sb
3 Tr H−1 Σ
≤
n − sb
≤

Next, we proceed to the fourth term:
2

−1
−1 −1
2 hH, HL HR Tb Σi
γb
−s
2
−1
−1
−1
=
2 hΣ, Tb HR HL Hi
γb nb − s
2
−1 −1
=
2 hΣ, Tb H i
n
γb b − s
2
−1
≤
2 hI, Tb Σi
γbλmin (H) nb − s

hH, T4 i =

≤

n
b

3 Tr H−1 Σ
n
b

γbλmin (H)

2
−s

which wraps up bounding hH, T4 i. We consider the final term:
hH, T5 i = kHk Tr (T5 )
√
≤ d kHk kT5 kF
≤
≤

n

−s+1

√

ρbb
d kHk ·
kΣkF
(n − bs)λmin (H)λmin (Tb )

√

2ρbb
d kHk ·
kΣkF
(n − bs)λmin (H)2

n

32

−s+1

(29)

n

−s+1
√
ρb
≤ 2 dκ2 · b
Tr H−1 Σ
(n − bs)

which wraps up bounding all terms to establish the proof of Theorem 2.
A.3.2

Proof of Theorem 3


n
− 1.
Proof. We first prove a recursive bound on the error R(we ) in epoch e using Theorem 2 for e ≤ log bt
2
−1
R
d
(H
+H
)
Σ
k
k
L
R
Since in epoch e, we use a minibatch size of be = 2e−1 b = 2e−1 kHk Tr(HL +HR )−1 Σ2 , number of initial iterations
2

def

s = t − 1 and total samples ne = be t, we obtain from Theorem 2 the following bound on the error of the
iterate we :
√

 2√
 s 
ne
36 dκ2
3κ d
∗
∗
{R(we−1 ) − R(w )} + exp −
R(we ) − R(w ) ≤ 
{R(we−1 ) − R(w∗ )}
2 exp −
ne
3κ
3b
κ
−s
ne
e
b
e
be − s
!
( nbee − s + 1)
2κ2 Tr H−1 Σ
9κ Tr H−1 Σ
3 Tr H−1 Σ
+ 
exp −
+
2 +
ne − be s
(ne − be s)
3κ
be nbee − s
√


√

3 dκ2 Tr H−1 Σ
t
≤ 54 dκ2 exp −
{R(we−1 ) − R(w∗ )} +
using be = 2e−1 b, ne = be t, s = t − 1
e−1
3κ
2 b
√ 2


t
3 dκ Tr H−1 Σ
≤ exp −
( since t ≥ 16κ log (dκ)) .
{R(we−1 ) − R(w∗ )} +
3κ log(dκ)
2e−1 b
By expanding out the recursion, we obtain:
√

e
X
te
1
3 dκ2 Tr H−1 Σ
∗
R(we ) − R(w ) ≤ exp −
· e−`
{R(w0 ) − R(w )} +
`−1
3κ log(dκ)
2 b
16
`=1
√


e
te
3σ 2 dκ2 Tr H−1 Σ X 1
= exp −
{R(w0 ) − R(w∗ )} +
3κ log(dκ)
2e−1 b
8e−`
`=1
√


te
3 dκ2 Tr H−1 Σ 8
∗
·
≤ exp −
{R(w0 ) − R(w )} +
3κ log(dκ)
2e−1 b
7
√


2
−1
te
16 dκ Tr H Σ
≤ exp −
{R(w0 ) − R(w∗ )} +
.
3κ log(dκ)
2e+1 b
∗



(30)

n
In the final step of Algorithm 2, we have a minibatch size of b = 2t
, number of samples n = n/2 and
t
number of initial iterations s = 2 . Plugging this into Theorem 2 again, we obtain the following bound on the
final error.
√
√




t
144κ2 d
t
6κ2 d
∗
∗
−
R(w) − R(w ) ≤
exp
−
{R(w
)
−
R(w
)}
+
exp
{R(we ) − R(w∗ )}
e
t2
6κ
t
3κ
!
√
( nb − s + 1)
3 Tr H−1 Σ 9κ Tr H−1 Σ 2κ2 d Tr H−1 Σ
+
+ 
exp −
2 +
3κ
n − bs
(n − bs)
b nb − s





√
√
t
4 Tr H−1 Σ
18κ
t
≤ 3 dκ2 exp −
{R(we ) − R(w∗ )} +
3+
+ 2κ2 d exp −
3κ
n
t
6κ


n
using b = , n = n/2, s = t/2

2t
√ 2
t
20 Tr H−1 Σ
≤ 3 dκ exp −
{R(we ) − R(w∗ )} +
( since t ≥ 16κ log (dκ)) .
(31)
3κ
n

33

n
Since in the final iteration, we have e = log bt
− 1, using (30), we have
√


te
16 dκ2 Tr H−1 Σ
∗
∗
R(we ) − R(w ) ≤ exp −
{R(w0 ) − R(w )} +
3κ log(dκ)
2e+1 b
√ 2
t

 3κ log(dκ)
2bt
16t dκ Tr H−1 Σ
{R(w0 ) − R(w∗ )} +
≤
.
n
n

Plugging the above inequality into (31), we obtain


√
t
20 Tr H−1 Σ
R(w) − R(w∗ ) ≤ 3 dκ2 exp −
{R(we ) − R(w∗ )} +
3κ
n
!
√
t

 
 3κ log(dκ)
√ 2
t
2bt
16t dκ2 Tr H−1 Σ
20 Tr H−1 Σ
∗
{R(w0 ) − R(w )} +
≤ 3 dκ exp −
+
3κ
n
n
n
t

 3κ log(dκ)
2bt
40 Tr H−1 Σ
≤
{R(w0 ) − R(w∗ )} +
( since t > 16κ log(dκ)) .
n
n
This proves the theorem.

A.4

Proof of Lemma 4

For this problem instance, we begin by noting that (HL + HR )−1 Σ is diagonal as well, with entries:
(
1/2
if i = 1
1
−1
−1
{(HL + HR ) Σ}ii = {H Σ}ii =
2
1/2(d − 1)
if i > 1
Let us consider the case with batch size b = 1. With the appropriate choice of step size γ that ensure
contracting operators, we require considering Tr[Tb −1 Σ] as in equation 29, which corresponds to bounding
the leading order term in the variance. We employ the taylor’s expansion (just as in claim 1 of lemma 15) to
expand the term of interest Tb −1 Σ:
Tb −1 Σ =

∞
X

i
γ(HL + HR )−1 M (HL + HR )−1 Σ

i=0

= (HL + HR )−1 Σ +

∞
X

i
γ(HL + HR )−1 M (HL + HR )−1 Σ

i=1
∞
X

⇒ Tr Tb −1 Σ = Tr(HL + HR )−1 Σ +

Tr

h

i
i
γ(HL + HR )−1 M (HL + HR )−1 Σ

i=1
∞
h
i
X
i
1
Tr γ(HL + HR )−1 M (HL + HR )−1 Σ
Tr Tb −1 Σ = Tr H−1 Σ +
2
i=1

We observe that the term corresponding to i = 0 works out regardless of the choice of stepsize γ; we then
switch our attention to the second term, i.e., the term corresponding to i = 1:

d+2
Tr γ(HL + HR )−1 M (HL + HR )−1 Σ =
· Tr Σ
4
We require that this term should be ≤ Tr(HL + HR )−1 Σ, implying,
γ<

4 Tr(HL + HR )−1 Σ
(d + 2) Tr Σ

For this example, we observe that this yields γ <

4
1 ,
(d+2)(1+ d
)

d
1 ,
(d+2)(1+ d
)

which clearly is off by a factor d compared to

the realizable case which requires γ <
thus indicating a clear separation between the step sizes
required by SGD for the realizable and agnostic noise cases.
34

