Title: PARALLEL STOCHASTIC GRADIENT DESCENT WITH SOUND COMBINERS

Abstract: Stochastic gradient descent (SGD) is a well-known method for regression and classification tasks. However, it is an inherently sequential algorithm -at each step, the processing of the current example depends on the parameters learned from the previous examples. Prior approaches to parallelizing SGD, such as HOG-WILD! and ALLREDUCE, do not honor these dependences across threads and thus can potentially suffer poor convergence rates and/or poor scalability. This paper proposes SYMSGD, a parallel SGD algorithm that retains the sequential semantics of SGD in expectation. Each thread in this approach learns a local model and a probabilistic model combiner that allows the local models to be combined to produce the same result as what a sequential SGD would have produced, in expectation. This SYMSGD approach is applicable to any linear learner whose update rule is linear. This paper evaluates SYMSGD's accuracy and performance on 9 datasets on a shared-memory machine shows up-to 13× speedup over our heavily optimized sequential baseline on 16 cores.

1. INTRODUCTION:

Stochastic Gradient Descent (SGD) is an effective method for many regression and classification tasks. It is a simple algorithm with few hyper-parameters and its convergence rates are well understood both theoretically and empirically. However, its performance scalability is severely limited by its inherently sequential computation. SGD iteratively processes its input dataset where the computation at each iteration depends on the model parameters learned from the previous iteration.

Current approaches for parallelizing SGD do not honor this inter-step dependence across threads. Each thread learns a local model independently and combine these models in ways that can break sequential behavior. For instance, threads in HOGWILD! Recht et al. (2011) racily update a shared global model without holding any locks. In parameter-server Li et al. (2014a) , each thread (or machine) periodically sends its model deltas to a server that applies them to a global model. In ALLRE-DUCE Agarwal et al. (2014) , threads periodically reach a barrier where they compute a weightedaverage of the local models. Although these asynchronous parallel approaches reach the optimal solution eventually, they can produce a model that is potentially different from what a sequential SGD would have produced after processing a certain number of examples. Our experiments indicate that this makes their convergence rate slower than sequential SGD in terms of total number of examples studied. Our experiments show that all these algorithms either do not scale or their accuracy on the same number of examples falls short of a sequential baseline.

To address this problem, this paper presents SYMSGD, a parallel SGD algorithm that seeks to retain its sequential semantics. The key idea is for each thread to generate a sound combiner that allows the local models to be combined into a model that is the same as the sequential model. This paper describes a method for generating sound combiners for a class of SGD algorithms in which the inter-step dependence is linear in the model parameters. This class includes linear regression, linear regression with L2 regularization, and polynomial regression. While logistic regression is not in this class, our experiments show that linear regression performs equally well in classification tasks as logistic regression for the datasets studied in this paper. Also, this approach works even if the SGD computation is non-linear on the input examples and other parameters such as the learning rate; only the dependence on model parameters has to be linear.

Generating sound combiners can be expensive. SYMSGD uses random projection techniques to reduce this overhead but still retaining sequential semantics in expectation. We call this approach probabilistically sound combiners. Even though SYMSGD is expected to produce the same answer as the sequential SGD, controlling the variance introduced by the random projection requires care -a large variance can result in reduced accuracy. This paper describes the factors that affect this variance and explores the ensuing design trade-offs.

The resulting algorithm is fast, scales well on multiple cores, and achieves the same accuracy as sequential SGD on sparse and dense datasets. When compared to our optimized sequential baseline, SYMSGD achieves a speedup of 3.5× to 13× on 16 cores, with the algorithm performing better for denser datasets. Moreover, the cost of computing combiners can be efficiently amortized in a multiclass regression as a single combiner is sufficient for all of the classes. Finally, SYMSGD (like ALLREDUCE) is deterministic, producing the same result for a given dataset, configuration, and random seed. Determinism greatly simplifies the task of debugging and optimizing learning.

2. SOUND AND PROBABILISTIC MODEL COMBINERS:

Stochastic gradient descent (SGD) is a robust method for finding the parameters of a model that minimize a given error function. Figure 1 shows an example of a (convex) error function over two dimensions x and y reaching the minimum at parameter w * . SGD starts from some, not necessarily optimal, parameter w g (as shown in Figure 1 ), and repeatedly modifies w by taking a step along the gradient of the error function for a randomly selected example at the current w. The magnitude of the step is called the learning rate and is usually denoted by α. The gradient computed from one example is not necessarily the true gradient at w. Nevertheless, SGD enjoys robust convergence behavior by moving along the "right" direction over a large number of steps. This is shown pictorially in Figure 1 , where SGD processes examples in dataset D 1 to reach w 1 from w g . Subsequently, SGD starts from w 1 and processes a different set D 2 to reach w h . There is a clear dependence between the processing of D 1 and the processing of D 2 -the latter starts from w 1 , which is only determined after processing D 1 . Our goal is to parallelize SGD despite this dependence.

State of the art parallelization techniques such as HOGWILD! and ALLREDUCE approach this problem by processing D 1 and D 2 starting from the same model w g (let us assume that there only two processors for now), and respectively reaching w 1 and w 2 . Then, they combine their local models into a global model, but do so in an ad-hoc manner. For instance, ALLREDUCE computes a weighted average of w 1 and w 2 , where the per-feature weights are chosen so as to prefer the processor that has larger update for that feature. This weighted average is depicted pictorially as w a . Similarly, in HOGWILD!, the two processors race to update the global model with their respective local model without any locking. (HOGWILD! performs this udpate after every example, thus the size of D 1 and D 2 is 1.) Both approaches do not necessarily reach w h , the model that a sequential SGD would have reached on D 1 and D 2 . While SGD is algorithmically robust to errors, such ad-hoc combinations can result in slower convergence or poor performance, as we demonstrate in Section 4.

Sound Combiner::

The goal of this paper is to soundly combine local models. Looking at Figure 1 , a sound combiner combines local models w 1 and w 2 , respectively generated from datasets D 1 and D 2 , into a global model w h that is guaranteed to be the same as the model achieved by the sequential SGD processing D 1 and then D 2 . In effect, a sound combiner allows us to parallelize the sequential computation without changing its semantics.

If we look at the second processor, it starts its computation at w g , while in a sequential execution it would have started at w 1 , the output of the first processor. To obtain sequential semantics, we need to "adjust" its computation from w g to w 1 . To do so, the second processor performs its computation starting from w g +∆w, where ∆w is an unknown symbolic vector. This allows the second processor to both compute a local model (resulting from the concrete part) and a sound combiner (resulting from the symbolic part) that accounts for changes in the initial state. Once both processors are done learning, second processor finds w h by setting ∆w to w 1 -w g where w 1 is computed by the first processor. This parallelization approach of SGD can be extended to multiple processors where all processor produce a local model and a combiner (except for the first processor) and the local models are combined sequentially using the combiners.

When the update to the model parameters is linear in a SGD computation, then the dependence on the unknown ∆w can be concisely represented by a combiner matrix, as formally described in Section 3. Many interesting machine learning algorithms, such as linear regression, linear regression with L2 regularization, and polynomial regression already have linear update to the model parameters (but not necessarily linear on the input example).

Probabilistically Sound Combiner: The main problem with generating a sound combiner is that the combiner matrix has as many rows and columns as the total number of features. Thus, it can be effectively generated only for datasets with modest number of features. Most interesting machine learning problems involve learning over tens of thousands to billions of features, for which maintaining a combiner matrix is clearly not feasible.

We solve this problem through dimensionality reduction. Johnson-Lindenstrauss (JL) lemma Johnson & Lindenstrauss (1984) allows us to project a set of vectors from a high-dimensional space to a random low-dimensional space while preserving distances. We use this property to reduce the size of the combiner matrix without losing the fidelity of the computation -our parallel algorithm produces the same result as the sequential SGD in expectation.

Of course, a randomized SGD algorithm that generates the exact result in expectation is only useful if the resulting variance is small enough to maintain accuracy and the rate of convergence. We observe that for the variance to be small, the combiner matrix should have small singular values. Interestingly, the combiner matrix resulting from SGD is dominated by the diagonal entries as the learning rate has to be small for effective learning. We use this property to perform the JL projection only after subtracting the identity matrix. Also, other factors that control the singular values are the learning rate, number of processors, and the frequency of combining local models. This paper explores this design space and demonstrates the feasibility of efficient parallelization of SGD that retains the convergence properties of sequential SGD while enjoying parallel scalability.

3. PARALLEL SYMSGD ALGORITHM:

Consider a training dataset (X n×f , y n×1 ), where f is the number of features, n is the number of examples in the dataset, the i th row of matrix X, X i , represents the features of the i th example, and y i is the dependent value (or label) of that example. A linear model seeks to find a w * = arg min

w∈R f n i=0 Q(X i • w, y i ) that minimizes an error function Q. For linear regression, Q(X i • w, y i ) = (X i • w -y i ) 2 . When (X i , y i

) is evident from the context, we will simply refer to the error function as Q i (w).

SGD iteratively finds w * by updating the current model w with a gradient of Q r (w) for a randomly selected example r. For the linear regression error function above, this amounts to the update

w i = w i-1 -α∇Q r (w i-1 ) = w i-1 -α(X r • w i-1 -y r )X T r (1)

Here, α is the learning rate that determines the magnitude of the update along the gradient. As it is clear from this equation, w i is dependent on w i-1 which creates a loop-carried dependence and consequently makes parallelization of SGD across iterations using naïve approaches impossible.

The complexity of SGD for each iteration is as follows. Assume that X r has z non-zeros. Therefore, the computation in Equation 1requires O(z) amount of time for the inner product computation, X r •w i-1 , and the same O(z) amount of time for scalar-vector multiplication, α(X r •w i-1 -y r )X T r . If the updates to the weight vector happen in-place meaning that w i and w i-1 share the same memory location, the computation in Equation 1takes O(z) amount of time.

3.1. SYMBOLIC STOCHASTIC GRADIENT DESCENT:

This section explains a new approach to parallelize the SGD algorithm despite its loop-carried dependence. As shown in Figure 1 , the basic idea is to start each processor (except the first) on a concrete model w along with a symbolic unknown ∆w that captures the fact that the starting model can change based on the output of the previous processor. If the dependence on ∆w is linear during an SGD update, which is indeed the case for linear regression, then the symbolic dependence on ∆w on the final output can be captured by an appropriate matrix M a→b that is a function of the input examples X a , . . . , X b processed (y a , . . . , y b do not affect this matrix). Specifically, as Lemma A.1 in the Appendix shows, this combiner matrix is given by

EQUATION

In effect, the combiner matrix above is the symbolic representation of how a ∆w change in the input will affect the output of a processor. M a→b is referred by M when the inputs are not evident.

The parallel SGD algorithm works as follows (see Figure 1 ). In the learning phase, each processor i starting from w 0 , computes both a local model l i and a combiner matrix M i . In a subsequent reduction phase, each processor in turn computes its true output using

EQUATION

Lemma A.1 ensures that this combination of local models will produce the same output as what these processors would have generated had they run sequentially. We call such combiners sound.

One can compute a sound model combiner for other SGD algorithms provided the loop-carried dependence on w is linear. In other words, there should exist a matrix A i and vector b i in iteration i such that w i = A i • w i-1 + b i . Note that A i and b i can be nonlinear in terms of input datasets.

3.2. DIMENSIONALITY REDUCTION OF A SOUND COMBINER:

The combiner matrix M generate above can be quite large and expensive to compute. The sequential SGD algorithm maintains and updates the weight vector w, and thus requires O(f ) space and time, where f is the number of features. In contrast, M is a f × f matrix and consequently, the space and time complexity of parallel SGD is O(f 2 ). In practice, this would mean that we would need O(f ) processors to see constant speedups, an infeasible proposition particularly for datasets that can have thousands if not millions of features.

SYMSGD resolves this issue by projecting M into a smaller space while maintaining its fidelity. This projection is inspired by the Johnson-Lindenstrauss (JL) lemma Johnson & Lindenstrauss (1984) and follows the treatment of Achlioptas Achlioptas (2001) . Lemma 3.1.1 Let A be a random f × k matrix with

a ij = d ij / √ k

where a ij is the element of A at the i th row and j th column and d ij is independently sampled from a random distribution D with I E[D] = 0 and Var

[D] = 1. Then I E[A • A T ] = I f ×f

The matrix A from Lemma 3.1 projects from R f → R k where k can be much smaller than f . This allows us to approximate Equation 3 as

EQUATION

Lemma 3.1 guarantees that the approximation above is unbiased.

I E[l i + M i • A • A T (w i-1 -w 0 )] = l i + M i • I E[A • A T ](w i-1 -w 0 ) = w i

This allows an efficient algorithm that only computes the projected version of the combiner matrix while still producing the same answer as the sequential algorithm in expectation. We call such combiners probabilistically sound.

Algorithm 1: SYMSGD learning a local model and a model combiner.

1 <vector,matrix,matrix> SymSGD( 2 float α, vector: w0, X1..Xn, 3 scalar: y1..yn) { 4 vector w = w0;

5 matrix A = 1 √ k random(D,f,k); 6 matrix MA = A; 7 for i in (1..n) { 8 w = w -α(Xi•w -yi)Xi T ; 9 MA = MA -α Xi•(Xi T MA); } 10

return <w,MA,A>; } Algorithm 2: SYMSGD combining local models using model combiners.

1 vector SymSGDCombine(vector w0, 2 vector w, vector l,

3 matrix MA, matrix A) { 4 parallel { 5 matrix NA = MA -A; 6 w = l+w-w0+NA•A T (w-w0); 7 } 8 return w; }

Algorithm 1 shows the resulting symbolic SGD learner. The random function in line 5 returns a f × k matrix with elements chosen independently from the random distribution D according to Lemma 3.1. When compared to the sequential SGD, the additional work is the computation of M A in Line 9. It is important to note that this algorithm maintains the invariant that M A = M • A at every step. This projection incurs a space and time overhead of O(z × k) where z is the number of non-zeros in X i . This overhead is acceptable for small k and infact in our experiments in Section 4, k is between 7 to 15 across all benchmarks. Most of the overhead for such a small k is hidden by utilizing SIMD hardware within a processor (SymSGD with one thread is only half as slow as the sequential SGD as discussed in Section 4.1). After learning a local model and a probabilistically sound combiner in each processor, Algorithm 2 combines the resulting local model using the combiners, but additionally employs the optimizations discussed in Section 3.3.

Note that the correctness and performance of SYMSGD do not depend on the sparsity of a dataset and as Section 4 demonstrates, it works for very sparse and completely dense datasets. Also, note that X 1 , . . . , X n may contain a subset of size f' of all f features. Our implementation of Algorithm 1 takes advantage of this property and allocates and initializes A for only the observed features. This optimization is omitted from the pseudo code in Algorithm 1 for the sake of simplicity.

3.3. CONTROLLING THE VARIANCE:

While the dimensionality reduction discussed above is expected to produce the right answer, this is useful only if the variance of the approximation is acceptably small. Computing the variance is involved and is discussed in the associated technical report SymSGDTR. But we discuss the main result that motivates the rest of the paper.

Consider the approximation of 

M •∆w with v = M •A•A T •∆w. Let C(v

∆w 2 2 k i σ 2 i (M ) ≤ tr(C(v)) ≤ ∆w 2 2 k ( i σ 2 i (M ) + σ 2 max (M ))

The covariance is small if k, the dimension of the projected space, is large. But increasing k proportionally increases the overhead of the parallel algorithm. Similarly, covariance is small if the projection happens on small ∆w. Looking at Equation 4, this means that w i-1 should be as close to w 0 as possible, implying that processors should communicate frequently enough such that their models are roughly in sync. Finally, the singular values of M should be as small as possible. The next section describes a crucial optimization that achieves this.

Taking Identity Off: Expanding Equation 2, we see that the combiner matrices are of the form

I -αR 1 + α 2 R 2 -α 3 R 3 + • • •

where R i matrices are formed from the sum of products of X j • X T j matrices. Since α is a small number, this sum is dominated by I. In fact, for a combiner matrix M generated from n examples, M -I has at most n non-zero singular values SymSGDTR. We use these observation to lower the variance of dimensionality reduction by projecting matrix N = M -I instead of M . Appendix A.3 empirically shows the impact of this optimization. Rewriting Equations 3 and 4, we have

EQUATION

Lemma 3.1 guarantees that the approximation above is unbiased. Algorithm 2 shows the pseudo code for the resulting probabilistically sound combination of local models. The function SymSGDCombine is called iteratively to combine the model of the first processor with the local models of the rest. Note that each model combination is executed in parallel (Line 4) by parallelizing the underlying linear algebra operations.

An important factor in controlling the singular values of N is the frequency of model combinations which is a tunable parameter in SYMSGD. As it is shown in Appendix A.3, the fewer the number of examples learned, the smaller the singular values of N and the less variance (error) in Equation 5.

Implementation For the implementation of SymSGD function, matrix M and weight vector w are stored next to each other. This enables better utilization of vector units in the processor and improves the performance of our approach significantly. Also, most of datasets are sparse and therefore, SGD and SymSGD only copy the observed features from w 0 to their learning model w. Moreover, for the implementation of matrix A, we used Achlioptas (2001) theorem to minimize the overhead of creating A. In this approach, each element of A is independently chosen from { 1 3 , -1 3 , 0} with probability { 1 6 , 1 6 , 2 3 }, respectively.

4. EVALUATION:

All experiments described in this section were performed on an Intel Xeon E5-2630 v3 machine clocked at 2.4 GHz with 256 GB of RAM. The machine has two sockets with 8 cores each, allowing us to study the scalability of the algorithms across sockets. We disabled hyper-threading and turbo boost. We also explicitly pinned threads to cores in a compact way which means that thread i + 1 was placed as close as possible to thread i. The machine runs Windows 10. All of our implementations were compiled with Intel C/C++ compiler 16.0 and relied heavily on OpenMP primitives for parallelization and MKL for efficient linear algebra computations. And, finally, to measure runtime, we use the average of five independent runs on an otherwise idle machine.

There are several algorithms and implementations that we used for our comparison: Vowpal Wabbit Langford et al. (2007) , a widely used public library, Baseline, a fast sequential implementation, HW-Paper, the implementation from Recht et al. ( 2011), HW-Release, an updated version, Hog-Wild, which runs Baseline in multiple threads without any synchronization, and ALLREDUCE, the implementation from Agarwal et al. (2014) . Each of these algorithms have different parameters and settings and we slightly modified to ensure a fair comparison; see Appendix A.4 for more details.

When studying the scalability of a parallel algorithm, it is important to compare the algorithms against an efficient baseline Bailey (1991); McSherry et al. (2015) . Otherwise, it is empirically not possible to differentiate between the scalability achieved from the parallelization of the inefficiencies and the scalability inherent in the algorithm. We spent a significant effort to implement a well-tuned sequential algorithm which we call Baseline in our comparisons. Baseline is between 1.97 to 7.62 (3.64 on average) times faster than Vowpal Wabbit and it is used for all speedup graphs in this paper. Table 1 : Datasets used for evaluation with their settings, maximum accuracies using logistic and linear regression and maximum speedup using SYMSGD and HOGWILD!. Red speedup numbers represent the only cases where HOGWILD! performs faster than SYMSGD.

Datasets:

Wabbit with the configurations discussed in Appendix A.4 to measure the maximum accuracies that can be achieved using linear and logistic regression and the result is presented in columns 8 and 9 of Table 1 . In the case of aloi dataset, even after 500 passes (the default for our evaluation was 100 passes) the accuracies did not saturate to the maximum possible and we reported that both linear and logistic achieved at least 80% accuracy. The last two columns show the maximum speedup of SYMSGD and HOGWILD! over the baseline. 

4.1. RESULTS:

Figure 2 shows the accuracy and speedup measurements on three benchmarks: rcv1.binary, a sparse binary dataset, rcv1.multiclass, a sparse multiclass dataset, and epsilon, a dense binary dataset. The results for the other six benchmarks are presented in Appendix A.5.

Sparse Binary, rcv1.binary: Figure 2a compares the scalability of all the algorithms studied in this paper. HW-Paper is around six times slower than HW-Release. While this could potentially be a result of us running HW-Release on a Ubuntu VM, our primary aim of this comparison was to ensure that HogWild is a competitive implementation of HOGWILD!. Thus, we remove HW-Paper and HW-Release in our subsequent comparisons.

SYMSGD is half as slow as the Baseline on one thread as it performs lot more computation, but scales to a 3.5× speedup to 16 cores. Note, this represents a roughly 7× strong-scaling speedup with respect to its own performance on one thread. Analysis of the hardware performance counters shows the current limit to SYMSGD's scalability arises from load-imbalance across barrier synchronization, which provides an opportunity for future improvement.

Figure 2d shows the accuracy as a function of the number of examples processed by different algorithms. SYMSGD "stutters" at the beginning, but it too matches the accuracy of Baseline. The initial stuttering happens because the magnitude of the local models on each processor are large during the first set of examples. This directly affect the variance of the combiner matrix approximation. However, as more examples are given to SYMSGD, the magnitude of the local models are smaller and thus SYMSGD better matches the Baseline accuracy. One way to avoid this stuttering is to combine models more frequently (lower variance) or running single threaded for the first few iterations.

HogWild does approach sequential accuracy, however, it does so at the cost of scalablity (i.e., see Figure 2a (a)). Likewise, ALLREDUCE scales slightly better but does so at the cost of accuracy.

Sparse Multiclass, rcv1.multiclass: Figure 2b shows the scalability on rcv1.multiclass. Since this is a multiclass dataset, SYMSGD is competitive with the baseline on one thread as it is able to amortize the combiner matrix computation across all of the classes (M is the same across different classes). Thus, it enjoys much better scalability of 7× when compared to rcv1.binary. HogWild scales similar to SYMSGD up-to 8 threads but suffers when 16 threads across multiple sockets are used. Figure 2e shows that SYMSGD meets the sequential accuracy after an initial stutter. ALLREDUCE suffers from accuracy.

Dense Binary, epsilon: Figure 2c in Appendix A.5 shows that SYMSGD achieves a 7× speedup over the baseline on 16 cores. This represents a 14× strong scaling speedup over SYMSGD on one thread. As HOGWILD! is not designed for dense workloads, its speedup suffers when 16 cores across multiple sockets are used. This shows that SYMSGD scales to both sparse and dense datasets.

Similarly, ALLREDUCE suffers from accuracy.

5. RELATED WORK:

Most schemes for parallelizing SGD learn local models independently and communicate to update the global model. The algorithms differ in how and how often the update is performed. These choices determine the applicability of the algorithm to shared-memory or distributed systems.

To the best of our knowledge, our approach is the only one that retain the semantics of the sequential SGD algorithm. While some prior work provides theoretical analysis of the convergence rates that justify a specific parallelization, convergence properties of SYMSGD simply follow from the sequential SGD algorithm. On the other hand, SYMSGD is currently restricted to class of SGD computations where the inter-step dependence is linear in the model parameters.

Given a tight coupling of the processing units, Langford et al. Langford et al. (2009) suggest on a round-robin scheme to update the global model allowing for some staleness. However, as the SGD computation per example is usually much smaller when compared to the locking overhead, HOGWILD! Recht et al. ( 2011) improves on this approach to perform the update in a "racy" manner. While HOGWILD! is theoretically proven to achieve good convergence rates provided the dataset is sparse enough and the processors update the global model fast enough, our experiments show that the generated cache-coherence traffic limits its scalability particularly across multiple sockets. Moreover, as HOGWILD! does not update the model atomically, it potentially loses correlation among more frequent features resulting in loss of accuracy. Lastly, unlike SYMSGD, which works for both sparse and dense datasets, HOGWILD! is expclitly designed for sparse data. Recently, Sallinen et al. (2016) proposed applying lock-free HOGWILD! approach to mini-batch. However, mini-batch converges slower than SGD and also they did not study multi-socket scaling.

Zinkevich et al. Zinkevich et al. (2010) propose a MapReduce-friendly framework for SGD. The basic idea is for each machine/thread to run a sequential SGD on its local data. At the end, the global model is obtained by averaging these local models. Alekh et al. Agarwal et al. (2014) extend this approach by using MPI AllReduce operation. Additionally, they use the adagrad Duchi et al. (2011) approach for the learning rates at each node and use weighted averaging to combine local models with processors that processed a feature more frequently having a larger weight. Our experiments on our datasets and implementation shows that it does not achieve the sequential accuracy.

Several distributed frameworks for machine learning are based on parameter server Li et al. (2014b; a) where clients perform local learning and periodically send the changes to a central parameter server that applies the changes. For additional parallelism, the models themselves can be split across multiple servers and clients only contact a subset of the servers to perform their updates.

6. CONCLUSION:

With terabytes of memory available on multicore machines today, our current implementation has the capability of learning from large datasets without incurring the communication overheads of a distributed system. That said, we believe the ideas in this paper apply to distributed SGD algorithms and how to pursue in future work.

Many machine learning SGD algorithms require a nonlinear dependence on the parameter models. While SYMSGD does not directly apply to such algorithms, it is an interesting open problem to devise linear approximations (say using Taylor expansion) to these problems and subsequently parallelize with probabilistically sound combiners. This is an interesting study for future work. 

M a→b = a i=b (I -αX T i • X i )

Proof. The proof follows from a simple induction. Starting from w s , let the models computed by SGD after processing (X a , y a ), (X a+1 , y a+1 ), . . . , (X b , y b ) respectively be w a , w a+1 , . . . w b . Consider the base case of processing of (X a , y a ). Starting from w s + ∆w, SGD computes the model w a using Equation 1 (reminder:

w i = w i-1 -α(X i • w i-1 -y i )X T i ): w a = w s + ∆w -α(X a • (w s + ∆w) -y a )X T a = w s + ∆w -α(X a • w s -y a )X T a -α(X a • ∆w)X T a = w s -α(X a • w s -y a )X T a + ∆w -α(X a • ∆w)X T a = w a + ∆w -α(X a • ∆w)X T a (6) = w a + ∆w -αX T a (X a • ∆w) (7) = w a + ∆w -α(X T a • X a ) • ∆w (8) = w a + (I -αX T a • X a ) • ∆w

Step 6 uses Equation 1, Step 7 uses the fact that X a • ∆w is a scalar (allowing it to be rearranged), and Step 8 follows from associativity property of matrix multiplication.

The induction is very similar and follows from replacing ∆w with M a→i-1 ∆w and the property that A.2 PROOF OF LEMMA 3.1 

EQUATION

Proof. Let's call B = A • A T . Then b ij , the element of B at row i and column j, is s a is a js . Therefore, I E[b ij ] = k s=1 I E[a is a js ] = ( 1 √ k ) 2 k s=1 I E[d is d js ] = 1 k k s=1 I E[d is d js ]. For i = j, I E[b ij ] = 1 k k s=1 I E[d is ]I E[d js ] because d ij are chosen independently. Since I E[D] = 0 and d is , d js ∈ D, I E[d is ] = I E[d js ] = 0 and consequently, I E[b ij ] = 0. For i = j, I E[b ii ] = 1 k s I E[d is d is ] = 1 k s I E[d 2 is ]. Since I E[D 2 ] = 1 and d is ∈ I E[d 2 is ] = 1. As a result, I E[b ii ] = 1 k k s=1 I E[d 2 is ] = 1 k k s=1 1 = 1. A.3 EMPIRICAL EVALUATING SINGULAR VALUES OF M

A.4 ALGORITHM DETAILS AND SETTINGS:

This section provides details of all algorithms we used in this paper. Each algorithm required slight modification to ensure fair comparison.

Vowpal Wabbit: Vowpal Wabbit Langford et al. ( 2007) is one of the widely used public libraries for machine learning algorithms. We used this application as a baseline for accuracy of different datasets and as a comparison of logistic and linear regression and also an independent validation of the learners without any of our implementation bias. Vowpal Wabbit applies accuracy-improving optimizations such as adaptive and individual learning steps or per feature normalize updates. While all of these optimizations are applicable to SYMSGD, we avoided them since the focus of this paper is the running time performance of our learner. The non-default flags that we used are: --sgd, --power t 0, --holdout off, --oaa nc for multiclass datasets where nc is the number of classes, --loss function func where func is squared or logistic. For learning rate, we searched for α, the learning rate, in the set of {.1, .5, .01, .05, .001, .005, . . . } and used --learning rate α. We went through dataset 100 times for each dataset (--passes 100) and saved the learned model after each pass (--save per pass). At the end, for linear and logistic regressions, we reported the maximum accuracies achieved among different passes and different learning rates.

Baseline: Baseline uses a mixture of MKL Intel and manually vectorized implementations of linear algebra primitives in order to deliver the fastest performance. Baseline processes up-to 3.20 billion features per second at 6.4 GFLOPS.

HOGWILD!: HOGWILD! Recht et al. ( 2011) is a lock-free approach to parallelize SGD where multiple thread apply Equation 1 simultaneously. Although this approach may have race condition when two threads process instances with a shared feature but the authors discuss that this does not hurt the accuracy significantly for sparse datasets. There are multiple implementations of this approach that we studied and evaluated in this section. Below is a description of each:

• HW-Paper: This is the implementation used to report the measurements in Recht et al. (2011) which is publicly available Hogwild. This code implements SVM algorithm. Therefore, we modified the update rule to linear regression. The modified code was compiled and run on our Windows machine described above using an Ubuntu VM since the code is configured for Linux systems.

• HW-Release: This is an optimized implementation that the authors built after the HOGWILD! paper Recht et al. ( 2011) was published. Similar to HW-Paper, we changed the update rule accordingly and executed it on the VM.

• HogWild: We implemented this version which runs Baseline by multiple threads without any synchronizations. This code runs natively on Windows and enjoys all the optimizations applied to our Baseline such as call to MKL library and manual vectorization of linear algebra primitives.

ALLREDUCE: ALLREDUCE Agarwal et al. ( 2014) is an approach where each thread makes a copy from the global model and applies the SGD update rule to the local model for certain number of instances. Along with the local model, another vector g is computed which indicates the confidence in an update for the weight of a feature in the local model. After the learning phase, the local weight vectors are averaged based on the confidence vectors from each thread. We implemented this approach similarly using MKL calls and manual vectorization.

A.5 SPEEUPS ON REMAINING BENCHMARKS 

:

See proof in Appendix A.2.

Figures and Tables:

Figure 1: Convex error function for a two-dimensional feature space.

) be the covariance matrix of v. The trace of the covariance matrix tr(C(v)) is the sum of the variance of individual elements of v. Let λ i (M ) by the ith eigenvalue of M and σ i (M ) = λ i (M T M ) the ith singular value of M . Let σ max (M ) be the maximum singular value of M . Then the following holds SymSGDTR:

Figure 2: Speedup and accuracy comparison

Figure 3: Distribution of singular values of M for rcv1 binary dataset for α = 0.01 and α = 0.001.

Figure3empirically demonstrates the benefit of taking identity off. This figure plots the singular values of M for the rcv1.binary dataset (described in Section 4) after processing64, 128, 192, 256   examples for two different learning rates. As it can be seen, the singular values are close to 1. However, the singular values of N = M -I are roughly the same as those of M minus 1 and consequently, are small. Finally, the smaller α, the closer the singular values of M are to 1 and the singular values of N are to 0. Also, note that the singular values of M decrease as the numbers of examples increase and therefore, the singular values of N increase. As a result, the more frequent the models are combined, the less variance (and error) is introduced into Equation 5.

Figure 4: Speedups on remaining benchmarks

C. Cortes, N. D. Lawrence, and K. Q. Weinberger (eds.), Advances in Neural Information Processing Systems 27, pp. 19-27. Curran Associates, Inc., 2014b. URL http: //papers.nips.cc/paper/5597-communication-efficient-distributedmachine-learning-with-the-parameter-server.pdf. If the SGD algorithm for linear regression processes examples (X a , y a ), (X a+1 , y a+1 ), . . . , (X b , y b ) starting from model w s to obtain w b , then its outcome starting on model w s + ∆w is given by w b + M a→b • ∆w where the combiner matrix M a→b is given by

Table content: <table><tr><td>Frank McSherry, Michael Isard, and Derek G. Murray. Scalability! but at what cost? In 15th</td></tr><tr><td>Workshop on Hot Topics in Operating Systems (HotOS XV), Kartause Ittingen, Switzerland, May</td></tr><tr><td>2015. USENIX Association. URL https://www.usenix.org/conference/hotos15/</td></tr><tr><td>workshop-program/presentation/mcsherry.</td></tr><tr><td>Benjamin Recht, Christopher Re, Stephen Wright, and Feng Niu. Hogwild: A lock-free approach to</td></tr><tr><td>parallelizing stochastic gradient descent. In Advances in Neural Information Processing Systems,</td></tr><tr><td>pp. 693-701, 2011.</td></tr><tr><td>S. Sallinen, N. Satish, M. Smelyanskiy, S. S. Sury, and C. R. High performance parallel stochastic</td></tr><tr><td>gradient descent in shared memory. In 2016 IEEE International Parallel and Distributed Pro-</td></tr><tr><td>cessing Symposium (IPDPS), pp. 873-882, May 2016. doi: 10.1109/IPDPS.2016.107.</td></tr><tr><td>SymSGDTR. Symsgd technical report. https://www.microsoft.com/en-us/research/</td></tr><tr><td>publication/parallelizing-sgd-symbolically/, 2016. [accessed Sep 2016].</td></tr><tr><td>Martin Zinkevich, Markus Weimer, Lihong Li, and Alex J Smola. Parallelized stochastic gradient</td></tr><tr><td>descent. In Advances in neural information processing systems, pp. 2595-2603, 2010.</td></tr><tr><td>A APPENDIX</td></tr><tr><td>A.1 COMBINER MATRIX</td></tr><tr><td>Lemma A.1.</td></tr></table>

Peer review (iclr2017):
This paper propose a parallel mechanism for stochastic gradient descent method (SGD) in case of gradient can be computed via linear operations (including least square linear regression and polynomial regression problems). The motivation is to recover the same effect compared with sequential SGD, by using a proposed sound combiner. To make such combiner more efficient, the authors also use a randomized projection matrix to do dimension reduction. Experiments shows the proposed method has better speedup than previous methods like Hogwild! and Allreduce. I feel that there might be some fundamental misunderstanding on SGD. **The combiner matrixM generate above can be quite large and expensive to compute. The sequential SGD algorithm maintains and updates the weight vector w , and thus requires O(f) space and time, where f is the number of features. In contrast,M is a f f matrix and consequently, the space and time complexity of parallel SGD is O(f^2) . In practice, this would mean that we would need O(f) processors to see constant speedups, an infeasible proposition particularly for datasets that can have thousands if not millions of features.* I do not think one needs O(f^2) space and complexity for updating M_i * v, where v is an f-dimensional vector. Note that M_i is a low rank matrix in the form of (I - a_i a_i*). The complexity and space can be reduced to O(f) if compute it by O(v - a_i (a_i* v)) equivalently. If M_i is defined in the form of the product of n number of rank 1 matrices. The complexity and space complexity is O(fn). In the context of this paper, n should be much smaller than f. I seriously doubt that all author*s assumptions, experiments, and strategies in this paper are based on this incorrect assumption on space and complexity of SGD. Why one can have speedup is unclear for me. It is unclear what computations are in parallel and why this sequential algorithms can bring speedup if M_i*v is computed in the most efficient way. I suggest authors to make the following changes to make this paper more clear and theoretically solid - provide computational complexity per step of the proposed algorithm - convergence rate analysis (convergence analysis is not enough): we would like to see how the dimension reduction can affect the complexity.