A Hybrid-Order Distributed SGD Method for Non-Convex Optimization to Balance Communication Overhead, Computational Complexity, and Convergence Rate
Naeimeh Omidvar Mohammad Ali Maddah-Ali Hamed Mahdavi
cs.LG cs.DC cs.IT math.IT math.OC stat.ML G.1.6


[







A Hybrid-Order Distributed SGD Method for Non-Convex Optimization to Balance Communication Overhead, Computational Complexity, and Convergence Rate
     














equal*



Naeimeh Omidvaripm
Mohammad Ali Maddah-Alinokia
Hamed Mahdavice_sharif






ipmSchool of Computer Science, Institute for Research in Fundamental Sciences (IPM), Tehran, Iran (email: nomidvar@connect.ust.hk).


nokiaNokia Bell Labs (email: mohammad.maddah-ali@nokia-bell-labs.com).
ce_sharifDepartment of Computer Engineering, Sharif University of Technology, Tehran, Iran (email: hmdmahdavi@ce.sharif.edu)







Machine Learning, ICML

0.3in

]









  





In this paper, we propose a method of distributed stochastic gradient descent (SGD), with low communication load and computational complexity, and still fast convergence. To reduce the communication load, at each iteration of the algorithm, the worker nodes calculate and communicate some scalers, that are the directional derivatives of the sample functions in some pre-shared directions. However, to maintain accuracy, after every specific number of iterations, they communicate the vectors of stochastic gradients. To reduce the computational complexity in each iteration, the worker nodes approximate the directional derivatives with zeroth-order stochastic gradient estimation, by performing just two function evaluations rather than computing a first-order gradient vector. The proposed method highly improves the convergence rate of the zeroth-order methods, guaranteeing order-wise faster convergence. Moreover, compared to the famous communication-efficient methods of model averaging (that perform local model updates and periodic communication of the gradients to synchronize the local models), we prove that for the general class of non-convex stochastic problems and with reasonable choice of parameters, the proposed method guarantees the same orders of communication load and convergence rate, while having order-wise less computational complexity. Experimental results on various learning problems in neural networks applications demonstrate the effectiveness of the proposed approach compared to various state-of-the-art distributed SGD methods.


Better long version for journal:
In this paper, we propose a method of distributed stochastic gradient descent (SGD), with low communication load and computational complexity, and still fast convergence. To reduce the communication load, at each iteration of the algorithm, the worker nodes calculate and communicate some scalers, that are the directional derivatives of the sample functions in some pre-shared directions. However, to maintain accuracy, after every τ∈ℕ iterations, they communicate the vectors of stochastic gradients. To reduce the computational complexity in each iteration, the worker nodes approximate the directional derivatives with zeroth-order stochastic gradient estimation, by performing just two function evaluations rather than computing a first-order gradient vector. 
For a problem of size d ∈ℕ,  the communication load of the proposed method is 1+(τ-1)/d times the communication load of the famous communication-efficient methods of model averaging (that perform local model updates and periodic communication of the gradients to synchronize the local models). Moreover, its computational load is 𝒪( 1/τ+1/d ) times the computational load of model averaging schemes. Finally, by theoretical analyses, we show that for the general class of non-convex stochastic problems, the proposed scheme converges to a stationary point with the rate of 𝒪( d/√(mN)), which highly outperforms the convergence rate of the existing zeroth-order methods, guaranteeing order-wise faster convergence. Moreover, the convergence rate of the fastest model averaging scheme is 𝒪(  τ/√(mN)). Therefore, by a reasonable choice of τ = 𝒪( d ), we can guarantee the same-order convergence rate and communication load as in the fastest converging communication-efficient methods, with order-wise less computational complexity. Experimental results on various learning problems in neural networks applications demonstrate the effectiveness of the proposed approach compared to various state-of-the-art distributed SGD methods.




Dr. Maddah-Ali's version:
In this paper, we propose a method of distributed stochastic gradient descent (SGD), with low communication overhead and computational complexity, and still converges fast. To reduce the communication load, at each iteration of the algorithm, the worker nodes calculate and communicate some scalers, that are the directional derivatives of the sample functions in some pre-shared directions. However, to maintain accuracy,  after every τ∈ℕ iterations, they communicate the vectors of stochastic gradients.  To reduce the computational complexity in each iteration, the worker nodes approximate the directional derivatives with zeroth-order stochastic gradient estimation, by performing just two function evaluations rather than computing a first-order gradient vector. For a problem of size d ∈ℕ,  the communication load of the proposed method is 1+(τ-1)/d times the communication load of the famous communication-efficient methods of model averaging, that perform local model updates and periodic communication of the gradients to synchronize the local models. Moreover, its computational load is 𝒪( 1/τ+1/d ) times the computational load of model averaging schemes. Finally, by theoretical analyses, we show that for the general class of non-convex stochastic problems, the proposed scheme converges to a stationary point with the rate of 𝒪( d/√(mN)), which highly improves/outperforms the rate of zeroth-order methods. Moreover, the convergence rate of the fastest model averaging scheme is 𝒪(  τ/√(mN)). Thus, by a reasonable choice of τ = 𝒪( d ), we can guarantee the same-order communication load and convergence rate as in the fastest converging communication-efficient methods, with order-wise less computational complexity.  Experimental results on various learning problems on neural networks applications demonstrate the effectiveness of the proposed approach compared to various state-of-the-art methods. 






§ INTRODUCTION



Stochastic gradient descent (SGD) is an optimization tool 
 which is widely used  for solving many machine learning problems, 
due to its simplicity 
 and acceptable empirical performance 
  <cit.>. 
 It iteratively updates the model parameters at the opposite direction of the stochastic gradient of the cost function. 
With the emergence of large-scale machine learning models such as deep neural networks (DNNs), however, centralized deployment of SGD has become intractable in its memory and time requirements. 
As such, and accelerated with the recent advances in multi-core parallel processing technology, distributed deployment of SGD 
 has become a trend, which trains machine learning models 
 on multiple computation nodes (a.k.a., worker nodes) in parallel, to enhance the scalability of the training procedure   <cit.>. 
Under a distributed SGD method, at each iteration, each 
worker node evaluates a stochastic estimation of the gradient of the objective 
function over a randomly chosen/arrived data sample. The workers  
then broadcast their gradient updates 
to their peers, and aggregate 
 the gradients to update the model in parallel. 

Some of the major factors in designing distributed SGD algorithms, which are in conflict with each other and need to be addressed carefully, can be identified as follows:

 
1- Communication overhead: 
Distributed implementation of SGD algorithms mainly suffers from high 
level of communication overhead due to exchanging  
the stochastic gradient vectors among the worker nodes at each iteration, especially when the dimension of the model is large.  
In fact, the communication overhead has been observed to be 
the major bottleneck in scaling distributed SGD, and the time required to exchange this information will increase the overall time required for the algorithm to converge <cit.>.  
As such, alleviating the communication cost has recently gained lots of attention from different research communities. 

To address the issue of communication overhead, two main approaches can be identified 
in the literature: 1- 
Reducing the number of rounds of communication among the worker nodes by 
allowing them to perform local model updates at each iteration, and limiting their synchronization 
to a periodic exchange of local models after some 
iterations. This approach is often referred to as model averaging <cit.>, 2- Reducing the number of bits communicated among the worker nodes 
by compression or quantization of the gradient vectors <cit.>. This approach  propagates the quantization error through the iterations, reduces the convergence rate, and would add to the computational complexity of each iteration. 


2- Computation load: In many practical applications, the worker nodes are generally commodity devices with limited computational resources. 
Therefore, in order to guarantee a scalable solution, 
the distributed algorithm should essentially 
impose  as low computational cost 
per worker node as possible. 

Note that in general, one of the major derives for the computational load in SGD algorithms is the calculation of the first-order gradients at each iteration, which is very expensive, if not impossible, to obtain in many real-world problems (especially for learning large-scale models) <cit.>. In many applications, calculating  
the stochastic gradient requires extensive simulation for each worker node, where the complexity of each simulation may require significant computational time <cit.>. 
Moreover, in many scenarios of training 
DNNs,  
there is a highly complicated relationship between the 
parameters of the model and the objective function, so that deriving explicit form of the gradient is impossible and computing the gradient at each iteration costs significantly high computational complexity for the worker node <cit.>.[It is noted that although 
under fast differentiation techniques, 
the gradient of the sample function in neural network can be derived with less complexity, there exist some restrictions for applying such technique to a general neural network. In particular, such techniques require to store the results of all the intermediate computations, which is 
impractical 
for many scenarios due to memory limitations <cit.>.] 

We note that a promising approach to reduce the computational load of SGD methods is utilizing zeroth-order (ZO) gradient estimations instead of deriving the first-order gradients. 
In fact, since calculating a ZO gradient requires just two function evaluations <cit.>, it highly reduces the computational cost compared to  
the calculation of a first-order gradient (which imposes 𝒪(d) times higher  computational complexity in general <cit.>).  
Furthermore, by using pre-shared seeds for the generation of the random directions involved in zeroth-order gradients calculation, the worker nodes will not need to communicate the whole ZO gradient vectors, and can send just the scalar values of the  computed 
derivative approximations, as will be shown in Section <ref>. As such, the communication overhead is significantly reduced as compared to the communication overhead imposed by sending the whole gradient vectors, 
especially when the dimension of the problem is large. 
However, loosing accuracy in calculating the gradients directly reflects in the convergence rates of zeroth-order SGD methods, especially for non-convex problems (e.g., see the convergence rates comparison in Table <ref>). This is one of the main 
challenges that we overcome in this work. 

3- Convergence rate: 
Finally achieving fast convergence to the solution of a problem is essentially desired, especially for large-scale problems dealing with huge datasets. 



 §.§ Main Contributions


In this work, we propose a distributed 
optimization method that 
strikes a better balance between communication efficiency, computation efficiency and accuracy, compared to various distributed methods. The proposed methods enjoys the low computational complexity and communication overhead of the zeroth-order gradient calculation, and at the same time, benefits from periodic first-order gradient calculation and model updating.  
We theoretically prove that by proper combining of these two, we can achieve 
very good 
convergence rate as well. 
The main contributions of this work can be summarized as follows: 






 








  * We develop a new distributed SGD method, with low communication overhead and computational complexity, and yet fast convergence. To reduce the communication load, at each iteration of the algorithm, the worker nodes calculate and communicate some scalers, that are the directional derivatives of the sample functions in some pre-shared directions. To reduce the computational complexity, the worker nodes approximate the directional derivatives with zeroth-order stochastic gradient estimation, by performing just two function evaluations.  
Finally, to alleviate the  approximation error of the zeroh-order stochastic gradient estimations, 
after every τ∈ℕ iterations, the worker nodes compute and communicate the first-order stochastic gradient vectors.


  * We provide theoretical analyses for the convergence rate guarantee of the proposed method. 
In particular, we prove that for the general class of non-convex stochastic problems, the proposed scheme converges to a stationary point with the rate of 𝒪 ( d/√(mN) ), which highly outperforms the convergence rates  of zeroth-order methods (e.g., see <cit.>). 
Moreover, compared to the fastest first-order communication-efficient algorithms (e.g., 
the model averaging scheme in <cit.> with the convergence rate of 𝒪 ( τ/√(m N)  )), the proposed algorithm exhibits the same convergence rate in terms of both the number of iterations and the number of worker nodes.  Finally, by a reasonable choice of  τ = 𝒪( d ), we can guarantee 
the same-order communication load and convergence rate as in the fastest converging 
model averaging methods, with order-wise less computational complexity (see Table <ref>). 
 


  * Due to utilizing zeroth-order stochastic gradient updates, the proposed algorithm exhibits a sufficiently low level of computational complexity, which is 
 much lower than the one in the communication-efficient methods including the model averaging schemes (e.g., <cit.>) or the schemes with gradients compression (e.g., <cit.>), and is comparable to the one in the existing zeroth-order stochastic optimization algorithms (e.g., <cit.>). 
As a baseline, computational complexity of the proposed method is 𝒪( 1/τ+1/d ) times the computational complexity of the model averaging schemes. 



  * Due to utilizing pre-shared seeds for the generation of random directions, at the iterations with zeroth-order gradient updates, the worker nodes just need to communicate a scalar rather than a d-dimensional vector, 
which drastically reduces the number of bits required for communication among the worker nodes.  
In particular,  
the proposed algorithm exhibits a communication load equal to sending  
d + τ -1 scalar values by each worker node per  τ iterations. This communication load is much lower than that of 
fully synchronous SGD (syncSGD)  <cit.>, and is comparable to that of the communication-efficient schemes with model averaging, with a ratio of ( 1+τ/d ).




  * Using numerical experiments, we empirically demonstrate the accuracy and convergence properties of the proposed method compared to various  state-of-the-art baselines. 


 
Table <ref> summarizes  
the comparison of the proposed method to the most related works in the literature,  
in terms of convergence rate, communication overhead per iteration, and computational load per iteration normalized to the computational complexity of  
computing a first-order stochastic gradient.




§ RELATED WORK


As mentioned before, communication bottlenecks in distributed SGD originate from two sources that can be significant; first, the number of bits communicated at each communication round, and second, the number of rounds of communication. 
To tackle the first barrier, the existing works mainly try to quantize and compress the gradient vectors before communicating them <cit.>. 
Although those methods are generally effective in reducing the number of bits communicated at each iteration, 
their quantization error increases the error variance of the communicated gradient vectors, leading to a slow convergence. 
To tackle the second  
source of communication overhead, 
the idea of model averaging has been proposed, where  the worker nodes perform local updates at each iteration, and communicate their updated 
models after every τ iterations 
 to periodically synchronize them
 <cit.>. In consequence, the number of rounds of communication is significantly reduced. 


There exist various model averaging schemes in the literature. 
<cit.>  
investigated periodic averaging SGD (PA-SGD), where the models are averaged across the nodes afters every τ local updates. 
<cit.> 
 demonstrated that the convergence error of PA-SGD grows linearly in the number of local updates τ. <cit.> provided some theoretical studies on 
why SGD with model averaging works well.  
Moreover, <cit.> 
proposed elastic averaging SGD (EASGD), where 
a more complicated averaging is done by encompassing 
 an auxiliary variable 
to allow some slackness between the models. 
They empirically validated the effectiveness of EASGD, without providing any convergence analysis. 
Furthermore,  <cit.> 
 proposed consensus-based distributed SGD (D-PSGD), in which worker nodes 
 synchronize their local models only with their neighboring nodes. 
By incorporating extra memory as well as high storage, <cit.> infuse redundancy to the training data to further reduce the residual error in local averaging, and improve the convergence rate. 
However, this method needs more storage per worker node, and the data at different worker nodes are overlapping. 

Finally, note that a key feature for both the aforementioned 
communication-efficient approaches is that they require  
stochastic first-order gradient information at each iteration (via subsequent calls to a stochastic first-order oracle (SFO)), in order to guarantee the convergence. However, as illustrated before, 
 in many real-world applications and scenarios, 
it may be computationally costly 
 for the worker nodes to obtain such information at every iteration <cit.>. 
Therefore, although the approaches mentioned above 
provide communication-efficient solutions to distributed learning, 
their computational complexity may not be 
tolerable for general commodity worker nodes, 
and restricts the applicability of such methods, in practice. 

To reduce the computational load of SGD methods, the idea of zeroth-order gradient estimations can be utilized.  
This idea has been widely used for gradient-free optimization where the explicit expressions of gradients of the objective function are expensive or infeasible to obtain, and only function evaluations are accessible <cit.>. As earlier discussed in the Introduction section, ZO gradient estimation can highly reduce the communication overhead as well as the computational complexity. However, the aforementioned benefits of zeroth-order SGD methods   
comes at the cost of significantly deteriorated convergence rates. 
This is because a zeroth-order gradient is in fact a biased estimation of the true gradient, and the involved approximation error leads to a high residual error and consequently,  
inferior convergence rates in zeroth-order SGD algorithms 
<cit.>. 
For example, for general non-convex problems, the (centralized) zeroth-order SGD algorithm proposed in  <cit.> has a 
convergence rate of 𝒪 (   d^1/3/ m^1/3 N^1/4 ) (where  m  denotes the number of sampled directions at each iteration of the algorithm, and hence, can be considered as equivalent to the number of worker nodes for a distributed setting). 

By incorporating multiple restarts and extra memory, <cit.> proposed a zeroth-order extension of stochastic variance reduced gradient method (SVRG) and proved that their proposed zeroth-order SVRG (ZO-SVRG-Ave) achieves a convergence rate of 
𝒪  (d/N+1/min{d,m}  ) 
for non-convex problems, in which the second term 
highly deteriorates the convergence performance of the algorithm. 
To eliminate this error term, they proposed a coordinate-wise version of ZO-SVRG, but it costs 𝒪 (d) times more function evaluations 
per iteration, which causes high computational complexity, especially for high-dimensional problems. 
Moreover, both of their proposed ZO-SVRG methods required full dataset storage, which is not affordable for distributed deployment on commodity devices. 
Furthermore, <cit.> proposed a zeroth-order version of the stochastic alternating direction method of multipliers with a guaranteed  convergence rate of 𝒪 (1/√(N)) for convex problems but cannot be applied to general non-convex problems. 
We note that the slow convergence of zeroth-order stochastic gradient methods, especially for non-convex problems, is one of the main challenges that we address in this paper. 






§ COMMUNICATION-EFFICIENT DISTRIBUTED STOCHASTIC OPTIMIZATION ALGORITHM





 §.§ Problem Definition


Consider a setting with m distributed worker nodes that are interested in solving the following non-convex stochastic optimization problem in a distributed manner:

    min_x∈ℝ^d  f ( x) = min_x∈ℝ^d𝔼_ζ[ F ( x , ζ) ],

where f is a generic non-convex loss function, and ζ is a 
random variable with unknown distribution. 
As a special case of the above 
formulation,  
consider the following non-convex  
sum 
problem  
which appears in a variety of machine learning applications, ranging from generalized linear models to deep neural networks: 

    min_x∈ℝ^d  f ( x) = min_x∈ℝ^d1K∑_k=1^K F ( x , ζ_k ) ,


in which F ( x , ζ_k ) is the  
training loss over the
sample ζ_k, 
 and K is the total number of samples. 
Note that for generality, in the rest of this paper, we focus on the problem formulation (<ref>). 
However, our analyses apply equally to both objectives. 





 §.§ The Proposed Algorithm


In the proposed algorithm, named as Hybrid-Order Distributed SGD (HO-SGD), at each iteration t, each worker node i 
samples or receives data sample ζ_t+1,i identically and independently from the dataset. 
Then, the data sample is used to evaluate 
a stochastic zeroth-order gradient 
approximation 
via computing a finite difference using two 
function queries 
as follows <cit.>. 

    G̃_t,i = 
    dμ[  F ( x^t  + μv_t+1,i , ζ_t+1,i)  -  F ( x^t , ζ_t+1,i) ] v_t+1,i,

where d is 
the dimension of the model to be learned, μ > 0 is a smoothing parameter, and v_t+1,i is a  
direction randomly drawn from a uniform distribution over a unit sphere. 

Next, each worker node i 
communicates its computed zeroth-order gradient to the other nodes. For this purpose, note that  
the 
direction of the derivatives, i.e., the vectors v_t+1,i ,  ∀ i, are some randomly generated directions, where the seeds are pre-shared among the nodes before optimization. As such, each worker node i does not need to send the vector and just needs to send  the value of the finite-difference approximated directional derivative, i.e., the scalar dμ[ F( x^t + μv_t+1,i , ζ_t+1,i) - F( x^t , ζ_t+1,i) ].  
Therefore, instead of communicating a d-dimensional gradient vector, each node  communicates just a single scalar. As a result, the communication load reduces to d times less than the communication load of exchanging the 
stochastic gradient vector. 
The model is then updated by each worker node in parallel, using the average of the local zeroth-order gradients of all the nodes, as shown by (<ref>)-(<ref>).

Finally, after every τ -1 iterations, the worker nodes perform one iteration of first-order stochastic gradient computation and communication, and update the model accordingly. 
The pseudo-code of the proposed algorithm is shown in Algorithm <ref>. 
We note that the introduced algorithm does not necessarily assume 
that each worker node has access to 
the entire data. 
Rather, as long as each data sample is assigned to each worker node uniformly at random, Algorithm <ref> works, and all the  results apply. 




 §.§ Discussion on the Proposed Algorithm and Comparison to the Related Works
  

In the following, we review some remarks regarding  Algorithm <ref>. 
First, note that if we choose   τ=1, Algorithm <ref> reduces to fully synchronous 
distributed SGD method  <cit.>, where the workers perform first-order gradient computation and communication at all iterations. Moreover, if we consider τ≥ N, the workers always perform zeroth-order gradient updates, thereby  
the algorithm reduces to distributed zeroth-order stochastic gradient method. 
Therefore, as the two ends of its spectrum, the proposed algorithm encompasses both zeroth-order and first-order distributed SGD algorithms as special cases. 





To address the challenge of slow convergence of zeroth-order SGD iterations, the proposed method employs periodic rounds of first-order stochastic gradient updates. This  significantly reduces the residual error of the zeroth-order stochastic gradient updates, as compared to the existing zeroth-order methods. 
Therefore, our algorithm can be viewed as a zero-order stochastic optimization with periodic rounds of first-order updates, 
which reduces the communication overhead and  computational complexity, and at the same time guarantees  improved convergence rate (as will be seen in the next section). 
In particular, the main advantages of the proposed algorithm can be identified as follows:
  


Low computational complexity: By employing 
zeroth-order stochastic gradient estimations, each worker node 
performs just two function evaluations,  rather than a  
gradient computation. 
 This contributes to a significant reduction in the computational complexity of the proposed algorithm, as compared to the previous first-order distributed methods. 
In particular, it is estimated that in general, computing a zeroth-order gradient estimation costs 𝒪( d ) times less computational load than computing a first-order gradient estimation <cit.>. 
Therefore, considering τ -1 iterations of zeroth-order update and one iteration of first-order update at each period of τ iterations of  Algorithm <ref>, the computational complexity of the proposed algorithm is extremely lower than that of the first-order 
 communication-efficient methods, 
with a ratio of 𝒪( 1/d+1/τ), and is comparable to the one in the distributed zeroth-order methods.  

 
Low communication overhead: As aforementioned, due to utilizing pre-shared seeds for the generation of random directions, at the iterations with zeroth-order gradient updates, the worker nodes just need to communicate a scalar instead of a d-dimensional vector, 
which drastically reduces the number of bits required for communication. 
In particular, 
the communication overhead of the proposed scheme is ( 1 + τ-1/d) times the communication overhead of the fastest first-order communication-efficient methods (i.e., model-averaging schemes such as <cit.>), and is comparable to the communication overhead in the zeroth-order methods. 
 
Fast convergence rate: Finally, 
the proposed algorithm guarantees a fast convergence rate for a general class of non-convex problems. As will be shown in the next section, in terms of the number of iterations and the number of workers, the proposed method guarantees the same convergence rate as those of the fastest first-order communication-efficient methods, and order-wisely better convergence rate than those of the zeroth-order methods. 

 



§ CONVERGENCE ANALYSIS
 

In this section, we present the convergence analysis of the proposed algorithm for a general class of non-convex problems.  
Prior to that, we first state the  
main assumptions and definitions used for the convergence analysis. 
 



 §.§ Assumptions and Definitions


Our convergence analysis is based on the following assumptions, which are all standard and 
widely used in the context of non-convex optimization <cit.>. 

[Unbiased and finite variance first-order stochastic gradient estimation]
The stochastic gradient evaluated on each data sample ζ by the first-order oracle is an unbiased estimator of the full (exact) gradient, i.e., 

    𝔼[ ∇ F ( x,ζ) ] = ∇ f ( x),   ∀x∈ℝ^d,
 
with a finite variance, i.e., there exists a constant σ≤ 0 such that 

    𝔼[ ∥∇ F ( x,ζ) - ∇ f ( x) ∥^2  ] ≤σ^2 ,   ∀x∈ℝ^d.
 




[Lipschitz continuous and bounded gradient]
The objective function isf ( x) is differentiable and L-smooth, i.e., its gradient ∇ f is L-Lipschitz continuous: 

    ∥∇ f ( x) - ∇ f ( y) ∥≤ L ∥x - y∥,    ∀x,y∈ℝ^d .
 
Moreover, the norm of the gradient of the objective function is bounded, i.e., there exists a constant M ≤ 0 such that 

    ∥∇ f ( x) ∥≤ M ,   ∀x∈ℝ^d.
 


[Bounded bellow objective function value]
The objective function value is bounded below by a scalar f^∗.








 §.§ Main Results


First, note that since  f ( x) is non-convex, we need a proper measure to show  
 the gap between the the output of the algorithm 
and the set of stationary solutions. 
As such and similar to the previous works, we consider the expected gradient norm as an indicator of convergence, 
and state the algorithm achieves an ϵ-suboptimal solution if 
<cit.>:  

    𝔼[ 1N∑_t=1^N ∥∇ f ( x^t ) ∥^2 ] ≤ϵ.
 
Noted that this condition guarantees convergence of the algorithm to a stationary point <cit.>. 

The following theorem presents the main convergence result of the proposed algorithm. 
Note that all the intermediate results and analytical proofs have been deferred to Appendix <ref>. 


In the proposed Algorithm 1, under Assumptions <ref>-<ref>, if the step-size (a.k.a., the learning rate) and the smoothing parameter are chosen 
such that α_t = √( B m)/L √(N),  ∀ t and μ≤1/√(dN), respectively, 
and the total iterations N is sufficiently large, i.e., 
N > 16 ( d+ B m-1 )^2/B m, 
then the average-squared gradient norm after N iterations is bounded by: 

    1N   ∑_t=0^N-1𝔼[ ∥∇ f ( x^t) ∥^2 ]   
       ≤4 L ( f ( x^0 ) - f^∗) /√( B m N)  
    +   2 σ^2 /√( B mN)τ
          + 1 (τ > 1 )  ( 4 L^2/ d^2 √( B m N)τ 
    + 4 L^2/ d^2 N √( B m N)
          + L^2/√( B mN)τ-1/τ 
    +  L^2 /N √( B mN)τ
    + 4 d σ^2/√( B m N)τ-1/τ
            + 4 d σ^2/ N √( B m N)τ
    + L^2 /√( B mN)τ-1τ 
    +  L^2 / N √( B mN)τ),
 
where 1 (·) is the indicator function. 


 
 

The above theorem indicates that 
for any τ > 1, the convergence rate is of the order of 
(d/√(mN)τ - 1/τ) = 𝒪(d/√(mN))
with respect to the parameters of the problem, which significantly outperforms the zeroth-order stochastic gradient methods (see Table <ref> for more details). 
Moreover, by setting the period of the first-order gradient updates such that τ = 𝒪( d ), we can guarantee the same convergence rate as in the fastest converging communication-efficient methods such as <cit.>, which has the convergence rate of 𝒪(τ/√(mN)).  
Also note that in the special case of  τ≥ N, i.e., when the workers only use stochastic zeroth-order oracle for gradient estimation, the proposed algorithm recovers zeroth-order SGD algorithm. 
In this case also, our obtained convergence rate significantly improves the previous results on the convergence rate of the zeroth-order SGD algorithms  <cit.>. 
Finally, when 
τ = 1, i.e., when the workers only use stochastic first-order oracle for gradient estimation, the above theorem guarantees a convergence rate of 𝒪 (1/√(mN) ), 
which is consistent with 
the convergence rate of fully synchronous SGD <cit.>, <cit.>.  




It is noted that the error upper bound (<ref>) is decomposed into two parts. The first part, which contains the first two terms,  
is similar to  
 the error bound in fully synchronous SGD <cit.>, and is resulted from the periodic steps of first-order stochastic gradient updates. The second part, containing the remaining terms,  
is a measure of the approximation error in zeroth-order stochastic gradient steps, and vanishes when τ =1, i.e., when there is no zeroth-order gradient updates. 


 
It should be noted that the upper bound (<ref>)  
grows 
very slowly with the period of first-order updates τ, with an order of 𝒪( 1 ). 
As a result, 
for a fixed number of 
 iterations,  our algorithm needs 
 few rounds of first-order stochastic gradient vectors computation and  communication among the nodes to reach a certain error bound.   
This significantly contributes to the communication efficiency and computation efficiency of the proposed algorithm. 
Moreover, such growth rate significantly outperforms the existing results in the model averaging schemes, where the error upper bound grows quadratically or linearly with τ, which is a result of high model discrepancies due to local model updates <cit.>. 





§ EXPERIMENTS



In this section, we experimentally evaluate the performance of the proposed algorithm and compare it to various 
algorithms, including model averaging schemes of PA-SGD <cit.>, RI-SGD <cit.>, 
fully synchronous SGD 
<cit.>, 
zeroth-order stochastic gradient method ZO-SGD <cit.>,  
and the zeroth-order stochastic variance reduced gradient method 
ZO-SVRG-Ave <cit.>. 
We evaluate the performance of the proposed algorithm on two applications, as follows:  
The first application is generation of adversarial examples from DNNs, 
which arises in testing the robustness of a deployed DNN to adversarial attacks 
<cit.>, and is performed on a system with an Nvidia Tesla  K80 GPU. 
The second application is training a 2-layer fully connected DNN with the dataset of COVTYPE7. This experiment is performed on a system with 8 Cores of CPU and 4 Nvidia Tesla K80 GPUs. 
More details on the second experiment can be found in the Appendix.


In this section, we experimentally evaluate the performance of the proposed algorithm and compare it to various state-of-the-art  distributed 
algorithms, including the model averaging scheme  
RI-SGD <cit.>, 

fully synchronous SGD 
(syncSGD) <cit.>, 
zeroth-order stochastic gradient method (ZO-SGD) <cit.>, 
and the zeroth-order stochastic variance reduced gradient method 
(ZO-SVRG-Ave) <cit.>.

We evaluate the performance of the proposed algorithm on on the application of generation of adversarial examples from DNNs, 

which arises in testing the robustness of a deployed DNN to adversarial attacks 
<cit.>.
We evaluate the performance of the proposed algorithm on two different applications, as follows.   
*****

The first application is generation of adversarial examples from DNNs, 

which arises in testing the robustness of a deployed DNN to adversarial attacks 
<cit.>. , and is performed on a system with a on a system with one Nvidia Tesla  K80 GPU. 
The second application is a multi-class classification task with multiple workers. 








 §.§ Generation of Adversarial Examples from DNNs


The first application is generation of adversarial examples from DNNs, which arises in testing the robustness of a de- ployed DNN to adversarial attacks. 
In the context of image classification, adversarial examples are carefully crafted perturbed images that are barely noticable and visually imperceptible, but when added to natural images, can fool the target model to mis-classify <cit.>. 
In many applications dealing with mission-critical information, the robustness of a deployed DNN to adversarial attacks is highly critical for reliability of the model, e.g., traffic sign identification for autonomous driving. 
The task of generating a universal adversarial perturbation to K natural images can be regarded as an optimization problem of the form (<ref>). 
More details on the problem formulation of  generating adversarial examples can be found in Appendix <ref>. 
Note that in general, the attacker can utilize the model evaluations and the parameters of the model to acquire its gradients <cit.>.  


Experimental Setup: 

Similar to  <cit.>,  
we apply the proposed algorithm and the baselines to generate adversarial examples to attack a 
well-trained DNN7 on the MNIST handwritten digit classification task, which achieves 99.4 test accuracy on natural examples. [<https://github.com/carlini/nn_robust_attacks>.] 
In our experiments, performed on a system with an Nvidia Tesla K80 GPU, we choose n = 10 examples from the same class, and set the 
 batch size and the number of workers to B = 5 and  m=5, respectively, for all the methods. We also use a constant step-size of 30/d, where d=900  is the image dimension, and the smoothing parameter follows μ = 𝒪(1/√(dN)), where N is the number of iterations. 
 

Experimental Results: 
Fig. <ref> depicts the attack loss versus the number of iterations, and Table <ref> shows the least l_2 distortion of the successful (universal) adversarial examples. 
It can be verified that compared to the zeroth-order methods of ZO-SGD and ZO-SVRG-Ave, the proposed method achieves significantly faster convergence and lower loss. Moreover, its convergence speed and attained loss is comparable to those of the fastest first-order methods, i.e., RI-SGD and Synchronous SGD. 
In terms of  l_2 distortion, 
the proposed method suggests better visual quality of the resulting adversarial examples than those in the previous zeroth-order methods, and its visual quality is similar to those of the first-order methods. 








Table <ref> shows the original natural images along with the resulted adversarial examples generated by different methods, 

 under the elaborated experimental setup. 








 §.§ Multi-Class Classification Tasks with Multiple Workers









The second experiment includes various multi-class classification tasks performed by multiple worker nodes in a distributed environment. We use four  different famous datasets, including COVTYPE, 
SensIT Vehicle (both ACOUSTIC and SEISMIC),  
and SENSORLESS. 
[All the datasets 
are available online at <https://www.csie.ntu.edu.tw/c̃jlin/libsvmtools/datasets/multiclass.html>. ] Details and the description of each dataset can be found in Table <ref>.  







Experimental Setup: As the base model for training, we choose a high-dimensional fully connected two-layer neural network with more than 1.69 M parameters (1.3 K and 1.3 K hidden neurons, consecutively), i.e., d> 1.69 × 10^6. 

We use PyTorch <cit.> to develop the proposed algorithm and the baselines in a distributed environment, and conduct different types of experiments on 
a system with 8 Cores of CPU and 4 Nvidia Tesla K80 GPUs. 


The number of worker nodes and the batch size are set to m=4 and B= 64, respectively, for all the methods, and the period of the first-order gradient exchanges is set to τ = 8 for the proposed method and the periodic model averaging method of RI-SGD. Moreover, a redundancy factor of μ= 0.25 is considered for the RI-SGD method. All the methods are run from the same initial points.  
For all the methods, the batch size is set to B= 64, and their 
learning rates are optimized, for each dataset. 
Furthermore, for the proposed method and the model averaging method of RI-SGD, the period of the first-order gradient exchanges is set to τ = 8. 
Finally, it should be noted that for each dataset, we have optimized the learning rates of all the methods, in order to 
have a fair comparison. 





Experimental Results: The performance of different methods in distributed training of the considered model for various datasets is compared in Fig. <ref>. It shows the training loss versus iterations, the training loss versus wall-clock time (in seconds), and the testing accuracy versus wall-clock time, for different  datasets as aforementioned. As can be verified from the figures, the proposed method significantly outperforms the zeroth-order method ZO-SGD, in terms of convergence speed, wall-clock time, and testing accuracy. 
Moreover, 
despite the high dimension of the model, the performance of the proposed method is still comparable to the first-order methods of sync-SGD and RI-SGD (which is a model averaging method with the complementary help of infused redundancy), while the proposed method benefits from lower computational complexity as discussed before. 
These experimental findings comply with our theoretical results discussed in Section <ref>. 





















 §.§ Multi-class Classification with Multiple Workers


In this experiment, we implemented our method on a single machine with 4 GPU's workers to solve a multi-class classification task. We used the COVTYPE7 dataset <cit.> and a two layer neural network with more than 800000 parameters (1000 and 800 hidden neurons, consecutively). More details on this experiment can be found in the Appendix, and are omitted here due to space limitations. Here, we  demonstrate some  results of comparing our method to the baselines.

Fig.s  <ref> and <ref> depict  the training loss versus iterations and the training loss versus wall-clock time, respectively. As can be verified from these figure, the performance of the proposed method is significantly better than the ZO-SGD method, and it is comparable to the first-order methods. These results demonstrate that despite the high dimension of the model, the performance of the proposed method is still comparable to the first-order method, in terms of the convergence speed and the wall clock time., and significantly outperforms the zeroth-order methods. 














§ CONCLUSION


In this paper, we proposed a hybrid zeroth/first-order 
distributed SGD method for solving non-convex unconstrained stochastic optimization problems 
that benefits from low communication overhead and computational complexity, and yet converges fast. 
that strikes a very good balance 
between communication overhead, computational complexity, and convergence rate. 
By theoretical analyses, we showed that the proposed algorithm reaches the same iteration complexity as the first-order distributed SGD algorithms, while it enjoys order-wisely lower computational complexity and comparable communication overhead.   Moreover, the proposed algorithm significantly outperforms the convergence rates of the existing zeroth-order algorithms, guaranteeing a convergence rate of the order of 𝒪( d / √(m N)), with a comparable computational complexity. 
Experimental results 
demonstrate the effectiveness of the proposed approach compared to various state-of-the-art methods.





icml2020























Appendix section

section.subsection



§ PROBLEM FORMULATION OF GENERATING ADVERSARIAL EXAMPLES FROM DNN






Here, we elaborate the details of  
the problem formulation for the  task of 
generating adversarial examples considered 
in Section <ref>.  


The task of generating a universal adversarial perturbation to K natural images can be regarded as an optimization problem of the form (<ref>), in which F ( x , ζ_k ) is the designed attack loss function of the k^th image, and is defined as <cit.>  

    F ( x , ζ_k ) =     c. max{ 0, f_y_k(  0.5 tanh( tanh^-1 2 a_k + x) ) - max_j ≠ y_j f_j(  0.5 tanh( tanh^-1 2 a_k + x) )  }
       + ∥ 0.5 tanh( tanh^-1 2 a_k + x) - a_k  ∥^2_2,

where ζ_k = ( a_k, y_k ) in which, a_k ∈[ -0.5,0.5 ]^d and y_k denote 
the k^th natural image and its original class label, respectively. Furthermore, for a generated adversarial example z and for each image class i=1,⋯,I, the function 
f_i ( z) outputs prediction score of the model (e.g., log-probability) in that class. Note that in the above formulation, the tanh operation is used to keep the generated adversarial example, i.e., z = 0.5 tanh( tanh^-1 2 a_k + x), in the valid image space [ -0.5,0.5 ]^d. 
Finally,  c is a regularization parameter to balance a trade-off between the adversarial attack success and the l_2 distortion of the adversarial example. Note that lower values for c result in lower l_2 distortion which suggests better visual quality of the resulting adversarial examples, while higher values for c results in higher success rate of the adversarial attack (i.e., it can perturb more natural images to be misclassified). 







ICML2020:
Table <ref> shows the original natural images along with their adversarial examples generated by different methods, under the experimental setup elaborated in Section <ref>. 







 §.§ Multi-Class Classification Tasks with Multiple Workers









The second experiment includes various multi-class classification tasks performed by multiple worker nodes in a distributed environment. We use four  different famous datasets, including COVTYPE, 
SensIT Vehicle (both ACOUSTIC and SEISMIC),  
and SENSORLESS. 
[All the datasets 
are available online at <https://www.csie.ntu.edu.tw/c̃jlin/libsvmtools/datasets/multiclass.html>. ] Details and the description of each dataset can be found in Table <ref>.  






  §.§.§ Experimental Setup


As the base model for training, we choose a high-dimensional fully connected two-layer neural network with more than 1.69 M parameters (1.3 K and 1.3 K hidden neurons, consecutively), i.e., d> 1.69 × 10^6. 

To develop the proposed algorithm as well as the baselines in a distributed environment, we use PyTorch <cit.>, and conduct different types of experiments on clusters with 4 GPUs. 
The number of worker nodes and the batch size are set to m=4 and B= 64, respectively, for all the methods, and the period of the first-order gradient exchanges is set to τ = 8 for the proposed method and the periodic model averaging method of RI-SGD. Moreover, a redundancy factor of μ= 0.25 is considered for the RI-SGD method. All the methods are run from the same initial points.  
For all the methods, the batch size is set to B= 64, and their 
learning rates are optimized, for each dataset. 
Furthermore, for the proposed method and the model averaging method of RI-SGD, the period of the first-order gradient exchanges is set to τ = 8. 
Finally, it should be noted that for each dataset, we have optimized the learning rates of all the methods, in order to 
have a fair comparison. 




  §.§.§ Experimental Results



The performance of different methods in distributed training of the considered model for various datasets is compared in Fig. <ref>. It shows the training loss versus iterations, the training loss versus wall-clock time (in seconds), and the testing accuracy versus wall-clock time, for different  datasets as aforementioned. As can be verified from the figures, the proposed method significantly outperforms the zeroth-order method ZO-SGD, in terms of convergence speed, wall-clock time, and testing accuracy. 
Moreover, 
despite the high dimension of the model, the performance of the proposed method is still comparable to the first-order methods of sync-SGD and RI-SGD (which is a model averaging method with the complementary help of infused redundancy), while the proposed method benefits from lower computational complexity as discussed before. 
These experimental findings comply with our theoretical results discussed in Section <ref>. 




















§ CONVERGENCE PROOFS
 










In order to prove 
 Theorem <ref>, we first define some auxiliary variables and functions, then we present and prove some useful lemmas which will be utilized for proving Theorem <ref>. Note that unless stated otherwise, the expectations are with respect to all the random variables.



 §.§ Preliminary Definitions


We use the following smoothing scheme 
to approximate the first-order information of a given  function f. 

For any smoothing parameter μ > 0, the co-called smoothing function of any original function f ( x) is defined as <cit.>:  

    f_μ( x) ≜𝔼_u∼ U_b[ f ( x + μu) ],

where 
U_b is the uniform distribution over the Euclidean d-dimensional  unit  ball. 


The above smoothing scheme 
exhibits several interesting features 
 that we shall use for the analysis of the convergence of the proposed method. 
 For example, as shown by Lemma 4.1(b) in <cit.>,  
for any function f ∈ C_L^1 ( ℝ^d ) (this condition is satisfied under (<ref>) in Assumption <ref>), we have

    | f_μ( x) - f ( x) | ≤μ^2 L2,

and

    ∥∇ f_μ( x) - ∇ f ( x) ∥≤μ L d2.







Here, we define three functions that will be extensively used in the analysis of the proposed algorithm.






  (a) The 
function f̃_t ( x) indicates the 
considered function at each iteration of the proposed algorithm (that is being minimized by using its stochastic gradient direction): 


    f̃_t ( x) ≜{[     f ( x) If(t,τ)=0,;           ;    f_μ( x) Otherwise. ].




  (b) Moreover, G̃_t indicates the update direction, or equivalently, the gradient estimation of the objective function, utilized at each iteration of the proposed algorithm: 

    G̃_t ( x) ≜{[          1Bm∑_i=1^m ∑_b=1^B ∇ F ( x , ζ_t+1,i,b)                                       If(t,τ)=0,;                                                 ; 1Bm∑_i=1^m ∑_b=1^B G_μ( x , ζ_t+1,i,b , v_t+1,i)                                       Otherwise, ].
 
where G_μ( x , ζ_t+1,i,b , v_t+1,i) is the zeroth-order gradient estimation and is defined as 

    G_μ(  x , ζ_t+1,i,b , v_t+1,i) ≜dμ[ F( x^t + μv_t+1,i , ζ_t+1,i,b) - F( x^t , ζ_t+1,i,b) ] v_t+1,i.




  (c) Finally, δ̃_t 
 indicates the difference between the utilized gradient estimation and the true gradient of the function considered at each iteration:  

    δ̃_t ( x) ≜∇f̃_t ( x)  -  G̃_t ( x).









 §.§ Preliminary Lemmas


First, note that according to the definition of G̃_t in (<ref>) and Assumption <ref>, using 
Lemma 4.2 in <cit.>, 
it can be shown that for all t ≥ 0, 

    𝔼_ζ_t+1 , v_t+1[ δ̃_t ( x) ] = 0,   ∀x∈ℝ^d.
 



In the following, some useful lemmas are in order.


For any L-smooth function f ( x):  
∀ t ≥ 0,  ∀x∈ℝ^d,

    f̃_t+1( x)  ≤{[          f̃_t( x) + μ^2 L2 If(t+1,τ)=0   XOR (t,τ)=0,;                           ;                   f̃_t( x)                 Otherwise. ].





If one and only one of t and t+1 divides τ, it can be followed that exactly one of the two functions f̃_t+1( .) and f̃_t(. ) is f ( .) and the other one is f_μ( .). Therefore, according to (<ref>), it can be concluded that the 
difference of these two functions at any point x∈ℝ^d is upper-bounded by μ^2 L2. Otherwise (i.e., if both t and t+1 divide τ or both t and t+1 do not divide τ), then according to (<ref>), both the associated functions are the same, and hence, their difference 
is zero at any point x∈ℝ^d.




When running Algorithm <ref> for N iterations (t=0,…,N-1), it can be shown that 
if τ > 1, the first case of the above lemma occurs exactly 2 ⌊N-1τ⌋ - 1 times. Otherwise (i.e., if τ = 1), the first case never occurs. This will be used later on in upper-bounding the error terms in the proof of the main theorem, i.e., Theorem <ref>.



For any L-smooth function f ( x): 
 ∀ t ≥ 0,  ∀x∈ℝ^d,

    ∥∇f̃_t ( x) ∥^2 ≥{[                     ∥∇ f ( x^t ) ∥^2                           If(t,τ)=0,;                                     ; 12∥∇ f ( x^t ) ∥^2   - μ^2 d^2 L^2 4                           Otherwise. ].





For the case when t satisfies (t, τ)=0, the result is obvious as we have  f̃_t ( x) = f ( x) (due to the definition in (<ref>)), and hence, ∇f̃_t ( x) = ∇ f ( x). 
For the case when (t, τ) ≠ 0, starting from (<ref>) and using the fact that ∥a∥^2 ≤ 2 ∥b∥^2 + 2 ∥a-b∥^2 ,  ∀a, b∈ℝ^d, we have

    ∥∇ f ( x) ∥^2    ≤ 2 ∥∇f̃_t ( x) ∥^2 + ∥∇f̃_t ( x) - ∇ f ( x) ∥^2, 
       ≤ 2 ∥∇f̃_t ( x) ∥^2 + μ^2 d^2 L^22,

and consequently,

    ∥∇f̃_t ( x) ∥^2 ≥12∥∇ f ( x) ∥^2 - μ^2 d^2 L^24,

which is the desired result.





Under Assumptions <ref>-<ref>, in the proposed Algorithm <ref>, the second moment of the update direction (<ref>) can be bounded as: ∀ t ≥ 0, 

    𝔼[ ∥G̃_t ( x^t  )  ∥^2 ] ≤{[                                      𝔼[ ∥∇ f ( x^t ) ∥^2 ]  + σ^2  B m                                                             If(t,τ)=0,;                                                                       ; 2 ( d+B m-1 ) B m 𝔼[ ∥∇ f ( x^t ) ∥^2 ]   +  2 d σ^2B m +μ^2 L^2 d^2 2                                                             Otherwise, ].




 




First note that according to (<ref>) and (<ref>), we have 

    𝔼[ ∥G̃_t ( x^t  ) ∥^2 ]      = 𝔼_Ω_t[ 𝔼_ζ_t+1 , v_t+1[ ∥G̃_t ( x^t  ) ∥^2 | Ω_t . ]  ] , 
        = 𝔼_Ω_t[ 𝔼_ζ_t+1 , v_t+1[ ∥G̃_t ( x^t  ) - 𝔼_ζ_t+1 , v_t+1[ G̃_t ( x^t  ) ] ∥^2  + ∥𝔼_ζ_t+1 , v_t+1[ G̃_t ( x^t  ) ] ∥^2 | Ω_t . ] ],  
        = 𝔼_Ω_t[ 𝔼_ζ_t+1 , v_t+1[ ∥G̃_t ( x^t  ) - ∇f̃_t ( x^t ) ∥^2  | Ω_t .  ]  + ∥∇f̃_t ( x^t ) ∥^2  ],  
        = 𝔼[ ∥G̃_t ( x^t  ) - ∇f̃_t ( x^t )  ∥^2 ] + 𝔼[ ∥∇f̃_t ( x^t )  ∥^2 ] .

 


Now, note that at each iteration t, one of the following two cases occurs. In the following, we expand the last equality above for each of these two cases, separately.  

Case 1: If (t, τ)=0, then 
substituting (<ref>) in (<ref>), we have 

    𝔼[ ∥G̃_t ( x^t  )  ∥^2 ]  
        =  𝔼[ ∥1Bm∑_i=1^m ∑_b=1^B ( ∇ F ( x^t , ζ_t+1,i,b) - ∇ f ( x^t ) ) ∥^2 ]  + 𝔼[ ∥∇ f ( x^t )  ∥^2 ] , 
       (a)=1B^2 m^2∑_i=1^m ∑_b=1^B 𝔼[ ∥∇ F ( x^t , ζ_t+1,i,b) - ∇ f ( x^t ) ∥^2 ] + 𝔼[ ∥∇ f ( x^t )  ∥^2 ] , 
       (b)≤σ^2Bm + 𝔼[ ∥∇ f ( x^t )  ∥^2 ] ,
 
where equation (a) results from independence of observed random samples 
ζ_t+1,i,b,  ∀ i=1,…,m, ∀ b=1,…,B and (<ref>) in Assumption <ref>, and inequality (b) is due to (<ref>) in Assumption <ref>.
 


Case 2: If (t, τ) ≠ 0,  then 
substituting (<ref>) in (<ref>) results in 

    𝔼[ ∥G̃_t ( x^t  )  ∥^2 ]      = 𝔼[ ∥1Bm∑_i=1^m ∑_b=1^B( G_μ( x^t , ζ_t+1,i,b , v_t+1,i)  - ∇f̃_μ( x^t ) ) ∥^2 ]   + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] , 
       (a)=1B^2m^2∑_i=1^m ∑_b=1^B  𝔼[ ∥G_μ( x^t , ζ_t+1,i,b , v_t+1,i) - ∇f̃_μ( x^t ) ∥^2 ] + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] 
    , 
        = 1B^2m^2∑_i=1^m 𝔼_Ω_t[ 𝔼_ζ_t+1 , v_t+1[ ∥G_μ( x^t , ζ_t+1,i,b , v_t+1,i)    - ∇f̃_μ( x^t )  ∥^2 | Ω_t . ]  ] 
          + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ]  , 
       (b)=1B^2m^2∑_i=1^m ∑_b=1^B 𝔼_Ω_t[ 𝔼_ζ_t+1 , v_t+1[ ∥G_μ( x^t , ζ_t+1,i,b , v_t+1,i) ∥^2 ] - ∥∇f̃_μ( x^t )  ∥^2  ] 
          + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] ,   
        = 1B^2m^2∑_i=1^m ∑_b=1^B 𝔼_Ω_t[ 𝔼_ζ_t+1[ 𝔼_v_t+1[ ∥G_μ( x^t , ζ_t+1,i,b , v_t+1,i) ∥^2 ] ]    - ∥∇f̃_μ( x^t )  ∥^2  ] 
          + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] , 
    
       (c)≤1B^2m^2∑_i=1^m  ∑_b=1^B 𝔼_Ω_t[ 𝔼_ζ_t+1[ 2 d  ∥∇ F ( x^t , ζ_t+1,i,b) ∥^2 + μ^2 L^2 d^22]   - ∥∇f̃_μ( x^t )  ∥^2  ] 
          + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] , 
       (d)=1B^2m^2∑_i=1^m ∑_b=1^B 𝔼_Ω_t[ 2 d  𝔼_ζ_t+1[  ∥∇ F ( x^t , ζ_t+1,i,b) -  ∇ f ( x^t )  ∥^2 ]    + 2 d   ∥∇ f ( x^t  ) ∥^2  + μ^2 L^2 d^22
        - ∥∇f̃_μ( x^t )  ∥^2   ]   + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] , 
    
       (e)≤1B^2 m^2∑_i=1^m  ∑_b=1^B 𝔼_Ω_t[ 2 d  σ^2 + 2 d   ∥∇ f ( x^t  ) ∥^2  + μ^2 L^2 d^22 - ∥∇f̃_μ( x^t )  ∥^2  ]  + 𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] , 
    
        = 2dBm 𝔼[ ∥∇ f ( x^t  ) ∥^2 ] + 2 d  σ^2Bm + μ^2 L^2 d^22Bm + Bm-1Bm𝔼[ ∥∇f̃_μ( x^t )  ∥^2 ] , 
    
       (f)≤ 2 ( d + Bm - 1 ) Bm𝔼[ ∥∇ f ( x^t  ) ∥^2 ] + 2 d  σ^2Bm + μ^2 L^2 d^22 ,
 
where ζ_t+1≜{ζ_t+1,i,b}_i = 1,…, m, b=1,…,B, v_t+1≜{v_t+1,i}_i = 1,…, m, and Ω_t≜{ζ_j,i,b, v_j,i}_j=0,…,t ,   i = 1,…, m,   b=1,…,B.
Note that equation (a) results from independence of random vectors  {ζ_t+1,i,b , v_t+1,i  , ∀ i=1,…,m, ∀ b=1,…, B } 

 and equation (<ref>). Moreover, equation (b) comes from (<ref>) and (<ref>). Inequality (c) comes from the definition of G_μ along with  Proposition 6.6. in <cit.>.  
Furthermore, equation (d) and inequality (e) come from (<ref>) and (<ref>)  in Assumption <ref>, respectively. Finally, inequality (f) comes from (<ref>) and the fact that ∥a∥^2 ≤ 2 ∥b∥^2 + 2 ∥a-b∥^2 ,  ∀a, b∈ℝ^d.








 §.§ Proof of Theorem <ref>


Having the above 
lemmas, we are now ready to prove the main theorem. 
First note that since the objective function f ( x) is L-smooth, 
its smoothing function f_μ( x) is also L-smooth <cit.>.  
Therefore, it is concluded from definition (<ref>) that the function f̃_t ( x),  ∀ t ≥ 0 is also L-smooth. As such, from the properties of smooth functions <cit.>, it is concluded that for any t ≥ 0,

    f̃_t ( x^t+1)    ≤f̃_t ( x^t) + ⟨∇f̃_t ( x^t) ,  x^t+1 - x^t⟩    + L2∥x^t+1 - x^t∥^2, 
        = f̃_t ( x^t) -  α_t ⟨∇f̃_t ( x^t) , G̃_t ⟩ + L  α_t^2 2∥G̃_t ∥^2 ,

where the  equality is due to the update equation (<ref>) in Algorithm <ref>. Moreover, combining (<ref>) and (<ref>) results in

    f̃_t+1( x^t+1)  ≤f̃_t ( x^t+1) + L μ^22.

Therefore, summing up the above two inequalities results in


    f̃_t+1( x^t+1)    ≤f̃_t ( x^t) -  α_t ⟨∇f̃_t ( x^t) , G̃_t ⟩ + L  α_t^2 2∥G̃_t ∥^2.


Now, by applying expectation to both sides of the above inequality, it is concluded that

    𝔼[  f̃_t+1( x^t+1) ] 
       ≤𝔼[f̃_t ( x^t) ] -  α_t 𝔼[⟨∇f̃_t ( x^t) , G̃_t ⟩]   + L  α_t^2 2𝔼[ ∥G̃_t ∥^2 ], 
        = 𝔼[f̃_t ( x^t) ] + L  α_t^2 2𝔼[ ∥G̃_t ∥^2 ]  -  α_t 𝔼_Ω_t[ 𝔼_ζ_t+1,v_t+1[⟨∇f̃_t ( x^t) , G̃_t ⟩| Ω_t] ] , 
       (a)=𝔼[f̃_t ( x^t) ] + L  α_t^2 2𝔼[ ∥G̃_t ∥^2 ]  -  α_t 𝔼_Ω_t[ ⟨∇f̃_t ( x^t) , 𝔼_ζ_t+1,v_t+1[ G̃_t | Ω_t] ⟩] , 
       (b)=𝔼[f̃_t ( x^t) ]  + L  α_t^2 2𝔼[ ∥G̃_t ∥^2 ]  -  α_t 𝔼_Ω_t[ ∥∇f̃_t ( x^t) ∥^2  ],

where equality (a) is due to independence of ∇f̃_t ( x^t) to the random vectors ζ_t+1 and v_t+1, and equality (b) follows from  combining (<ref>)  and (<ref>). 
Using Lemmas <ref> and <ref>, we can further bound the second and third terms in the right hand side of (<ref>), as follows.




    𝔼   [ f̃_t( x^t+1) ] 
      ≤
       {[                                                       𝔼[f̃_t ( x^t) ] -  ( α_t - L2α_t^2 ) 𝔼[ ∥∇ f ( x^t) ∥^2  ] +  L σ^2 2 B mα_t^2                                                                                                                           If(t,τ)=0,;                                                                                                                                     ; 𝔼[f̃_t ( x^t) ] -  ( 12α_t - ( d+ B m-1)L B mα_t^2 ) 𝔼[ ∥∇ f ( x^t) ∥^2  ] + μ^2 d^2 L^2 4α_t  + ( d L σ^2 B m + μ^2 d^2 L^34) α_t^2                                                                                                                           Otherwise. ].




Now, summing up the inequality (<ref>) for all t=0,…,N-1, applying
 Lemma <ref>, and using the telescopic rule result in

    𝔼[ f̃_N( x^N) ]    ≤𝔼[ f̃_0( x^0) ] - ∑_t=0^N-1β_t 𝔼[ ∥∇ f ( x^t) ∥^2 ]   + 1(τ > 1 )  ( 2 ⌊N-1τ⌋ - 1 ) μ^2 L2  + ∑_t=0^N-1 A_t, 
    
       (a)≤𝔼[ f̃_0( x^0) ] - ∑_t=0^N-1β_t 𝔼[ ∥∇ f ( x^t) ∥^2 ]  + 1(τ > 1 ) μ^2 L N τ  + ∑_t=0^N-1 A_t

where 1(τ > 1 ) is an indicator function which takes non-zero value of 1 only when its input argument is true; otherwise, it takes 0. Therefore, the third term in the right hand side of the above inequality does not exist if the period τ is set to 1.[Note that this corresponds to the case where the worker nodes use SFO estimation in all the iterations of Algorithm <ref>. Therefore, the proposed Algorithm <ref> 
reduces to the method of 
 synchronous distributed  SGD described in <cit.>.] Moreover, inequality (a) is due to the basic property of the floor function. 
Finally, for any t ≥ 0, the parameters β_t and A_t in the above inequality are defined as  

    β_t ≜{[                 α_t - L2α_t^2                    If(t,τ)=0,;                              ; 12α_t - ( d+ B m-1)L B mα_t^2                    Otherwise, ].

 and


    A_t ≜{[                                       L σ^2 2 B mα_t^2                                             If(t,τ)=0,;                                                       ; μ^2 d^2 L^2 4α_t + ( d L σ^2 B m + μ^2 d^2 L^34) α_t^2                                             Otherwise, ].


respectively. 

Using Definition (<ref>), we can write inequality (<ref>) for t=0 and t=N as follows:  

    f̃_0 ( x^0 ) - f ( x^0 ) ≤μ^2 L2,

and

    f ( x^N ) - f̃_N ( x^N ) ≤μ^2 L2,

respectively. 
Moreover, note that in the special case of τ=1, we have ( t,τ)=0,  ∀ t, and hence, according to Definition <ref>, we have f̃_t ( x) - f ( x) = 0. As such, in this case, we can further tighten the upper-bounds in (<ref>) and (<ref>) and write them as equal to  0. Consequently, combining these upper-bounds and then using Assumption <ref> results in 


    f̃_0 ( x^0 ) - f̃_N ( x^N )    ≤ f ( x^0 ) - f ( x^N ) + 1( τ > 1 ) μ^2 L, 
       ≤ f ( x^0 ) - f^∗ + 1( τ > 1 ) μ^2 L .

Taking expectation from both sides of (<ref>) and then combining it with (<ref>)  concludes that 

    ∑_t=0^N-1β_t  𝔼[ ∥∇ f ( x^t) ∥^2 ]  
    
        ≤ f ( x^0 )  -  f^∗ + 1(τ > 1 ) ( μ^2 L N τ   + μ^2 L )  + ∑_t=0^N-1 A_t


 
Now, choosing α_t ←√( B m)L √(N) , ∀ t ≥ 0 as stated in Theorem <ref>, and defining β≜min_tβ_t, it can be derived from (<ref>) that 

    β    = 12α_t - ( d+ B m-1)L B mα_t^2 |_α_t =√( B m)L √(N), 
        = √( B mN) - 2 ( d+ B m-1 )  2 L N , 
    
        > √( B m)4L √(N),

where the inequality is due to the assumption on the number of iterations stated in Theorem <ref>. 
Therefore, we have

    1β≤4 L N√( B mN).


Multiplying 1β N on both sides of inequality (<ref>), it follows that

    1N∑_t=0^N-1𝔼[ ∥∇ f ( x^t) ∥^2 ]  ≤   f ( x^0 ) - f^∗β N + 1(τ > 1 )  ( μ^2 Lβτ + μ^2 Lβ N)  + 1β N∑_t=0^N-1 A_t,

where the left hand side is the desired error term we aim to upper-bound. For this purpose, it suffices to upper-bound the last term on the right hand side. Using (<ref>), we can write

    ∑_t=0^N-1 A_t     =  L σ^2 2 B m∑_t=0, ( t,τ)=0^N-1α_t^2  + μ^2 d^2 L^2 4∑_t=0, ( t,τ) ≠ 0^N-1α_t  + ( d L σ^2 B m + μ^2 d^2 L^34) ∑_t=0, ( t,τ) ≠ 0^N-1α_t^2 , 
       (a)= L σ^2 2 L^2 N∑_t=0, ( t,τ)=0^N-1 1  + μ^2 d^2 L √( B m)4 √(N)∑_t=0, ( t,τ) ≠ 0^N-1 1     + ( d σ^2L N + μ^2 d^2 L  B m4 N) ∑_t=0, ( t,τ) ≠ 0^N-1 1 ,
 
where equality (a) comes from substituting the value of α_t. Moreover, note that the summation terms in 
the right hand side of this equality 
 can be computed 
as follows: 

    ∑_t=0, ( t,τ)=0^N-1 1 = ⌊N-1τ⌋≤Nτ,


    ∑_t=0, ( t,τ) ≠ 0^N-1 1 = ( N -1 - ⌊N-1τ⌋) ≤ N - N-1τ ,

where the inequalities above 
come from the basic properties of the floor function. 
Substituting (<ref>) and (<ref>) in (<ref>) concludes that

    ∑_t=0^N-1 A_t     =  L σ^2 2 L^2 N⌊N-1τ⌋+ μ^2 d^2 L √( B m)4 √(N)( N - 1 - ⌊N-1τ⌋)+ ( d σ^2L N + μ^2 d^2 L  B m4 N) ( N - 1 - ⌊N-1τ⌋) , 
       
     ≤ L σ^2 2 L^2 NNτ   + 1( τ > 1 ) ( N -  N-1τ) (μ^2 d^2 L √( B m)4 √(N) + d σ^2L N + μ^2 d^2 L  B m4 N)  , 
        = σ^2 2 L τ  + 1( τ > 1 ) ( N τ-1τ + 1τ) (μ^2 d^2 L √( B m)4 √(N) + d σ^2L N + μ^2 d^2 L  B m4 N) , 
        = σ^2 2 L τ  
    + 1( τ > 1 )  ( μ^2 d^2 L √( B m N)4τ-1τ
    + μ^2 d^2 L √( B m)4 √(N)τ
    + d σ^2Lτ-1τ 
    + d σ^2L N τ 
    + μ^2 d^2 L  B m4τ-1τ
       
    + μ^2 d^2 L  B m4 N τ) .
  

Now,  
combining (<ref>) and (<ref>) with (<ref>),  
it follows that 

    1N∑_t=0^N-1𝔼[ ∥∇ f ( x^t) ∥^2 ] ≤   4 L ( f ( x^0 ) - f^∗) √( B m N)  
    +   2 σ^2 √( B mN)τ
        + 1( τ > 1 )  ( 4 μ^2 L^2 N√( B m N)τ 
    + 4 μ^2 L^2√( B m N) + μ^2 d^2 L^2 τ-1τ
    + μ^2 d^2 L^2 N τ
    + 4 d σ^2√( B m N)τ-1τ
        + 4 d σ^2 N √( B m N)τ 
    + μ^2 d^2 L^2 √( B m)√(N)τ-1τ 
    + μ^2 d^2 L^2 √( B m) N √(N)τ),


Finally, choosing the smoothing parameter to be any value satisfying μ≤1√(dN) (and consequently, μ^2 ≤1d N≤1d √( B mN), where the last inequality comes from the fact that N >  B m, 
resulted from the assumption N > 16 ( d+ B m-1 )^2 B m 
 as stated in Theorem <ref>), the above inequality can be further bounded as 

    1N∑_t=0^N-1𝔼[ ∥∇ f ( x^t) ∥^2 ]  ≤   4 L ( f ( x^0 ) - f^∗) √( B m N)  
    +   2 σ^2 √( B mN)τ
        + 1(τ > 1 )  ( 4 L^2 d √( B m N)τ 
    + 4 L^2 d N √( B m N) 
    + d L^2√( B mN)τ-1τ
    +  d L^2 N √( B mN)τ
       
    + 4 d σ^2√( B m N)τ-1τ  
    + 4 d σ^2 N √( B m N)τ
    + d L^2 √( B mN)τ-1τ 
    + d L^2  N √( B mN)τ),

which completes the proof of Theorem <ref>. ▪



