Stability and Generalization of the Decentralized Stochastic Gradient Descent^††thanks: This work is sponsored in part by National Key R&D Program of China (2018YFB0204300), and the National Science Foundation of China (No. 61932001 and 61906200).

In this paper, we establish the first theoretical result on the stability and generalization of the D-SGD. We elaborate on our contributions below.

1. 1.

   Stability of D-SGD: We provide the uniform stability of D-SGD in the general convex, strongly convex, and nonconvex cases. Our theory shows that besides the learning rate, data size, and iteration number, the stability and generalization of D-SGD are also dependent on the connected graph structure. To the best of our knowledge, our result is the first theoretical stability guarantee for D-SGD. In the general convex setting, we also present the stability of D-SGD in terms of the ergodic average instead of the last iteration for the excess generalization analysis.

2. 2.

   Computational errors for D-SGD with convexity and projection: We consider more general schemes of D-SGD, that is, D-SGD with projection. In the previous work (Ram, Nedić, and Veeravalli 2010b), to get the convergence rate, the authors need to make additional assumptions on the graph ([Assumptions 2 and 3, (Ram, Nedić, and Veeravalli 2010b)]). In this paper, we remove these assumptions, and we present the computational errors of D-SGD with projections in the strongly convex setting.

3. 3.

   Generalization bounds for D-SGD with convexity: We derive (excess) generalization bounds for convex D-SGD. The excess generalization is controlled by the computational error and the generalization bound, which can be directly obtained from the stability.

4. 4.

   Numerical results: We numerically verify our theoretical results by using various benchmark machine learning models, ranging from strongly convex and convex to nonconvex settings, in different decentralized settings.

The population risk minimization is an important model in machine learning and statistics, whose mathematical formulation reads as

where $f({\bf x};\xi)$ denotes the loss of the model associated with data $\xi$ and $\mathcal{D}$ is the data distribution. Due to the fact that $\mathcal{D}$ is usually unknown or very complicated, we consider the following surrogate ERM

where $S:=\{\xi_{1},\xi_{2},\ldots,\xi_{N}\}$ and $~{}\xi_{i}\sim\mathcal{D}$ is a given data.

For a specific stochastic algorithm $\mathcal{A}$ act on $S$ with output $\mathcal{A}(S)$, the generalization error of $\mathcal{A}$ is defined as $\epsilon_{\textrm{gen}}:=\mathbb{E}_{S,\mathcal{A}}[R(\mathcal{A}(S))-R_{S}(\mathcal{A}(S))].$ Here, the expectation is taken over the algorithm and the data. The generalization bound reflects the joint effects caused by the data $S$ and the algorithm $\mathcal{A}$. We are also interested in the excess generalization error, which is defined as $\epsilon_{\textrm{ex-gen}}:=\mathbb{E}_{S,\mathcal{A}}[R(\mathcal{A}(S))-R({\bf x}^{*})]$, where ${\bf x}^{*}$ is the minimizer of $R$. Let $\overline{{\bf x}}$ be the minimizer of $R_{S}$. Due to the unbiased expectation of the data $S$, we have $\mathbb{E}_{S}[R_{S}({\bf x}^{*})]=\mathbb{E}[R({\bf x}^{*})]$. Thus, Bottou and Bousquet (2008) point out $\epsilon_{\textrm{ex-gen}}$ can be decomposed as follows

Notice that $R_{S}(\overline{{\bf x}})\leq R_{S}({\bf x}^{*})$, therefore

The uniform stability is used to bound the generalization error of a given algorithm $\mathcal{A}$ (Hardt, Recht, and Singer 2016; Elisseeff, Evgeniou, and Pontil 2005).

It has been proved that the uniform stability directly implies the generalization bound.

Thus, to get the generalization bound of a random algorithm, we just need to compute the uniform stability bound $\epsilon$.

Notation: We use the following notations throughout the paper. We denote the $\ell_{2}$ norm of ${\bf x}\in\mathbb{R}^{d}$ as $\|{\bf x}\|$. For a matrix ${\bf A}$, ${\bf A}^{\top}$ denotes its transpose, we denote the spectral norm of ${\bf A}$ as $\|{\bf A}\|_{\textrm{op}}$. Given another matrix ${\bf B}$, ${\bf A}\succ{\bf B}$ means that ${\bf A}-{\bf B}$ is positive define; and ${\bf A}\succeq{\bf B}$ means ${\bf A}-{\bf B}$ is positive semidefinite. The identity matrix is defined as $\mathbb{I}$. We use $\mathbb{E}[\cdot]$ to denote the expectation of $\cdot$ with respect to the underlying probability space. For two positive constants $a$ and $b$, we denote $a=\mathcal{O}(b)$ if there exists $C>0$ such that $a\leq Cb$, and $\widetilde{\mathcal{O}}(b)$ hides a logarithmic factor of $b$.

Note that ${\bf W}$ is a doubly stochastic matrix (Marshall, Olkin, and Arnold 1979), and the mixing matrix is non-unique for a given graph. Several common examples for ${\bf W}$ include the Laplacian matrix and the maximum-degree matrix (Boyd, Diaconis, and Xiao 2004). A crucial constant that characterizes the mixing matrix is

where $\lambda_{i}$ denotes the $i$th largest eigenvalue of ${\bf W}\in\mathbb{R}^{m\times m}$. The definition of the mixing matrix implies that $0\leq\lambda<1$.

Note the fact that the stationary distribution of an irreducible aperiodic finite Markov chain is uniform if and only if its transition matrix is doubly stochastic. Thus, ${\bf W}$ corresponds to some Markov chain’s transition matrix, and the parameter $0\leq\lambda<1$ characterizes the speed of convergence to the stationary state.

We consider a general decentralized stochastic gradient descent with projection, which carries out in the following manner: in the $t$-th iteration, 1) client $i$ applies an approximate copy ${\bf x}^{t}(i)\in\mathbb{R}^{d}$ to calculate a unbiased gradient estimate $\nabla f({\bf x}^{t}(i);\xi_{j_{t}(i)})$, where $j_{t}(i)\in\mathbb{Z}^{+}$ is the local random index; 2) client $i$ replaces its local parameters with the weighted average of its neighbors, i.e.,

3) client $i$ updates its parameters as

with learning rate $\alpha_{t}>0$, and ${\bf Proj}_{V}(\cdot)$ stands for projecting the quantity $\cdot$ into the space $V$. We stress that, in practice, we do not need to compute the average ${\bf x}^{t}=\frac{1}{m}\sum_{i=1}^{m}{\bf x}^{t}(i)$ in each iteration, and we take the average only in the last iteration.

In the following, we draw necessary assumptions, which are all common and widely used in the nonconvex analysis community.

Assumption 1 implies that $|f({\bf x};\xi)-f({\bf y};\xi)|\leq B\|{\bf x}-{\bf y}\|,$ for all ${\bf x},{\bf y}\in V$ and any $\xi\sim\mathcal{D}$.

Compared with the scheme presented in (Lian et al. 2017), our algorithm accommodates a projection after each update in each client. When $\nabla f({\bf x};\xi)$ is non-strongly convex, $V$ can be set as the full space and Algorithm 1 reduces to the scheme given in (Lian et al. 2017), whose convergence has been well studied. Such a projection is more general and is necessary for the strongly convex analysis; we explain this necessary claim as follows: if $f$ is $\nu$-strongly convex, then $\|\nabla f({\bf x})\|^{2}\geq\nu\|{\bf x}-{\bf x}^{*}\|^{2}$ with ${\bf x}^{*}$ being the minimizer of $f({\bf x})$ (Karimi, Nutini, and Schmidt 2016). Thus, when ${\bf x}$ is far from ${\bf x}^{*}$, the gradient is unbounded, which breaks Assumption 1. However, with the projection procedure, D-SGD (Algorithm 1) actually minimizes function (1) over the set $V$. The strong convexity gives us $f({\bf x})-f({\bf x}^{*})\geq\frac{\nu}{2}\|{\bf x}-{\bf x}^{*}\|^{2}$, which indicates ${\bf x}^{*}\in\textbf{B}(\textbf{0},\sqrt{{2(f({\bf x}^{0})-\min f)}/{\nu}})$. Thus, when the radius of $V$ is large enough, the projection does not change the output of D-SGD.

This part contains the stability result of D-SGD when $f(\cdot;\xi)$ is generally convex.

Compared to the results of minimizing (1) by using centralized SGD with step sizes $(\alpha_{t})_{t\geq 1}$ [Theorem 3.8, (Hardt, Recht, and Singer 2016)], which yields the uniformly stable bound as ${2B^{2}\sum_{t=1}^{T-1}\alpha_{t}}/{(mn)}$. Theorem 1 shows that D-SGD suffers from an additional term $4B^{2}\sum_{t=1}^{T-1}(1+\alpha_{t}B)\sum_{j=0}^{t-1}\alpha_{j}\lambda^{t-1-j}$, which does not vanish when $\lambda>0$.

If we set $\alpha_{t}={1}/{(t+1)}$, with Lemma 3, it is easy to check that $\epsilon_{\textrm{stab}}=\mathcal{O}(\frac{\ln T}{mn}+C_{\lambda}\ln T)$; However, if we use a constant learning rate, (i.e., $\alpha_{t}\equiv\alpha$), when $0<\lambda<1$, we have $4B^{2}\sum_{t=1}^{T-1}(1+\alpha_{t}B)\sum_{j=0}^{t-1}\alpha_{j}\lambda^{t-1-j}=\mathcal{O}(\frac{\alpha T}{1-\lambda})$ and $\epsilon_{\textrm{stab}}=\mathcal{O}(\frac{\alpha T}{1-\lambda}+\frac{\alpha T}{mn})$. The result indicates that although decentralization reduces the busiest node’s communication, it hurts the stability.

Theorem 1 provides the uniform stability for the last-iterate of D-SGD. However, the computational error of D-SGD in general convexity case uses the following average

Such a mismatch leads to the difficulty in characterizing the excess generalization bound. It is thus necessary describe to the uniform stability in terms of $\textrm{ave}({\bf x}^{T})$. To this end, we consider that D-SGD outputs $\textrm{ave}({\bf x}^{t})$ instead of ${\bf x}^{t}$ in the $t$-th iteration. The uniform stability, in this case, is defined as $\epsilon_{\textrm{ave-stab}}$, and we have the following result.

Unlike the uniform stability for ${\bf x}^{T}$, the average turns out to be a very complicated one. We thus just present two classical kinds of step size.

In the convex setting, for a fixed iteration number $T$, as the data size $mn$ increases and $\lambda$ decreases, $\epsilon_{\textrm{stab}}$ gets smaller for both diminishing and constant learning rates. However, similar to SGD, D-SGD also fails to have $\epsilon_{\textrm{stab}}$ under control when $T$ increases. This drawback does not exist in the strongly convex setting.

Strongly convex loss functions appear in the $\ell_{2}$ regularized machine learning models. As mentioned in Section 2, to guarantee the bounded gradient, the set $V$ should be restricted to a closed ball. We formulate the uniform stability results in this case in Theorem 2.

The uniformly stability bound for SGD with strong convexity is ${2B^{2}}/{(mn\nu)}$ (Hardt, Recht, and Singer 2016), which is smaller than the one of D-SGD. From Theorem 2, we see that the uniform stability bound of D-SGD is independent on the iterative number $T$. Moreover, D-SGD enjoys a smaller uniformly stable bound when the data size $mn$ is larger and $\lambda$ is smaller.

We now present the stability result for nonconvex loss functions.

Without the convexity assumption, the uniform stable bound of D-SGD deteriorated. Theorem 3 shows that $\epsilon_{\textrm{stab}}=\mathcal{O}((1+C_{\lambda}){T^{\frac{cL}{1+cL}}}/{(mn)})$, which is much larger than the bounds in the convex case $\left(\mathcal{O}({\ln T}/{(mn)}+C_{\lambda}\ln T)\right)$.

This part consists of optimization errors of D-SGD for convex and strongly convex settings. Assume ${\bf x}^{*}$ is the minimizer of $f({\bf x})$ over the set $V$, i.e., $f({\bf x}^{*})=\min_{{\bf x}\in V}f({\bf x})$.

It is worth mentioning that the optimization error is established on the average point $\textrm{ave}({\bf x}^{T})$ for technical reasons.

In the following, we provide the results for the strongly convex setting.

The result shows that D-SGD with projection converges sublinearly in the strongly convex case. To reach an $\epsilon>0$ error, we shall set the iteration number as $\widetilde{\mathcal{O}}({1}/{\epsilon})$. Our result coincides with the existing convergence results of SGD with strong convexity (Rakhlin, Shamir, and Sridharan 2012). What is different is that D-SGD is affected by the parameter $\lambda$, which is determined by the structure of the connected graph²²2For the strongly convex case, we avoid showing the result under general step size due to the complicated form..

Notice that the computational error of D-SGD, in this case, is described by $\textrm{ave}({\bf x}^{T})$. Thus, we need to estimate the generalization bound about $\textrm{ave}({\bf x}^{T})$.

Now, we present the excess generalization of D-SGD under strong convexity.

We consider the following optimization problem

which arises from a simple regression problem. Here, we use the Body Fat dataset (Johnson 1996) which contains 252 samples. We run D-SGD on two subsets of the Body Fat dataset, and both of size 200. Let ${\bf x}^{k}$ and $\widehat{{\bf x}}^{k}$ be the outputs of the D-SGD on the two different subsets. We define the absolute loss difference as $|\Phi({\bf x}^{k})-\Phi(\widehat{{\bf x}}^{k})|.$ For the above six graphs, we record the absolute difference in the value of function $\Phi$ for a set of learning rate, namely, $\{0.001,0.004,0.016,0.064\}$ in Figure 3. In the second test, we use the learning rate $0.001$ and compare the absolute loss difference with different graphs in Figure 2 (a). Our results show that the smaller learning rate usually yields a smaller loss difference, and the complete graph can achieve the smallest bound. These observations are consistent with our theoretical results for the convex D-SGD.

To verify our theory on the strongly convex case, we consider the regularized logistic regression model as follows

We use the benchmark ijcnn1 dataset(Rennie and Rifkin 2001) and set $\lambda=10^{-4}$. Two 8000-sample sub-datasets with 1000 different samples are used as the test set. We conduct experiments on the two datasets with the same set of learning rates that are used in the last subsection. The absolute loss difference under different learning rates is plotted in Figure 4, and the performance under different graphs is reported in Figure 2 (b). The results of D-SGD in the strongly convex case is similar to the convex case. Also, note that the absolute loss difference increases as the learning rates grow.

We test ResNet-20 (He et al. 2016) for CIFAR10 classification (Krizhevsky 2009). We adopt two different 40000-sample subsets. The loss values are built on the test set. The absolute loss difference with the learning rate set $\{0.0001,0.0004,0.0016,0.0064\}$ versus the epochs is presented in Figure 5, and the absolute loss difference with different graphs are shown in Figure 5 (c). 100 epochs are used in the nonconvex test. The results show that the nonconvex D-SGD is much more unstable than the convex ones, which matches our theoretical findings.