Deep learning with Elastic Averaging SGD

Sixin Zhang, Anna Choromanska, Yann LeCun

Introduction

One of the most challenging problems in large-scale machine learning is how to parallelize the training of large models that use a form of stochastic gradient descent (SGD) . There have been attempts to parallelize SGD-based training for large-scale deep learning models on large number of CPUs, including the Google’s Distbelief system . But practical image recognition systems consist of large-scale convolutional neural networks trained on few GPU cards sitting in a single computer . The main challenge is to devise parallel SGD algorithms to train large-scale deep learning models that yield a significant speedup when run on multiple GPU cards.

In this paper we introduce the Elastic Averaging SGD method (EASGD) and its variants. EASGD is motivated by quadratic penalty method , but is re-interpreted as a parallelized extension of the averaging SGD algorithm . The basic idea is to let each worker maintain its own local parameter, and the communication and coordination of work among the local workers is based on an elastic force which links the parameters they compute with a center variable stored by the master. The center variable is updated as a moving average where the average is taken in time and also in space over the parameters computed by local workers. The main contribution of this paper is a new algorithm that provides fast convergent minimization while outperforming DOWNPOUR method and other baseline approaches in practice. Simultaneously it reduces the communication overhead between the master and the local workers while at the same time it maintains high-quality performance measured by the test error. The new algorithm applies to deep learning settings such as parallelized training of convolutional neural networks.

The article is organized as follows. Section 2 explains the problem setting, Section 3 presents the synchronous EASGD algorithm and its asynchronous and momentum-based variants, Section 4 provides stability analysis of EASGD and ADMM in the round-robin scheme, Section 5 shows experimental results and Section 6 concludes. The Supplement contains additional material including additional theoretical analysis.

Problem setting

EASGD update rule

The communication period τ\tau controls the frequency of the communication between every local worker and the master, and thus the trade-off between exploration and exploitation.

2 Momentum EASGD

The momentum EASGD (EAMSGD) is a variant of our Algorithm 1 and is captured in Algorithm 2. It is based on the Nesterov’s momentum scheme , where the update of the local worker of the form captured in Equation 3 is replaced by the following update

where δ\delta is the momentum term. Note that when δ=0\delta=0 we recover the original EASGD algorithm.

As we are interested in reducing the communication overhead in the parallel computing environment where the parameter vector is very large, we will be exploring in the experimental section the asynchronous EASGD algorithm and its momentum-based variant in the relatively large τ\tau regime (less frequent communication).

Stability analysis of EASGD and ADMM in the round-robin scheme

In this section we study the stability of the asynchronous EASGD and ADMM methods in the round-robin scheme . We first state the updates of both algorithms in this setting, and then we study their stability. We will show that in the one-dimensional quadratic case, ADMM algorithm can exhibit chaotic behavior, leading to exponential divergence. The analytic condition for the ADMM algorithm to be stable is still unknown, while for the EASGD algorithm it is very simpleThis condition resembles the stability condition for the synchronous EASGD algorithm (Condition 23 for p=1p=1) in the analysis in the Supplement..

The analysis of the synchronous EASGD algorithm, including its convergence rate, and its averaging property, in the quadratic and strongly convex case, is deferred to the Supplement.

In our setting, the ADMM method involves solving the following minimax problemThe convergence analysis in is based on the assumption that “At any master iteration, updates from the workers have the same probability of arriving at the master.”, which is not satisfied in the round-robin scheme.,

where λi\lambda^{i}’s are the Lagrangian multipliers. The resulting updates of the ADMM algorithm in the round-robin scheme are given next. Let t0t\geq 0 be a global clock. At each tt, we linearize the function F(xi)F(x^{i}) with F(xti)+\scalprodL2F(xti)xixti+12η\normL2xixti2F(x^{i}_{t})+\scalprod{L2}{\nabla{F}(x^{i}_{t})}{x^{i}-x^{i}_{t}}+\frac{1}{2\eta}\norm{L2}{x^{i}-x^{i}_{t}}^{2} as in . The updates become

The EASGD algorithm in the round-robin scheme is defined similarly and is given below

At time tt, only the ii-th local worker (whose index i1i-1 equals tt modulo pp) is activated, and performs the update in Equations 18 which is followed by the master update given in Equation 19.

For each of the pp linear maps, it’s possible to find a simple condition such that each map, where the ithi^{\text{th}} map has the form F3iF2iF1iF^{i}_{3}\circ F^{i}_{2}\circ F^{i}_{1}, is stable (the absolute value of the eigenvalues of the map are smaller or equal to one). However, when these non-symmetric maps are composed one after another as follows F=F3pF2pF1pF31F21F11\mathcal{F}=F^{p}_{3}\circ F^{p}_{2}\circ F^{p}_{1}\circ\ldots\circ F^{1}_{3}\circ F^{1}_{2}\circ F^{1}_{1}, the resulting map F\mathcal{F} can become unstable! (more precisely, some eigenvalues of the map can sit outside the unit circle in the complex plane).

We now present the numerical conditions for which the ADMM algorithm becomes unstable in the round-robin scheme for p=3p=3 and p=8p=8, by computing the largest absolute eigenvalue of the map F\mathcal{F}. Figure 1 summarizes the obtained result.

For the composite map FpF1F^{p}\circ\ldots\circ F^{1} to be stable, the condition that needs to be satisfied is actually the same for each ii, and is furthermore independent of pp (since each linear map FiF^{i} is symmetric). It essentially involves the stability of the 2×22\times 2 matrix \left(\begin{array}[]{cc}1-\eta-\alpha&\alpha\\ \alpha&1-\alpha\end{array}\right), whose two (real) eigenvalues λ\lambda satisfy (1ηαλ)(1αλ)=α2(1-\eta-\alpha-\lambda)(1-\alpha-\lambda)=\alpha^{2}. The resulting stability condition (λ1|\lambda|\leq 1) is simple and given as 0η2,0α42η4η0\leq\eta\leq 2,0\leq\alpha\leq\frac{4-2\eta}{4-\eta}.

Experiments

In this section we compare the performance of EASGD and EAMSGD with the parallel method DOWNPOUR and the sequential method SGD, as well as their averaging and momentum variants.

All the parallel comparator methods are listed belowWe have compared asynchronous ADMM with EASGD in our setting as well, the performance is nearly the same. However, ADMM’s momentum variant is not as stable for large communication periods.:

DOWNPOUR , the pseudo-code of the implementation of DOWNPOUR used in this paper is enclosed in the Supplement.

Momentum DOWNPOUR (MDOWNPOUR), where the Nesterov’s momentum scheme is applied to the master’s update (note it is unclear how to apply it to the local workers or for the case when τ>1\tau>1). The pseudo-code is in the Supplement.

All the sequential comparator methods (p=1p=1) are listed below:

Momentum SGD (MSGD) with constant momentum δ\delta.

ASGD with moving rate αt+1=1t+1\alpha_{t+1}=\frac{1}{t+1}.

MVASGD with moving rate α\alpha set to a constant.

We perform experiments in a deep learning setting on two benchmark datasets: CIFAR-10 (we refer to it as CIFAR) Downloaded from http://www.cs.toronto.edu/~kriz/cifar.html. and ImageNet ILSVRC 2013 (we refer to it as ImageNet) Downloaded from http://image-net.org/challenges/LSVRC/2013.. We focus on the image classification task with deep convolutional neural networks. We next explain the experimental setup. The details of the data preprocessing and prefetching are deferred to the Supplement.

For all our experiments we use a GPU-cluster interconnected with InfiniBand. Each node has 44 Titan GPU processors where each local worker corresponds to one GPU processor. The center variable of the master is stored and updated on the centralized parameter server Our implementation is available at https://github.com/sixin-zh/mpiT..

To describe the architecture of the convolutional neural network, we will first introduce a notation. Let (c,y)(c,y) denotes the size of the input image to each layer, where cc is the number of color channels and yy is both the horizontal and the vertical dimension of the input. Let CC denotes the fully-connected convolutional operator and let PP denotes the max pooling operator, DD denotes the linear operator with dropout rate equal to 0.50.5 and SS denotes the linear operator with softmax output non-linearity. We use the cross-entropy loss and all inner layers use rectified linear units. For the ImageNet experiment we use the similar approach to with the following 1111-layer convolutional neural network (3,221)C(96,108)P(96,36)C(256,32)P(256,16)C(384,14) C(384,13)C(256,12)P(256,6)D(4096,1)D(4096,1)S(1000,1). For the CIFAR experiment we use the similar approach to with the following 77-layer convolutional neural network (3,28)C(64,24)P(64,12)C(128,8)P(128,4)C(64,2)D(256,1)S(10,1).

In our experiments all the methods we run use the same initial parameter chosen randomly, except that we set all the biases to zero for CIFAR case and to 0.1 for ImageNet case. This parameter is used to initialize the master and all the local workersOn the contrary, initializing the local workers and the master with different random seeds ’traps’ the algorithm in the symmetry breaking phase.. We add l2l_{2}-regularization λ2\normx2\frac{\lambda}{2}\norm{}{x}^{2} to the loss function F(x)F(x). For ImageNet we use λ=105\lambda=10^{-5} and for CIFAR we use λ=104\lambda=10^{-4}. We also compute the stochastic gradient using mini-batches of sample size 128128.

2 Experimental results

For all experiments in this section we use EASGD with β=0.9\beta=0.9Intuitively the ’effective β\beta’ is β/τ=pα=pηρ\beta/\tau=p\alpha=p\eta\rho (thus ρ=βτpη\rho=\frac{\beta}{\tau p\eta}) in the asynchronous setting. , for all momentum-based methods we set the momentum term δ=0.99\delta=0.99 and finally for MVADOWNPOUR we set the moving rate to α=0.001\alpha=0.001. We start with the experiment on CIFAR dataset with p=4p=4 local workers running on a single computing node. For all the methods, we examined the communication periods from the following set τ={1,4,16,64}\tau=\{1,4,16,64\}. For comparison we also report the performance of MSGD which outperformed SGD, ASGD and MVASGD as shown in Figure 6 in the Supplement. For each method we examined a wide range of learning rates (the learning rates explored in all experiments are summarized in Table 1, 2, 3 in the Supplement). The CIFAR experiment was run 33 times independently from the same initialization and for each method we report its best performance measured by the smallest achievable test error.

From the results in Figure 2, we conclude that all DOWNPOUR-based methods achieve their best performance (test error) for small τ\tau (τ{1,4}\tau\in\{1,4\}), and become highly unstable for τ{16,64}\tau\in\{16,64\}. While EAMSGD significantly outperforms comparator methods for all values of τ\tau by having faster convergence. It also finds better-quality solution measured by the test error and this advantage becomes more significant for τ{16,64}\tau\in\{16,64\}. Note that the tendency to achieve better test performance with larger τ\tau is also characteristic for the EASGD algorithm.

We next explore different number of local workers pp from the set p={4,8,16}p=\{4,8,16\} for the CIFAR experiment, and p={4,8}p=\{4,8\} for the ImageNet experimentFor the ImageNet experiment, the training loss is measured on a subset of the training data of size 50,000.. For the ImageNet experiment we report the results of one run with the best setting we have found. EASGD and EAMSGD were run with τ=10\tau=10 whereas DOWNPOUR and MDOWNPOUR were run with τ=1\tau=1. The results are in Figure 3 and 4. For the CIFAR experiment, it’s noticeable that the lowest achievable test error by either EASGD or EAMSGD decreases with larger pp. This can potentially be explained by the fact that larger pp allows for more exploration of the parameter space. In the Supplement, we discuss further the trade-off between exploration and exploitation as a function of the learning rate (section 9.5) and the communication period (section 9.6). Finally, the results obtained for the ImageNet experiment also shows the advantage of EAMSGD over the competitor methods.

Conclusion

In this paper we describe a new algorithm called EASGD and its variants for training deep neural networks in the stochastic setting when the computations are parallelized over multiple GPUs. Experiments demonstrate that this new algorithm quickly achieves improvement in test error compared to more common baseline approaches such as DOWNPOUR and its variants. We show that our approach is very stable and plausible under communication constraints. We provide the stability analysis of the asynchronous EASGD in the round-robin scheme, and show the theoretical advantage of the method over ADMM. The different behavior of the EASGD algorithm from its momentum-based variant EAMSGD is intriguing and will be studied in future works.

The authors thank R. Power, J. Li for implementation guidance, J. Bruna, O. Henaff, C. Farabet, A. Szlam, Y. Bakhtin for helpful discussion, P. L. Combettes, S. Bengio and the referees for valuable feedback.

References

Additional theoretical results and proofs

We provide here the convergence analysis of the synchronous EASGD algorithm with constant learning rate. The analysis is focused on the convergence of the center variable to the local optimum. We discuss one-dimensional quadratic case first, then the generalization to multi-dimensional setting (Lemma 7.3) and finally to the strongly convex case (Theorem 7.1).

It follows from Lemma 7.1 that for the center variable to be stable the following has to hold

It can be verified that ϕ\phi and γ\gamma are the two zero-roots of the polynomial in λ\lambda: λ2(2a)λ+(1a+c2)\lambda^{2}-(2-a)\lambda+(1-a+c^{2}). Recall that ϕ\phi and λ\lambda are the functions of η\eta and α\alpha. Thus (see proof in Section 7.1.2)

γ<1\gamma<1 iff c2>0c^{2}>0 (i.e. η>0\eta>0 and α>0\alpha>0).

ϕ>1\phi>-1 iff (2ηh)(2pα)>2α(2-\eta h)(2-p\alpha)>2\alpha and (2ηh)+(2pα)>α(2-\eta h)+(2-p\alpha)>\alpha.

ϕ=γ\phi=\gamma iff a2=4c2a^{2}=4c^{2} (i.e. ηh=α=0\eta h=\alpha=0).

The proof the above Lemma is based on the diagonalization of the linear gradient map (this map is symmetric due to the relation β=pα\beta=p\alpha). The stability analysis of the asynchronous EASGD algorithm in the round-robin scheme is similar due to this elastic symmetry.

Substituting the gradient from Equation 20 into the update rule used by each local worker in the synchronous EASGD algorithm (Equation 5 and 6) we obtain

where η\eta is the learning rate, and α\alpha is the moving rate. Recall that α=ηρ\alpha=\eta\rho and A=hA=h.

and the (diffusion) vector bt=(ηξt1,,ηξtp,0)Tb_{t}=(\eta\xi^{1}_{t},\ldots,\eta\xi^{p}_{t},0)^{T}.

By combining Equation 25 and 26 as follows

where the last step results from the following relations: pα21pαϕ=1αηhϕ\frac{p\alpha^{2}}{1-p\alpha-\phi}=1-\alpha-\eta h-\phi and ϕ+γ=1αηh+1pα\phi+\gamma=1-\alpha-\eta h+1-p\alpha. Thus we obtained

Substituting u0,u1,,utu_{0},u_{1},\dots,u_{t}, each given through Equation 28, into Equation 29 we obtain

To be more specific, the Equation 30 is obtained by integrating by parts,

Since the random variables ξl\xi_{l} are i.i.d, we may sum the variance term by term as follows

1.2 Condition in Equation 23

γ<1\gamma<1 iff c2>0c^{2}>0 (i.e. η>0\eta>0 and β>0\beta>0).

ϕ>1\phi>-1 iff (2ηh)(2β)>2β/p(2-\eta h)(2-\beta)>2\beta/p and (2ηh)+(2β)>β/p(2-\eta h)+(2-\beta)>\beta/p.

ϕ=γ\phi=\gamma iff a2=4c2a^{2}=4c^{2} (i.e. ηh=β=0\eta h=\beta=0).

Recall that a=ηh+(p+1)αa=\eta h+(p+1)\alpha, c2=ηhpαc^{2}=\eta hp\alpha, γ=1aa24c22\gamma=1-\frac{a-\sqrt{a^{2}-4c^{2}}}{2}, ϕ=1a+a24c22\phi=1-\frac{a+\sqrt{a^{2}-4c^{2}}}{2}, and β=pα\beta=p\alpha. We have

γ<1aa24c22>0a>a24c2a2>a24c2c2>0\gamma<1\Leftrightarrow\frac{a-\sqrt{a^{2}-4c^{2}}}{2}>0\Leftrightarrow a>\sqrt{a^{2}-4c^{2}}\Leftrightarrow a^{2}>a^{2}-4c^{2}\Leftrightarrow c^{2}>0.

ϕ>12>a+a24c224a>a24c24a>0,(4a)2>a24c24a>0,42a+c2>04>ηh+β+α,42(ηh+β+α)+ηhβ>0\phi>-1\Leftrightarrow 2>\frac{a+\sqrt{a^{2}-4c^{2}}}{2}\Leftrightarrow 4-a>\sqrt{a^{2}-4c^{2}}\Leftrightarrow 4-a>0,(4-a)^{2}>a^{2}-4c^{2}\Leftrightarrow 4-a>0,4-2a+c^{2}>0\Leftrightarrow 4>\eta h+\beta+\alpha,4-2(\eta h+\beta+\alpha)+\eta h\beta>0.

ϕ=γa24c2=0a2=4c2\phi=\gamma\Leftrightarrow\sqrt{a^{2}-4c^{2}}=0\Leftrightarrow a^{2}=4c^{2}.

The next corollary is a consequence of Lemma 7.1. As the number of workers pp grows, the averaging property of the EASGD can be characterized as follows

Let the Elastic Averaging relation β=pα\beta=p\alpha and the condition 23 hold, then

Note that when β\beta is fixed, limpa=ηh+β\lim_{p\rightarrow\infty}a=\eta h+\beta and c2=ηhβc^{2}=\eta h\beta. Then limpϕ=min(1β,1ηh)\lim_{p\rightarrow\infty}\phi=\min(1-\beta,1-\eta h) and limpγ=max(1β,1ηh)\lim_{p\rightarrow\infty}\gamma=\max(1-\beta,1-\eta h). Also note that using Lemma 7.1 we obtain

Corollary 7.1 is obtained by plugining in the limiting values of ϕ\phi and γ\gamma. ∎

The crucial point of Corollary 7.1 is that the MSE in the limit tt\rightarrow\infty is in the order of 1/p1/p which implies that as the number of processors pp grows, the MSE will decrease for the EASGD algorithm. Also note that the smaller the β\beta is (recall that β=pα=pηρ\beta=p\alpha=p\eta\rho), the more exploration is allowed (small ρ\rho) and simultaneously the smaller the MSE is.

2 Generalization to multidimensional case

The next lemma (Lemma 7.2) shows that EASGD algorithm achieves the highest possible rate of convergence when we consider the double averaging sequence (similarly to ) {z1,z2,}\{z_{1},z_{2},\dots\} defined as below

If the condition in Equation 23 holds, then the normalized double averaging sequence defined in Equation 32 converges weakly to the normal distribution with zero mean and variance σ2/ph2\sigma^{2}/ph^{2},

Also recall that {ξti}\{\xi^{i}_{t}\}’s are i.i.d. random variables (noise) with zero mean and the same covariance Σ0\Sigma\succ 0. We are interested in the asymptotic behavior of the double averaging sequence {z1,z2,}\{z_{1},z_{2},\dots\} defined as

Recall the Equation 30 from the proof of Lemma 7.1 (for the convenience it is provided below):

where ξt=i=1pξti\xi_{t}=\sum_{i=1}^{p}\xi_{t}^{i}. Therefore

Therefore the expression in Equation 35 is asymptotically normal with zero mean and variance σ2/ph2\sigma^{2}/ph^{2}. ∎

The asymptotic variance in the Lemma 7.2 is optimal with any fixed η\eta and β\beta for which Equation 23 holds. The next lemma (Lemma 7.3) extends the result in Lemma 7.2 to the multi-dimensional setting.

Let hh denotes the largest eigenvalue of AA. If (2ηh)(2β)>2β/p(2-\eta h)(2-\beta)>2\beta/p, (2ηh)+(2β)>β/p(2-\eta h)+(2-\beta)>\beta/p, η>0\eta>0 and β>0\beta>0, then the normalized double averaging sequence converges weakly to the normal distribution with zero mean and the covariance matrix V=A1Σ(A1)TV=A^{-1}\Sigma(A^{-1})^{T},

Let the spatial average of the local parameters at time tt be denoted as yty_{t} where yt=1pi=1pxtiy_{t}=\frac{1}{p}\sum_{i=1}^{p}x_{t}^{i}, and let the average noise be denoted as ξt\xi_{t}, where ξt=1pi=1pξti\xi_{t}=\frac{1}{p}\sum_{i=1}^{p}\xi_{t}^{i}. Equations 24 and 25 can then be reduced to the following

From Equations 37 and 38 it follows that Ut+1=MUt+ΞtU_{t+1}=MU_{t}+\Xi_{t}. Note that this linear system has a degenerate noise Ξt\Xi_{t} which prevents us from directly applying results of . Expanding this recursive relation and summing by parts, we have

Note that the only non-vanishing term (in weak convergence) of 1tk=0tUk\frac{1}{\sqrt{t}}\sum_{k=0}^{t}U_{k} is 1t(ηL)1k=1tΞk1\frac{1}{\sqrt{t}}(\eta L)^{-1}\sum_{k=1}^{t}\Xi_{k-1} thus we have

The eigenvalue λ\lambda of MM and the (non-zero) eigenvector (y,z)(y,z) of MM satisfy

Since (y,z)(y,z) is assumed to be non-zero, we can write z=βy/(λ+β1)z=\beta y/(\lambda+\beta-1). Then the Equation 50 can be reduced to

Thus yy is the eigenvector of AA. Let λA\lambda_{A} be the eigenvalue of matrix AA such that Ay=λAyAy=\lambda_{A}y. Thus based on Equation 51 it follows that

where a=ηλA+(p+1)αa=\eta\lambda_{A}+(p+1)\alpha, c2=ηλApαc^{2}=\eta\lambda_{A}p\alpha. It follows from the condition in Equation 23 that 1<λ<1-1<\lambda<1 iff η>0\eta>0, β>0\beta>0, (2ηλA)(2β)>2β/p(2-\eta\lambda_{A})(2-\beta)>2\beta/p and (2ηλA)+(2β)>β/p(2-\eta\lambda_{A})+(2-\beta)>\beta/p. Let hh denote the maximum eigenvalue of AA and note that 2ηλA2ηh2-\eta\lambda_{A}\geq 2-\eta h. This implies that the condition of our lemma is sufficient. ∎

As in Lemma 7.2, the asymptotic covariance in the Lemma 7.3 is optimal, i.e. meets the Fisher information lower-bound. The fact that this asymptotic covariance matrix VV does not contain any term involving ρ\rho is quite remarkable, since the penalty term ρ\rho does have an impact on the condition number of the Hessian in Equation 2.

3 Strongly convex case

We have thus the update for the spatial average,

From Equation 54 the following relation holds,

By the cosine rule (2\scalprodL2abcd=\normL2ad2\normL2ac2+\normL2cb2\normL2db22\scalprod{L2}{a-b}{c-d}=\norm{L2}{a-d}^{2}-\norm{L2}{a-c}^{2}+\norm{L2}{c-b}^{2}-\norm{L2}{d-b}^{2}), we have

By the Cauchy-Schwarz inequality, we have

Combining the above estimates in Equations 57, 58, 59, 60, we obtain

Now we apply similar idea to estimate \normL2ytx2\norm{L2}{y_{t}-x^{\ast}}^{2}. From Equation 56 the following relation holds,

By \scalprodL21pi=1pai1pj=1pbj=1pi=1p\scalprodL2aibi1p2i>j\scalprodL2aiajbibj\scalprod{L2}{\frac{1}{p}\sum_{i=1}^{p}a_{i}}{\frac{1}{p}\sum_{j=1}^{p}b_{j}}=\frac{1}{p}\sum_{i=1}^{p}\scalprod{L2}{a_{i}}{b_{i}}-\frac{1}{p^{2}}\sum_{i>j}\scalprod{L2}{a_{i}-a_{j}}{b_{i}-b_{j}}, we have

Denote ξt=1pi=1pξti\xi_{t}=\frac{1}{p}\sum_{i=1}^{p}\xi^{i}_{t}, we can rewrite Equation 65 as

By combining the above Equations 66, 67 with 68, we obtain

Thus it follows from Equation 57 and 70 that

Recall yt=1pi=1pxtiy_{t}=\frac{1}{p}\sum_{i=1}^{p}x^{i}_{t}, we have the following bias-variance relation,

By the Cauchy-Schwarz inequality, we have

Combining the above estimates in Equations 71, 72, 73, we obtain

Since 2μLμ+L1\frac{2\sqrt{\mu L}}{\mu+L}\leq 1, we need also bound the non-linear term \scalprodL2ftiftjxtixtjL\normL2xtixtj2\scalprod{L2}{\nabla f^{i}_{t}-\nabla f^{j}_{t}}{x_{t}^{i}-x_{t}^{j}}\leq L\norm{L2}{x_{t}^{i}-x_{t}^{j}}^{2}. Recall the bias-variance relation 1pi=1p\normL2xtix2=1p2i>j\normL2xtixtj2+\normL2ytx2\frac{1}{p}\sum_{i=1}^{p}\norm{L2}{x^{i}_{t}-x^{\ast}}^{2}=\frac{1}{p^{2}}\sum_{i>j}\norm{L2}{x^{i}_{t}-x^{j}_{t}}^{2}+\norm{L2}{y_{t}-x^{\ast}}^{2}. The key observation is that if 1pi=1p\normL2xtix2\frac{1}{p}\sum_{i=1}^{p}\norm{L2}{x^{i}_{t}-x^{\ast}}^{2} remains bounded, then larger variance i>j\normL2xtixtj2\sum_{i>j}\norm{L2}{x^{i}_{t}-x^{j}_{t}}^{2} implies smaller bias \normL2ytx2\norm{L2}{y_{t}-x^{\ast}}^{2}. Thus this non-linear term can be compensated.

Again choose η\eta small enough such that η2+η2α1α2ημ+L0\eta^{2}+\frac{\eta^{2}\alpha}{1-\alpha}-\frac{2\eta}{\mu+L}\leq 0 and take expectation in Equation 75,

As for the center variable in Equation 55, we apply simply the convexity of the norm \normL22\norm{L2}{\cdot}^{2} to obtain

as long as 0β10\leq\beta\leq 1, 0α<10\leq\alpha<1 and η2+η2α1α2ημ+L0\eta^{2}+\frac{\eta^{2}\alpha}{1-\alpha}-\frac{2\eta}{\mu+L}\leq 0, i.e. 0η2μ+L(1α).0\leq\eta\leq\frac{2}{\mu+L}(1-\alpha).

Additional pseudo-codes of the algorithms

Algorithm 3 captures the pseudo-code of the implementation of the DOWNPOUR used in this paper.

2 MDOWNPOUR pseudo-code

Algorithms 4 and 5 capture the pseudo-codes of the implementation of momentum DOWNPOUR (MDOWNPOUR) used in this paper. Algorithm 4 shows the behavior of each local worker and Algorithm 5 shows the behavior of the master.

Experiments - additional material

For the ImageNet experiment, we re-size each RGB image so that the smallest dimension is 256256 pixels. We also re-scale each pixel value to the interval $.Wethenextractrandomcrops(andtheirhorizontalflips)ofsize. We then extract random crops (and their horizontal flips) of size3\times 221\times 221pixelsandpresentthesetothenetworkinminibatchesofsizepixels and present these to the network in mini-batches of size128$.

For the CIFAR experiment, we use the original RGB image of size 3×32×323\times 32\times 32. As before, we re-scale each pixel value to the interval $.Wethenextractrandomcrops(andtheirhorizontalflips)ofsize. We then extract random crops (and their horizontal flips) of size3\times 28\times 28pixelsandpresentthesetothenetworkinminibatchesofsizepixels and present these to the network in mini-batches of size128$.

The training and test loss and the test error are only computed from the center patch (3×28×283\times 28\times 28) for the CIFAR experiment and the center patch (3×221×2213\times 221\times 221) for the ImageNet experiment.

2 Data prefetching (Sampling the dataset by the local workers)

We will now explain precisely how the dataset is sampled by each local worker as uniformly and efficiently as possible. The general parallel data loading scheme on a single machine is as follows: we use kk CPUs, where k=8k=8, to load the data in parallel. Each data loader reads from the memory-mapped (mmap) file a chunk of cc raw images (preprocessing was described in the previous subsection) and their labels (for CIFAR c=512c=512 and for ImageNet c=64c=64). For the CIFAR, the mmap file of each data loader contains the entire dataset whereas for ImageNet, each mmap file of each data loader contains different 1/k1/k fractions of the entire dataset. A chunk of data is always sent by one of the data loaders to the first worker who requests the data. The next worker requesting the data from the same data loader will get the next chunk. Each worker requests in total kk data chunks from kk different data loaders and then process them before asking for new data chunks. Notice that each data loader cycles through the data in the mmap file, sending consecutive chunks to the workers in order in which it receives requests from them. When the data loader reaches the end of the mmap file, it selects the address in memory uniformly at random from the interval [0,s][0,s], where s=(number of images in the mmap filemodulomini-batch size)s=(\textsf{number of images in the mmap file}\text{\>\>modulo\>\>}\textsf{mini-batch size}), and uses this address to start cycling again through the data in the mmap file. After the local worker receives the kk data chunks from the data loaders, it shuffles them and divides it into mini-batches of size 128128.

3 Learning rates

In Table 1 we summarize the learning rates η\eta (we used constant learning rates) explored for each method shown in Figure 2. For all values of τ\tau the same set of learning rates was explored for each method.

In Table 2 we summarize the learning rates η\eta (we used constant learning rates) explored for each method shown in Figure 3. For all values of pp the same set of learning rates was explored for each method.

In Table 3 we summarize the initial learning rates η\eta we use for each method shown in Figure 4. For all values of pp the same set of learning rates was explored for each method. We also used the rule of the thumb to decrease the initial learning rate twice, first time we divided it by 55 and the second time by 22, when we observed that the decrease of the online predictive (training) loss saturates.

4 Comparison of SGD, ASGD, MVASGD and MSGD

Figure 6 shows the convergence of the training and test loss (negative log-likelihood) and the test error computed for the center variable as a function of wallclock time for SGD, ASGD, MVASGD and MSGD (p=1p=1) on the CIFAR experiment. For all CIFAR experiments we always start the averaging for the ADOWNPOURADOWNPOUR and ASGDASGD methods from the very beginning of each experiment. For all ImageNet experiments we start the averaging for the ASGDASGD at the same time when we first reduce the initial learning rate.

Figure 7 shows the convergence of the training and test loss (negative log-likelihood) and the test error computed for the center variable as a function of wallclock time for SGD, ASGD, MVASGD and MSGD (p=1p=1) on the ImageNet experiment.

5 Dependence of the learning rate

This section discusses the dependence of the trade-off between exploration and exploitation on the learning rate. We compare the performance of respectively EAMSGD and EASGD for different learning rates η\eta when p=16p=16 and τ=10\tau=10 on the CIFAR experiment. We observe in Figure 8 that higher learning rates η\eta lead to better test performance for the EAMSGD algorithm which potentially can be justified by the fact that they sustain higher fluctuations of the local workers. We conjecture that higher fluctuations lead to more exploration and simultaneously they also impose higher regularization. This picture however seems to be opposite for the EASGD algorithm for which larger learning rates hurt the performance of the method and lead to overfitting. Interestingly in this experiment for both EASGD and EAMSGD algorithm, the learning rate for which the best training performance was achieved simultaneously led to the worst test performance.

6 Dependence of the communication period

This section discusses the dependence of the trade-off between exploration and exploitation on the communication period. We have observed from the CIFAR experiment that EASGD algorithm exhibits very similar convergence behavior when τ=1\tau=1 up to even τ=1000\tau=1000, whereas EAMSGD can get trapped at worse energy (loss) level for τ=100\tau=100. This behavior of EAMSGD is most likely due to the non-convexity of the objective function. Luckily, it can be avoided by gradually decreasing the learning rate, i.e. increasing the penalty term ρ\rho (recall α=ηρ\alpha=\eta\rho), as shown in Figure 9. In contrast, the EASGD algorithm does not seem to get trapped at all along its trajectory. The performance of EASGD is less sensitive to increasing the communication period compared to EAMSGD, whereas for the EAMSGD the careful choice of the learning rate for large communication periods seems crucial.

Compared to all earlier results, the experiment in this section is re-run three times with a new randomTo clarify, the random initialization we use is by default in Torch’s implementation. seed and with faster cuDNNhttps://developer.nvidia.com/cuDNN packagehttps://github.com/soumith/cudnn.torch. All our methods are implemented in Torchhttp://torch.ch. The Message Passing Interface implementation MVAPICH2http://mvapich.cse.ohio-state.edu is used for the GPU-CPU communication.

7 Breakdown of the wallclock time

In addition, we report in Table 4 the breakdown of the total running time for EASGD when τ=10\tau=10 (the time breakdown for EAMSGD is almost identical) and DOWNPOUR when τ=1\tau=1 into computation time, data loading time and parameter communication time. For the CIFAR experiment the reported time corresponds to processing 400×128400\times 128 data samples whereas for the ImageNet experiment it corresponds to processing 1024×1281024\times 128 data samples. For τ=1\tau=1 and p{8,16}p\in\{8,16\} we observe that the communication time accounts for significant portion of the total running time whereas for τ=10\tau=10 the communication time becomes negligible compared to the total running time (recall that based on previous results EASGD and EAMSGD achieve best performance with larger τ\tau which is ideal in the setting when communication is time-consuming).

8 Time speed-up

In Figure 10 and 11, we summarize the wall clock time needed to achieve the same level of the test error for all the methods in the CIFAR and ImageNet experiment as a function of the number of local workers pp. For the CIFAR (Figure 10) we examined the following levels: {21%,20%,19%,18%}\{21\%,20\%,19\%,18\%\} and for the ImageNet (Figure 11) we examined: {49%,47%,45%,43%}\{49\%,47\%,45\%,43\%\}. If some method does not appear on the figure for a given test error level, it indicates that this method never achieved this level. For the CIFAR experiment we observe that from among EASGD, DOWNPOUR and MDOWNPOUR methods, the EASGD method needs less time to achieve a particular level of test error. We observe that with higher pp each of these methods does not necessarily need less time to achieve the same level of test error. This seems counter intuitive though recall that the learning rate for the methods is selected based on the smallest achievable test error. For larger pp smaller learning rates were selected than for smaller pp which explains our results. Meanwhile, the EAMSGD method achieves significant speed-up over other methods for all the test error levels. For the ImageNet experiment we observe that all methods outperform MSGD and furthermore with p=4p=4 or p=8p=8 each of these methods requires less time to achieve the same level of test error. The EAMSGD consistently needs less time than any other method, in particular DOWNPOUR, to achieve any of the test error levels.