Hogwild!: A Lock-Free Approach to Parallelizing Stochastic
Gradient DescentFeng Niu, Benjamin Recht, Christopher Ré and Stephen J. WrightComputer Sciences Department, University of Wisconsin-Madison
1210 W Dayton St, Madison, WI 53706
June 2011
Abstract
Stochastic Gradient Descent (SGD) is a popular algorithm that can achieve state-of-the-art
performance on a variety of machine learning tasks. Several researchers have recently proposed schemes to parallelize SGD, but all require performance-destroying memory locking andsynchronization. This work aims to show using novel theoretical analysis, algorithms, and implementation that SGD can be implemented without any locking. We present an update scheme
called Hogwild! which allows processors access to shared memory with the possibility of overwriting each other’s work. We show that when the associated optimization problem is sparse,
meaning most gradient updates only modify small parts of the decision variable, then Hogwild!
achieves a nearly optimal rate of convergence. We demonstrate experimentally that Hogwild!
outperforms alternative schemes that use locking by an order of magnitude.

Keywords. Incremental gradient methods, Machine learning, Parallel computing, Multicore

1

Introduction

With its small memory footprint, robustness against noise, and rapid learning rates, Stochastic Gradient Descent (SGD) has proved to be well suited to data-intensive machine learning tasks [3,5,26].
However, SGD’s scalability is limited by its inherently sequential nature; it is difficult to parallelize. Nevertheless, the recent emergence of inexpensive multicore processors and mammoth,web-scale data sets has motivated researchers to develop several clever parallelization schemes for
SGD [4, 10, 12, 16, 30]. As many large data sets are currently pre-processed in a MapReduce-likeparallel-processing framework, much of the recent work on parallel SGD has focused naturally on
MapReduce implementations. MapReduce is a powerful tool developed at Google for extracting
information from huge logs (e.g., “find all the urls from a 100TB of Web data”) that was designed
to ensure fault tolerance and to simplify the maintenance and programming of large clusters ofmachines [9]. But MapReduce is not ideally suited for online, numerically intensive data analysis. Iterative computation is difficult to express in MapReduce, and the overhead to ensure fault
tolerance can result in dismal throughput. Indeed, even Google researchers themselves suggestthat other systems, for example Dremel, are more appropriate than MapReduce for data analysis
tasks [21].
For some data sets, the sheer size of the data dictates that one use a cluster of machines. However, there are a host of problems in which, after appropriate preprocessing, the data necessary for1statistical analysis may consist of a few terabytes or less. For such problems, one can use a single
inexpensive work station as opposed to a hundred thousand dollar cluster. Multicore systems havesignificant performance advantages, including (1) low latency and high throughput shared main
memory (a processor in such a system can write and read the shared physical memory at over
12GB/s with latency in the tens of nanoseconds); and (2) high bandwidth off multiple disks (athousand-dollar RAID can pump data into main memory at over 1GB/s). In contrast, a typical
MapReduce setup will read incoming data at rates less than tens of MB/s due to frequent checkpointing for fault tolerance. The high rates achievable by multicore systems move the bottlenecks
in parallel computation to synchronization (or locking) amongst the processors [2, 13]. Thus, toenable scalable data analysis on a multicore machine, any performant solution must minimize theoverhead of locking.
In this work, we propose a simple strategy for eliminating the overhead associated with locking:run SGD in parallel without locks, a strategy that we call Hogwild!. In Hogwild!, processors areallowed equal access to shared memory and are able to update individual components of memory
at will. Such a lock-free scheme might appear doomed to fail as processors could overwrite eachother’s progress. However, when the data access is sparse, meaning that individual SGD steps only
modify a small part of the decision variable, we show that memory overwrites are rare and thatthey introduce barely any error into the computation when they do occur. We demonstrate boththeoretically and experimentally a near linear speedup with the number of processors on commonlyoccurring sparse learning problems.
In Section 2, we formalize a notion of sparsity that is sufficient to guarantee such a speedup
and provide canonical examples of sparse machine learning problems in classification, collaborative
filtering, and graph cuts. Our notion of sparsity allows us to provide theoretical guarantees oflinear speedups in Section 4. As a by-product of our analysis, we also derive rates of convergence
for algorithms with constant stepsizes. We demonstrate that robust 1/k convergence rates are
possible with constant stepsize schemes that implement an exponential back-off in the √constantover time. This result is interesting in of itself and shows that one need not settle for 1/ k rates
to ensure robustness in SGD algorithms.
In practice, we find that computational performance of a lock-free procedure exceeds even ourtheoretical guarantees. We experimentally compare lock-free SGD to several recently proposed
methods. We show that all methods that propose memory locking are significantly slower thantheir respective lock-free counterparts on a variety of machine learning applications.

2

Sparse Separable Cost Functions

Our goal throughout is to minimize a function f : X ⊆ Rn → R of the formXf (x) =fe (xe ) .

(2.1)

e∈E

Here e denotes a small subset of {1, . . . , n} and xe denotes the values of the vector x on thecoordinates indexed by e. The key observation that underlies our lock-free approach is that thenatural cost functions associated with many machine learning problems of interest are sparse in thesense that |E| and n are both very large but each individual fe acts only on a very small numberof components of x. That is, each subvector xe contains just a few components of x.

2

(a)

(b)

(c)

Figure 1: Example graphs induced by cost function. (a) A sparse SVM induces a hypergraphwhere each hyperedge corresponds to one example. (b) A matrix completion example induces
a bipartite graph between the rows and columns with an edge between two nodes if an entry is
revealed. (c) The induced hypergraph in a graph-cut problem is simply the graph whose cuts
we aim to find.

The cost function (2.1) induces a hypergraph G = (V, E) whose nodes are the individual components of x. Each subvector xe induces an edge in the graph e ∈ E consisting of some subset ofnodes. A few examples illustrate this concept.
Sparse SVM. Suppose our goal is to fit a support vector machine to some data pairs E ={(z1 , y1 ), . . . , (z|E| , y|E| )} where z ∈ Rn and y is a label for each (z, y) ∈ E.
X
minimizex
max(1 − yα xT zα , 0) + λkxk22 ,
(2.2)
α∈E

and we know a priori that the examples zα are very sparse (see for example [14]). To write thiscost function in the form of (2.1), let eα denote the components which are non-zero in zα and let
du denote the number of training examples which are non-zero in component u (u = 1, 2, . . . , n).
Then we can rewrite (2.2) as
!
X
X x2
u
T
minimizex
max(1 − yα x zα , 0) + λ
.
(2.3)
du
u∈e
α∈E

α

Each term in the sum (2.3) depends only on the components of x indexed by the set eα .
Matrix Completion. In the matrix completion problem, we are provided entries of a low-rank,
nr × nc matrix Z from the index set E. Such problems arise in collaborative filtering, Euclidean
distance estimation, and clustering [8,17,24]. Our goal is to reconstruct Z from this sparse samplingof data. A popular heuristic recovers the estimate of Z as a product LR∗ of factors obtained fromthe following minimization:X
minimize(L,R)(Lu Rv∗ − Zuv )2 + µ2 kLk2F + µ2 kRk2F ,(2.4)(u,v)∈E

where L is nr × r, R is nc × r and Lu (resp. Rv ) denotes the uth (resp. vth) row of L (resp. R)[17, 24, 27]. To put this problem in sparse form, i.e., as (2.1), we write (2.4) asoX nminimize(L,R)(Lu Rv∗ − Zuv )2 + 2|Eµu− | kLu k2F + 2|Eµ−v | kRv k2F(u,v)∈E

where Eu− = {v : (u, v) ∈ E} and E−v = {u : (u, v) ∈ E}.
3

Graph Cuts. Problems involving minimum cuts in graphs frequently arise in machine learning(see [6] for a comprehensive survey). In such problems, we are given a sparse, nonnegative matrix W
which indexes similarity between entities. Our goal is to find a partition of the index set {1, . . . , n}that best conforms to this similarity matrix. Here the graph structure is explicitly determined by
the similarity matrix W ; arcs correspond to nonzero entries in W . We want to match each string
to some list of D entities. P
Each node is associated with a vector xi in the D-dimensional simplex
D
D
SD = {ζ ∈ R : ζv ≥ 0
v=1 ζv = 1}. Here, two-way cuts use D = 2, but multiway-cuts withtens of thousands of classes also arise in entity resolution problems [18]. For example, we may havea list of n strings, and Wuv might index the similarity of each string. Several authors (e.g., [7])propose to minimize the cost function
X
minimizex
wuv kxu − xv k1 subject to xv ∈ SD for v = 1, . . . , n .
(2.5)(u,v)∈E

In all three of the preceding examples, the number of components involved in a particular termfe is a small fraction of the total number of entries. We formalize this notion by defining thefollowing statistics of the hypergraph G:Ω := max |e|, ∆ :=e∈E

max1≤v≤n |{e ∈ E : v ∈ e}|maxe∈E |{ê ∈ E : ê ∩ e 6= ∅}|, ρ :=
.
|E|
|E|

(2.6)

The quantity Ω simply quantifies the size of the hyper edges. ρ determines the maximum fraction ofedges that intersect any given edge. ∆ determines the maximum fraction of edges that intersect any
variable. ρ is a measure of the sparsity of the hypergraph, while ∆ measures the node-regularity.
For our examples, we can make the following observations about ρ and ∆.
1. Sparse SVM. ∆ is simply the maximum frequency that any feature appears in an example,
while ρ measures how clustered the hypergraph is. If some features are very common acrossthe data set, then ρ will be close to one.
2. Matrix Completion. If we assume that the provided examples are sampled uniformly atr)r)
random and we see more than nc log(nc ) of them, then ∆ ≈ log(nand ρ ≈ 2 log(n. Thisnrnrfollows from a coupon collector argument [8].
3. Graph Cuts. ∆ is the maximum degree divided by |E|, and ρ is at most 2∆.
We now describe a simple protocol that achieves a linear speedup in the number of processors
when Ω, ∆, and ρ are relatively small.

3

The Hogwild! Algorithm

Here we discuss the parallel processing setup. We assume a shared memory model with p processors. The decision variable x is accessible to all processors. Each processor can read x, and cancontribute an update vector to x. The vector x is stored in shared memory, and we assume thatthe componentwise addition operation is atomic, that isxv ← xv + a4

Algorithm 1 Hogwild! update for individual processors
1: loop2:
Sample e uniformly at random from E
3:
Read current state xe and evaluate Ge (x)4:for v ∈ e do xv ← xv − γbTv Ge (x)5: end loop
can be performed atomically by any processor for a scalar a and v ∈ {1, . . . , n}. This operation
does not require a separate locking structure on most modern hardware: such an operation isa single atomic instruction on GPUs and DSPs, and it can be implemented via a compare-andexchange operation on a general purpose multicore processor like the Intel Nehalem. In contrast,
the operation of updating many components at once requires an auxiliary locking structure.
Each processor then follows the procedure in Algorithm 1. To fully describe the algorithm, letbv denote one of the standard basis elements in Rn , with v ranging from 1, . . . , n. That is, bv isequal to 1 on the vth component and 0 otherwise. Let Pv denote the Euclidean projection matrix
onto the vth coordinate, i.e., Pv = bv bTv . Pv is a diagonal matrix equal to 1 on the vth diagonal andzeros elsewhere. Let Ge (x) ∈ Rn denote a gradient or subgradient of the function fe multiplied by|E|. That is, we extend fe from a function on the coordinates of e to all of Rn simply by ignoring
the components in ¬e (i.e., not in e). Then
|E|−1 Ge (x) ∈ ∂fe (x).
Here, Ge is equal to zero on the components in ¬e. Using a sparse representation, we can calculate
Ge (x), only knowing the values of x in the components indexed by e. Note that as a consequenceof the uniform random sampling of e from E, we have
E[Ge (xe )] ∈ ∂f (x) .
In Algorithm 1, each processor samples an term e ∈ E uniformly at random, computes thegradient of fe at xe , and then writes
xv ← xv − γbTv Ge (x),for each v ∈ e(3.1)

Importantly, note that the processor modifies only the variables indexed by e, leaving all of thecomponents in ¬e (i.e., not in e) alone. We assume that the stepsize γ is a fixed constant. Even
though the processors have no knowledge as to whether any of the other processors have modified x,
we define xj to be the state of the decision variable x after j updates have been performed1 . Since
two processors can write to x at the same time, we need to be a bit careful with this definition, but
we simply break ties at random. Note that xj is generally updated with a stale gradient, which
is based on a value of x read many clock cycles earlier. We use xk(j) to denote the value of thedecision variable used to compute the gradient or subgradient that yields the state xj .
In what follows, we provide conditions under which this asynchronous, incremental gradientalgorithm converges. Moreover, we show that if the hypergraph induced by f is isotropic andsparse, then this algorithm converges in nearly the same number of gradient steps as its serial
counterpart. Since we are running in parallel and without locks, this means that we get a nearlylinear speedup in terms of the number of processors.
1

Our notation overloads subscripts of x. For clarity throughout, subscripts i, j, and k refer to iteration counts,
and v and e refer to components or subsets of components.

5

4

Fast Rates for Lock-Free Parallelism

We now turn to our theoretical analysis of Hogwild! protocols. To make the analysis tractable, we
assume that we update with the following “with replacement” procedure: each processor samples
an edge e uniformly at random and computes a subgradient of fe at the current value of the decisionvariable. Then it chooses an v ∈ e uniformly at random and updates
xv ← xv − γ|e|bTv Ge (x)Note that the stepsize is a factor |e| larger than the step in (3.1). Also note that this update iscompletely equivalent tox ← x − γ|e|PvT Ge (x) .
(4.1)
This notation will be more convenient for the subsequent analysis.
This with replacement scheme assumes that a gradient is computed and then only one of its
components is used to update the decision variable. Such a scheme is computationally wasteful as
the rest of the components of the gradient carry information for decreasing the cost. Consequently,
in practice and in our experiments, we perform a modification of this procedure. We partition out
the edges without replacement to all of the processors at the beginning of each epoch. The processorsthen perform full updates of all of the components of each edge in their respective queues. However,
we emphasize again that we do not implement any locking mechanisms on any of the variables. We
do not analyze this “without replacement” procedure because no one has achieved tractable analyses
for SGD in any without replacement sampling models. Indeed, to our knowledge, all analysis of
without-replacement sampling yields rates that are comparable to a standard subgradient descentalgorithm which takes steps along the full gradient of (2.1) (see, for example [22]). That is, theseanalyses suggest that without-replacement sampling should require a factor of |E| more steps than
with-replacement sampling. In practice, this worst case behavior is never observed. In fact, it isconventional wisdom in machine learning that without-replacement sampling in stochastic gradientdescent actually outperforms the with-replacement variants on which all of the analysis is based.
To state our theoretical results, we must describe several quantities that important in theanalysis of our parallel stochastic gradient descent scheme. We follow the notation and assumptionsof Nemirovski et al [23]. To simplify the analysis, we will assume that each fe in (2.1) is a convex
function. We assume Lipschitz continuous differentiability of f with Lipschitz constant L:k∇f (x0 ) − ∇f (x)k ≤ Lkx0 − xk, ∀ x0 , x ∈ X.

(4.2)

We also assume f is strongly convex with modulus c. By this we mean thatcf (x0 ) ≥ f (x) + (x0 − x)T ∇f (x) + kx0 − xk2 , for all x0 , x ∈ X.
2

(4.3)

When f is strongly convex, there exists a unique minimizer x? and we denote f? = f (x? ). We
additionally assume that there exists a constant M such thatkGe (xe )k2 ≤ M almost surely for all x ∈ X .

(4.4)

We assume throughout that γc < 1. (Indeed, when γc > 1, even the ordinary gradient descentalgorithms will diverge.)
Our main results are summarized by the following
6Proposition 4.1 Suppose in Algorithm 1 that the lag between when a gradient is computed andwhen it is used in step j — namely, j − k(j) — is always less than or equal to τ , and γ is defined
to be
ϑc
.
γ=
(4.5)
2
2LM Ω 1 + 6ρτ + 4τ 2 Ω∆1/2
for some  > 0 and ϑ ∈ (0, 1). Define D0 := kx0 − x? k2 and let k be an integer satisfying
2LM 2 Ω 1 + 6τ ρ + 6τ 2 Ω∆1/2 log(LD0 /)
k≥
.
c2 ϑ

(4.6)

Then after k component updates of x, we have E[f (xk ) − f? ] ≤ .
In the case that τ = 0, this reduces to precisely the rate achieved by the serial SGD protocol.
A similar rate is achieved if τ = o(n1/4 ) as ρ and ∆ are typically both o(1/n). In our setting, τis proportional to the number of processors, and hence as long as the number of processors is lessn1/4 , we get nearly the same recursion as in the linear rate.
We prove Proposition 4.1 in two steps in the Appendix. First, we demonstrate that the sequence
aj = 21 E[kxj − x? k2 ] satisfies a recursion of the form aj ≤ (1 − cr γ)(aj+1 − a∞ ) + a∞ for someconstant a∞ that depends on many of the algorithm parameters but not on the state, and someconstant cr < c. This cr is an “effective curvature” for the problem which is smaller that the true
curvature c because of the errors introduced by our update rule. Using the fact that cr γ < 1, we
will show in Section 5 how to determine an upper bound on k for which ak ≤ /L. Proposition 4.1
then follows because E[f (xk ) − f (x? )] ≤ Lak since the gradient of f is Lipschitz. A full proof isprovided in the appendix.
Note that up to the log(1/) term in (4.6), our analysis nearly provides a 1/k rate of convergence
for a constant stepsize SGD scheme, both in the serial and parallel cases. Moreover, note that ourrate of convergence is fairly robust to error in the value of c; we pay linearly for our underestimateof the curvature of f . In contrast, Nemirovski et al demonstrate that when the stepsize is inversely
proportional to the iteration counter, an overestimate of c can result in exponential slow-down [23]!
We now turn to demonstrating that we can eliminate the log term from (4.6) by a slightly more
complicated protocol where the stepsize is slowly decreased after a large number of iterations.

5

Robust 1/k rates.

Suppose we run Algorithm 1 for a fixed number of gradient updates K with stepsize γ < 1/c. Then,we wait for the threads to coalesce, reduce γ by a constant factor β ∈ (0, 1), and run for β −1 K
iterations. In some sense, this piecewise constant stepsize protocol approximates a 1/k diminishingstepsize. The main difference with the following analysis from previous work is that our stepsizes
are always less than 1/c in contrast to beginning with very large stepsizes. Always working withsmall stepsizes allows us to avoid the possible exponential slow-downs that occur with standard
diminishing stepsize schemes.
To be precise, suppose ak is any sequence of real numbers satisfying
ak+1 ≤ (1 − cr γ)(ak − a∞ (γ)) + a∞ (γ)
where a∞ is some non-negative function of γ satisfyinga∞ (γ) ≤ γB
7

(5.1)

and cr and B are constants. This recursion underlies many convergence proofs for SGD where ak
denotes the distance to the optimal solution after k iterations. We will derive appropriate constantsfor Hogwild! in the Appendix. We will also discuss below what these constants are for standard
stochastic gradient descent algorithms.
Factoring out the dependence on γ will be useful in what follows. Unwrapping (5.1) we haveak ≤ (1 − cr γ)k (a0 − a∞ (γ)) + a∞ (γ) .
Suppose we want this quantity to be less than . It is sufficient that both terms are less than /2.
For the second term, this means that it is sufficient to set
γ≤


.
2B

(5.2)

For the first term, we then need
(1 − γcr )k a0 ≤ /2
which holds ifk≥log(2a0 /)
.
γcr

(5.3)

ϑ
for ϑ ∈ (0, 1]. Combining this with (5.3) tells us that afterBy (5.2), we should pick γ = 2Bk≥

2B log(2a0 /)ϑcr

iterations we will have ak ≤ . This right off the bat almost gives us a 1/k rate, modulo the log(1/)
factor.
To eliminate the log factor, we can implement a backoff scheme where we reduce the stepsize by
a constant factor after several iterations. This backoff scheme will have two phases: the first phasewill consist of converging to the ball about x? of squared radius less than 2B
cr at an exponential
rate. Then we will converge to x? by shrinking the stepsize.
To calculate the number of iterates required to get inside a ball of squared radius 2B
cr , supposeϑthe initial stepsize is chosen as γ = cr (0 < ϑ < 1). This choice of stepsize guarantees thatthe ak converge to a∞ . We use the parameter ϑ to demonstrate that we do not suffer much forunderestimating the optimal stepsize (i.e., ϑ = 1) in our algorithms. Using (5.3) we find thata c 
0 r
k ≥ ϑ−1 log(5.4)
ϑB
iterations are sufficient to converge to this ball. Note that this is a linear rate of convergence.
Now assume that a0 < 2ϑB
cr . Let’s reduce the stepsize by a factor of β each epoch. This reduces
the achieved  by a factor of β. Thus, after logβ (a0 /) epochs, we will be at accuracy . The total
number of iterations required is then the sum of terms with the form (5.3), with a0 set to be theradius achieved by the previous epoch and  set to be β times this a0 . Hence, for epoch number ν,
the initial distance is β ν−1 a0 and the final radius is β ν . Summing over all of the epochs (except for8the initial phase) gives
logβ (a0 /)X
k=1

log(2/β)log(2/β)=ϑϑβ klogβ (a0 /)

X

β −k

k=1

log(2/β) β −1 (a0 /) − 1=
ϑ
β −1 − 1a0 log(2/β)≤ϑ 1 − β2B log(2/β)≤.
cr  1 − β(5.5)

This expression is minimized by selecting a backoff parameter ≈ 0.37. Also, note that when we
reduce the stepsize by β, we need to run for β −1 more iterations.
Combining (5.4) and (5.5), we estimate a total number of iterations equal tok ≥ ϑ−1 log

a c 
0 r

ϑB+

2B log(2/β)
cr  1 − β

are sufficient to guarantee that ak ≤ .
Rearranging terms, the following two expressions give  in terms of all of the algorithm parameters:
2 log(2/β) B
1
.
≤
(5.6)
·
·
0 cr
1−β
cr k − ϑ−1 log aϑB

5.1

Consequences for serial SGD

Let us compare the results of this constant step-size protocol to one where the stepsize at iteration k
is set to be γ0 /k for some initial step size γ for the standard (serial) incremental gradient algorithm
applied to (2.1). Nemirovski et al [23] show that the expected squared distance to the optimal
solution, ak , satisfies
ak+1 ≤ (1 − 2cγk )ak + 12 γk2 M 2 .
22γMWe can put this recursion in the form (5.1) by setting γk = γ, cr = 2c, B = M
4c , and a∞ = 4c .
ΘThe authors of [23] demonstrate that a large step size: γk = 2ck with Θ > 1 yields a bound
 2

1
M
Θ2
ak ≤ max·, D0
k
c2 4Θ − 4On the other hand, a constant step size protocol achievesak ≤

log(2/β) M 2
1

.
·
·
4(1 − β) c2 k − ϑ−1 log 4D0 c2
ϑM 2

This bound is obtained by plugging the algorithm parameters into (5.6) and letting D0 = 2a0 .
Note that both bounds have asymptotically the same dependence on M , c, and k. The expressionlog(2/β)4(1 − β)9is minimized when β ≈ 0.37 and is equal to 1.34. The expression
Θ24Θ − 4is minimized when Θ = 2 and is equal to 1 at this minimum. So the leading constant is slightly
worse in the constant stepsize protocol when all of the parameters are set optimally. However, if
D0 ≥ M 2 /c2 , the 1/k protocol has error proportional to D0 , but our constant stepsize protocol still
has only a logarithmic dependence on the initial distance. Moreover, the constant stepsize scheme
is much more robust to overestimates of the curvature parameter c. For the 1/k protocols, if one
overestimates the curvature (corresponding to a small value of Θ), one can get arbitrarily slow
rates of convergence. An simple, one dimensional example in [23] shows that Θ = 0.2 can yield aconvergence rate of k −1/5 . In our scheme, ϑ = 0.2 simply increases the number of iterations by afactor of 5.
The proposed fix in [23] for√ the sensitivity to curvature estimates results in asymptoticallyslower convergence rates of 1/ k. It is important to note that we need not settle for these slower
rates and can still achieve robust convergence at 1/k rates.

5.2

Parallel Implementation of a Backoff Scheme

The scheme described about results in a 1/k rate of convergence for Hogwild! with the only
synchronization overhead occurring at the end of each “round” or “epoch” of iteration. When
implementing a backoff scheme for Hogwild!, the processors have to agree on when to reduce thestepsize. One simple scheme for this is to run all of the processors for a fixed number of iterations,wait for all of the threads to complete, and then globally reduce the stepsize in a master thread. We
note that one can eliminate the need for the threads to coalesce by sending out-of-band messages
to the processors to signal when to reduce γ. This complicates the theoretical analysis as there
may be times when different processors are running with different stepsizes, but in practice could
allow one to avoid synchronization costs. We do not implement this scheme, and so do not analyze
this idea further.

6

Related Work

Most schemes for parallelizing stochastic gradient descent are variants of ideas presented in theseminal text by Bertsekas and Tsitsiklis [4]. For instance, in this text, they describe using stale
gradient updates computed across many computers in a master-worker setting and describe settings
where different processors control access to particular components of the decision variable. They
prove global convergence of these approaches, but do not provide rates of convergence (This is one
way in which our work extends this prior research). These authors also show that SGD convergenceis robust to a variety of models of delay in computation and communication in [29].
We also note that constant stepsize protocols with backoff procedures are canonical in SGD
practice, but perhaps not in theory. Some theoretical work which has at least demonstrated convergence of these protocols can be found in [20, 28]. These works do not establish the 1/k rates
which we provided above.
Recently, a variety of parallel schemes have been proposed in a variety of contexts. In MapReduce settings, Zinkevich et al proposed running many instances of stochastic gradient descent on
10

type
SVM
MC
Cuts

data
set
RCV1
Netflix
KDD
Jumbo
DBLifeAbdomen

size
(GB)
0.9
1.5
3.9
30
3e-3
18

ρ

∆

0.44
2.5e-33.0e-32.6e-78.6e-39.2e-41.0
2.3e-31.8e-3
1.4e-7
4.3e-39.2e-4

Hogwild!
time
train
test
(s)
error error
9.5 0.297 0.339
301.0 0.754 0.928
877.5
19.5
22.6
9453.5 0.031 0.013
230.0
10.6 N/A
1181.43.99 N/A

Round Robin
time
train
test
(s)
error error
61.8 0.297 0.339
2569.1 0.754 0.927
7139.019.5
22.6
N/A N/A N/A
413.5
10.5 N/A
7467.25
3.99 N/A

Figure 2: Comparison of wall clock time across of Hogwild! and RR. Each algorithm is run
for 20 epochs and parallelized over 10 cores.

different machines and averaging their output [30]. Though the authors claim this method can
reduce both the variance of their estimate and the overall bias, we show in our experiments thatfor the sorts of problems we are concerned with, this method does not outperform a serial scheme.
Schemes involving the averaging of gradients via a distributed protocol have also been proposedby several authors [10, 12]. While these methods do achieve linear speedups, they are difficult
to implement efficiently on multicore machines as they require massive communication overhead.
Distributed averaging of gradients requires message passing between the cores, and the cores needto synchronize frequently in order to compute reasonable gradient averages.
The work most closely related to our own is a round-robin scheme proposed by Langford etal [16]. In this scheme, the processors are ordered and each update the decision variable in order.
When the time required to lock memory for writing is dwarfed by the gradient computation time,
this method results in a linear speedup, as the errors induced by the lag in the gradients are not
too severe. However, we note that in many applications of interest in machine learning, gradientcomputation time is incredibly fast, and we now demonstrate that in a variety of applications,
Hogwild! outperforms such a round-robin approach by an order of magnitude.

7

Experiments

We ran numerical experiments on a variety of machine learning tasks, and compared against around-robin approach proposed in [16] and implemented in Vowpal Wabbit [15]. We refer to thisapproach as RR. To be as fair as possible to prior art, we hand coded RR to be nearly identicalto the Hogwild! approach, with the only difference being the schedule for how the gradients are
updated. One notable change in RR from the Vowpal Wabbit software release is that we optimized
RR’s locking and signaling mechanisms to use spinlocks and busy waits (there is no need for genericsignaling to implement round robin). We verified that this optimization results in nearly an orderof magnitude increase in wall clock time for all problems that we discuss.
We also compare against a model which we call AIG which can be seen as a middle groundbetween RR and Hogwild!. AIG runs a protocol identical to Hogwild! except that it locks allof the variables in e in before and after the for loop on line 4 of Algorithm 1. Our experiments
demonstrate that even this fine-grained locking induces undesirable slow-downs.
All of the experiments were coded in C++ are run on an identical configuration: a dual Xeon
X650 CPUs (6 cores each x 2 hyperthreading) machine with 24GB of RAM and a software RAID-0

11

3

4

SpeedupSpeedup4

5

Hogwild
AIG
RR

2
1
0
0

3

10

Hogwild
AIG
RR

8

Speedup

5

2
1

(a)
2

4

6

8

Number of Splits

10

0
0

6

Hogwild
AIG
RR

4
2

(b)
2

4

6

8

Number of Splits

10

0
0

(c)
2

4

6

8

10

Number of Splits

Figure 3: Total CPU time versus number of threads for (a) RCV1, (b) Abdomen, and (c)DBLife.

over 7 2TB Seagate Constellation 7200RPM disks. The kernel is Linux 2.6.18-128. We never use
more than 2GB of memory. All training data is stored on a seven-disk raid 0. We implemented acustom file scanner to demonstrate the speed of reading data sets of disk into small shared memory.
This allows us to read data from the raid at a rate of nearly 1GB/s.
All of the experiments use a constant stepsize γ which is diminished by a factor β at the end ofeach pass over the training set. We run all experiments for 20 such passes, even though less epochs
are often sufficient for convergence. We show results for the largest value of the learning rate γ
which converges and we use β = 0.9 throughout. We note that the results look the same across alarge range of (γ, β) pairs and that all three parallelization schemes achieve train and test errors
within a few percent of one another. We present experiments on the classes of problems describedin Section 2.
Sparse SVM. We tested our sparse SVM implementation on the Reuters RCV1 data set on
the binary text classification task CCAT [19]. There are 804,414 examples split into 23,149 training
and 781,265 test examples, and there are 47,236 features. We swapped the training set and thetest set for our experiments to demonstrate the scalability of the parallel multicore algorithms.
In this example, ρ = 0.44 and ∆ = 1.0—large values that suggest a bad case for Hogwild!.
Nevertheless, in Figure 3(a), we see that Hogwild! is able to achieve a factor of 3 speedup withwhile RR gets worse as more threads are added. Indeed, for fast gradients, RR is worse than aserial implementation.
For this data set, we also implemented the approach in [30] which runs multiple SGD runs inparallel and averages their output. In Figure 5(b), we display at the train error of the ensemble
average across parallel threads at the end of each pass over the data. We note that the threads
only communicate at the very end of the computation, but we want to demonstrate the effect ofparallelization on train error. Each of the parallel threads touches every data example in each pass.
Thus, the 10 thread run does 10x more gradient computations than the serial version. Here, theerror is the same whether we run in serial or with ten instances. We conclude that on this problem,
there is no advantage to running in parallel with this averaging scheme.
Matrix Completion. We ran Hogwild! on three very large matrix completion problems.
The Netflix Prize data set has 17,770 rows, 480,189 columns, and 100,198,805 revealed entries. The
KDD Cup 2011 (task 2) data set has 624,961 rows, 1,000,990, columns and 252,800,275 revealed
entries. We also synthesized a low-rank matrix with rank 10, 1e7 rows and columns, and 2e9
revealed entries. We refer to this instance as “Jumbo.” In this synthetic example, ρ and ∆ areboth around 1e-7. These values contrast sharply with the real data sets where ρ and ∆ are both
12

4

5

SpeedupSpeedup5

6

Hogwild
AIG
RR

3
2

4

8

Hogwild
AIG
RR

6

Speedup

6

3
2

4

2
1
0
0

Hogwild
AIG

1
2

4

6

8

Number of Splits

10

0
0

2

4

6

8

Number of Splits

10

0
0

2

4

6

8

10

Number of Splits

Figure 4: Total CPU time versus number of threads for the matrix completion problems (a)Netflix Prize, (b) KDD Cup 2011, and (c) the synthetic Jumbo experiment.

on the order of 1e-3.
Figure 5(a) shows the speedups for these three data sets using Hogwild!. Note that theJumbo and KDD examples do not fit in our allotted memory, but even when reading data off disk,
Hogwild! attains a near linear speedup. The Jumbo problem takes just over two and a half
hours to complete. Speedup graphs comparing Hogwild! to AIG and RR on the three matrix
completion experiments are provided in Figure 4. Similar to the other experiments with quickly
computable gradients, RR does not show any improvement over a serial approach. In fact, with
10 threads, RR is 12% slower than serial on KDD Cup and 62% slower on Netflix. In fact, it is
too slow to complete the Jumbo experiment in any reasonable amount of time, while the 10-way
parallel Hogwild! implementation solves this problem in under three hours.
Graph Cuts. Our first cut problem was a standard image segmentation by graph cuts problem
popular in computer vision. We computed a two-way cut of the abdomen data set [1]. This data
set consists of a volumetric scan of a human abdomen, and the goal is to segment the image into
organs. The image has 512 × 512 × 551 voxels, and the associated graph is 6-connected withmaximum capacity 10. Both ρ and ∆ are equal to 9.2e-4 We see that Hogwild! speeds up the cut
problem by more than a factor of 4 with 10 threads, while RR is twice as slow as the serial version.
Our second graph cut problem sought a mulit-way cut to determine entity recognition in alarge database of web data. We created a data set of clean entity lists from the DBLife website
and of entity mentions from the DBLife Web Crawl [11]. The data set consists of 18,167 entities
and 180,110 mentions and similarities given by string similarity. In this problem each stochastic
gradient step must compute a Euclidean projection onto a simplex of dimension 18,167. As aresult, the individual stochastic gradient steps are quite slow. Nonetheless, the problem is still
very sparse with ρ=8.6e-3 and ∆=4.2e-3. Consequently, in Figure 3, we see the that Hogwild!
achieves a ninefold speedup with 10 cores. Since the gradients are slow, RR is able to achieve aparallel speedup for this problem, however the speedup with ten processors is only by a factor of 5.
That is, even in this case where the gradient computations are very slow, Hogwild! outperforms
a round-robin scheme.
What if the gradients are slow? As we saw with the DBLIFE data set, the RR method
does get a nearly linear speedup when the gradient computation is slow. This raises the question
whether RR ever outperforms Hogwild! for slow gradients. To answer this question, we ran the
RCV1 experiment again and introduced an artificial delay at the end of each gradient computation
to simulate a slow gradient. In Figure 5(c), we plot the wall clock time required to solve the SVM
problem as we vary the delay for both the RR and Hogwild! approaches.
13

0.34

10

1 Thread
3 Threads
10 Threads

0.335

Train Error

Speedup6

Jumbo
Netflix
KDD

4

2

0.33
0.325
0.32

0
0

2

4

6

8

Number of Splits

10

0.31
0

6

Hogwild
AIG
RR

4
2

0.315

(a)

8

Speedup

8

(b)
5

10

Epoch

15

20

0 0
10

(c)
2

10

4

10

6

10

Gradient Delay (ns)Figure 5: (a) Speedup for the three matrix completion problems with Hogwild!. In all threecases, massive speedup is achieved via parallelism. (b) The training error at the end of eachepoch of SVM training on RCV1 for the averaging algorithm [30]. (c) Speedup achieved overserial method for various levels of delays (measured in nanoseconds).

Notice that Hogwild! achieves a greater decrease in computation time across the board. The
speedups for both methods are the same when the delay is few milliseconds. That is, if a gradienttakes longer than one millisecond to compute, RR is on par with Hogwild! (but not better). At
this rate, one is only able to compute about a million stochastic gradients per hour, so the gradientcomputations must be very labor intensive in order for the RR method to be competitive.

8

Conclusions

Our proposed Hogwild! algorithm takes advantage of sparsity in machine learning problems
to enable near linear speedups on a variety of applications. Empirically, our implementations
outperform our theoretical analysis. For instance, ρ is quite large in the RCV1 SVM problem, yet
we still obtain significant speedups. Moreover, our algorithms allow parallel speedup even when
the gradients are computationally intensive.
Our Hogwild! schemes can be generalized to problems where some of the variables occur quite
frequently as well. We could choose to not update certain variables that would be in particularlyhigh contention. For instance, we might want to add a bias term to our Support Vector Machine,
and we could still run a Hogwild! scheme, updating the bias only every thousand iterations orso.
For future work, it would be of interest to enumerate structures that allow for parallel gradientcomputations with no collisions at all. That is, it may be possible to bias the SGD iterations tocompletely avoid memory contention between processors. For example, recent work proposed abiased ordering of the stochastic gradients in matrix completion problems that completely avoids
memory contention between processors [25]. An investigation into how to generalize this approach toother structures and problems would enable even faster computation of machine learning problems.

Acknowledgements
BR is generously supported by ONR award N00014-11-1-0723 and NSF award CCF-1139953. CRis generously supported by the Air Force Research Laboratory (AFRL) under prime contract no.
FA8750-09-C-0181, the NSF CAREER award under IIS-1054009, ONR award N000141210041, andgifts or research awards from Google, LogicBlox, and Johnson Controls, Inc. SJW is generously
14supported by NSF awards DMS-0914524 and DMS-0906818 and DOE award DE-SC0002283. Any
opinions, findings, and conclusion or recommendations expressed in this work are those of theauthors and do not necessarily reflect the views of any of the above sponsors including DARPA,
AFRL, or the US government.

References[1] Max-flow problem instances in vision. From http://vision.csd.uwo.ca/data/maxflow/.
[2] K. Asanovic and et al. The landscape of parallel computing research: A view from berkeley. Technical
Report UCB/EECS-2006-183, Electrical Engineering and Computer Sciences, University of California
at Berkeley, 2006.
[3] D. P. Bertsekas. Nonlinear Programming. Athena Scientific, Belmont, MA, 2nd edition, 1999.
[4] D. P. Bertsekas and J. N. Tsitsiklis. Parallel and Distributed Computation: Numerical Methods. Athena
Scientific, Belmont, MA, 1997.
[5] L. Bottou and O. Bousquet. The tradeoffs of large scale learning. In Advances in Neural Information
Processing Systems, 2008.
[6] Y. Boykov and V. Kolmogorov. An experimental comparison of min-cut/max-flow algorithms for energy
minimization in vision. IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(9):1124–
1137, 2004.
[7] G. Călinescu, H. Karloff, and Y. Rabani. An improved approximation algorithm for multiway cut. In
Proceedings of the thirtieth annual ACM Symposium on Theory of Computing, pages 48–52, 1998.
[8] E. Candès and B. Recht. Exact matrix completion via convex optimization. Foundations of Computational Mathematics, 9(6):717–772, 2009.
[9] J. Dean and S. Ghemawat. MapReduce: simplified data processing on large clusters. Communicationsof the ACM, 51(1):107–113, 2008.
[10] O. Dekel, R. Gilad-Bachrach, O. Shamir, and L. Xiao. Optimal distributed online prediction usingmini-batches. Technical report, Microsoft Research, 2011.
[11] A. Doan. http://dblife.cs.wisc.edu.
[12] J. Duchi, A. Agarwal, and M. J. Wainwright. Distributed dual averaging in networks. In Advances in
Neural Information Processing Systems, 2010.
[13] S. H. Fuller and L. I. Millett, editors. The Future of Computing Performance: Game Over or Next
Level. Committee on Sustaining Growth in Computing Performance. The National Academies Press,
Washington, D.C., 2011.
[14] T. Joachims. Training linear svms in linear time. In Proceedings of the ACM Conference on Knowledge
Discovery and Data Mining (KDD), 2006.
[15] J. Langford. https://github.com/JohnLangford/vowpal_wabbit/wiki.
[16] J. Langford, A. J. Smola, and M. Zinkevich. Slow learners are fast. In Advances in Neural Information
Processing Systems, 2009.
[17] J. Lee, , B. Recht, N. Srebro, R. R. Salakhutdinov, and J. A. Tropp. Practical large-scale optimization
for max-norm regularization. In Advances in Neural Information Processing Systems, 2010.
[18] T. Lee, Z. Wang, H. Wang, and S. Hwang. Web scale entity resolution using relational evidence.
Technical report, Microsoft Research, 2011. Available at http://research.microsoft.com/apps/
pubs/default.aspx?id=145839.

15[19] D. Lewis, Y. Yang, T. Rose, and F. Li. RCV1: A new benchmark collection for text categorization
research. Journal of Machine Learning Research, 5:361–397, 2004.
[20] Z. Q. Luo and P. Tseng. Analysis of an approximate gradient projection method with applications tothe backpropagation algorithm. Optimization Methods and Software, 4:85–101, 1994.
[21] S. Melnik, A. Gubarev, J. J. Long, G. Romer, S. Shivakumar, M. Tolton, and T. Vassilakis. Dremel:Interactive analysis of web-scale datasets. In Proceedings of VLDB, 2010.
[22] A. Nedic and D. P. Bertsekas. Convergence rate of incremental subgradient algorithms. In S. Uryasev
and P. M. Pardalos, editors, Stochastic Optimization: Algorithms and Applications, pages 263–304.
Kluwer Academic Publishers, 2000.
[23] A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro. Robust stochastic approximation approach tostochastic programming. SIAM Journal on Optimization, 19(4):1574–1609, 2009.
[24] B. Recht, M. Fazel, and P. Parrilo. Guaranteed minimum rank solutions of matrix equations via nuclear
norm minimization. SIAM Review, 52(3):471–501, 2010.
[25] B. Recht and C. Ré. Parallel stochastic gradient algorithms for large-scale matrix completion. Submitted
for publication. Preprint available at http://pages.cs.wisc.edu/~brecht/publications.html, 2011.
[26] S. Shalev-Shwartz and N. Srebro. SVM Optimization: Inverse dependence on training set size. In
Proceedings of the 25th Internation Conference on Machine Learning, 2008.
[27] N. Srebro, J. Rennie, and T. Jaakkola. Maximum margin matrix factorization. In Advances in Neural
Information Processing Systems, 2004.
[28] P. Tseng. An incremental gradient(-projection) method with momentum term and adaptive stepsize
rule. SIAM Joural on Optimization, 8(2):506–531, 1998.
[29] J. Tsitsiklis, D. P. Bertsekas, and M. Athans. Distributed asynchronous deterministic and stochastic
gradient optimization algorithms. IEEE Transactions on Automatic Control, 31(9):803–812, 1986.
[30] M. Zinkevich, M. Weimer, A. Smola, and L. Li. Parallelized stochastic gradient descent. Advances in
Neural Information Processing Systems, 2010.

A

Analysis of Hogwild!

It follows by rearrangement of (4.3) that1(x − x0 )T ∇f (x) ≥ f (x) − f (x0 ) + ckx − x0 k2 , for all x ∈ X.
2

(A.1)

In particular, by setting x0 = x? (the minimizer) we have1(x − x? )T ∇f (x) ≥ ckx − x? k2 , for all x ∈ X.
2
We will make use of these identities frequently in what follows.

16

(A.2)

A.1

Developing the recursion

For 1 ≤ v ≤ n, let Pv denote the projection onto the vth component of Rn . (Pv is a diagonal matrix
with a 1 in the vth diagonal position and zeros elsewhere.) Similarly, for e ∈ E, we let Pe denotethe projection onto the components indexed by e — a diagonal matrix with ones in the positionsindexed by e and zeros elsewhere.
We start with the update formula (4.1). Recall that k(j) is the state of the decision variable’s
counter when the update to xj was read. We havexj+1 = xj − γ|ej |Pvj Gej (xk(j) ) .
By subtracting x? from both sides, taking norms, we have11
kxj+1 − x? k22 = kxj − x? k22 − γ|ej |(xj − x? )T Pvj Gej (xk(j) )2
2
1+ γ 2 |ej |2 kPvj Gej (xk(j) )k2
2
1= kxj − x? k22 − γ|ej |(xj − xk(j) )T Pvj Gej (xj )2− γ|ej |(xj − xk(j) )T Pvj (Gej (xk(j) ) − Gej (xj ))1− γ|ej |(xk(j) − x? )T Pvj Gej (xk(j) ) + γ 2 |ej |2 kPvj Gej (xk(j) )k2
2
Let aj = 21 E[kxj − x? k22 ]. By taking expectations of both sides and using the bound (4.4), we obtain
aj+1 ≤ aj − γE[(xj − xk(j) )T Gej (xj )] − γE[(xj − xk(j) )T (Gej (xk(j) ) − Gej (xj ))]1− γE[(xk(j) − x? )T Gej (xk(j) )] + γ 2 ΩM 2 .
2

(A.3)

where we recall that Ω = maxe∈E |e|. Here, in several places, we used the useful identity: for any
function ζ of x1 , . . . , xj and any i ≤ j, we have E[|ej |Pvj Gej (xi )T ζ(x1 , . . . , xj )] = E E |ej |Pvj Gej (xi )T ζ(x1 , . . . , xj ) | e1 , . . . , ej , v1 , . . . , vj−1
= E Gej (xi )T ζ(x1 , . . . , xj ) .
We will perform many similar calculations throughout, so, before proceeding, we denotee[i] := (e1 , e2 , . . . , ei , v1 , v2 , . . . , vi ),to be the tuple of all edges and vertices selected in updates 1 through i. Note that xl depends
on e[l−1] but not on ej or vj for any j ≥ l. We next consider the three expectation terms in thisexpression.
Let’s first bound the third expectation term in (A.3). Since xk(j) is independent of ej we haveE (xk(j) − x? )T Gej (xk(j) ) = E E (xk(j) − x? )T Gej (xk(j) )|e[k(j)−1]= E (xk(j) − x? )T E Gej (xk(j) )|e[k(j)−1]= E (xk(j) − x? )T ∇f (xk(j) ) ,17

where ∇f denotes an element of ∂f . It follows from (A.2) thatE (xk(j) − x? )T ∇f (xk(j) ) ≥ cak(j) .

(A.4)

The first expectation can be treated similarly:E[(xj − xk(j) )T Gej (xj )] = E E[(xj − xk(j) )T Gej (xj ) | e[j−1] ]= E (xj − xk(j) )T E[Gej (xj ) | e[j−1] ]= E[(xj − xk(j) )T ∇f (xj )]c≥ E[f (xj ) − f (xk(j) )] + E[kxj − xk(j) k2 ]
2

(A.5)

where the final inequality is from (A.1). Moreover, we can estimate the difference between f (xj )and f (xk(j) ) asE[f (xk(j) ) − f (xj )] =j−1XE[f (xi ) − f (xi+1 )]i=k(j)=j−1 X
X

E[fe (xi ) − fe (xi+1 )]i=k(j) e∈E
j−1

γ X X
≤
E[Ge (xi )T Gei (xi )]|E|
i=k(j) e∈E≤ γτ ρM 2 .

(A.6)

Here we use the inequalityfe (xi ) − fe (xi+1 ) ≤1γGe (xi )T (xi − xi+1 ) =Ge (xi )T Gei (xi )|E||E|

which follows because fe is convex. By combining (A.5) and (A.6), we obtainc
E[(xj − xk(j) )T Gej (xj )] ≥ −γτ ρM 2 + E[kxj − xk(j) k2 ] .
2

18

(A.7)

We turn now to the second expectation term in (A.3). We haveE (xj − xk(j) )T (Gej (xk(j) ) − Gej (xj ))


j−1
X
(xi+1 − xi )T (Gej (xk(j) ) − Gej (xj ))
= E
i=k(j)


= E

j−1
X


γ|ei |Gei (xk(i) )T (Gej (xk(j) ) − Gej (xj ))

i=k(j)





j−1
 X

= E
γ|ei |Gei (xk(i) )T (Gej (xk(j) ) − Gej (xj ))

i=k(j)ei ∩ej 6=∅≥ −E 



j−1
X


γ|ei |kGei (xk(i) )k kGej (xk(j) ) − Gej (xj )k

i=k(j)ei ∩ej 6=∅





j−1
 X

≥ −E 
2ΩM 2 γ 

i=k(j)ei ∩ej 6=∅≥ −2ΩM 2 γρτ(A.8)

where ρ is defined by (2.6). Here, the third line follows from our definition of the gradient update.
The fourth line is tautological: only the edges where ei and ej intersect nontrivially factor into thesum. The subsequent inequality is Cauchy-Schwarz, and the following line follows from (4.4).
By substituting (A.4), (A.7), and (A.8) into (A.3), we obtain the following bound: M 2γ2(Ω + 2τ ρ + 4Ωρτ ) .
aj+1 ≤ aj − cγ ak(j) + 12 E[kxj − xk(j) k2 ] +
2

(A.9)

To complete the argument, we need to bound the remaining expectation in (A.9). We expand
out the expression multiplied by cγ in (A.9) to findak(j) + 12 E[kxj − xk(j) k2 = aj − E (xj − xk(j) )T (xk(j) − x? )


j−1
X
= aj − E (xi+1 − xi )T (xk(j) − x? )
i=k(j)

= aj − E 

j−1
X


γ|ei |Gei (xk(i) )T Pvi (xk(j) − x? ) .

i=k(j)

Let e[¬i] denote the set of all sampled edges and vertices except for ei and vi . Since ei and vi

19are both independent of xk(j) , we can proceed to bound


j−1
X
γ|ei |Gei (xk(i) )T Pvi (xk(j) − x? )
E
i=k(j)


≤E 

j−1
X


γΩM kPvi (xk(j) − x? )k2 

i=k(j)

=γΩM

j−1
X




E Eei ,vi kPei ,vi (xk(j) − x? )k2 | e[¬i]i=k(j)

≤γΩM

j−1
X

h


1/2 iEei ,vi (xk(j) − x? )T Pvi (xk(j) − x? ) | e[¬i]h(xk(j) − x? )T Ee,v [Pv ](xk(j) − x? )E

i=k(j)

=γΩM

j−1
X

E

1/2 ii=k(j)≤τ γΩM Eh1/2 i(xk(j) − x? )T Ee,v [Pv ](xk(j) − x? )≤τ γΩM ∆1/2 E[kxk(j) − x? k2 ]≤τ γΩM ∆1/2 (E[kxj − x? k2 ] + τ γΩM )√ 1/2≤τ γΩM ∆1/2 ( 2aj + τ γΩM ) ,where ∆ is defined in (2.6). The first inequality is Cauchy-Schwartz. The next inequality is Jensen.
The second to last inequality follows from our definition of xj , and the final inequality is Jensen
again.
Plugging the last two expressions into (A.9), we obtain

√
1
1/2
2cΩM τ ∆1/2 aj + M 2 γ 2 Q
aj+1 ≤ (1 − cγ)aj + γ 2(A.10)2
whereQ = Ω + 2τ ρ + 4Ωρτ + 2τ 2 Ω2 ∆1/2 .
Here we use the fact that cγ < 1 to get a simplified form for Q. This recursion only involves
constants involved with the structure of f , and the nonnegative sequence aj . To complete theanalysis, we will perform a linearization to put this recursion in a more manageable form.
To find the steady state, we must solve the equation
√

M 2γ2
1/2
2cΩM τ ∆1/2 a∞a∞ = (1 − cγ)a∞ + γ 2+
Q.
(A.11)
2This yields the fixed point
s
!2
Q
M 2γ2
a∞ =
Ωτ ∆1/2 + Ω2 τ 2 ∆+
2
cγ

2
p
M 2γ
M 2γ
1/2
2
2
≤
Ωτ ∆ + Ω τ ∆ + Q = C(τ, ρ, ∆, Ω)
.
2c
2c
20

(A.12)

Note that for ρ and ∆ sufficiently small, C(τ, ρ, ∆, Ω) ≈ 1.
Since the square root is concave, we can linearize (A.10) about the fixed point a∞ to yield√γ22cΩM τ ∆1/2
(aj − a∞ )aj+1 ≤ (1 − cγ)(aj − a∞ ) +
√
2 a∞
√

1 2 2
+ (1 − cγ)a∞ + γ 22cΩM τ ∆1/2 a1/2
∞ + M γ Q
2
= (1 − cγ(1 − δ)) (aj − a∞ ) + a∞ .
Here
δ=

1
1
q
q
≤
Q
1 + 1 + cγΩ2 τ 2 ∆1 + 1 + Ω2Qτ 2∆

To summarize, we have shown that the sequence aj of squared distances satisfies
aj+1 ≤ (1 − cγ(1 − δ(τ, ρ, ∆, Ω)))(aj − a∞ ) + a∞

(A.13)

2

with a∞ ≤ C(τ, ρ, ∆, Ω) M2cγ . In the case that τ = 0 (the serial case), C(τ, ρ, ∆, Ω) = Ω and√δ(τ, ρ, ∆, Ω) = 0. Note that if τ is non-zero, but ρ and ∆ are o(1/n) and o(1/ n) respectively,
then as long as τ = o(n1/4 ), C(τ, ρ, ∆, Ω) = O(1). In our setting, τ is proportional to the numberof processors, and hence as long as the number of processors is less n1/4 , we get nearly the same
recursion as in the linear rate.
In the next section, we show that (A.13) is sufficient to yield a 1/k convergence rate. Since we
can run p times faster in our parallel setting, we get a linear speedup.

A.2

Proof of Proposition 4.1: Final Steps

Since ∇f is Lipschitz, we havef (x) ≤ f (x0 ) + ∇f (x0 )T (x − x0 ) +Lkx − x0 k2
2

for all x, x0 ∈ Rn .

Setting x0 = x? gives f (x) − f (x? ) ≤ L2 kx − x? k2 . Hence,E[f (xk ) − f (x? )] ≤ Lak
for all k. To ensure the left hand side is less than , it suffices to guarantee that ak ≤ /L.
To complete the proof of Proposition 4.1, we use the results of Section 4. We wish to achieve atarget accuracy of /L. To apply (5.3), choose a∞ as in (A.12) and the values
cr = c(1 − δ(τ, ρ, ∆, Ω)
B = C(τ, ρ, ∆, Ω)

M 2γ
.
2c

By (A.12), we have a∞ ≤ γB.
Choose γ satisfying (4.5). With this choice, we automatically have
γ≤

c


2 ≤
q
2LB
LM 2 Ω2 τ 2 ∆ 1 + 1 + Ω2Qτ 2∆21because (1 +√1 + x)2 ≤ 4 + 2x for all x ≥ 0. Substituting this value of γ into (5.3), we see thatk≥

C(τ, ρ, ∆, Ω)
LM 2 log(LD0 /)·2c
1 − δ(τ, ρ, ∆, Ω)

iterations suffice to achieve ak ≤ /L. Now observe thatC(τ, ρ, ∆, Ω)= Ω2 τ 2 ∆ ·1 − δ(τ, ρ, ∆, Ω)

2

q
Q
1 + 1 + Ω2 τ 2 ∆
1−

q 1
1+ 1+

Q
Ω2 τ 2 ∆2
qQ
1 + 1 + Ω2 τ 2 ∆2 2q=Ω τ ∆·
1 + Ω2Qτ 2∆≤ 8Ω2 τ 2 ∆ + 2Q
= 8Ω2 τ 2 ∆ + 2Ω + 4τ ρ + 8Ωρτ + 4τ 2 Ω2 ∆1/2
≤ 2Ω(1 + 6τ ρ + 6τ 2 Ω∆1/2 ) .
Here,the second to last inequality follows because
3√
1+ 1+x√≤ 8 + 2x
1+x
for all x ≥ 0. Plugging this bound into (A.14) completes the proof.

22(A.14)