Decentralized and Parallel Primal and Dual Accelerated Methods for Stochastic Convex Programming Problems

TL;DR

This paper presents decentralized primal and dual stochastic gradient methods that are optimal in communication steps and can be parallelized at each node, improving efficiency for large-scale convex optimization problems.

Contribution

It introduces new primal and dual stochastic gradient methods for decentralized convex optimization that are optimal in communication and nearly optimal in oracle calls, with parallelization capabilities.

Findings

Methods are optimal in communication steps.

Algorithms can be parallelized at each node.

Applicable to various data science and inverse problems.

Abstract

We introduce primal and dual stochastic gradient oracle methods for decentralized convex optimization problems. Both for primal and dual oracles, the proposed methods are optimal in terms of the number of communication steps. However, for all classes of the objective, the optimality in terms of the number of oracle calls per node takes place only up to a logarithmic factor and the notion of smoothness. By using mini-batching technique, we show that the proposed methods with stochastic oracle can be additionally parallelized at each node. The considered algorithms can be applied to many data science problems and inverse problems.

Tables

Table 1. Table 1: The optimal number of first-order oracle calls (number of iterations N 𝑁 N )

	$\mu$-strongly convex and $L$-smooth	$L$-smooth	$\mu$-strongly convex
#iterations	$\sqrt{\frac{L}{\mu }}\ln \left(\frac{\mu R^2}{\epsilon }\right)$	$\sqrt{\frac{L{R}^2}{\epsilon }}$	$\frac{M^2}{\mu \epsilon }$	$\frac{M^2{R}^2}{{\epsilon }^2}$
#oracle calls of $\nabla f\left(x\right)$	$\sqrt{\frac{L}{\mu }}\ln \left(\frac{\mu R^2}{\epsilon }\right)$	$\sqrt{\frac{L{R}^2}{\epsilon }}$	$\frac{M^2}{\mu \epsilon }$	$\frac{M^2{R}^2}{{\epsilon }^2}$

Table 2. Table 2: The optimal number of stochastic (unbiased) first-order oracle calls

	$\mu$-strongly convex and $L$-smooth	$L$-smooth	$\mu$-strongly convex
#iterations	$\sqrt{\frac{L}{\mu }}\ln \left(\frac{\mu R^2}{\epsilon }\right)$	$\sqrt{\frac{L{R}^2}{\epsilon }}$	$\frac{M^2}{\mu \epsilon }$	$\frac{M^2{R}^2}{{\epsilon }^2}$
#oracle calls of $\nabla f(x,\xi )$	$\max \{\frac{{\sigma }^2}{\mu \epsilon },$ $\sqrt{\frac{L}{\mu }}\ln \left(\frac{\mu R^2}{\epsilon }\right)\}$	$\max \{\frac{{\sigma }^2{R}^2}{{\epsilon }^2},\sqrt{\frac{L{R}^2}{\epsilon }}\}$	$\frac{M^2}{\mu \epsilon }$	$\frac{M^2{R}^2}{{\epsilon }^2}$

Table 3. Table 3: The optimal bounds for primal deterministic oracle

	$f_k$ is $\mu$-strongly convex and $L$-smooth	$f_k$ is $L$-smooth	$f_k$ is $\mu$-strongly convex
#communication rounds	$\tilde{O}\left(\sqrt{\frac{L}{\mu }\chi }\right)$	$\tilde{O}\left(\sqrt{\frac{L{R}^2}{\epsilon }\chi }\right)$	$O\left(\sqrt{\frac{M^2}{\mu \epsilon }\chi }\right)$	$O\left(\sqrt{\frac{M^2{R}^2}{{\epsilon }^2}\chi }\right)$
#oracle calls of $\nabla f_k\left(x_k\right)$ per node $k$	$\tilde{O}\left(\sqrt{\frac{L}{\mu }}\right)$	$O\left(\sqrt{\frac{L{R}^2}{\epsilon }}\right)$	$O\left(\frac{M^2}{\mu \epsilon }\right)$	$O\left(\frac{M^2{R}^2}{{\epsilon }^2}\right)$
Algorithm	PSTM, $Q={ℝ}^n$	PSTM, $Q={ℝ}^n$	R-Sliding	Sliding

Table 4. Table 4: The optimal bounds for primal stochastic (unbiased) oracle

	$f_k$ is $\mu$-strongly convex and $L$-smooth	$f_k$ is $L$-smooth	$f_k$ is $\mu$-strongly convex
#communication rounds	$\tilde{O}\left(\sqrt{\frac{L}{\mu }\chi }\right)$	$\tilde{O}\left(\sqrt{\frac{L{R}^2}{\epsilon }\chi }\right)$	$O\left(\sqrt{\frac{M^2}{\mu \epsilon }\chi }\right)$	$O\left(\sqrt{\frac{M^2{R}^2}{{\epsilon }^2}\chi }\right)$
#oracle calls of $\nabla f_k(x_k,{\xi }_k)$ per node $k$	$\tilde{O}\left(\max \{\frac{{\sigma }^2}{\mu \epsilon },\sqrt{\frac{L}{\mu }}\}\right)$	$O\left(\max \{\frac{{\sigma }^2{R}^2}{{\epsilon }^2},\sqrt{\frac{L{R}^2}{\epsilon }}\}\right)$	$O\left(\frac{M^2+{\sigma }^2}{\mu \epsilon }\right)$	$O\left(\frac{\left(M^2+{\sigma }^2\right)R^2}{{\epsilon }^2}\right)$
Algorithm	PBSTM, $Q={ℝ}^n$	PBSTM, $Q={ℝ}^n$	Stochastic R-Sliding	Stochastic Sliding

Table 5. Table 5: The optimal bounds for dual deterministic oracle

	$f_k$ is $\mu$-strongly convex and $L$-smooth	$f_k$ is $\mu$-strongly convex
#communication rounds	$\tilde{O}\left(\sqrt{\frac{L}{\mu }\chi }\right)$	$O\left(\sqrt{\frac{M^2}{\mu \epsilon }\chi }\right)$
#oracle calls of $\nabla {\phi }_k\left({\lambda }_k\right)$ per node $k$	$\tilde{O}\left(\sqrt{\frac{L}{\mu }\chi }\right)$	$O\left(\sqrt{\frac{M^2}{\mu \epsilon }\chi }\right)$
Algorithm	ROGM-G or STM, $Q={ℝ}^n$	OGM-G or PDSTM

Table 6. Table 6: The optimal bounds for dual stochastic (unbiased) oracle

	$f_k$ is $\mu$-strongly convex and $L$-smooth	$f_k$ is $\mu$-strongly convex
#communication rounds	$\tilde{O}\left(\sqrt{\frac{L}{\mu }\chi }\right)$	$O\left(\sqrt{\frac{M^2}{\mu \epsilon }\chi }\right)$
#oracle calls of $\nabla {\phi }_k({\lambda }_k,{\xi }_k)$ per node $k$	$\tilde{O}\left(\max \{\frac{M^2{\sigma }_{\phi }^2}{{\epsilon }^2}\chi ,\sqrt{\frac{L}{\mu }\chi }\}\right)$	$O\left(\max \{\frac{M^2{\sigma }_{\phi }^2}{{\epsilon }^2}\chi ,\sqrt{\frac{M^2}{\mu \epsilon }\chi }\}\right)$
Algorithm	R-RRMA+AC-SA2, $Q={ℝ}^n$	SPDSTM

Equations

$$\underset{x\in Q\subseteq {\mathbb{R}}^n}{\min }f(x):=\mathbb{E}[f(x,\xi )].$$

$$\min \left\{O\left(\frac{M^2{R}^2}{{\epsilon }^2}\right),O\left(\frac{M^2}{\mu \epsilon }\right)\right\}$$

$$\min \left\{O\left(\sqrt{\frac{L{R}^2}{\epsilon }}\right)+O\left(\frac{{\sigma }^2{R}^2}{{\epsilon }^2}\right),O\left(\sqrt{\frac{L}{\mu }}\ln \left(\frac{\mu R^2}{\epsilon }\right)\right)+O\left(\frac{{\sigma }^2}{\mu \epsilon }\right)\right\}$$

$$O\left(\frac{{\sigma }^2{R}^2/{\epsilon }^2}{\sqrt{L{R}^2/\epsilon }}\right)\text{ or }O\left(\frac{{\sigma }^2/(\mu \epsilon )}{\sqrt{L/\mu }\ln \left(\mu R^2/\epsilon \right)}\right)$$

$$\underset{x\in Q\subseteq {\mathbb{R}}^n}{\min }\tilde{f}(x):=\frac{1}{m}\sum _{k=1}^{m}f(x,{\xi }^k)$$

$$m=\min \left\{\tilde{O}\left(\frac{n{M}^2{R}^2}{{\epsilon }^2}\right),\tilde{O}\left(\frac{M^2}{\mu \epsilon }\right)\right\}.$$

$$\underset{x\in Q\subseteq {\mathbb{R}}^n}{\min }\tilde{f}(x):=\frac{1}{q}\sum _{k=1}^q\mathbb{E}[f_k(x,{\eta }^k)],$$

$$\underset{x\in Q}{\min }f(x):=\frac{1}{m}\sum _{k=1}^m{f}_k(x).$$

$$\underset{x\in Q\subseteq {\mathbb{R}}^n}{\min }f(x),$$

$$\forall x\in Q\Vert \nabla f(x){\Vert }_2\le M.$$

$$\forall x,y\in Q\Vert \nabla f(y)-\nabla f(x){\Vert }_2\le L\Vert y-x{\Vert }_2.$$

$$\forall x,y\in Qf(y)\ge f(x)+⟨\nabla f(x),y-x⟩+\frac{\mu }{2}\Vert y-x{\Vert }_2^2.$$

$$z^{k+1}=\arg \underset{z\in Q}{\min }\{\sum _{l=0}^{k+1}{\alpha }_l\left(⟨\nabla f({\tilde{x}}^l),z-{\tilde{x}}^l⟩+h(z)+\frac{\mu }{2}\Vert z-{\tilde{x}}^l{\Vert }_2^2\right)$$

$$+\frac{1}{2}\Vert z-{\tilde{x}}^0{\Vert }_2^2\}.$$

$$O\left(\frac{({\alpha }_{k+1}\epsilon {)}^2(A_{k+1}\mu +1)}{(A_{k+1}{L}_{h}R{)}^2}\right)\ge O\left(\frac{{\epsilon }^3}{L{L}_h^2{R}^4}\right),$$

$${\Vert \mathbb{E}[\nabla f(x,\xi )]-\nabla f(x)\Vert }_2\le \delta =O(\epsilon /R)$$

$$\mathbb{E}\left[\exp \left(\Vert \nabla f(x,\xi )-\mathbb{E}[\nabla f(x,\xi )]{\Vert }_2^2{\sigma }^{-2}\right)\right]\le \exp (1),$$

$$N=\min \left\{O\left(\sqrt{\frac{L{R}^2}{\epsilon }}\right),O\left(\sqrt{\frac{L}{\mu }}\ln \left(\frac{L{R}^2}{\epsilon }\right)\right)\right\}$$

$${\nabla }^{r_{k+1}}f({\tilde{x}}^{k+1},\{{\xi }_i^{k+1}{\}}_{i=1}^{r_{k+1}})=\frac{1}{r_{k+1}}\sum _{i=1}^{r_{k+1}}\nabla f({\tilde{x}}^{k+1},{\xi }_i^{k+1}),$$

$$r_{k+1}=O\left(\frac{{\sigma }^2{\alpha }_{k+1}\ln (N/\beta )}{(1+A_{k+1}\mu )\epsilon }\right).$$

$$\sum _{k=0}^N{r}_k=O(N)+\min$$

$$O\left(\frac{{\sigma }^2}{\mu \epsilon }\ln \left(\frac{L{R}^2}{\epsilon }\right)\ln \left(\frac{\sqrt{L/\mu }}{\beta }\right)\right)\}.$$

$$f(x)+⟨\tilde{\nabla }f(x),y-x⟩-{\delta }_1\Vert y-x{\Vert }_2\le f(y)$$

$$\le f(x)+⟨\tilde{\nabla }f(x),y-x⟩+\frac{L}{2}\Vert y-x{\Vert }_2^2+{\delta }_2,$$

$$f(x^N)-f(x^{\ast })=O\left(\frac{L{R}^2}{N^2}+{\delta }_1\tilde{R}+N{\delta }_2\right),$$

$$⟨\tilde{\nabla }f(x)-\nabla f(x),y-x⟩\le \frac{1}{2L}\Vert \tilde{\nabla }f(x)-\nabla f(x){\Vert }_2^2+\frac{L}{2}\Vert y-x{\Vert }_2^2,$$

$$\underset{\begin{array}{c}Ax=0,\\ x\in Q\end{array}}{\min }f(x),$$

$$R_y^2\le \frac{\Vert \nabla f(x^{\ast }){\Vert }_2^2}{{\lambda }_{\min }^{+}(A^{T}A)}.$$

$$F(x)=f(x)+\frac{R_y^2}{\epsilon }\Vert Ax{\Vert }_2^2\to \underset{x\in Q}{\min }.$$

$$F(x^N)-\underset{x\in Q}{\min }F(x)\le \epsilon ,$$

$$f(x^N)-\underset{x\in Q,Ax=0}{\min }f(x)\le \epsilon ,\Vert A{x}^N|{|}_2\le \frac{(1+\sqrt{5})\epsilon }{2R_y}.$$

$$\tilde{O}\left(\sqrt{{\lambda }_{\max }(A^{T}A)/{\lambda }_{\min }^{+}(A^{T}A)}\right)$$

Full text

1 Introduction

We consider the stochastic convex optimization problem

[TABLE]

Such kind of problems arise in many applications of data science [Shalev-Shwartz and Ben-David, 2014, Shapiro et al., 2009] and mathematical statistics [Spokoiny et al., 2012]. To solve this problem with the average precision $\varepsilon$ in the function value (i.e., to find such $x^{N}$ that $\mathbb{E}f(x^{N})-\min\limits_{x\in Q}f(x)\leq\varepsilon$), one can use stochastic gradient (mirror) descent [Juditsky and Nemirovski, 2012] with

[TABLE]

number of calculations of unbiased stochastic subgradients $\nabla f(x,\xi)$. Here $\mu$ is the constant of strong convexity of $f$, $R=\|x^{0}-x_{*}\|_{2}$ is the Euclidean distance between starting point $x^{0}$ and the solution $x^{*}$ of (1) that corresponds to the minimum of this norm. We also used $\mathbb{E}[\|\nabla f(x,\xi)\|_{2}^{2}]\leq M^{2}$. Generally, we can parallelize (2) on no more than $\tilde{O}(1)$ processors by using batch-parallelization [Dvurechensky et al., 2018]. If we additionally assume that $f$ has $L$-Lipschitz (continuous) gradient and $\mathbb{E}[\|\nabla f(x,\xi)-\nabla f(x)\|_{2}^{2}]\leq\sigma^{2}$, then (2) is replaced by

[TABLE]

number of calculations of unbiased stochastic subgradients $\nabla f(x,\xi)$ by using batch parallelization [Devolder, 2013, Dvurechensky and Gasnikov, 2016, Gasnikov and Nesterov, 2018, Ghadimi and Lan, 2013]. In this case we can parallelize subgradients calculations on no more than

[TABLE]

processors (depending on where the minimum in (3) is reached). Notice that this is much better than in previous case. Since this result cannot be improved [Woodworth et al., 2018], it is the best possible way (in general) to solve (1) by using parallel architecture in online context [Shalev-Shwartz et al., 2009].

For many reasons, in some situations in practice, it can be impossible to organize model-based request111For desired $x$ and independently generated $\xi$, a request returns $\nabla f(x,\xi)$. This allows not to keep the set of functions $\{f(x,\xi^{k})\}_{k}$ for different $k$ in the memory. for calculating stochastic gradient $\nabla f(x^{k},\xi^{k})$ in online regime. Typically, in machine learning applications [Hastie et al., 2001, Shalev-Shwartz and Ben-David, 2014], instead of online access to $\left\{\nabla f(x^{k},\xi^{k})\right\}_{k=1}^{m}$ we have offline access. This means that the set of functions $\left\{f(x,\xi^{k})\right\}_{k=1}^{m}$ are stored in the memory and to use them in algorithms, we need to request corresponding function and then calculate its gradient. This may significantly change the complexity of the problem. Indeed, it is known from [Guigues et al., 2017, Shalev-Shwartz et al., 2009, Shapiro et al., 2009] that with high probability the exact solution of problem

[TABLE]

is an $\varepsilon$-solution (in the function value) of problem (1) if

[TABLE]

If $\mu=0$ or it is small enough one may use a regularization technique (see e.g., [Gasnikov, 2017, Shalev-Shwartz et al., 2009]). This allows to reduce the first part of the estimate from $\tilde{O}\left(nM^{2}R^{2}/\varepsilon^{2}\right)$ to $\tilde{O}\left(M^{2}R^{2}/\varepsilon^{2}\right).$ Moreover, we cannot typically find the exact solution of (4) but in the $\mu$-strongly convex (or regularized) case it suffices to solve (4) with accuracy $O(\mu\varepsilon^{2}/M^{2})$ (see [Shalev-Shwartz et al., 2009]).

To solve (4) in offline context we have to store $\{f(x,\xi^{k})\}_{k=1}^{m}$ in the memory. Since $m$ can be large, centralized distributed architecture is often more preferable in this context [Bertsekas and Tsitsiklis, 1989]. In the general case, centralized architecture is based on communication network and it can be obtained by building a spanning tree of a given network [Scaman et al., 2017]. For $m$-node centralized distributed architecture, the number of gradient oracle calls per each node is defined by (3) with $\sigma^{2}=0$, $L$ and $\mu$ corresponding to $\tilde{f}$. The number of communication steps will be $d$ times more, where $d$ is the distance between the origin (root) and farthest node. If we have only $q\ll m$ nodes, then we divide the data $\{f(x,\xi^{k})\}_{k=1}^{m}$ into $q$ blocks with $l=m/q$ terms in each block. If $l$ is too large by itself one can reformulate (4) as follows [McMahan et al., 2016]

[TABLE]

where $f_{k}(x,\eta^{k})=f(x,\xi^{kl+\eta^{k}})$ and $\eta^{k}=i$ with probability $1/l$, $i=1,...,l$. Representation (5) allows to use bound (3) in the stochastic case in a parallel manner at each node. The number of oracle calls per node also corresponds (in general) to (3) and the number of communication steps is also $d$ times more than (3) with $\sigma^{2}=0$.

Unfortunately, centralized architecture has a synchronization drawback and a high requirement for the master node [Scaman et al., 2017]. To address these disadvantages to some extent, a decentralized distributed architecture should be used [Bertsekas and Tsitsiklis, 1989, Kibardin, 1979]. This architecture relies on two basic principles [Nedic, 2020]: every node communicates only with its neighbors, and all communications are performed simultaneously. The main difference here is a simple strategy of communications: each node communicates only with all available direct neighbors. This architecture is more robust. In particular, it can be applied to time-varying (wireless) communication networks [Rogozin et al., 2018].

As problems (4) and (5) have a definite structure of the sum type, they can be solved much faster on one machine. For instance, using some incremental algorithms [Allen-Zhu, 2017, Lan and Zhou, 2017, Lin et al., 2015, Woodworth and Srebro, 2016], one can solve (4) $\sqrt{m}$ times cheaper in terms of the number of oracle calls, but not in terms of the number of iterations ( =communication steps). Unfortunately, this result prohibits parallelization. Note, that for this problem in asynchronized mode ( at each step only two randomly chosen nodes can communicate ), one can obtain such ($\sim\sqrt{m}$) an acceleration for the star-type communication network [Lan and Zhou, 2018]. Moreover, a ‘dual’ analogue of this acceleration has recently been proposed for (4) [Hendrikx et al., 2018] and (5) [Hendrikx et al., 2019a, Hendrikx et al., 2019b] with arbitrary communication networks.

Before stating the contribution we introduce the notions of condition number ${\chi}$ for Laplacian communication matrix of some network and the height of spanning tree denoted by $d$. Note, that $\sqrt{\chi}\geq d$ and typically $\sqrt{\chi}\leq nd$ (see [Nedić et al., 2018]). The last bound corresponds to a star topology [Gasnikov, 2017] (the most simple centralized type architecture). In many interesting cases $\sqrt{\chi}=\tilde{O}(d)$ (see [Nedić et al., 2018, Scaman et al., 2017]).

1.1 Contribution

• •

We justify the transition from the optimal centralized distributed complexity bounds for problems (4) and (5) in the smooth case to decentralized ones222 In the deterministic case this was partially done in [Li et al., 2018]. by replacing $d$ with $\sqrt{\chi}$, the average $L$ with the worse one and variance of $f$ with the variance of $f_{k}$, that can be $m$ times more.333For instance, this takes place in the case when we have independent noise at each $f_{k}$. Note also that in this case in decentralized distributed optimization one can improve the variance dependence and eliminate factor $m$ (see [Olshevsky et al., 2019a, Olshevsky et al., 2019b]). But, this is possible due to the worse estimate for the number of communication steps. Here and everywhere below we keep $\mu$ at the average level (without loss of generality, we can assume that each $f_{k}$ has the same $\mu$) by using a trick from [Scaman et al., 2017]. The announced results are also not improvable in terms of communication steps (rounds) [Arjevani and Shamir, 2015, Scaman et al., 2017].

By using different smoothing techniques [Allen-Zhu and Hazan, 2016, Nesterov, 2005, Scaman et al., 2018], we may lead the non-smooth case to the smooth one with $L\sim 1/\varepsilon$. This allows to reduce the complexity estimate (2) by using (3). However, in general, this reduction makes the cost of oracle calls more expensive. Thus, we can only improve the communication steps (rounds) bound that corresponds to (3) (up to a $\sqrt{\chi}$ factor) with $L\sim 1/\varepsilon$ and $\sigma^{2}=0$.444Note that such tricks sometimes allow to obtain the optimal (in terms of dependence on $\varepsilon$) communication round estimates [Arjevani and Shamir, 2015, Scaman et al., 2018]. Can we preserve the bound (2) for standard conception of oracle calls (primal oracle that gives $\nabla f_{k}$) per node in decentralized approach by improving the number of communication steps? The answer is positive up to the replacement of the average $M$ to the worst one [Lan et al., 2017, Scaman et al., 2018]. In our paper, we simplify the approaches proposed in these articles to prove this result.

As a motivation for the second main result, we consider the problem of type

[TABLE]

where $f_{k}(x)$ has Fenchel–Legendre representation $f_{k}(x)=\max\limits_{y}\left\{\langle x,y\rangle-\varphi_{k}(y)\right\}$ with convex $\varphi_{k}(y)$. Such type of optimization problems arise, for instance, in different applications of inverse problems, especially for linear problems in Hilbert space where we have to use discretization and sum-type functions naturally arise [Byrne, 2014, Gao and Blumensath, 2017, Gasnikov et al., 2017, Vogel, 2002, Ye et al., 2019]. The case of dual-friendly $f_{k}$ also arises in the problem of Wasserstein barycenter calculation [Dvinskikh et al., 2019, Dvurechenskii et al., 2018, Uribe et al., 2018]. Suppose that $\nabla\varphi_{k}(y)$ is available but $f_{k}(x),\nabla f_{k}(x)$ are unavailable. Moreover, sometimes $\varphi_{k}(y)=\mathbb{E}[\varphi_{k}(y,\xi)]$ and it is worth to use $\nabla\varphi_{k}(y,\xi)$ instead of $\nabla\varphi_{k}(y)$ [Dvinskikh et al., 2019, Dvurechenskii et al., 2018]. Considering this example as one of the motivation for using dual oracle instead of primal one, we provide the second main result.

• •

We develop optimal decentralized distributed algorithms with dual (stochastic) oracle for strongly convex objective in (6). The approach is based on dual reformulation of (6) [Scaman et al., 2017]. An optimal algorithm for non-strongly convex dual function with stochastic oracle was recently proposed in [Dvinskikh et al., 2019]. To propose an optimal method with stochastic dual oracle for strongly convex primal objective, we use recent work [Foster et al., 2019]. We notice a rather unexpected result: we cannot improve (up to a logarithmic factor) the bound for the number of dual stochastic gradient calculations in comparison with non-strongly convex dual objective.

We also notice that initially we were motivated by the study of the dual oracle not only as an application from [Dvinskikh et al., 2019, Dvurechenskii et al., 2018, Uribe et al., 2018]. We also tried to find a simple explanation for the optimal communication step bounds [Arjevani and Shamir, 2015, Scaman et al., 2018] in non-smooth case. One of the ways to do it is Nesterov’s dual smoothing technique [Nesterov, 2005] that builds a bridge to the notion of dual oracle. This plan was partially (in the deterministic case) implemented in [Scaman et al., 2017, Uribe et al., 2020, Uribe et al., 2020]. Here we generalize the results of these works for the stochastic dual oracle.

1.2 Paper organization

The paper is organized as follows. In Section 2, we propose optimal stochastic (parallelized) accelerated gradient methods for stochastic convex optimization problems. In Sections 3 and 4, we apply the results of Section 2 to stochastic convex optimization problems with affine type of constraints (of type $Ax=0$). We describe the modern stochastic (parallelized) accelerated gradient methods which are optimal both in terms of (stochastic) oracle calls and matrix-vector multiplications $Ax$. In Sections 3, we are focusing on primal methods, in Section 4, we present dual ones. Section 5, describes the distributed primal and dual formulation of the finite-sum minimization problem, and presents distributed algorithms. In Section 6, we incorporate the proposed distributed decentralized method to get the optimal bounds for the finite-sum minimization problem using primal or dual oracle. Finally, we discuss future work and possible extensions. We notice that all proposed methods are optimal in terms of communication steps and in many cases in terms of (parallel stochastic) primal/dual oracle calls.

2 Stochastic convex optimization

First-order methods for the optimization problem of minimizing a convex function $f$ on a simple convex set $Q$, e.g.,

[TABLE]

play a fundamental role in modern problems arising in machine learning and statistics. The complexity of these methods is measured by the number of iterations or (and) the number of oracle calls. For a deterministic oracle, this concept can be identified. By the first-order oracle, we mean a black-box model that for a given input $x\in Q$, returns the vector $\nabla f(x)$.

We say that a function $f$ is $M$-Lipschitz continuous if 555Here and below in such type of assumptions (especially in the case when $Q$ is unbounded) instead of $\forall x\in Q$ we may write $\forall x\in Q:\|x-x^{*}\|_{2}\leq 2R$ [Gasnikov, 2017] (analogously for $y$).

[TABLE]

We say that function $f$ is $L$-smooth or has $L$-Lipschitz continuous gradient if

[TABLE]

We also say that function $f$ is $\mu$-strongly convex if

[TABLE]

Accelerated gradient methods (e.g., Algorithm 1 (STM) [Gasnikov and Nesterov, 2018, Nesterov, 2018b, Lan, 2019]) allow to obtain the optimal number of iterations and number of gradient oracle calls for problem (7) as described in Table 1, where $R=\|x^{0}-x^{*}\|_{2}$ is the Euclidean distance from the starting point $x^{0}$ to the solution $x^{*}$ of (7) that corresponds to the minimum of this norm, and $\varepsilon$ is the desired precision in function value.

Remark 1

For a composite optimization problem with composite term $h(x)$, step 4 of Algorithm 1 is replaced by the more general operator [Gasnikov and Nesterov, 2018, Nesterov, 2018b]

[TABLE]

If $h(x)$ has $L_{h}$-Lipschitz gradient in the $\ell_{2}$-norm then due to Theorem 9 [Stonyakin et al., 2019b] and Theorem 19 [Stonyakin et al., 2019a] it suffices to solve auxiliary problem with accuracy (in terms of function value)

[TABLE]

*where $\varepsilon$ is desired accuracy (in function value) for initial problem (7).

If $\mu=0$ we can also generalize this step for the non-Euclidean case and using restarts [Gasnikov and Nesterov, 2018] generalize such a method on $\mu>0$. Note, that by using restarts with STM($L$,0,$x^{0}$) one can eliminate the gap from $\ln(LR^{2}/\varepsilon)$ to $\ln(\mu R^{2}/\varepsilon)$ between lower bounds and the bounds for STM($L$,$\mu$,$x^{0}$) without restarts [Gasnikov and Nesterov, 2018]. The same remains true for the stochastic oracle.*

Generally, iteration complexity is determined by the complexity of calculating the gradient, which can be computationally expensive. Thus, stochastic approximations of the true gradient can be used instead. In this case, or when the true gradient is unavailable (if e.g., function $f$ is given in the form of expectation $f(x):=\mathbb{E}[f(x,\xi)$]) we denote the inexact (or noise-corrupted) first-order oracle as $\nabla f(x,\xi)$, given by a blackbox model with stochasticity (noise) $\xi$ corrupting the true gradient. Assume that

[TABLE]

and

[TABLE]

then with probability at least $1-\beta$, we have $f(x^{N})-f(x^{*})\leq\varepsilon$ after

[TABLE]

iterations of STM using the approximated gradient (instead of the real one $\nabla f(\tilde{x}^{k+1})$)

[TABLE]

where $\xi_{1}^{k+1},\dots,\xi_{r_{k+1}}^{k+1}$ are i.i.d from the same distribution as $\xi$ and the batch size is

[TABLE]

Moreover, the total number of oracle calls666 Oracle calls can be easily and fully parallelized (on $r_{k}$ processors) at each iteration. Note, that for $\nabla^{r_{k}}f(x,\{\xi_{i}\}_{i=1}^{r_{k}})$ we can reduce the variance $\sigma^{2}:=O(\sigma^{2}/r_{k})$. is (this bound is optimal up to logarithmic factors)

[TABLE]

We refer to such a variant of STM as BSTM($L$,$\mu$,$\sigma^{2}$,$x^{0}$) (batched STM($L$,$\mu$,$x^{0}$)).

Thus, using minibatches for constructing an approximation of the true gradient allows us to keep the optimal number of iterations for stochastic methods, as presented in Table 1, where we skip high probability logarithmic multipliers. The number of stochastic oracle calls for this case is shown in Table 2.

Remark 2

*We notice that in assumption (8), $R=\|x^{0}-x^{*}\|_{2}$. Generally, in such type of assumptions, $R$ is the diameter of $Q$ (see [Cohen et al., 2018, d’Aspremont, 2008]) (it is not a compact set when $Q=\mathbb{R}^{n}$). To obtain such a generalization we have to use the advanced recurrent technique to bound $\|z^{k}-x^{*}\|_{2}$ from [Dvinskikh et al., 2019, Gorbunov et al., 2018] and [Gasnikov, 2017, Chapter 2]. Further we provide the sketches how to get this result (for simplicity $\sigma=0,\mu=0$).

1. For an inexact gradient $\tilde{\nabla}f(x)$ satisfying for all $x,y$*

[TABLE]

STM* outputs777Note, that according to [Poljak, 1981] even for the last point of gradient descent on a simple quadratic optimization problem we cannot guarantee convergence without proper stopping rule. With proper stopping rule in (8) it is required (see [Polyak, 1987, Theorem 7, item 6.1.3]) $\delta\sim\varepsilon^{2}$ that is worse than what we have $\delta\sim\varepsilon$. But we can guarantee standard convergence of noisy gradient descent under (8) in (Cesaro) average [Gasnikov, 2017] (not for the last point). The results below generalize [Gasnikov, 2017] on a proper accelerated method (STM). $x^{N}$ such that [Devolder et al., 2014, Dvinskikh et al., 2020]*

[TABLE]

*where $\max\left\{\|\tilde{x}^{k}-x^{*}\|_{2},\|z^{k}-x^{*}\|_{2},\|x^{k}-x^{*}\|_{2}\right\}\leq\tilde{R}$.

2. Since*

[TABLE]

*one can consider $\delta_{2}:=\delta^{2}/(2L)$ and $L:=2L$ in (2).

3. In the deterministic case, $\tilde{R}=R$ with the proper stopping rule of the algorithm (see [Gasnikov, 2017, formulas (2.17), (2.18) in Chapter 2] and [Gasnikov and Nesterov, 2018]). In the stochastic case (with high probability) $\tilde{R}=O(R)$ (see [Dvinskikh et al., 2019, Gorbunov et al., 2018]) for STM with the proper batch-size (10).*

In particular, for the case of non-smooth objective, the stochastic oracle does not yield gains compared to its deterministic counterpart.

In Table 2, in the appropriate cells, we assumed that the following inequalities hold: $\mathbb{E}\|\nabla f(x,\xi)-\nabla f(x)\|_{2}^{2}\leq\sigma^{2}$ and $\mathbb{E}\|\nabla f(x,\xi)\|^{2}_{2}\leq M^{2}$.

Both in Table 1 and in Table 2, the last two columns can be obtained from the corresponding first columns by choosing $L=M^{2}/(2\delta)$, where $\delta=\varepsilon/N$ (see [Gasnikov and Nesterov, 2018]). This is the idea of universal accelerated methods [Nesterov, 2015], but with predefined $L$. Here and in all further tables we skip numerical constants.

3 Primal methods for stochastic convex optimization with affine constraints

To build the complete theory of distributed primal and dual method we need to generalize the result of Tables 1 and 2 for the convex optimization problem888In decentralized optimization $A$ is taken to be $\sqrt{W}$ (square root of the Laplacian matrix of the communication network).

[TABLE]

where $A{\succeq}0$ and $\text{Ker}A\neq\varnothing$. The purpose of this section is to develop such algorithms for (12) that are optimal in terms of the number of $\nabla f(x)$ calculations and the number of $A^{T}Ax$ calculations. In this section we use Euclidean proximal setup [Ben-Tal and Nemirovski, 2001]. This is the only section where we significantly rely on Euclidean prox-structure.

Denote by $R_{y}=\|y^{*}\|_{2}$ the $\ell_{2}$-norm of the smallest solution $y^{*}$ of dual (up to a sign) problem (4). Solution $y^{*}$ is not unique since $\text{Ker}A\neq\varnothing$. From [Lan et al., 2017] we have such a bound

[TABLE]

Using the penalty method we rewrite (12) as follows

[TABLE]

Next we use [Gasnikov, 2017, Remark 4.2] and get if the following holds

[TABLE]

then

[TABLE]

We start with the smooth case and assume that $Q=\mathbb{R}^{n}$. If $f$ has $L$-Lipschitz continuous gradient then we can solve (14) by STM (or BSTM in the stochastic case) considering the second term to be composite [Gasnikov and Nesterov, 2018, Nesterov, 2013]. In this case, we obtain the optimal number of $\nabla f(x)$ (or $\nabla f(x,\xi)$) calculations, see Tables 1, 2. But the total number of $A^{T}Ax$ calculations is

[TABLE]

times more since $\text{Im}A=\text{Im}A^{T}=(\text{Ker}A)^{\perp}$ and $Q=\mathbb{R}^{n}$ and as a consequence of these facts the auxiliary problem can be divided into two subproblems: minimization of quadratic form with matrix of the form $(R_{y}^{2}/\varepsilon){A^{T}A}+cI$ ($c$ is some positive constant and $I$ is identity matrix) on $(\text{Ker}A)^{\perp}$ and minimization of quadratic form with matrix of the form $cI$ on $\text{Ker}A$. Linear terms do not play any role in complexity. The complexity of the auxiliary problem is determined by the worst (reduced on corresponding subspace) conditional number of these two subproblems. Obviously, the first one is worse. The reduced the conditional number is

[TABLE]

This factor arises because of the complexity of the auxiliary problem. We refer to these approaches as PSTM and PBSTM (Penalty STM and BSTM). Here and below we skip arguments of the algorithms if they are obvious from the context.

In non-smooth case ($f$ is $M$-Lipschitz), we use the Sliding algorithm [Lan, 2016], [Lan, 2019] . If $\mu=0$ according to [Lan, 2016] this algorithm requires (see Tables 1, 2 for comparison)

$O\left(\sqrt{\frac{\lambda_{\max}(A^{T}A)R_{y}^{2}R_{x}^{2}}{\varepsilon^{2}}}\right)$

calculations of $A^{T}Ax$ and

$O\left(\frac{M^{2}R_{x}^{2}}{\varepsilon^{2}}\right)$

calculations of $\nabla f(x)$, where $R_{x}=\|x^{0}-x^{*}\|_{2}$.

If we have unbiased $\nabla f(x,\xi)$ with $\sigma^{2}$-sub-Gaussian variance [Jin et al., 2019] instead of $\nabla f(x)$, i.e.,

[TABLE]

with $\sigma^{2}=O(M^{2})$ (for compact notation999In general $M^{2}$ is replaced by $M^{2}+\sigma^{2}$ in the stochastic case.), then the bound for calculations of $A^{T}A$ does not change and the bound for calculations of $\nabla f(x,\xi)$ is the same as it was for the number of calculations of $\nabla f(x)$ in deterministic case (up to a logarithmic high-probability deviations factor).

By using a restart technique [Uribe et al., 2020] we can generalize this method for $\mu$-strongly convex $f$:

[TABLE]

calculations of $A^{T}Ax$ and

[TABLE]

calculations of $\nabla f(x)$ ($\nabla f(x,\xi)$).

We call this approach by R-Sliding (Restart Sliding).

4 Dual methods for stochastic convex optimization with affine constraints

Now we assume that we can build a dual problem for

[TABLE]

where $\text{Ker}A\neq\varnothing$.

Remark 3

We notice that turning to a dual problem does not oblige us using dual oracle. Instead, we can use a primal oracle and the Moreau theorem [Rockafellar, 2015] with Fenchel-Legendre representation. This maximization problem can be solved using the first-order oracle for the function $f$. But such an approach does not allow to obtain the optimal bounds on the number of primal first-order oracle calls. Note that typically in decentralized optimization $A$ in (15) is taken as the square root of the Laplacian matrix $W$ of the communication network [Scaman et al., 2017]. But in the asynchronized case the square root $\sqrt{W}$ replaced by incidence matrix $M$ [Hendrikx et al., 2018] ($W=M^{T}M$). Then in asynchronized case instead of accelerated methods for (4) one should use an accelerated (block) coordinate descent method [Dvurechensky et al., 2017, Gasnikov, 2017, Hendrikx et al., 2018, Shalev-Shwartz and Zhang, 2014].

The dual problem (up to a sign) is following

[TABLE]

If $f$ is $\mu$-strongly convex in the $\ell_{2}$-norm, then $\psi$ has $L_{\psi}=\frac{\lambda_{\max}(A^{T}A)}{\mu}$–Lipschitz continuous gradient in the $\ell_{2}$-norm101010Here and below we can also consider other norms (see [Uribe et al., 2020] for details). [Kakade et al., 2009, Rockafellar, 2015]. In this case we can apply STM($L{{{}_{\psi}}}$,0,0) to (4). Note that due to Demyanov–Danskin’s theorem $\nabla\psi(y)=Ax(A^{T}y)$ [Rockafellar, 2015]. Similarly to [Anikin et al., 2017, Chernov et al., 2016] one can prove that

[TABLE]

where $R_{y}=\|y^{*}\|_{2}$ is the radius of solution of (4) which is the smallest in the $\ell_{2}$-norm, see (13). We call this approach by PDSTM (Primal-Dual STM).

If we have only a stochastic (randomized) unbiased model $\nabla\varphi(\lambda,\xi)|_{\lambda=A^{T}y}=x(A^{T}y,\xi)$ with $\sigma^{2}_{\varphi}$-sub-Gaussian variance, i.e.

[TABLE]

then for BSTM($L_{\psi}$,0,$\sigma^{2}_{\psi}$,0) where $\sigma^{2}_{\psi}=\lambda_{\max}(A^{T}A)\sigma^{2}_{\varphi}$ with probability $\geq 1-\beta$ (4) holds true [Dvinskikh et al., 2019]. We refer to this algorithm as SPDSTM (Stochastic PDSTM).

In the case when $\psi$ from (4) is additionally $\mu_{\psi}$-strongly convex in the $\ell_{2}$-norm in111111Since $\text{Im}A=(\text{Ker}A^{T})^{\perp}$ we will have that all the points $\tilde{y}^{k},z^{k},y^{k},$ generated by STM and methods based on STM, belong to $y^{0}+(\text{Ker}A^{T})^{\perp}$. That is, from the point of view of estimates this means, that we can consider $\psi$ to be $\mu_{\psi}$-strongly convex everywhere. $y^{0}+(\text{Ker}A^{T})^{\perp}$ (if $f$ has $L$-Lipschitz gradient in the $\ell_{2}$-norm and $Q=\mathbb{R}^{n}$ then $\mu_{\psi}=\lambda_{\min}^{+}(A^{T}A)/L$ [Kakade et al., 2009, Rockafellar, 2015], where $\lambda_{\min}^{+}(A^{T}A)$ is the minimal positive eigenvalue of $A^{T}A$) we need to use another approach. Because of primal-duality [Nesterov, 2009, Nemirovski et al., 2010] we have to put $\mu_{\psi}=0$ in STM and in methods based on STM (STM($L_{\psi}$,$\mu_{\psi}$,$y^{0}$) is not primal-dual method when $\mu_{\psi}>0$). The restart technique (see, e.g. [Gasnikov, 2017]) also does not work here because in (4) we have to use $R_{y}=\|y^{0}\|_{2}+\|y^{0}-y^{*}\|_{2}$ in general. That is why we take here $y^{0}=0$. So the main trick here is the following relation [Allen-Zhu, 2018, Anikin et al., 2017, Nesterov, 2012]

[TABLE]

From (18), to satisfy

[TABLE]

it is sufficient to find such $y^{N}$ ($\|y^{N}\|_{2}\leq 2R_{y}$) that

[TABLE]

Recently, there appear accelerated methods with the proper rate of convergence in terms of the norm of the gradient OGM-G [Gasnikov, 2017, Kim and Fessler, 2018]:

[TABLE]

After $\bar{N}=O\left(\sqrt{\frac{L_{\psi}}{\mu_{\psi}}}\right)$ iterations of OGM-G we have

[TABLE]

So after $l=\log_{2}\left({\|}\nabla\psi(y^{0})\|_{2}\frac{R_{y}}{\varepsilon}\right)$ restarts ($y^{0}:=y^{\bar{N}}$) we have (19). We denote such an approach by ROGM-G (Restart OGM-G). This approach requires

[TABLE]

of $\nabla\psi(y)$ (that is $Ax(A^{T}y)$) calculations. The key inequality to prove this fact is

[TABLE]

This holds due to

[TABLE]

The same result with the replacement

[TABLE]

can be obtained by using STM($L_{\psi}$,$\mu_{\psi}$,0) with bound ${\sqrt{\frac{L_{\psi}}{\mu_{\psi}}}}\ln\left(\frac{L_{\psi}R_{y}^{2}}{\varepsilon^{\prime}}\right)$ (see [Nesterov, 2010]) and desired accuracy $\varepsilon^{\prime}=\frac{\varepsilon^{2}}{2L_{\psi}R^{2}_{y}}$. This follows from

[TABLE]

Now we consider RRMA+AC-SA2 [Foster et al., 2019] (see also [Allen-Zhu, 2018] in the non-accelerate, but composite case). This algorithm converges as follows (for simplicity we skip polylogarithmic factors and high probability terminology)

[TABLE]

We assume that at each iteration, $\nabla\psi(y,\xi)$ with sub-Gaussian variance $\sigma^{2}_{\psi}$ is available [Jin et al., 2019] (see also above). If we use restarts with the size of each restart $\bar{N}=\tilde{O}\left(\sqrt{\frac{L_{\psi}}{\mu_{\psi}}}\right)$ (see above) and use batched gradient (9) with batch size (at $k$-th restart; $\bar{y}^{k}$ is the output point from the previous restart)

[TABLE]

then $\|\nabla\psi(\bar{y}^{l})\|_{2}\leq\varepsilon/R_{y}$ after $l=O\left(\log_{2}\left(\|\nabla\psi(y^{0})\|_{2}R_{y}/\varepsilon\right)\right)$ restarts. Therefore, the total number of oracle calls is

[TABLE]

Note that the same bound takes place in the non-strongly convex case ($\mu_{\psi}=0$). From [Allen-Zhu, 2018, Jin et al., 2019] it is known that this bound cannot be improved. However, we may expect that this bound can be reduced to $\tilde{O}(\sigma^{2}_{\psi}/(\mu_{\psi}\varepsilon))$ ( see Table 2 for stochastic primal oracle). For the stochastic dual oracle, such a reduction is probably impossible. We call this approach by R-RRMA+AC-SA2 (Restart RRMA+AC-SA2).

5 Decentralized distributed optimization

Now we show how to present (6) in a decentralized distributed manner

[TABLE]

This particular representation of the objective in (P1) allows involving distributed methods which are particularly necessary for large-scale problems handling the large quantities of data and which are based on the idea of agents’ cooperative solution of the global problem [Bertsekas and Tsitsiklis, 1989]. For a given multi-agent network system, we privately assign each function $f_{k}$ to the agent $k$ and suppose that agents can exchange the information with their neighbors (e.g., send and receive vectors). We define this system through the Laplacian matrix $\bar{W}\in\mathbb{R}^{m\times m}$ of some graph (communication network) $G=(V,E)$ with the set $V$ of $m$ vertices and the set of edges $E=\{(i,j):i,j\in V\}$ as follows

[TABLE]

where ${\rm deg}(i)$ is the degree of vertex $i$ (i.e., the number of neighboring nodes).

From the definition of the matrix $\bar{W}$ it can be easily seen that $\bar{W}$ establishes the communication of agents and allows only the communication between neighboring nodes. Moreover, due to connectivity of graph $G$ the vector $\bm{1}_{m}=(1,...,1)^{T}\in\mathbb{R}^{m}$ is the unique (up to a scaling factor) eigenvector of $\bar{W}$ associated with the eigenvalue $\lambda=0$, which allows us to compactly rewrite the consensus agreement $x_{1}=...=x_{m}\in\mathbb{R}^{n}$ as $W\mathbf{x}=0$, moreover, as $\sqrt{W}{\mathbf{x}}=0$ (see [Scaman et al., 2017]), where $W{=\bar{W}}\otimes I_{n}$ is the Kronecker product of the Laplacian matrix $\bar{W}\in\mathbb{R}^{m}$ and the identity (unit) matrix $I_{n}$ and $\mathbf{x}=[x_{1}^{T},...,x_{m}^{T}]^{T}\in\mathbb{R}^{mn}$.

To present problem (P1) in a distributed fashion we rewrite it with introducing the artificial consensus equality constraints and then change these constraints to one affine constraint with the communication matrix $W$ as follows

[TABLE]

or in another form

[TABLE]

where all $f_{k}$ are $M$-Lipschitz, $L$-smooth and $\mu$-strongly convex (it is possible that, $L=\infty$ or (and) $\mu=0$).

We also consider the stochastic version of problem (P2), where $f_{k}(x_{k})=\mathbb{E}[f_{k}(x_{k},\xi_{k})]$. We consider the unbiased stochastic primal oracle that returns $\nabla f_{k}(x_{k},\xi_{k})$ (where $\xi=\{\xi_{k}\}_{k=1}^{m}$ are independent) under the following $\sigma^{2}$-sub-Gaussian variance condition (for all $k=1,...,m$)

[TABLE]

Problem (P2) can be considered to be a particular case of problem (12) with the following replacements

1:•    $A=\sqrt{W}$

2:•    $L_{F}={L}/{m}$

3:•    $\mu_{F}={\mu}/{m}$

4:•    $\|\nabla F(\mathbf{x})\|_{2}^{2}\leq M_{F}^{2}={M^{2}}/{m}$

5:•    $\sigma_{F}^{2}=O\left({\sigma^{2}}/{m}\right)$

6:•    $R_{\mathbf{x}}^{2}=\|\mathbf{x}^{0}-\mathbf{x}^{*}\|_{2}^{2}=m\|x^{0}-x^{*}\|_{2}^{2}=mR^{2}$

7:•    $R_{\mathbf{y}}^{2}=\|\mathbf{y}^{*}\|_{2}^{2}\leq{\|\nabla F(\mathbf{x}^{*})\|_{2}^{2}}/{\lambda_{\min}^{+}(W)}\leq{M^{2}}/{(m\lambda_{\min}^{+}(W))}$

The main observation in the primal approach (see Section 3) is as follows [Scaman et al., 2017]:

$A^{T}Ax=Wx$ (calculated in a decentralized distributed manner)

If each function $f_{k}$ is a dual-friendly (dual function is available calculated by the Fenchel–Legendre transform [Uribe et al., 2020]) then we can construct the dual problem to problem (P2) with dual Lagrangian variables $\mathbf{y}=[y_{1}^{T}\in\mathbb{R}^{n},\cdots,y_{m}^{T}\in\mathbb{R}^{n}]^{T}\in\mathbb{R}^{mn}$

[TABLE]

where $\varphi_{k}(\lambda_{k})=\max\limits_{x_{k}\in Q\subseteq\mathbb{R}^{n}}\{\langle\lambda_{k},x_{k}\rangle-f_{k}(x_{k})\}$ and the vector $[\sqrt{W}\mathbf{x}]_{k}$ represents the $k$-th $n$-dimensional block of $\sqrt{W}\mathbf{x}$.

We also consider the stochastic version of problem (D2), where $\varphi_{k}({\lambda}_{k})=\mathbb{E}[\varphi_{k}({\lambda}_{k},\xi_{k})]$. We consider the unbiased stochastic dual oracle returns $\nabla\varphi_{k}({\lambda}_{k},\xi_{k})$ (where $\xi=\{\xi_{k}\}_{k=1}^{m}$ are independent) under the following $\sigma_{\varphi}^{2}$-sub-Gaussian variance condition (for all $k=1,...,m$)

[TABLE]

Problem (D2) can be considered as a particular case of problem (4) with

1:•    $A=\sqrt{W}$

2:•    $\sigma_{{\Psi}}^{2}=O\left({\lambda_{\max}(W)m\sigma_{\varphi}^{2}}\right)$

The main observation in the dual approach (see Section 4) is as follows [Scaman et al., 2017]: since $x(A^{T}y)=\mathbf{x}(\sqrt{W}\mathbf{y})$ we should change the variables as follows

1:    $\tilde{\mathbf{y}}:=\sqrt{W}\tilde{\mathbf{y}}$

2:    $\mathbf{z}:=\sqrt{W}\mathbf{z}$

3:    $\mathbf{y}:=\sqrt{W}\mathbf{y}$

It is obvious that input, output and steps 3–5 of Algorithm 2 are changed such that they can be performed in a decentralized distributed manner. For that we just multiply the corresponding steps by $\sqrt{W}$.

6 Main Results

In this section, we present the rates of convergence for problems (P1) and (D2) (and their stochastic counterparts) in terms of the number of iterations (communication steps) and the number of (parallelized) oracle calls. For the primal problem, we present the results to achieve $\varepsilon$-precision in objective residuals, and for the dual problem we seek to achieve $\varepsilon$-precision in duality gap or primal objective residuals (in smooth strongly convex case). Feasibility constrains are smaller than $\varepsilon/R_{\mathbf{y}}$.

For brevity, we introduce the condition number of the Laplacian matrix $W$ as follows

[TABLE]

where $\lambda_{\min}^{+}(W)$ is the minimal positive eigenvalue of $W$, and $\lambda_{\max}(W)$ is the maximal eigenvalue of $W$. Now we are ready to present our main results incorporated in multiple tables: Tables 3– 6. These results are obtained by direct substitution of constants from Section 5 to the problems from Sections 3 and 4.

Note that the bounds on communication steps (rounds) are optimal (up to a logarithmic factor) due to [Arjevani and Shamir, 2015, Scaman et al., 2017, Scaman et al., 2018]. Bounds for the oracle calls per node are probably optimal in the class of methods with optimal number of communication steps (up to a logarithmic factor) in the deterministic case [Allen-Zhu, 2018, Foster et al., 2019, Woodworth et al., 2018] and optimal for the non-smooth stochastic primal oracle and stochastic dual oracle for parallel architecture.121212In parallel architecture the bounds on stochastic oracle calls per node of type $\max\{B,D\}$ can be parallel up to $B/D$ processors. For stochastic oracle the bounds hold in terms of high probability deviations (we skip the corresponding logarithmic factor).

We emphasize that the difference between centralized (or parallel) estimates and obtained decentralized ones is not only in the replacement of $d$ by $\tilde{O}(\sqrt{\chi})$ in the smooth cases for the primal oracle and the meaning of $L$. In the stochastic smooth strongly convex case (one can also consider the convex case) we know that the total number of primal oracle calls [Kulunchakov and Mairal, 2019a, Kulunchakov and Mairal, 2019b, Kulunchakov and Mairal, 2019c, Lan and Zhou, 2018] is

[TABLE]

This bound is optimal but it uses an incremental oracle and does not imply full parallelization. The best known way to parallelize it is described in [Lan and Zhou, 2018]. For full parallelization one should use a standard accelerated scheme without variance reduction and incremental oracle [Woodworth et al., 2018]. In this case, another bound for the total number of oracle calls occurs, that is

[TABLE]

But this bound assumes the natural way of parallelization or centralized distribution of calculations. In the last case for a graph of diameter $d$ with $m$ nodes we have the following number of oracle calls per node

[TABLE]

and the following number of communication steps

[TABLE]

For decentralized architecture (see Table 4) the number of oracle calls per node and the number of communication steps are

[TABLE]

respectively. Unfortunately, the factor $m$ is no longer presented in $\sigma^{2}$ in the decentralized case. It is interesting to note, that it is possible to propose such a decentralized distributed algorithm that requires

[TABLE]

oracle calls per node (stochastic gradients calculations) [Olshevsky et al., 2019a, Olshevsky et al., 2019b]. However, this algorithm is not optimal in terms of communication steps. Moreover, to the best of our knowledge, it is an open question whether is (21) optimal bound in terms of oracle calls (per node) in the class of methods with the optimal number of communication steps.

Note also that the blue bound in Table 6 seems to be rather unexpected at first sight for us. But we hypothesize that this bound is optimal not only in terms of the number of communication steps but also in terms of the number of oracle calls (per node) in the class of methods with the optimal number of communication steps.

The detailed proofs of the statements collected in this paper takes more than 90 pages. These can be found in the arXiv preprint [Gorbunov et al., 2019]:

https://arxiv.org/pdf/1911.07363.pdf.

7 Discussion

Below we outline various areas for further work.

• •

We can expect that the results can be improved by replacing the first-order methods with primal and dual deterministic oracle by tensor methods ($p=2,3$) from [Nesterov, 2018a]. However, for the moment we do not know any such results. For the dual approach we also do not know how to use the trick $A=\sqrt{W}$ (see [Scaman et al., 2017]). Here we should take $A=W$, which (with additional increased complexity of auxiliary problem) makes the bounds on communication steps worse. The basic fact in the dual approach is the following. To solve auxiliary problems we have to calculate the values of the form $\nabla^{p}_{y}\varphi(Wy)$ on different vectors [Carmon and Duchi, 2016, Nesterov, 2018a, Nesterov, 2018b]. This can be done multiplying $W$ by vectors (communications) and multiplying corresponding (block) diagonal tensor ( $\nabla^{p}_{\lambda}\varphi(\lambda)|_{\lambda=Wy}$) by vectors (can be distributed among nodes)

• •

The primal approach in the smooth case can be generalized for the (stochastic inexact) gradient-free oracle. The number of communication steps remains the same. The number of oracle calls becomes $\sim n$ times larger [Gorbunov et al., 2018, Dvurechensky et al., 2017]. In the non-smooth case gradient-free (stochastic) decentralized distributed algorithm was developed in [Beznosikov et al., 2019]

• •

In [Hendrikx et al., 2019a], [Hendrikx et al., 2019b] asynchronized distributed optimization was considered via dual accelerated (block) coordinate descent algorithms. The primal approach proposed above allows asynchronized generalizations in the smooth case. For that we should use the (block) coordinate version of STM [Dvurechensky et al., 2017] and additional randomization of sum type when $[Wx]_{i}$ is calculated. This will increase the number of communication steps $\sim\sqrt{n}\div n$ times

• •

Most of the results of this paper can be generalized to composite problems [Nesterov, 2013]. Perhaps, it is possible to make the next step and try to generalize these results to more general types of models [Stonyakin et al., 2019b, Stonyakin et al., 2019a]

• •

For smooth convex centralized distributed optimization problems there exists a universal way to accelerate non-accelerated (stochastic, asynchronized etc.) algorithm Catalyst [Lin et al., 2015]. The basic idea is using a non-accelerated centralized distributed algorithm for the inner problem arising at each step of the Catalyst procedure

• •

Perhaps, it is possible to generalize the primal approach described above on time-varying graphs [Rogozin and Gasnikov, 2019]. Moreover, these generalizations can be done also for the smooth stochastic case

• •

It seems the result of [Rogozin et al., 2018] can be improved by using mixed communication: many decentralized steps alternate with centralized ones. In this case, one can use non-accelerated distributed decentralized algorithms, which are robust on time-varying graphs [Rogozin et al., 2018], and then accelerate them by using the Catalyst technique [Lin et al., 2015] and centralization. Since the graph is changing we should recalculate spanning tree whenever we apply a centralized step. We expect that this mixed communication will be useful also for tensor schemes in decentralized optimization

• •

The main scheme in the primal approach is based on the result formulated directly after (14). This result does not depend on convexity of the objective. So it would be interesting to apply this scheme for non-convex distributed optimization problems [Sun and Hong, 2018]

Since the first version of this paper was submitted on arXiv, there appeared alternative explanations (for smooth problems) of the results for the primal deterministic oracle [Fallah et al., 2019, Kovalev et al., 2020, Li and Lin, 2020, Xu et al., 2019, Hendrikx et al., 2020a] and the primal stochastic oracle (strongly convex case) [Fallah et al., 2019]. Moreover, in [Rogozin et al., 2020, Ye et al., 2020] (see also [Rogozin and Gasnikov, 2019] for the non-accelerated case and [Scaman et al., 2018] for lower bounds) for the primal deterministic oracle (strongly convex case) it was shown that $L:=\max_{k}L_{k}$ used in this paper can be improved to $L:=L_{f}$, where $L_{f}$ is the Lipschitz gradient constant of (6) (which can be much smaller [Tang et al., 2019]). The same holds true for $\mu$. More interestingly, we expect that combinations of [Rogozin and Gasnikov, 2019, Rogozin et al., 2020, Ye et al., 2020] allows to develop accelerated decentralized distributed primal algorithms on time-varying graphs (accelerated in $L$ and $\varepsilon$, but not in $\chi$). Moreover, based on [Rogozin and Gasnikov, 2019, Rogozin et al., 2020, Ye et al., 2020] we may expect that the answer to the open problem posed in the Section 6 is negative. The bound can be decreased by a factor $m$ for a stochastic primal oracle. The reasons for that are almost the same that we have for the $\max_{k}L_{k}\to L_{f}$ improvement.

Based on the recent works [Koloskova et al., 2020, Woodworth et al., 2020b, Woodworth et al., 2020a] we describe below a hypothesis about the optimal bounds in a more general situation. We consider the general sum-type problem (for simplicity we consider the case $Q=\mathbb{R}^{n}$, but it can be naturally generalized)

[TABLE]

We assume that all $f_{k}(x,\xi)$ in (22) satisfy

[TABLE]

Also we introduce131313Note that nowadays it’s quite popular to obtain estimates on rate of convergence that depend on the constants defined at the solution point $x^{*}$ [Stonyakin et al., 2020].

[TABLE]

and

[TABLE]

Assume now that at each iteration $t$ we may call one time an oracle (that returns an independent realization of $\nabla f_{k}(x^{t},\xi^{t,k})$) at each node and make no more than one communication step with (in general, random) communication matrix $\bar{W}_{k}$. Moreover, we assume that [Koloskova et al., 2020]141414In [Koloskova et al., 2020] it was also assumed that symmetric matrix $\bar{W}$ determined by doubly stochastic matrix $P$: $\bar{W}=P-I$, where $I$ – unit matrix. So do we. But below, when we generalize the results of [Koloskova et al., 2020] on a large class of algorithms that assumes more than one communication on one iteration, we may consider $\bar{W}$ to be the same as in section 5 and $\chi$ is determined in (20) style.

[TABLE]

[TABLE]

For the class of algorithms described above the required number of iterations $N$ to achieve an accuracy $\mathbb{E}[f(x^{N})]-f(x^{*})\leq\varepsilon$ for the best known (for the moment) algorithms is

[TABLE]

where $\beta\in[0,1]$, $\alpha=1\text{ or }1/2,$ here $I$ is a function such that $I[\rm true]=1$, $I[\rm false]=0$, and when $f_{k}$ is $\mu$-strongly convex

[TABLE]

This bound with $\beta=1,{\alpha=1}$ (non-accelerated case) was obtained in [Koloskova et al., 2020] for simple (local) decentralized SGD. Roughly speaking, $\beta=0$ corresponds to the lower bound for this algorithm [Koloskova et al., 2020, Yuan and Ma, 2020, Woodworth et al., 2020b, Woodworth et al., 2020a]. In [Karimireddy et al., 2019] by using a variance reduction trick in federated learning architecture ($\bar{W}_{p\tilde{\tau}+q}\equiv 0$ for $q=1,...,\tilde{\tau}-1$ and $\bar{W}_{p\tilde{\tau}}$ is a full communication matrix, hence $\chi=1$, $\tau=\tilde{\tau}$), the SCAFFOLD algorithm was proposed with $\alpha=1$ and (under some additional assumptions) without middle (blue) terms. Also in the federated learning architecture in case $f_{1}=...=f_{m}$ it seems that middle (blue) terms and the factor $\tau$ in the last terms of the described bounds can be eliminated under some natural additional assumptions [Woodworth et al., 2020b]. Here we assume at least one communication. If there is exactly one communication it will take place at the very end. This result means that in a federated learning setup under $f_{1}=...=f_{m}$ the frequency of communications $\tau$ does not play any significant role, which was previously mentioned in [Godichon-Baggioni and Saadane, 2020]. In centralized architecture it seems that it is possible to obtain additional acceleration of the algorithms mentioned above ($\alpha=1/2$, but the transition $\chi\to\sqrt{\chi}$ in the last terms requires additional assumptions) by using proper accelerated envelopes like Catalyst [Mathieu et al., 2020, Dvinskikh et al., 2020, Ivanova et al., 2019, Kulunchakov and Mairal, 2019c] with the precision $\varepsilon^{\prime}=\tilde{O}\left(\varepsilon\sqrt{\max\{\varepsilon/(LR^{2}),\mu/L}\}\right)$ (in the function) of the solution of auxiliary problem at each outer iteration. It is highly likely that this is true for more general architectures, see examples in [Hendrikx et al., 2020a, Li and Lin, 2020].

If we remove the condition that at each iteration we can only call stochastic oracle one time and make no more than one communication, then the blue term can be eliminated with $\alpha=1/2$. More precisely, the first two terms correspond to the total number of oracle calls per node and the last term – to the number of communications steps. Up to a factor $m$ in the denominator we have developed these results in the paper, but with $\tau=1$ and fixed $W$.151515For time-varying communication graphs in the eneral case for the moment it seems that we should put the factor $\chi$ instead of $\sqrt{\chi}$ in the last term [Rogozin and Gasnikov,]. We also note that it is an open problem whether it is possible in the general (accelerated) situation ( which differs from the [Koloskova et al., 2020]) to determine $\bar{\zeta}^{2}$ and $\bar{\sigma}^{2}$ only at point $x^{*}$? It seems, that for the current moment of time we have a positive answer only with respect to $\bar{\zeta}^{2}$.

The following generalizations are related with the case

[TABLE]

In this case with additional assumptions about proximal and dual friendly $f_{k}$ it is possible to reduce worth case constant $L$ to the ¡¡average¿¿ one [Hendrikx et al., 2020b]. In [Hendrikx et al., 2020b] this looks like a variance reduction acceleration, but the nature of the effect (also explained in [Hendrikx et al., 2020b]) based on coordinate descent acceleration [Nesterov and Stich, 2017] for the dual problem formulation. In [Hendrikx et al., 2020a] was proposed a dual-free generalization of [Hendrikx et al., 2020b] with a bit worse oracle complexity estimate. The main idea is to apply a non-accelerated coordinate descent for the dual problem with Bregman divergence determined by the dual function itself. In [Li et al., 2020] an optimal algorithm both for communications steps and oracle calls per node was developed.

Another way to obtain a better result for the required number of communication steps is possible if $\{f_{kj}\}$ have i.i.d. nature (that is typically for data science applications). In this case $f_{k}$ is statistically similar to $f$. Based on this fact in centralized architecture, we may use on a master node statistical preconditioned algorithms that can significantly reduce the required number of communications with slaves [Hendrikx et al., 2020c].

Acknowledgment

The work of D. Dvinskikh in Sections 1 and 5 was supported by RFBR 19-31-51001 and in Section 6 was funded by the Russian Science Foundation (project 18-71-10108). The work of A. Gasnikov was supported by the Ministry of Science and Higher Education of the Russian Federation (Goszadaniye) no. 075-00337-20-03, project No. 0714-2020-0005.

We would like to thank F. Bach, P. Dvurechensky, E. Gorbunov, A. Koloskova, A. Kulunchakov, J. Mairal, A. Nemirovski, A. Olshevsky, S. Parsegov, B. Polyak, N. Srebro A. Taylor and C. Uribe for useful discussions.