Convergence of Distributed Stochastic Variance Reduced Methods without Sampling Extra Data

TL;DR

This paper analyzes the convergence of distributed stochastic variance reduced methods without extra data sampling, providing theoretical guarantees and insights into their behavior with data heterogeneity and nonconvexity.

Contribution

It offers a comprehensive convergence analysis of distributed variance reduced algorithms, including accelerated and recursive variants, under data heterogeneity and nonconvex settings.

Findings

Linear convergence for strongly convex losses.

Convergence behavior depends on data homogeneity and regularization.

Extension to nonconvex loss functions.

Abstract

Stochastic variance reduced methods have gained a lot of interest recently for empirical risk minimization due to its appealing run time complexity. When the data size is large and disjointly stored on different machines, it becomes imperative to distribute the implementation of such variance reduced methods. In this paper, we consider a general framework that directly distributes popular stochastic variance reduced methods in the master/slave model, by assigning outer loops to the parameter server, and inner loops to worker machines. This framework is natural and friendly to implement, but its theoretical convergence is not well understood. We obtain a comprehensive understanding of algorithmic convergence with respect to data homogeneity by measuring the smoothness of the discrepancy between the local and global loss functions. We establish the linear convergence of distributed versions of a family of stochastic variance reduced algorithms, including those using accelerated and recursive gradient updates, for minimizing strongly convex losses. Our theory captures how the convergence of distributed algorithms behaves as the number of machines and the size of local data vary. Furthermore, we show that when the data are less balanced, regularization can be used to ensure convergence at a slower rate. We also demonstrate that our analysis can be further extended to handle nonconvex loss functions.

Figures

Tables

Table 1. Table 1: Communication rounds and runtime of the proposed and existing algorithms for strongly convex losses (ignoring logarithmic factors in κ 𝜅 \kappa ) to reach ϵ italic-ϵ \epsilon -accuracy. Algorithms with an asterisk are proposed/analyzed in this paper. Here, N 𝑁 N is the total number of input data points, n 𝑛 n is the number of worker machines, and κ 𝜅 \kappa is the condition number of the global loss function f 𝑓 f .

Algorithm	Communication Rounds	Runtime	Assumptions
DSVRG [LLMY17]	$\left(1+\kappa /\left(N/n\right)\right)\log \left(1/ϵ\right)$	$\left(N/n+\kappa \right)\log \left(1/ϵ\right)$	extra data
DASVRG [LLMY17]	$\left(1+\sqrt{\kappa /\left(N/n\right)}\right)\log \left(1/ϵ\right)$	$\left(N/n+\sqrt{\kappa N/n}\right)\log \left(1/ϵ\right)$	extra data
Dist. AGD	$\sqrt{\kappa }\log \left(1/ϵ\right)$	$\left(N/n\right)\sqrt{\kappa }\log \left(1/ϵ\right)$	none
ADMM	$\kappa \log \left(1/ϵ\right)$	$\left(N/n\right)\kappa \log \left(1/ϵ\right)$	none
SCOPE [ZXS+17]	$\kappa \log \left(1/ϵ\right)$	$\left(N/n+\kappa \log \kappa \right)\kappa \log \left(1/ϵ\right)$	uniform regularization
pSCOPE [ZZLL18]	$\log \left(1/ϵ\right)$	$\left(N/n+\kappa \right)\log \left(1/ϵ\right)$	good partition
D-SVRG*	$\log \left(1/ϵ\right)$	$\left(N/n+\kappa \right)\log \left(1/ϵ\right)$	distributed restricted smoothness
D-SARAH*	$\log \left(1/ϵ\right)$	$\left(N/n+\kappa \right)\log \left(1/ϵ\right)$	distributed restricted smoothness
D-MiG*	$\left(1+\sqrt{\kappa /\left(N/n\right)}\right)\log \left(1/ϵ\right)$	$\left(N/n+\sqrt{\kappa N/n}\right)\log \left(1/ϵ\right)$	distributed smoothness
D-RSVRG*	$\kappa \log \left(1/ϵ\right)$	$\left(N/n\right)\kappa \log \left(1/ϵ\right)$	large regularization
D-RSVRG*	$\log \left(1/ϵ\right)$	$\left(N/n+\kappa \right)\log \left(1/ϵ\right)$	small regularization

Equations

$$\underset{x\in {\mathbb{R}}^d}{\min }f(x):=\frac{1}{N}\sum _{z\in \mathcal{M}}\ell (x;z),$$

$$f(x):=\frac{1}{n}\sum _{k=1}^n{f}_k(x),$$

$$x^{\ast }:=\underset{x\in {\mathbb{R}}^d}{\mathrm{argmin}}f(x),f^{\ast }:=f(x^{\ast }).$$

$$\mathbb{E}\left[f({\tilde{x}}^{t+1})-f^{\ast }\right]<\frac{\sigma -2c}{2(\sigma -3c)}\mathbb{E}\left[f({\tilde{x}}^t)-f^{\ast }\right].$$

$$v_k^{t,s+1}=\{\begin{array}{ll}\nabla {\ell }_z(y_k^{t,s+1})-\nabla {\ell }_z({\tilde{x}}^t)+\nabla f({\tilde{x}}^t)& \text{SVRG}\\ \nabla {\ell }_z(y_k^{t,s+1})-\nabla {\ell }_z(y_k^{t,s})+v_k^{t,s}& \text{SARAH}\end{array}$$

$$\mathbb{E}\left[{\Vert \nabla f({\tilde{x}}^{t+1})\Vert }^2\right]<\frac{1}{2-8c^2/{\sigma }^2}\mathbb{E}\left[{\Vert \nabla f({\tilde{x}}^t)\Vert }^2\right].$$

$$x_k^{t,s+1}=\underset{x}{\mathrm{argmin}}\left\{\frac{1}{2\eta }{\Vert x-x_k^{t,s}\Vert }^2+⟨v_k^{t,s},x⟩+g(x)\right\}.$$

$$v_k^{t,s+1}=\nabla {\ell }_z(y_k^{t,s+1})-\nabla {\ell }_z({\tilde{x}}^t)+\nabla f({\tilde{x}}^t)+{\mu }_k(y_k^{t,s+1}-{\tilde{x}}^t),$$

$$\mathbb{E}\left[f({\tilde{x}}^{t+1})-f^{\ast }\right]\le$$

$$\cdot \mathbb{E}\left[f({\tilde{x}}^t)-f^{\ast }\right],$$

$$\mathcal{O}\left((N/n+{\zeta }^{-2}\bar{\kappa }){\zeta }^{-1}\min \left\{\kappa +\frac{\mu }{\sigma },\frac{1}{\max \{1-\mu /\sigma ,0\}}\right\}\log (1/ϵ)\right),$$

$$\frac{1}{Tm}\sum _{t=0}^{T-1}\sum _{s=0}^{m-1}\mathbb{E}\left[{\Vert \nabla f(y_{k(t)}^{t,s})\Vert }^2\right]\le \frac{2}{\eta Tm}\left(f({\tilde{x}}^0)-f^{\ast }\right),$$

$${\xi }^{s+1}={\xi }^s-\eta \left[\nabla {\ell }_z\left(y^s\right)-\nabla {\ell }_z\left(y^0\right)+\nabla f\left(y^0\right)\right],$$

$${\xi }^{s+1}=A{\xi }^s+B{w}^s,$$

$$w^s=\left[\begin{array}{c}\nabla {\ell }_z\left(y^s\right)-\nabla {\ell }_z\left(x^{\ast }\right)\\ \nabla {\ell }_z\left(x^{\ast }\right)-\nabla {\ell }_z\left(y^0\right)+\nabla f\left(y^0\right)\end{array}\right].$$

$$S_j(\xi ,w)=[{\xi }^{\top },w^{\top }]X_j[{\xi }^{\top },w^{\top }{]}^{\top },$$

$$V({\xi }^{s+1})\le {\rho }^{2}V({\xi }^k)+\sum _{j=1}^J{\lambda }_j{S}_j({\xi }^s,w^s)$$

$$\left[\begin{array}{cc}{A}^{\top }A-{\rho }^2{I}_d& A^{\top }B\\ B^{\top }A& B^{\top }B\end{array}\right]-\sum _{j=1}^J{\lambda }_j{X}_j⪯0.$$

$${\bar{X}}_1=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],{\bar{X}}_2=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right],{\bar{X}}_3=\left[\begin{array}{ccc}0& -1& -1\\ -1& 0& 0\\ -1& 0& 0\end{array}\right],$$

$$\mathbb{E}\left[S_1\right]$$

$$\mathbb{E}\left[S_2\right]$$

$$\mathbb{E}\left[S_3\right]$$

$$\{\begin{array}{l}\mathbb{E}\left[S_1\right]\le 2L\mathbb{E}\left[f(y^s)-f^{\ast }\right]+2cL\mathbb{E}\left[{\Vert y^s-x^{\ast }\Vert }^2\right];\\ \mathbb{E}\left[S_2\right]\le 4L\mathbb{E}\left[f(y^0)-f^{\ast }\right]+c(4L+2c)\mathbb{E}\left[{\Vert y^0-x^{\ast }\Vert }^2\right];\\ \mathbb{E}\left[S_3\right]\le -2\mathbb{E}\left[f(y^s)-f^{\ast }\right]+3c\mathbb{E}\left[{\Vert y^s-x^{\ast }\Vert }^2\right]+c\mathbb{E}\left[{\Vert y^0-x^{\ast }\Vert }^2\right].\end{array}$$

$${\lambda }_3-L{\lambda }_1-(2L{\lambda }_1+3{\lambda }_3)c/\sigma >0$$

$$\left[\begin{array}{ccc}0& {\lambda }_3-\eta & {\lambda }_3-\eta \\ {\lambda }_3-\eta & {\eta }^2-{\lambda }_1& {\eta }^2\\ {\lambda }_3-\eta & {\eta }^2& {\eta }^2-{\lambda }_2\end{array}\right]⪯0$$

$$\mathbb{E}\left[f(y^{+})-f^{\ast }\right]\le {\left({\lambda }_3-L{\lambda }_1-(2L{\lambda }_1+3{\lambda }_3)c/\sigma \right)}^{-1}\left[\frac{1}{\sigma m}+\frac{c}{\sigma }\left((4L+2c){\lambda }_2+{\lambda }_3\right)+2L{\lambda }_2\right]\mathbb{E}\left[f(y^0)-f^{\ast }\right],$$

$$\left[\begin{array}{ccc}0& 0& 0\\ 0& -{\eta }^2& {\eta }^2\\ 0& {\eta }^2& -{\eta }^2\end{array}\right]⪯0.$$

$$\mathbb{E}\left[f(y^{+})-f^{\ast }\right]\le$$

$$\mathbb{E}\left[f({\tilde{x}}^{t+1})-f^{\ast }\right]$$

$$\le \left(1-\frac{1}{2}\cdot \frac{\sigma -4c}{\sigma -3c}\right)\mathbb{E}\left[f({\tilde{x}}^t)-f^{\ast }\right].$$

$${\ell }_{{\mu }_k}^t(x;z)=\ell (x;z)+\frac{{\mu }_k}{2}{\Vert x-{\tilde{x}}^t\Vert }^2.$$

$$f_i^t(x)=\frac{1}{|{\mathcal{M}}_i|}\sum _{z\in {\mathcal{M}}_i}{\ell }_{{\mu }_k}^t(x;z)=f_i(x)+\frac{{\mu }_k}{2}{\Vert x-{\tilde{x}}^t\Vert }^2,$$

Full text

1 Introduction

Empirical risk minimization arises frequently in machine learning and signal processing, where the objective function is the average of losses computed at different data points. Due to the increasing size of data, distributed computing architectures, which assign the learning task over multiple computing nodes, are in great need to meet the scalability requirement in terms of both computation power and storage space. In addition, distributed frameworks are suitable for problems where there are privacy concerns to transmit and store all the data in a central location, a scenario related to the nascent field of federated learning [KMR15]. It is, therefore, necessary to develop distributed optimization frameworks that are tailored to solving large-scale empirical risk minimization problems with desirable communication-computation trade-offs, where the data are stored disjointly over different machines.

One prevalent model of distributed systems is the so-called master/slave model, where there is a central parameter server to coordinate the computation and information exchange across different worker machines. Due to the low per-iteration cost, a popular solution is distributed stochastic gradient descent (SGD) [RRWN11], where the parameter server aggregates gradients from each worker and does mini-batch gradient updates. However, distributed SGD is not communication-efficient and requires lots of communication rounds to converge, which partially diminishes the benefit of distribution. On the other hand, recent breakthroughs in developing stochastic variance reduced methods have made it possible to achieve fast convergence and small per-iteration cost at the same time, such as the notable SVRG [JZ13] algorithm. Yet, distributed schemes of such variance reduced methods that are both practical and theoretically sound are much less developed.

This paper focuses on a general framework of distributed stochastic variance reduced methods, to be presented in Alg. 1, which is natural and friendly to implement in a master/slave model. On a high level, SVRG-type algorithms [JZ13] contain inner loops for parameter updates via variance-reduced SGD, and outer loops for global gradient and parameter updates. Our general framework assigns outer loops to the parameter server, and inner loops to worker machines. The parameter server collects gradients from worker machines and then distributes the global gradient to each machine. Each worker machine then runs the inner loop independently in parallel using variance reduction techniques, which might be different when distributing different algorithms, and returns the updates to the parameter server at the end. Per iteration, two communication rounds are required: one communication round is used to average the parameter estimates, and the other is used to average the gradients, which is the same as distributed synchronous SGD. However, the premise is that by performing more efficient local computation using stochastic variance reduced methods, the algorithm converges in fewer iterations and is therefore more communication-efficient.

Due to the simplicity of this framework, similar methods have been implemented in several works [KMR15, DG16, RKR*+*16], and have achieved great empirical success. Surprisingly, a complete theoretical understanding of its convergence behavior is still missing at large. Moreover, distributed variants using accelerated variance reduction methods are not developed. The main analysis difficulty is that the variance-reduced gradient of each worker is no longer an unbiased gradient estimator when sampling from re-used local data.

To ease this difficulty, several variants of distributed SVRG, e.g. [LLMY17, Sha16, WWS17] have been proposed with performance guarantees, which try to bypass the biased gradient estimation issue by emulating the process of i.i.d. sampling from the global data using some complicated random data re-allocation protocol, which requires sampling extra data with or without replacement. These procedures lead to unnecessary data waste and potential privacy leakage, and can be cumbersome and difficult to implement in practice.

Consequently, a natural question arises: can we provide a mathematical analysis to the convergence of the natural framework of distributed stochastic variance reduced methods, under some simple and intuitive metric?

1.1 Contributions of This Paper

This paper provides a convergence analysis of a family of naturally distributed stochastic variance reduced methods under the framework described in Alg. 1, for both convex and nonconvex loss functions. By using different variance reduction schemes at the worker machines, we study distributed variants of three representative algorithms in this paper: SVRG [JZ13], SARAH employing recursive gradient updates [NLST17, NvDP*+*19], and MiG employing accelerated gradient updates [ZSC18]. Our methodology can be extended to study other variants in a similar fashion. The contributions of this paper are summarized below.

• •

We suggest a simple and intuitive metric called distributed smoothness to gauge data balancedness among workers, defined as the smoothness of the difference $f_{k}-f$ between the local loss function $f_{k}$ and the global loss function $f$, which is the average of the local loss functions. The metric is deterministic, easy-to-compute, applies for arbitrary dataset splitting, and is shown to play a critical role in the convergence analysis.

• •

We establish the linear convergence of distributed D-SVRG, D-SARAH, and D-MiG under strongly convex losses, as long as the distributed smoothness parameter is smaller than a constant fraction of the strong convexity parameter $\sigma$, e.g. $\sigma/4$, where the fraction might change for different algorithms. Our bounds capture the phenomenon that the convergence rate improves as the local loss functions become more similar to the global loss function, by reducing the distributed smoothness parameter. Furthermore, the run time complexity exhibits the so-called “linear speed-up” property in distributed computing, where the complexity depends on the local data size, instead of the global data size, which typically implies an improvement by a factor of $n$, where $n$ is the number of machines.

• •

When the local data are highly unbalanced, the distributed smoothness parameter becomes large, which implies that the algorithm might diverge. We suggest regularization as an effective way to handle this situation, and show that by adding larger regularization to machines that are less distributed smooth, one can still ensure linear convergence in a regularized version of D-SVRG, called D-RSVRG, though at a slower rate of convergence.

• •

More generally, the notion of distributed smoothness can also be used to establish the convergence under nonconvex losses. We demonstrate this through the convergence analysis of D-SARAH in the nonconvex setting.

1.2 Related Work

Distributed optimization is a classic topic [BPC*+*11, BT89], yet recent trends in data-intensive applications are calling for new developments with a focus on communication and computation efficiency. Examples of deterministic optimization methods include DANE [SSZ14], AIDE [RKR*+*16], DiSCo [ZL15], GIANT [WRKXM18], CoCoA [SFC*+*18], CEASE [FGW19], one-shot averaging [ZWLS10, ZWD12], etc.

Many stochastic variance reduced methods have been proposed recently, for example, SAG [SLRB17], SAGA [DBLJ14], SVRG [JZ13], SDCA [SSZ13], MiG [ZSC18], Katyusha [AZ17], Catalyst [LMH15], SCOPE [ZXS*+*17, ZZLL18], SARAH [NLST17], SPIDER [FLLZ18], SpiderBoost [WJZ*+*18], to name a few. Several previous works have studied distributed variants of SVRG. For example, the D-SVRG algorithm has been empirically studied before in [RKR*+*16, KMR15] without a theoretical convergence analysis. The pSCOPE algorithm [ZZLL18] is also a variant of distributed SVRG, and its convergence is studied under an assumption called good data partition in [ZZLL18], which is hard to interpret and verify in practice. The SCOPE algorithm [ZXS*+*17] is similar to the regularized variant D-RSVRG of D-SVRG under large regularization, however our analysis is much more refined by allowing different regularizations to different local workers with respect to the distributed smoothness of local data, and gracefully degenerates to the unregularized case when the distributed smoothness is benign. The general framework of distributed variance-reduced methods covering SARAH and MiG and various loss settings in this paper has not been studied before. For conciseness, Table 1 summarizes the communication and runtime complexities of the most relevant algorithms to the current paper.111Since in a master/slave model, all worker machines perform computations in parallel, the runtime is dictated by the worst-case computation complexity per worker.

There are also a lot of recent efforts on reducing the communication cost of distributed GD/SGD by gradient quantization [AGL*+*17, BWAA18, SFD*+*14, WXY*+*17], gradient compression and sparsification [AHJ*+*18, LZZ*+*17, LHM*+*18, WWLZ18, TLY*+*18]. In comparison, we communicate the exact gradient, and it is an interesting future direction to combine gradient compression schemes in distributed variance reduced stochastic gradient methods. Another line of works [YYLS17, MR16, LCCC20] seek to adopt variance reduction techniques in the decentralized setting, where each worker is only able to communicate with its neighbors over a network. In contrast, our work focuses on the master/slave model where each worker communicates with a parameter server.

1.3 Paper Organization

The rest of this paper is organized as follows. Section 2 presents the problem setup and a general framework of distributed stochastic optimization with variance-reduced local updates. Section 3 presents the convergence guarantees of D-SVRG, D-SARAH and D-MiG under appropriate distributed smoothness assumptions. Section 4 introduces regularization to D-SVRG to handle unbalanced data when distributed smoothness does not hold. Section 5 presents extensions to nonconvex losses for D-SARAH. Section 6 provides an outline to the analysis and the rest of the proofs are deferred to the appendix. Section 7 presents numerical experiments to corroborate the theoretical findings. Finally, we conclude in Section 8.

2 Problem Setup

Suppose we have a data set $\mathcal{M}=\{z_{1},\cdots,z_{N}\}$, where $z_{j}\in\mathbb{R}^{p}$ is the $j$th data point for $j=1,...,N$, and $N$ is the total number of data points. In particular, we do not make any assumptions on their statistical distribution. Consider the following empirical risk minimization problem

[TABLE]

where $x\in\mathbb{R}^{d}$ is the parameter to be optimized and $\ell:\mathbb{R}^{d}\times\mathbb{R}^{p}\mapsto\mathbb{R}$ is the sample loss function. For brevity, we use $\ell_{z}(x)$ to denote $\ell(x;z)$ throughout the paper.

In a distributed setting, where the data are distributed to $n$ machines or workers, we define a partition of the data set $\mathcal{M}$ as $\mathcal{M}=\bigcup_{k=1}^{n}\mathcal{M}_{k}$, where $\mathcal{M}_{i}\bigcap\mathcal{M}_{k}=\varnothing$, $\forall i\neq k$. The $k^{\text{th}}$ worker, correspondingly, is in possession of the data subset $\mathcal{M}_{k}$, $1\leq k\leq n$. We assume there is a parameter server (PS) that coordinates the parameter sharing among the workers. The sizes of data held by each worker machine is $N_{k}=|\mathcal{M}_{k}|$. When the data is split equally, we have $N_{k}=N/n$. The original problem (1) can be rewritten as minimizing the following objective function:

[TABLE]

where $f_{k}(x)=\frac{1}{(N/n)}\sum_{z\in\mathcal{M}_{k}}\ell_{z}(x)$ is the local loss function at the $k^{\text{th}}$ worker machine.222It is straightforward to state our results under unequal data splitting with proper rescaling.

2.1 A General Distributed Framework

Alg. 1 presents a general framework for distributed stochastic variance reduced methods, which assigns the outer loops of an SVRG-type algorithm [JZ13, NLST17, ZSC18] to PS and the inner loops to local workers. By using different variance reduction schemes at the worker machines (i.e. LocalUpdate), we obtain distributed variants of different algorithms. On a high level, the framework alternates between local computation by individual workers in parallel (i.e. Line 5), and global information sharing coordinated by the PS (i.e. Lines 8-12).

• •

The local worker conducts local computation LocalUpdate based on the current estimate $\tilde{x}^{t}$, the global gradient $\nabla f(\tilde{x}^{t})$, and its local data $f_{k}(\cdot)$; in this paper, we are primarily interested in local updates using stochastic variance-reduced gradients. A little additional information about the previous update is needed when employing acceleration, which will be specified in Alg. 3.

• •

The PS randomly selects from the local estimates $y_{k}^{t+}$ from all workers in each round, which is set as the global estimate $\tilde{x}^{t+1}$, then computes the global gradient $\nabla f(\tilde{x}^{t+1})$ by pulling the local gradient $\nabla f_{k}(\tilde{x}^{t+1})$. In total, each iteration requires two rounds of communications.

Remark 1**.**

A careful reader may suggest that there are better ways to set the global estimate $\tilde{x}^{t+1}$. One way is to let the PS randomly select a worker to compute a local estimate $y_{k}^{t+}$ and set other workers in an idle mode, which saves local computation. Another way is to set $\tilde{x}^{t+1}=\frac{1}{n}\sum_{k=1}^{n}y_{k}^{t+}$ as the average of all local updates, which may potentially improve the performance. The specific choice in Alg. 1 allows us to provide theoretical analysis for all the LocalUpdate rules considered in this paper.

2.2 Assumptions

Throughout, we invoke one or several of the following standard assumptions of the loss function in the convergence analysis.

Assumption 1** (Smoothness).**

The sample loss $\ell_{z}(\cdot)$ is $L$-smooth, i.e., the gradient of $\nabla\ell_{z}(\cdot)$ is $L$-Lipschitz for all $z\in\mathcal{M}$.

Assumption 2** (Convexity).**

The sample loss $\ell_{z}(\cdot)$ is convex for all $z\in\mathcal{M}$.

Assumption 3** (Strong Convexity).**

The empirical risk $f(\cdot)$ is $\sigma$-strongly convex.

When $f$ is strongly convex, the condition number of $f$ is defined as $\kappa:=L/\sigma$. Denote the unique minimizer and the optimal value of $f(x)$ as

[TABLE]

As it turns out, the smoothness of the deviation $f_{k}-f$ between the local loss function $f_{k}$ and the global loss function $f$ plays a key role in the convergence analysis, as it measures the balancedness between local data in a simple and intuitive manner. We refer to this as the “distributed smoothness”. In some cases, a weaker notion called restricted smoothness is sufficient, which is defined below.

Definition 1** (Restricted Smoothness).**

A differentiable function $f:\mathbb{R}^{d}\mapsto\mathbb{R}$ is called $c$-restricted smooth with regard to $x^{*}$ if $\left\lVert\nabla f(x^{*})-\nabla f(y)\right\rVert\leq c\left\lVert x^{*}-y\right\rVert$, for all $y\in\mathbb{R}^{d}$.

The restricted smoothness, compared to standard smoothness, fixes one of the arguments to $x^{*}$, and is therefore a much weaker requirement. The following assumption quantifies the distributed smoothness using either restricted smoothness or standard smoothness.

Assumption 4a** (Distributed Restricted Smoothness).**

The deviation $f-f_{k}$ is $c_{k}$-restricted smooth with regard to $x^{*}$ for all $1\leq k\leq n$.

Assumption 4b** (Distributed Smoothness).**

The deviation $f-f_{k}$ is $c_{k}$-smooth for all $1\leq k\leq n$.

It is straightforward to check that $c_{k}\leq L$ for all $1\leq k\leq n$. If all the data samples are generated following certain statistical distribution in an i.i.d. fashion, one can further link the distributed smoothness to the local sample size $N/n$, where $c_{k}$ decreases with the increase of $N/n$, see e.g. [SSZ14, FGW19] for further discussion.

Remark 2**.**

We provide a toy example to illustrate the difference between Assumptions 4a and 4b. Let $x\in\mathbb{R}$ and $\ell_{z}(x)=L_{\delta(z)}(x)+\frac{\sigma}{2}x^{2}$, where $L_{\delta(z)}(\cdot)$ denotes Huber loss with $\delta(z)>0$. We set $\delta(z)=0.99$ for half of $z\in\mathcal{M}$, and $1$ for the rest. Therefore, $f(x)=\frac{1}{2}(L_{0.99}(x)+L_{1}(x))+\frac{\sigma}{2}x^{2}$ with $x^{*}=0$. Now, consider some local loss function $f_{k}(x)=(\frac{1}{2}+\epsilon_{k})L_{0.99}(x)+(\frac{1}{2}-\epsilon_{k})L_{1}(x)+\frac{\sigma}{2}x^{2}$ for some $\epsilon_{k}$. We can verify that Assumption 4a holds with $c_{k}=|\epsilon_{k}|/100$ while Assumption 4b holds with $c_{k}=|\epsilon_{k}|$. Therefore, the distributed smoothness parameter can be much smaller when invoking the weaker Assumption 4a.

3 Convergence in the Strongly Convex Case

In this section, we describe three variance-reduced routines for LocalUpdate used in Alg. 1, namely SVRG [JZ13], SARAH [NLST17, NvDP*+*19], and MiG [ZSC18], and analyze their convergence when $f(\cdot)$ is strongly convex, respectively.

3.1 Distributed SVRG (D-SVRG)

The LocalUpdate routine of D-SVRG is described in Alg. 2. Theorem 1 provides the convergence guarantee of D-SVRG as long as the distributed restricted smoothness parameter is small enough.

Theorem 1** (D-SVRG).**

Suppose that Assumptions 1, 2 and 3 hold, and Assumption 4a holds with $c_{k}\leq c<\sigma/4$. With $m=O(\kappa(1-4c/\sigma)^{-2})$ and proper step size $\eta$, the iterates of D-SVRG satisfy

[TABLE]

The communication and runtime complexities of finding an $\epsilon$-optimal solution (in terms of function value) are $\mathcal{O}\left(\zeta^{-1}\log(1/\epsilon)\right)$ and $\mathcal{O}\left((N/n+\zeta^{-2}\kappa)\zeta^{-1}\log(1/\epsilon)\right)$ respectively, where $\zeta=1-4c/\sigma$.

Theorem 1 establishes the linear convergence of function values in expectation for D-SVRG, as long as the parameter $c$ is sufficiently small, e.g. $c<\sigma/4$. From the expressions of communication and runtime complexities, it can be seen that the smaller $c$, the faster D-SVRG converges – suggesting that the homogeneity of distributed data plays an important role in the efficiency of distributed optimization. When $c$ is set such that $c/\sigma$ is bounded above by a constant smaller than $1/4$, i.e. $\zeta=\mathcal{O}(1)$, the runtime complexity becomes $\mathcal{O}\left(\left(N/n+\kappa\right)\log(1/\epsilon)\right)$, which improves the counterpart of SVRG $\mathcal{O}\left(\left(N+\kappa\right)\log(1/\epsilon)\right)$ in the centralized setting.

Remark 3**.**

The above LocalUpdate routine corresponds to the so-called Option II (w.r.t. setting $y_{k}^{t,+}$ as uniformly at random selected from previous updates) specified in [JZ13]. Under similar assumptions, we also establish the convergence of D-SVRG using Option I, where the output $y_{k}^{t,+}$ is set as $y_{k}^{t,m}$. In addition, D-SVRG still converges linearly in the absence of Assumption 2. We leave these extensions in the appendix.

3.2 Distributed SARAH (D-SARAH)

The LocalUpdate of D-SARAH is also described in Alg. 2, which is different from SVRG in the update of stochastic gradient $v_{k}^{t,s}$, by using a recursive formula proposed in [NLST17]. Theorem 2 provides the convergence guarantee of D-SARAH as long as the distributed restricted smoothness parameter is small enough.

Theorem 2** (D-SARAH).**

Suppose that Assumptions 1, 2 and 3 hold, and Assumption 4a holds with $c_{k}\leq c<\sqrt{2}\sigma/4$. With $m=O(\kappa(1-2\sqrt{2}c/\sigma)^{-2})$ and proper step size $\eta$, the iterates of D-SARAH satisfy

[TABLE]

The communication and runtime complexities of finding an $\epsilon$-optimal solution (in terms of gradient norm) are $\mathcal{O}\left(\zeta^{-1}\log(1/\epsilon)\right)$ and $\mathcal{O}\left((N/n+\zeta^{-2}\kappa)\zeta^{-1}\log(1/\epsilon)\right)$ respectively, where $\zeta=1-2\sqrt{2}c/\sigma$.

Theorem 2 establishes the linear convergence of the gradient norm in expectation for D-SARAH, as long as the parameter $c$ is small enough. Similar to D-SVRG, a smaller $c$ leads to faster convergence of D-SARAH. When $c$ is set such that $c/\sigma$ is bounded above by a constant smaller than $\sqrt{2}/4$, the runtime complexity becomes $\mathcal{O}\left(\left(N/n+\kappa\right)\log(1/\epsilon)\right)$, which improves the counterpart of SARAH $\mathcal{O}\left(\left(N+\kappa\right)\log(1/\epsilon)\right)$ in the centralized setting. In particular, Theorem 2 suggests that D-SARAH may allow a larger $c$, compared with D-SVRG, to guarantee convergence.

3.3 Distributed MiG (D-MiG)

The LocalUpdate of D-MiG is described in Alg. 3, which is inspired by the inner loop of the MiG algorithm [ZSC18], a recently proposed accelerated variance-reduced algorithm. Compared with D-SVRG and D-SARAH, D-MiG uses additional information from the previous updates, according to Line 3-6 in Alg. 3. Theorem 3 provides the convergence guarantee of D-MiG, as long as the distributed smoothness parameter is small enough.

Theorem 3** (D-MiG).**

Suppose that Assumptions 1, 2 and 3 hold, and Assumption 4b holds with $c_{k}\leq c<\sigma/8$. Let $w=(1+\eta\sigma)/(1+3\eta c)$. With $m=O(N/n)$ and proper step size $\eta$, the iterates of D-MiG achieve an $\epsilon$-optimal solution within a communication complexity of $\mathcal{O}((1+\sqrt{\kappa/(N/n)})\log(1/\epsilon))$ and runtime complexity of $\mathcal{O}((N/n+\sqrt{\kappa N/n})\log(1/\epsilon))$.

Theorem 3 establishes the linear convergence of D-MiG under standard smoothness of $f-f_{k}$, in order to fully harness the power of acceleration. While we do not make it explicit in the theorem statement, the time complexity of D-MiG also decreases as $c$ gets smaller. Furthermore, the runtime complexity of D-MiG is smaller than that of D-SVRG/D-SARAH when $\kappa=\Omega(N/n)$.

Remark 4**.**

Theorem 3 continues to hold for regularized empirical risk minimization, where the loss function is given as $F(x)=f(x)+g(x)$, and $g(x)$ is a convex and non-smooth regularizer. In this case, line 11 in Alg. 3 is changed to

[TABLE]

4 Regularization Helps Unbalanced Data

So far, we have established the convergence when the distributed smoothness is not too large. While it may be reasonable in certain settings, e.g. in a data center where one has control over how to distribute data, it is increasingly harder to satisfy when the data are generated locally and heterogeneous across workers. However, when such conditions are violated, the algorithms might diverge. In this situation, adding a regularization term might ensure the convergence, at the cost of possibly slowing down the convergence rate. We consider regularizing the local gradient update of D-SVRG in Alg. 2 as

[TABLE]

where the last regularization term penalizes the proximity between the current iterates $y_{k}^{t,s+1}$ and the reference point $\tilde{x}^{t}$, where $\mu_{k}>0$ is the regularization parameter employed at the $k^{\text{th}}$ worker. We have the following theorem.

Theorem 4** (Distributed Regularized SVRG (D-RSVRG)).**

Suppose that Assumptions 1, 2 and 3 hold, and Assumption 4a holds with $c_{k}<(\sigma+\mu_{k})/4$. Let $\mu=\min_{1\leq k\leq n}\mu_{k}$. With proper $m$ and step size $\eta$, there exists some constant $0\leq\nu<1$ such that the iterates of D-RSVRG satisfy

[TABLE]

and the runtime complexity of finding an $\epsilon$-optimal solution is bounded by

[TABLE]

where $\zeta=1-4c/(\sigma+\mu)$ and $\bar{\kappa}=(L+\mu)/(\sigma+\mu)$.

Compared with Theorem 1, Theorem 4 relaxes the assumption $c_{k}<\sigma/4$ to $c_{k}<(\sigma+\mu_{k})/4$, which means that by inserting a larger regularization $\mu_{k}$ to local workers that are not distributed smooth, i.e. those with large $c_{k}$, one can still guarantee the convergence of D-RSVRG. However, increasing $\mu$ leads to a slower convergence rate: a large $\mu=8L$ leads to an iteration complexity $\mathcal{O}(\kappa\log(1/\epsilon))$, similar to gradient descent. Compared with SCOPE [ZXS*+*17] which requires a uniform regularization $\mu>L-\sigma$, our analysis applies tailored regularization to local workers, and potentially allows much smaller regularization to guarantee the convergence, since $c_{k}$’s can be much smaller than the smoothness parameter $L$.

5 Convergence in the Nonconvex Case

In this section, we extend the convergence analysis of D-SARAH to handle nonconvex loss functions, since SARAH-type algorithms are recently shown to achieve near-optimal performances for nonconvex problems [WJZ*+*18, NvDP*+*19, FLLZ18]. As a modification that eases the analysis, we make every worker return $y_{k}^{t+}=y_{k}^{t,m}$ in line 9 of Alg. 2. Our result is summarized in the theorem below.

Theorem 5** (D-SARAH for non-convex losses).**

Suppose that Assumption 1 and Assumption 4b hold with $c_{k}\leq c$. With the step size $\eta\leq\frac{2}{L\left(1+\sqrt{1+8m+4m(m-1)c^{2}/L^{2}}\right)}$, D-SARAH satisfies

[TABLE]

where $k(t)$ is the agent index selected in the $t^{\text{th}}$ round for parameter update, i.e. $\tilde{x}^{t+1}=y^{t+}_{k(t)}$ (c.f. line 8 of Alg. 1). To find an $\epsilon$-optimal solution, the communication complexity is $\mathcal{O}\left(1+(\sqrt{n/N}+c/L)L/\epsilon\right)$, and the runtime complexity is $\mathcal{O}(N/n+(\sqrt{N/n}+N/n\cdot c/L)L/\epsilon)$ by setting $m=\Theta(N/n)$.

Theorem 5 suggests that D-SARAH converges as long as the step size is small enough. Furthermore, a smaller $c$ allows a larger step size $\eta$, and hence faster convergence. To gain further insights, assuming i.i.d. data at each worker, by concentration inequalities it is known that $c/L=\mathcal{O}(\sqrt{\log(N/n)/(N/n)})$ under mild conditions [MBM18], and consequently, the runtime complexity of finding an $\epsilon$-accurate solution using D-SARAH is $\mathcal{O}(N/n+L\sqrt{\log(N/n)N/n}/\epsilon)$. This is comparable to the best known result $\mathcal{O}(N+L\sqrt{N}/\epsilon)$ for the centralized SARAH-type algorithms in the nonconvex setting [NvDP*+*19, FLLZ18, WJZ*+*18] up to logarithmic factors, where the data size is replaced from $N$ to the size of local data $N/n$ – demonstrating again the benefit of data distribution.

6 Proofs of Main Theorems

In this section, we outline the convergence proofs of D-SVRG (Theorem 1), D-RSVRG (Theorem 4), D-SARAH in the strongly convex (Theorem 2) and nonconvex (Theorem 5) settings, while leaving the details to Appendix D. The convergence proof of D-MiG (Theorem 3) is delegated to the supplemental materials due to space limits. Throughout this section, we simplify the notations $y_{k}^{t,s}$ and $v_{k}^{t,s}$ by dropping the $t$ superscript and the $k$ subscript, whenever the meaning is clear, since it is often sufficient to analyze the convergence of a specific worker $k$ during a single round.

6.1 D-SVRG (Theorem 1)

We generalize the analysis of SVRG using the dissipativity theory in [HWL18] to the analysis of D-SVRG, which might be of independent interest. Setting $\xi^{s}=y^{s}-x^{*}$, we can write the local update of D-SVRG via the following linear time-invariant system [HWL18]:

[TABLE]

where $z$ is selected uniformly at random from local data points $\mathcal{M}_{k}$, or equivalently,

[TABLE]

with $A=I_{d}$, $B=\begin{bmatrix}-\eta I_{d}&-\eta I_{d}\end{bmatrix}$, and

[TABLE]

Here, $I_{d}$ is the identity matrix of dimension $d$.

Dissipativity theory characterizes how “the inputs” $w^{s}$, $s=0,1,2,\ldots$ drive the internal energy stored in the “states” $\xi^{s}$, $s=0,1,2,\ldots$ via an energy function $V:\mathbb{R}^{d}\mapsto\mathbb{R}_{+}$ and a supply rate $S:\mathbb{R}^{d}\times\mathbb{R}^{2d}\mapsto\mathbb{R}$. For our purpose, it is sufficient to choose the energy function as $V(\xi)=\|\xi\|_{2}^{2}$. By setting $S(\xi,w)=\sum_{j=1}^{J}\lambda_{j}S_{j}(\xi,w)$ as the supply rate, where

[TABLE]

we have

[TABLE]

for some $\rho\in(0,1]$ as long as there exist non-negative scalars $\lambda_{j}$ such that

[TABLE]

In fact, by left multiplying $[\xi_{k}^{\top},w_{k}^{\top}]$ and right multiplying $[\xi_{k}^{\top},w_{k}^{\top}]^{\top}$ to (8), we recover (7).

To invoke the dissipativity theory for convergence analysis, we take two steps. First, we capture properties of the objective function such as strong convexity and co-coercivity using the supply rate functions by selecting proper matrices $X_{j}$ (c.f. Lemma 1). Next, we invoke (8) to determine a valid combination of supply rate functions to guarantee the dissipativity inequality (c.f. Lemma 2), which entails the rate of convergence.

To proceed, we start by defining the supply rate functions in (6). Since $\xi^{s}\in\mathbb{R}^{d}$ and $w^{s}\in\mathbb{R}^{2d}$, we will write $X_{j}$ as $X_{j}=\bar{X}_{j}\otimes I_{d}$, where $\bar{X}_{j}\in\mathbb{R}^{3\times 3}$. Following [HWL18], we consider the supply rates characterized by the following matrices:

[TABLE]

which correspond to the supply rates:

[TABLE]

The following lemma, proved in Appendix B.1, bounds the supply rates when $z$ is drawn from the local data $\mathcal{M}_{k}$.

Lemma 1**.**

Suppose that Assumption 1, 2, 3 and 4a hold with $c_{k}\leq c$. For the supply rates defined in (9), we have

[TABLE]

The next lemma, proved in Appendix B.2, details the dissipativity inequality when (8) holds when setting $\rho=1$.

Lemma 2**.**

Suppose that Assumption 1, 2, 3 and 4a hold with $c_{k}\leq c$. If there exist non-negative scalars $\lambda_{j},\ j=1,2,3$, such that

[TABLE]

and

[TABLE]

hold. Then D-SVRG satisfies

[TABLE]

where the final output $y^{+}$ is selected from $y^{1},\cdots,y^{m}$ uniformly at random.

We are now ready to prove Theorem 1. Set $\lambda_{1}=\lambda_{2}=2\eta^{2}$, and $\lambda_{3}=\eta$. We have (10) holds with the step size $\eta<\dfrac{\sigma-3c}{L(2\sigma+4c)}$, and (11) holds since

[TABLE]

Applying Lemma 2 with the above choice of parameters, (12) can be written as

[TABLE]

When $c<\sigma/4$, by choosing $\eta=(1-4c/\sigma)(40L)^{-1}$, $m=160\kappa(1-4c/\sigma)^{-2}$, a convergence rate no more than $\nu:=1-\frac{1}{2}\cdot\frac{\sigma-4c}{\sigma-3c}$ is obtained. Therefore, after line 8 of Alg. 1, we have for D-SVRG,

[TABLE]

To obtain an $\epsilon$-optimal solution in terms of function value, we need $\mathcal{O}(\xi^{-1}\log(1/\epsilon))$ communication rounds, where $\zeta=1-4c/\sigma$. Per round, the runtime complexity at each worker is $\mathcal{O}(N/n+m)=\mathcal{O}(N/n+\zeta^{-2}\kappa)$, where the first term corresponds to evaluating the batch gradient over the local data in parallel, and the second term corresponds to evaluating the stochastic gradients in the inner loop. Multiplying this with the communication rounds leads to the overall runtime.

6.2 D-RSVRG (Theorem 4)

Consider an auxiliary sample function at the $t^{\text{th}}$ round and the $k^{\text{th}}$ worker as

[TABLE]

This leads to the auxiliary local and global loss functions, respectively,

[TABLE]

for $1\leq i\leq n$, and

[TABLE]

Moreover, we have

[TABLE]

which means that D-RSVRG performs in exactly the same way as the unregularized D-SVRG with the auxiliary loss functions $\ell_{\mu_{k}}^{t}$ in the $t^{\text{th}}$ round. Note that $\ell_{\mu_{k}}^{t}$ is $(\mu_{k}+L)$-smooth and that $f^{t}$ is $(\mu_{k}+\sigma)$-strongly-convex, while the restricted smoothness of the $k^{\text{th}}$ worker, $f^{t}-f_{k}^{t}=f-f_{k}$ remains unchanged. Applying Theorem 1, we have

[TABLE]

where $f^{t*}$ is the optimal value of $f^{t}$, when Assumptions 1, 2 and 3 hold, and Assumption 4a holds, as long as $c_{k}<(\sigma+\mu_{k})/4$, along with proper $m$ and step size $\eta$.

However, the definition of the regularized loss functions $\ell^{t}_{\mu}$ and ${f}^{t}$ rely on $\tilde{x}^{t}$, which changes over different rounds. Our next step is to relate the descent of $f^{t}$ to $f$. To this end, we have

[TABLE]

where the second line uses (13), the third line uses (14), and the last line follows from $f(\tilde{x}^{t})=f^{t}(\tilde{x}^{t})$.

We can continue to bound $f(\tilde{x}^{t})-f^{t*}$ in two manners. First, $f(\tilde{x}^{t})-f^{t*}=f^{t}(\tilde{x}^{t})-f^{t*}>\frac{1}{2(L+\mu)}\left\lVert\nabla f(\tilde{x}^{t})\right\rVert^{2}>\frac{\sigma}{L+\mu}(f(\tilde{x}^{t})-f^{*})$. On the other hand, we have $f(\tilde{x}^{t})-f^{t*}\geq f(\tilde{x}^{t})-f^{t}(x^{*})=f(\tilde{x}^{t})-f(x^{*})-\frac{\mu}{2}\left\lVert\tilde{x}^{t}-x^{*}\right\rVert^{2}\geq(1-\mu/\sigma)(f(\tilde{x}^{t})-f^{*})$. Thus (4) follows immediately by combining the above two bounds.

6.3 D-SARAH in the Strongly Convex Case (Theorem 2)

We motivate the convergence analysis by citing the following lemma from [NLST17], which bounds the sum of gradient norm of iterations on a specific worker during $t^{\text{th}}$ round:

Lemma 3**.**

[NLST17*]**

Suppose that Assumption 1 holds, then*

[TABLE]

We note that the first term on RHS can be neglected when $m$ is large; the second term measures the effect of biased gradient estimator and can thus be controlled with the step size $\eta$ and distributed smoothness; the third term can be dropped when $\eta\leq 1/L$. A careful analysis leads to the following theorem, whose proof can be found in Appendix C.1.

Theorem 6**.**

Suppose that Assumption 1, 2, 3 and 4a hold with $c_{k}\leq c$. Then D-SARAH satisfies

[TABLE]

When $c<\frac{\sqrt{2}}{4}\sigma$, we can choose $\eta=\dfrac{2(1-8c^{2}/\sigma^{2})}{(9-8c^{2}/\sigma^{2})L}$ and $m=2\kappa\dfrac{9-8c^{2}/\sigma^{2}}{(1-8c^{2}/\sigma^{2})^{2}}$ in Theorem 6, leading to the following convergence rate:

[TABLE]

Consequently, following similar discussions as D-SVRG, the communication complexity of finding an $\epsilon$-optimal solution is $\mathcal{O}(\zeta^{-1}\log(1/\epsilon))$, and the runtime complexity is $\mathcal{O}\left((N/n+\zeta^{-2}\kappa)\zeta^{-1}\log(1/\epsilon)\right)$, where $\zeta=1-2\sqrt{2}c/\sigma$.

6.4 D-SARAH in the Nonconvex Case (Theorem 5)

The convergence analysis in the nonconvex case is also based upon Lemma 3. Due to lack of convexity, the tighter bound of the second term on RHS adopted in the analysis of Theorem 6 is not available, so a smaller step size $\eta$ is needed to cancel out the second term and the third term, i.e., to make

[TABLE]

Formally speaking, we have the following result, proved in Appendix C.2.

Theorem 7**.**

Suppose that Assumption 1 and 4b hold with $c_{k}\leq c$. By setting the step size

[TABLE]

For a single outer loop of D-SARAH, it satisfies:

[TABLE]

By setting $\tilde{x}^{t+1}=y^{m}$, from the above theorem, we have

[TABLE]

Hence, with $T$ outer loops, we have

[TABLE]

where $k(t)$ is the worker index selected in the $t^{\text{th}}$ round for parameter update. The communication complexity to achieve an $\epsilon$-optimal solution is

[TABLE]

with the choice $\eta=\Theta\left(\frac{1}{L\sqrt{m+m^{2}c^{2}/L^{2}}}\right)$. Per round, the runtime complexity at each worker is $\mathcal{O}(N/n+m)$. By choosing $m=\mathcal{O}(N/n)$, we achieve the runtime complexity

[TABLE]

7 Numerical Experiments

Thought the focus of this paper is theoretical, we illustrate the performance of the proposed distributed stochastic variance reduced algorithms in various settings as a proof-of-concept.

7.1 Logistic regression in the strongly convex setting

Consider $\ell_{2}$-regularized logistic regression, where the sample loss is defined as

[TABLE]

with the data $z_{i}=(a_{i},b_{i})\in\mathbb{R}^{d}\times\{\pm 1\}$. We evaluate the performance on the gisette dataset [GGBHD05] and the rcv1 dataset [LYRL04] by splitting the data equally to all workers. We scale the data according to $\max_{i\in[N]}\left\lVert a_{i}\right\rVert^{2}=1$, so that the smoothness parameter is estimated as $L=1/4+\lambda$. We choose $\lambda=N^{-0.5}$, $N^{-0.75}$ and $N^{-1}$ respectively to illustrate the performance under different condition numbers. We use the optimality gap, defined as $f(\tilde{x}^{t})-f^{*}$, to illustrate the convergence behavior.

For D-SVRG and D-SARAH, the step size is set as $\eta=1/(2L)$. For D-MiG, although the choice of $w$ in the theory requires knowledge of $c$, we simply ignore it and set $w=1+\eta\sigma$, $\theta=1/2$ and the step size $\eta=1/(3\theta L)$ to reflect the robustness of the practical performance to parameters. We further use $\tilde{x}^{t+1}=\frac{1}{n}\sum_{k=1}^{n}y_{k}^{t+}$ at the PS, which provides better empirical performance than the random selection rule in Alg. 1, as seen in Fig. 3. For D-AGD, the step size is set as $\eta=1/L$ and the momentum parameter is set as $\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$. Following [JZ13], which sets the number of inner loop iterations as $m=2N$, we set $m\approx 2N/n$ to ensure the same number of total inner iterations. We note that such parameters can be further tuned to achieve better trade-off between communication cost and computation cost in practice.

Fig. 1 illustrates the optimality gap of various algorithms with respect to the number of communication rounds with 4 local workers under different conditioning, and Fig. 2 shows the corresponding result with different numbers of local workers when $\lambda=N^{-1}$. The distributed stochastic variance-reduced algorithms outperform distributed AGD significantly. In addition, D-MiG outperforms D-SVRG and D-SARAH when the condition number is large. Fig. 3 shows the optimality gap after a fixed $20$ communication rounds with respect to the number of agents on the gisette dataset with $\lambda=N^{-1}$. We observe that the performance of distributed variance-reduced theorems degenerates as the number of agents grows, due to the shrinking size of local dataset, which leads to a larger distributed smoothness parameter $c$.

7.2 Dealing with Unbalanced data.

We justify the benefit of regularization by evaluating the proposed algorithms under unbalanced data allocation. We assign 50%, 30%, 19.9%, 0.1% percent of data to four workers, respectively, and set $\lambda=N^{-1}$ in the logistic regression loss (15). To deal with unbalanced data, we perform the regularized update, given in (3), on the worker with the least amount of data, and keep the update on the rest of the workers unchanged. A similar regularized update can be conceived for D-SARAH and D-MiG, resulting in regularized variants, D-RSARAH and D-RMiG. While our theory does not cover them, we still evaluate their numerical performance. We properly set $\mu$ according to the amount of data on this worker as $\mu=0.1/(0.1\%\cdot N)^{0.5}$. We set the number of iterations at workers $m=2N$ on all agents. Fig. 4 shows the optimality gap with respect to the number of communication rounds for all algorithms. It can be seen that all unregularized methods fail to converge, and the regularized algorithms still converge, verifying the role of regularization in addressing unbalanced data. It is also worth mentioning that the regularization can be flexibly imposed depending on the local data size, rather than homogeneously across all workers.

7.3 Performance of D-SARAH in the nonconvex setting

We follow the same setting as [WJZ*+*18] to evaluate D-SARAH and Distributed Gradient Descent (D-GD) on the gisette dataset with a nonconvex sample loss function:

[TABLE]

which consists of the logistic loss and a non-convex regularizer, where $x_{j}$ is the $j$th entry of $x$. The smoothness parameter of $\ell_{\mathrm{ncvx}}(x;z_{i})$ can be estimated as $L=1/4+2\lambda$. Fig. 5 plots the squared norm the gradient $\|\nabla f(\tilde{x}^{t})\|^{2}$ of D-SARAH and D-GD with respect to the number of communication rounds. It can be seen that D-SARAH achieves a much lower gradient norm than D-GD with the same number of communication rounds.

8 Conclusions

In this paper, we have developed a convergence theory for a family of distributed stochastic variance reduced methods without sampling extra data, under a mild distributed smoothness assumption that measures the discrepancy between the local and global loss functions. Convergence guarantees are obtained for distributed stochastic variance reduced methods using accelerations and recursive gradient updates, and for minimizing both strongly convex and nonconvex losses. We also suggest regularization as a means of ensuring convergence when the local data are unbalanced and heterogeneous. We believe the analysis framework is useful for studying distributed variants of other stochastic variance-reduced methods such as Katyusha [AZ17], and proximal variants such as [XZ14].

Acknowledgements

The results in this paper have been presented in part at the 2019 NeuriPS Workshop on Federated Learning for Data Privacy and Confidentiality. The work of S. Cen was partly done when visiting Microsoft Research Asia. The work of S. Cen and Y. Chi is supported in part by National Science Foundation under the grants CCF-1806154, CCF-1901199, and CCF-2007911, Office of Naval Research under the grants N00014-18-1-2142 and N00014-19-1-2404, and Army Research Office under the grant W911NF-18-1-0303.

Appendix A Preliminary

We first establish a lemma which will be useful later.

Lemma 4**.**

When Assumptions 1,2 and one of the distributed smoothness (Assumption 4a or 4b) hold, we have

[TABLE]

for any $\tilde{x}$, where the expectation is evaluated over $z$, and $D_{f}(x_{1},x_{2})=f(x_{1})-f(x_{2})-\left\langle\nabla f(x_{2}),x_{1}-x_{2}\right\rangle$.

Proof.

Given $f$ is $L$-smooth and convex, the Bregman divergence $D_{f}(x_{1},x_{2})$ is $L$-smooth and convex as a function of $x_{1}$. When Assumptions 1 and 2 hold, we have

[TABLE]

Averaging the above inequality over $z\in\mathcal{M}_{k}$ gives

[TABLE]

To further bound the left-hand side, Assumption 4a allows us to compare $D_{f}$ and $D_{f_{k}}$:

[TABLE]

Following similar arguments, using Assumption 4b we obtain a tighter bound by replacing $x^{*}$ with any $\tilde{x}$. In particular, setting $\tilde{x}=(x_{1}+x_{2})/2$ we have

[TABLE]

Combining the above estimates with (16) proves the lemma. ∎

Appendix B Proof for D-SVRG

B.1 Proof of Lemma 1

For $\mathbb{E}\left[S_{1}\right]$, we apply Lemma 4 directly:

[TABLE]

where the inequality follows from the definition of $D_{f}(y^{s},x^{*})$, and $\nabla f(x^{*})=0$. For $\mathbb{E}\left[S_{2}\right]$, we have

[TABLE]

where the first inequality is due to $\left\lVert a+b\right\rVert^{2}\leq 2\left\lVert a\right\rVert^{2}+2\left\lVert b\right\rVert^{2}$, the second inequality follows from evaluating the expectation and Assumption 4a, and the last step uses Lemma 4 again.

For $\mathbb{E}\left[S_{3}\right]$, we have

[TABLE]

where the first inequality is obtained by applying Assumption 2, Cauchy-Schwarz inequality and Assumption 4a.

B.2 Proof of Lemma 2

Setting $\rho=1$ in (8), we have that it becomes equivalent to (11), and consequently the dissipation inequality (7) holds. In view of Lemma 1, it can be written as

[TABLE]

Since $f$ is $\sigma$-strongly convex, $\left\lVert y^{s}-x^{*}\right\rVert^{2}<\frac{2}{\sigma}\left(f(y^{s})-f^{*}\right)$, then the above

[TABLE]

where $\gamma_{1}=(2\lambda_{3}-2L\lambda_{1}-(4L\lambda_{1}+6\lambda_{3})c/\sigma)$ and $\gamma_{2}=c\left[(4L+2c)\lambda_{2}+\lambda_{3}\right]$ are introduced as short-hand notations. In addition, $\gamma_{1}>0$ by assumption. Telescoping the above inequality by summing over $s=0,\ldots,m-1$, we have

[TABLE]

Note that the choice of $y^{+}$ implies

[TABLE]

Therefore,

[TABLE]

We obtain the final result by substituting $\mathbb{E}\left[\left\lVert y^{0}-x^{*}\right\rVert^{2}\right]\leq\frac{2}{\sigma}\mathbb{E}\left[f\left(y^{0}\right)-f^{*}\right]$ into the above inequality.

B.3 Convergence of D-SVRG without Assumption 2

When Assumption 2 does not hold, we can still use similar arguments as the proof of Theorem 1 and establish the convergence of D-SVRG, though at a slower rate. Using the same supply rates (9), Lemma 1 can be modified as below.

Lemma 5**.**

Suppose that Assumption 1, 3 and 4a hold. For the supply rates defined in (9), we have

[TABLE]

Proof.

With $L$-smoothness of $\ell_{z}$ we have the following estimate:

[TABLE]

So we have

[TABLE]

and

[TABLE]

The estimate of $\mathbb{E}\left[S_{3}\right]$ is identical to that in Lemma 1. ∎

Following the same process in the proof of Lemma 2 in Appendix B.2, we have the following inequality with proper choices of $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$:

[TABLE]

By summing the inequality and letting $\lambda_{1}=\lambda_{2}=2\eta^{2}$, $\lambda_{3}=\eta$, we have

[TABLE]

Therefore, with $1-2L^{2}\sigma^{-1}\eta-3c\sigma^{-1}>0$ the following convergence bound can be established:

[TABLE]

By choosing $\eta=\dfrac{1-4c\sigma^{-1}}{40\kappa L}$, $m=\dfrac{160\kappa^{2}}{(1-4c\sigma^{-1})^{2}}$, we get a convergence rate no more than $1-\dfrac{1}{2}\cdot\dfrac{\sigma-4c}{\sigma-3c}$. Hence the overall time complexity to find an $\epsilon$-optimal solution is

[TABLE]

where $\zeta=1-4c/\sigma$.

B.4 Convergence of D-SVRG with Option I

Another option for the output of the inner loops of SVRG, or the output of the local workers, is to output the last iterates, i.e. $y_{k}^{t+}=y_{k}^{t,m}$, which is called “Option 1” in [JZ13]. Here, we establish the convergence of D-SVRG using Option I in the following theorem.

Theorem 8** (D-SVRG with Option I).**

Suppose that Assumptions 1, 2 and 3 hold, and Assumption 4a holds with $c<\sigma/2$. With sufficiently large $m$ and sufficiently small step size $\eta$, there exists $0\leq\nu<1$ such that

[TABLE]

Theorem 8 indicates that the iterates $\tilde{x}^{t}$ of D-SVRG with Option I converge to the minimizer $x^{*}$ linearly in expectation as long as $c$ is sufficiently small. The proof is outlined in Section 6.1. By taking $m\to\infty$ and $\eta\to 0$, the rate approaches $\nu:=\frac{c}{2\sigma-3c}$, which suggests the algorithm admits faster convergence with the decrease of $c$, as expected.

Following [HWL18], we consider the following four supply rates:

[TABLE]

We have the following lemma which is proved at the end of this subsection.

Lemma 6**.**

Suppose that Assumptions 1, 2, 3 and 4a hold. For the supply rates defined in (17), we have

[TABLE]

Therefore, by choosing $P=I_{d}$, $\lambda_{1}=2\eta^{2}$, $\lambda_{2}=\eta-L\eta^{2}$, $\lambda_{3}=\eta^{2}$, $\lambda_{4}=L\eta^{2}$, $\rho^{2}=1-2\sigma(\eta-L\eta^{2})$, the condition (8) holds:

[TABLE]

This immediately leads to the inequality (7) which reads as:

[TABLE]

Let $\tilde{\rho}^{2}=(1-(2\sigma-3c)(\eta-L\eta^{2})+cL\eta^{2})$. Telescoping the inequality over $t=0,1,\cdots,m-1$ leads to

[TABLE]

Note that when $m\to\infty$ and $\eta\to 0$, the rate becomes $\nu:=\frac{c}{2\sigma-3c}$, hence we need $c<\sigma/2$ to get a rate $\nu<1$. We have

[TABLE]

Proof of Lemma 6.

The following inequalities can be viewed as combinations of standard inequalities in convex optimization (co-coercivity, etc) and the characterization of restricted smoothness.

[TABLE]

and

[TABLE]

$\mathbb{E}\left[S_{3}\right]\leq 0$ is simply the restatement of $L$-smoothness of $\ell(\cdot,z)$, $z\in\mathcal{M}$.

[TABLE]

∎

Appendix C Proof for D-SARAH

C.1 Proof of Theorem 6

To begin, we cite a supporting lemma from [NLST17].

Lemma 7**.**

[NLST17*]**

Suppose Assumption 1 and 2 hold and $\eta<2/L$, then*

[TABLE]

We also present a new lemma below with the proof given in Appendix C.3.

Lemma 8**.**

The update rule of D-SARAH satisfies

[TABLE]

Proof of Theorem 6.

By combining Lemmas 7 and 8, we have

[TABLE]

By Assumption 4a, we have

[TABLE]

Therefore, combining (19) and (C.1), we have

[TABLE]

Substituting it into Lemma 3 gives

[TABLE]

Since $f$ is $\sigma$-strongly convex, we have $4c^{2}\left\lVert y^{s}-x^{*}\right\rVert^{2}\leq\frac{4c^{2}}{\sigma^{2}}\left\lVert\nabla f(y^{s})\right\rVert^{2}$. Denote $y^{+}$ as the local update which is selected from $y^{0},\cdots,y^{m-1}$ uniformly at random. We have

[TABLE]

Since $\tilde{x}^{t+1}$ is randomly chosen from the local outputs $\{y_{k}^{t+},\ 1\leq k\leq n\}$, we have

[TABLE]

∎

C.2 Proof of Theorem 7

Recall Lemma 3. The theorem follows if

[TABLE]

The rest of this proof is thus dedicated to show the above inequality. Note that

[TABLE]

where the second inequality follows from Lemma 8 and Assumption 4b, and the third inequality follows from $y^{0}-y^{s}=\sum_{j=1}^{s}y^{j}-y^{j-1}$, and the last line follows from the definition. The $L$-smoothness of $\ell_{z}$ implies that

[TABLE]

So we have

[TABLE]

Therefore, with $0<\eta\leq\frac{2}{L\left(1+\sqrt{1+8(m-1)+4m(m-1)c^{2}/L^{2}}\right)}$, we have $m(m-1)c^{2}\eta^{2}+2L^{2}\eta^{2}(m-1)-(1-L\eta)\leq 0$ and the proof is finished.

C.3 Proof of Lemma 8

First, we write

[TABLE]

Let $\mathscr{F}_{s}$ denote the $\sigma$-algebra generated by all random sample selections in sub-iteration $0,\cdots,s-1$. We have

[TABLE]

where the second equality follows from

[TABLE]

Taking expectation over $\mathscr{F}_{s}$ gives

[TABLE]

Hence telescoping the above equality we obtain the claimed result.

Appendix D Proof for D-MiG (Theorem 3)

As earlier, we simplify the notations $y_{k}^{t,s}$, $x_{k}^{t,s}$ and $v_{k}^{t,s}$ by dropping the superscript $t$ and the subscript $k$. In this section we deal with the non-smooth target function $F(x)=f(x)+g(x)$ as mentioned in Remark 4, where $g$ is a convex and non-smooth function known to all agents. The analysis is done by carefully adapting the proof for the centralized algorithm (i.e. [ZSC18, Section B.1.]) with the distributed smoothness assumption. We impose the following constraint on the step size $\eta$:

[TABLE]

We restate the inequalities (8) and (9) [ZSC18] and change notations to match our context:

[TABLE]

where $\tilde{\nabla}=\nabla\ell_{z}(y^{s-1})-\nabla\ell_{z}(\tilde{x})+\nabla f(\tilde{x})$, with $z$ randomly selected from $\mathcal{M}_{k}$ and $u\in\mathbb{R}^{d}$ is an arbitrary vector. Following the $L$-smoothness argument in [ZSC18], (23) leads to

[TABLE]

By plugging in the constraint (21), we have

[TABLE]

By combining (22),(23),(24), [ZSC18, Lemma 3] and then taking expectation over the choice of random sample $z$, we have

[TABLE]

We further split the term $\mathbb{E}\left[\left\langle\nabla f(y^{s-1})-\tilde{\nabla},x^{s}-u\right\rangle\right]$ as

[TABLE]

where the first inequality is due to Cauchy-Schwarz inequality and Assumption 4b, with $\beta>0$ satisfying $\frac{1-\theta}{\theta}=\frac{L}{\beta}$, and the last inequality is obtained by combining

[TABLE]

with Lemma 4 under Assumption 4b, along with the inequality $ab\leq(a^{2}+b^{2})/2$. By substituting the inequality into (25) we have the following result.

[TABLE]

Note that the choice of $\beta$ cancels both $\frac{1-\theta}{\theta}\left\langle\nabla f(y^{s-1}),\tilde{x}-y^{s-1}\right\rangle$ and $\left(\frac{L\theta}{2(1-\theta)}+\frac{c}{2}\right)\mathbb{E}\left[\left\lVert x^{s}-x^{s-1}\right\rVert^{2}\right]$. By rearranging the above inequality, we can cancel the term $f(y^{s-1})$. We further use $-g(x^{s})\leq\frac{1-\theta}{\theta}g(\tilde{x})-\frac{1}{\theta}g(y^{s})$, which leads to

[TABLE]

Note that $\left\lVert\tilde{x}-x^{s-1}\right\rVert^{2}\leq 2\left\lVert\tilde{x}-x^{*}\right\rVert^{2}+2\left\lVert x^{s-1}-x^{*}\right\rVert^{2}$. By setting $u=x^{*}$ and using $\sigma$-strongly-convexity of $F$ we have

[TABLE]

To simplify the analysis we impose another constraint $\theta\leq 1/2$:

[TABLE]

Let $w=\frac{1+\eta\sigma}{1+3c\eta}$. Multiplying (26) by $w^{s-1}$ and then summing over $s$, we have

[TABLE]

Adding the superscript $t$ and the subscript $k$ back and applying Jensen’s inequality to the definition of $y_{k}^{t+}$, we get

[TABLE]

Averaging the inequality over $k=1,\cdots,n$, we have

[TABLE]

Case I

$\frac{m(\sigma-3c)}{L}\leq 3/4$. By setting $\theta=\sqrt{\frac{m(\sigma-3c)}{3L}}\leq 1/2$ and $\eta=\frac{1}{3L\theta+c}$ we satisfy the step size constraint (21). We aim to show that

[TABLE]

We have $w=\frac{1+\eta\sigma}{1+3c\eta}\leq 1+\eta\sigma-3c\eta$. Note that $L\theta=\sqrt{\frac{Lm(\sigma-3c)}{3}}\geq\sqrt{\frac{\sigma(\sigma-3c)}{3}}>c$, so $\eta\leq\frac{1}{4L\theta}$. Moreover, with $\sigma>8c$ we have

[TABLE]

Hence,

[TABLE]

where $\zeta=\frac{m}{4L\theta}(\sigma-3c)$. With $\zeta\leq\theta-\frac{4c}{\sigma}\theta^{2}\leq 1/2$ we have $\left(1-\zeta\right)\left(1+\zeta/m\right)^{m}\leq 1$. So we have

[TABLE]

Therefore we have

[TABLE]

The convergence rate over $T$ rounds of communication is $w^{-Tm}=\left(1+\Theta\left(\frac{1}{\sqrt{3\kappa m}}\right)\right)^{-Tm}$, so the communication complexity is $T=\mathcal{O}\left(\sqrt{\kappa/m}\log(1/\epsilon)\right)$. With the choice of $m=\Theta(N/n)$ the runtime complexity is $\mathcal{O}((N/n+m)T)=\mathcal{O}(\sqrt{\kappa N/n}\log(1/\epsilon))$.

Case II

$\frac{m(\sigma-3c)}{L}>3/4$. We set $\theta=1/2$ and $\eta=\frac{1}{3L\theta+c}$. We have $\eta\sigma=\frac{1}{3\kappa/2+c/\sigma}\leq 2/3$, thus

[TABLE]

and $w^{m}\geq 1+\frac{3}{5}(\sigma-3c)\eta m\geq 1+\frac{3m(\sigma-3c)}{5(3L/2+L/8)}>1+1/4$. So we have

[TABLE]

Therefore,

[TABLE]

This implies that $\mathcal{O}(\log(1/\epsilon))$ rounds of communication is sufficient to find an $\epsilon$-accurate solution, and that the runtime complexity is $\mathcal{O}(N/n\log(1/\epsilon))$.

Appendix E Discussions on Distributed Smoothness

In this paper, we have established that distributed variance reduced methods admit simple convergence analysis under the distributed smoothness for all worker machines, as long as the parameter $c$ is smaller than a constant fraction of $\sigma$, the strong convexity parameter of the global loss function $f$. In this section, we will show that the distributed smoothness can be guaranteed for many practical loss functions as long as the local data size is sufficiently large and homogeneous across machines. This is as expected, since SVRG and SARAH rely heavily on exploiting data similarities to reduce the variance.

Since the distributed smoothness only examines the gradient information of $f_{k}-f$, it can be applied to loss functions with non-smooth gradients, e.g. Huber loss. However, for simplicity of exposition, we limit our focus to the case when the sample loss $\ell_{z}(\cdot)$ is second-order differentiable and demonstrate the smoothness of $f-f_{k}$ via uniform concentration of Hessian matrices.

For simplicity, we consider the quadratic loss case, which allows us to compare with existing result for the DANE algorithm [SSZ14], which is a communication-efficient approximate Newton-type algorithm. Assume $\ell(\cdot,z)$ is quadratic for all $z$. Recall the following result on the concentration of Hessian matrices from [SSZ14].

Lemma 9**.**

[SSZ14*]**

If $0\preceq\nabla^{2}\ell_{z}(x)\preceq L$ holds for all $z$, then with probability at least $1-\delta$ over the samples, for all $x$,*

[TABLE]

Moreover, the iteration complexity of DANE is given by the theorem below.

Theorem 9**.**

[SSZ14*]**

If $0\preceq\nabla^{2}\ell_{z}(x)\preceq L$ holds for all $z$ and $\sigma\preceq\nabla^{2}f(x)\preceq L$, then with probability exceeding $1-\delta$, DANE needs*

[TABLE]

iterations to find an $\epsilon$-optimal solution.

By Theorem 9, [SSZ14] claims that when the local data size of every machine is sufficiently large, namely $N/n=\Omega\left(\kappa^{2}\log(dn)\right)$, DANE can find a desired $\epsilon$-optimal with $\mathcal{O}(\log(1/\epsilon))$ iterations and thus communication-efficient. Note that at this local data size, according to Lemma 9, it is sufficient to establish $c=\mathcal{O}(\sigma)$, which satisfies the convergence requirement of D-SVRG and D-SARAH. Consequently, the proposed D-SVRG and D-SARAH converges at the same iteration complexity as DANE, that is $\mathcal{O}(\log(1/\epsilon))$. Recall that DANE requires its local routines to be solved exactly. In contrast, our results formally justifies that SVRG and SARAH can be safely used as an inexact local solver for DANE without losing its performance guarantees, thus answering an open question in [RKR*+*16].