Optimal Decentralized Distributed Algorithms for Stochastic Convex Optimization

Introduction

In this paper, we are interested in the convex optimization problem

where $\xi$ is a random variable. Problems of this type play a central role in a bunch of applications of machine learning and mathematical statistics . Typically $x$ represents the feature vector defining the model, only samples of $\xi$ are available and the distribution of $\xi$ is unknown. One possible way to minimize generalization error (2) is to solve empirical risk minimization or finite-sum minimization problem instead, i.e. solve (1) with the objective

where $m$ should be sufficiently large to approximate the initial problem. Indeed, if $f(x,\xi)$ is convex and $M$-Lipschitz continuous for all $\xi$, $Q$ has finite diameter $D$ and $\hat{x}=\mathop{\arg\!\min}_{x\in Q}\hat{f}(x)$, then (see ) with probability at least $1-\beta$

and if additionally $f(x,\xi)$ is $\mu$-strongly convex for all $\xi$, then (see ) with probability at least $1-\beta$

In other words, to solve (1)+(2) with $\varepsilon$ functional accuracy via minimization of empirical risk (3) it is needed to have $m=\widetilde{\Omega}\left(\nicefrac{{M^{2}D^{2}n}}{{\varepsilon^{2}}}\right)$ in the convex case and $m=\widetilde{\Omega}\left(\max\left\{\nicefrac{{M^{2}D^{2}}}{{\mu\varepsilon}},\nicefrac{{M^{2}D^{2}}}{{\varepsilon^{2}}}\right\}\right)$ in the $\mu$-strongly convex case where $\widetilde{\Omega}(\cdot)$ hides a constant factor, a logarithmic factor of $\nicefrac{{1}}{{\beta}}$ and a polylogarithmic factor of $\nicefrac{{1}}{{\varepsilon}}$.

Stochastic first-order methods such as Stochastic Gradient Descent (SGD) or its accelerated variants like AC-SA or Similar Triangles Method (STM) are very popular choice to solve either (1)+(2) or (1)+(3). In contrast with their cheap iterations in terms of computational cost, these methods converge only to the neighbourhood of the solution, i.e. to the ball centered at the optimality and radius proportional to the standard deviation of the stochastic estimator. For the particular case of finite-sum minimization problem one can solve this issue via variance-reduction trick and its accelerated variants . Unfortunately, this technique is not applicable in general for the problems of type (1)+(2). Another possible way to reduce the variance is mini-batching. When the objective function is $L$-smooth one can accelerate the computations of batches using parallelization , and it is one of the examples where centralized distributed optimization appears naturally .

In other words, in some situations, e.g. when the number of samples $m$ is too big, it is preferable in practice to split the data into $q$ blocks, assign each block to the separate worker, e.g. processor, and organize computation of the gradient or stochastic gradient in the parallel or distributed manner. Moreover, in view of (4)-(5) sometimes to solve an expectation minimization problem it is needed to have such a big number of samples that corresponding information (e.g. some objects like images, videos and etc.) cannot be stored on $1$ machine because of the memory limitations (see Section 8 for the detailed example of such a situation). Then, we can rewrite the objective function in the following form

Here $f_{i}$ corresponds to the loss on the $i$-th data block and could be also represented as an expectation or a finite sum. So, the general idea for parallel optimization is to compute gradients or stochastic gradients by each worker, then aggregate the results by the master node and broadcast new iterate or needed information to obtain the new iterate back to the workers.

The visual simplicity of the parallel scheme hides synchronization drawback and high requirement to master node . The big line of works is aimed to solve this issue via periodical synchronization , error-compensation , quantization or combination of these techniques .

However, in this paper we mainly focus on another approach to deal with aforementioned drawbacks — decentralized distributed optimization . It is based on two basic principles: every node communicates only with its neighbours and communications are performed simultaneously. Moreover, this architecture is more robust, e.g. it can be applied to time-varying (wireless) communication networks .

One can consider this paper as a continuation of work where authors mentioned the key ideas that form a basis of this work. However, in this paper we provide formal proofs of some results announced in together with couple of new results that were not mentioned. Our contributions include:

Accelerated primal-dual method with biased stochastic dual oracle for convex and smooth dual problem. We extent the result from the recent work to the case when we have an access to the biased stochastic gradients. We emphasize that our analysis works for the minimization on whole space and we do not assume that the sequence generated by the method is bounded. It creates extra difficulties in the analysis, but we handle it via advanced technique for estimating recurrences (see also ).

Two accelerated methods with stochastic dual oracle for strongly convex and smooth dual problem. For the case when the dual function is strongly convex with Lipschitz continuous gradient we analyze two methods: one is R-RRMA-AC-SA2 and another is SSTM_sc. The first one was described in , but in this paper we formally state the method and prove high probability bounds for its convergence rate. The second method is also well-known, but to the best of our knowledge there were no convergence results for it in such generality that we handle. That is, we consider SSTM_sc with biased stochastic oracle applied to the unconstrained smooth and strongly convex minimization problem and prove high probability bounds for its convergence rate together with the bound for the noise level. As for the convex case, we also do not assume that the sequence generated by the method is bounded. Then we show how it can be applied to solve stochastic optimization problem with affine constraints using dual oracle.

Analysis of STM applied to convex smooth minimization problem with smooth convex composite term and inexact proximal step for unconstrained minimization. Surprisingly, but before this paper there were no analysis for STM in this case. The closest work to ours in this topic is , but in authors considered optimization problems on bounded sets.

2 Outline of the Paper

In Section 2, we introduce the notation and main definitions used in the paper. Then, we discuss optimal bounds for stochastic convex optimization in Section 3. In Section 4, we present the stochastic optimization problems with affine constraints and the state-of-the-art methods that solve the specific penalized unconstrained problem instead of the original one together with the novel approach which we call STP_IPS that aims to solve convex smooth unconstrained minimization problems with the smooth convex composite term and inexact proximal step. Next, we consider the same type of problems but using a dual approach and develop three different accelerated methods for this case together with the convergence analysis for each of them (Section 5). The first one is Stochastic Primal-Dual STM (SPDSTM), and it uses a biased stochastic dual oracle to solve primal and dual problems simultaneously for the case when the primal problem is $\mu$-strongly convex and Lipschitz continuous on some ball centered at zero. The next two methods are R-RRMA-AC-SA2 and SSTM_sc, and they solve the same problem when the primal functional is additionally $L$-smooth using a stochastic dual oracle. The difference between them is that R-RRMA-AC-SA2 uses special restarts technique and works with unbiased stochastic oracle, while SSTM_sc is directly accelerated and able to work with biased stochastic gradients. Then we show how to apply the results of the previous sections to the decentralized distributed optimization problems and derive the bounds for the proposed methods in Section 6. Finally, in Section 7, we compare bounds for the convergence rate in parallel and decentralized optimization, discuss the optimality of the obtained results, and present possible directions for future work. To illustrate how our theory works, we consider the problem of calculation of population Wasserstein barycenter in Section 8. We leave long proofs, auxiliary and technical results, and the whole section about STP_IPS in the appendix.

Notation and Definitions

Below we list some classical definitions for optimization (see, for example, for the details).

If $\mu>0$ then there exists unique minimizer of $f$ on $Q$ which we denote by $x^{*}$, except the situations when we explicitly specify $x^{*}$ in a different way. In the case when $\mu=0$, i.e. $f$ is convex, we assume that there exists at least one minimizer $x^{*}$ of $f$ on $Q$ and in the case when the set of minimizers of $f$ on the set $Q$ is not a singleton we choose $x^{*}$ to be either arbitrary or closest to the starting point of a method. When we consider some optimization method with a starting point $x^{0}$ we use $R$ or $R_{0}$ to denote the Euclidean distance between $x^{0}$ and $x^{*}$.

Optimal Bounds for Stochastic Convex Optimization

In this section our goal is to present the overview of the optimal methods and their convergence rates for the stochastic convex optimization problem (1)+(2) in the case when the gradient of the objective function is available only through (possibly biased) stochastic estimators with “light tails” or, equivalently, with $\sigma^{2}$-sub-Gaussian variance. That is, we are interested in the situation when for an arbitrary $x\in Q$ one can get such stochastic gradient $\nabla f(x,\xi)$ that

In this paper we are mainly focus on smooth optimization problems and use different modifications of Similar Triangles Method (STM) since it gives optimal rates in this case and it is easy enough to analyze at least in the deterministic case. For convenience, we state the method in this section as Algorithm 1.

Interestingly, if we run STM with $\mu>0$ to solve (1) with $\mu$-strongly convex and $L$-smooth objective, it will return $x^{N}$ such that $f(x^{N})-f(x^{*})\leq\varepsilon$ after $N=O\left(\sqrt{\nicefrac{{L}}{{\mu}}}\ln\left(\nicefrac{{LR^{2}}}{{\varepsilon}}\right)\right)$ iterations which is not optimal, seeIn some places we put references not to the first work where this bound was shown but to the works where this complexity bound was shown for either more convenient or more relevant to our work method. Table 1. To match the optimal bound in this case one should use classical restart of STM which is run with $\mu=0$ .

We notice that another highly widespread in machine learning applications type of problems is regularized or composite optimization problem

where $h$ is a convex proximable function. For this case STM can be generalized via modifying the update rule in the following way :

We address such problems with $L_{h}$-smooth composite term in the Appendix, see Section E for the details.

Next, we go back to the problem (1)+(2) and consider more general case when $\delta=0$ and $\sigma^{2}>0$. In this case one can construct unbiased estimator

where $\xi_{1},\ldots,\xi_{r}$ are i.i.d. samples and $\nabla f(x,\{\xi_{i}\}_{i=1}^{r})$ has $r$ times smaller variance than $\nabla f(x,\xi_{i})$:

Then in order to get such a point $x^{N}$ that $f(x^{N})-f(x^{*})\leq\varepsilon$ with probability at least $1-\beta$ where $\beta\in(0,1)$ and $f$ is $\mu$-strongly convex ($\mu\geq 0$) and $L$-smooth one can run STM for

which is optimal up to logarithmic factors. We call this modification Stochastic STM (SSTM). As for the deterministic case we summarize the state-of-the-art results for this case in Table 2.

Stochastic Convex Optimization with Affine Constraints: Primal Approach

Now, we are going to make a step towards decentralized distributed optimization and consider convex optimization problem with affine constraints:

where $A\succeq 0$ and $\text{Ker}A\neq\{0\}$. Up to a sign we can define the dual problem in the following way

However, in this section we are interested only in primal approaches to solve (17) and, in particular, the main goal of this section is to present first-order methods that are optimal both in terms of $\nabla f(x)$ and $A^{\top}Ax$ calculations. Before we start our analysis let us notice that typically in decentralized optimization matrix $A$ from (17) is chosen as a square root of Laplacian matrix $W$ of communication network (see Section 6 for the details). In asynchronous case the square root $\sqrt{W}$ is replaced by incidence matrix $M$ ($W=M^{\top}M$). Then in asynchronous case instead of accelerated methods for (18) one should use accelerated block-coordinate descent methods .

To solve problem (17) we use the following trick : instead of (17) we consider penalized problem

where $\varepsilon>0$ is the desired accuracy of the solution in terms of $f(x)$ that we want to achieve. The motivation behind this trick is revealed in the following theorem.

We start with the analysis of the case when $f$ is $L$-smooth and convex.

to produce point $x^{N}$ such that (22) holds.

That is, number of $A^{\top}Ax$ calculations matches the optimal bound for deterministic convex and $L$-smooth problems of type (1) multiplied by $\sqrt{\chi(A^{\top}A)}$ up to logarithmic factors (see Table 1).

We conjecture that the same technique in the case when $f$ is $\mu$-strongly convex and $L$-smooth gives the method that requires such number of $A^{\top}Ax$ calculations that matches the second rows of Tables 1 and 2 in the corresponding cases with additional factor $\sqrt{\chi(A^{\top}A)}$ and logarithmic factors. Recently such bounds were shown in for the distributed version of Multistage Accelerated Stochastic Gradient method from . However, this bounds were shown for the case when the stochastic gradient is unbiased.

Next, we assume that $Q$ is closed and convex and $f$ is $\mu$-strongly convex, but possibly non-smooth function with bounded gradients: $\|\nabla f(x)\|_{2}\leq M$ for all $x\in Q$. Let us start with the case $\mu=0$. Then, to achieve (22) one can run Sliding method from considering $f(x)$ as a composite term. In this case Sliding requires

In the case when $Q$ is a compact set and $\nabla f(x)$ is not available and unbiased stochastic gradient $\nabla f(x,\xi)$ is used instead (see inequalities (10)-(11) with $\delta=0$) one can show that Stochastic Sliding (S-Sliding) method can achieve (22) with probability at least $1-\beta$, $\beta\in(0,1)$, and it requires the same number of calculations of $A^{\top}Ax$ as in (26) up to logarithmic factors and

When $\mu>0$ one can apply restarts technique on top of S-Sliding (RS-Sliding) and get that to guarantee (22) with probability at least $1-\beta$, $\beta\in(0,1)$ RS-Sliding requires

We notice that bounds presented above for the non-smooth case are proved only for the case when $Q$ is bounded. For the case of unbounded $Q$ the convergence results with such rates were proved only in expectation. Moreover, it would be interesting to study S-Sliding and RS-Sliding in the case when $\delta>0$, i.e. stochastic gradient is biased, but we leave these questions for future works.

Stochastic Convex Optimization with Affine Constraints: Dual Approach

We notice that in this section we do not assume that $A$ is symmetric or positive semidefinite.

Assume additionally that $x(y,\xi)$ satisfies so-called “light-tails” inequality:

The size of the batch $r_{k}$ could always be restored from the context, so, we do not specify it here. Note that the batch version satisfies

Below we present the main convergence result of this section.

with probability at least $1-4\beta$. What is more, to guarantee (40) with probability at least $1-4\beta$ Algorithm 2 requires

2 Strongly Convex Dual Functions and Restarts Technique

In this section we assume that primal functional $f$ is additionally $L$-smooth. It implies that the dual function $\psi$ in (18) is additionally $\mu_{\psi}$-strongly convex in $y^{0}+(\text{Ker}A^{\top})^{\perp}$ where $\mu_{\psi}=\nicefrac{{\lambda_{\min}^{+}(A^{\top}A)}}{{L}}$ and $\lambda_{\min}^{+}(A^{\top}A)$ is the minimal positive eigenvalue of $A^{\top}A$.

From weak duality $-f(x^{*})\leq\psi(y^{*})$ and (20) we get the key relation of this section (see also )

This inequality implies the following theorem.

Consider function $f$ and its dual function $\psi$ defined in (20) such that problems (17) and (18) have solutions. Assume that $y^{N}$ is such that $\|\nabla\psi(y^{N})\|_{2}\leq\nicefrac{{\varepsilon}}{{R_{y}}}$ and $y^{N}\leq 2R_{y}$, where $\varepsilon>0$ is some positive number and $R_{y}=\|y^{*}\|_{2}$ where $y^{*}$ is any minimizer of $\psi$. Then for $x^{N}=x(A^{\top}y^{N})$ following relations hold:

Applying Cauchy-Schwarz inequality to (42) we get

The second part (43) immediately follows from $\|\nabla\psi(y^{N})\|_{2}\leq\nicefrac{{\varepsilon}}{{R_{y}}}$ and Demyanov-Danskin theorem which implies $\nabla\psi(y^{N})=Ax^{N}$. ∎

That is why, in this section we mainly focus on the methods that provides optimal convergence rates for the gradient norm. In particular, we consider Recursive Regularization Meta-Algorithm from (see Algorithm 3) with AC-SA2 (see Algorithm 5) as a subroutine (i.e. RRMA-AC-SA2) which is based on AC-SA algorithm (see Algorithm 4) from . We notice that RRMA-AC-SA2 is applied for a regularized dual function

In this section we consider the same oracle as in Section 5, but we additionally assume that $\delta=0$, i.e. stochastic first-order oracle is unbiased. To define batched version of the stochastic gradient we will use the following notation:

As before in the cases when the batch-size $r_{t}$ can be restored from the context, we will use simplified notation $\nabla\Psi(y,\boldsymbol{\xi}^{t})$ and $x(y,\boldsymbol{\xi}^{t})$.

In the AC-SA algorithm we use batched stochastic gradients of functions $\psi_{k}$ which are defined as follows:

The following theorem states the main result for RRMA-AC-SA2 that we need in the section.

Let $\psi$ be $L_{\psi}$-smooth and $\mu_{\psi}$-strongly convex function and $\lambda=\Theta\left(\nicefrac{{(L_{\psi}\ln^{2}N)}}{{N^{2}}}\right)$ for some $N>1$. If the Algorithm 3 performs $N$ iterations in totalThe overall number of performed iterations during the calls of AC-SA2 equals $N$. with batch size $r$ for all iterations, then it will provide such a point $\hat{y}$ that

where $C>0$ is some positive constant and $y^{*}$ is a solution of the dual problem (18).

Let us show that w.l.o.g. we can assume in this section that function $\psi$ defined in (20) is $\mu_{\psi}$-strongly convex everywhere with $\mu_{\psi}=\nicefrac{{\lambda_{\min}^{+}(A^{\top}A)}}{{L}}$. In fact, from $L$-smoothness of $f$ we have only that $\psi$ is $\mu_{\psi}$-strongly convex in $y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$ (see for the details). However, the structure of the considered here methods is such that all points generated by the RRMA-AC-SA2 and, in particular, AC-SA lie in $y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$.

We prove the statement of the theorem by induction. For $t=0$ the statement is trivial, since $y^{0}_{md}=y^{0}_{ag}=z^{0}\in y_{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. Assume that $y_{md}^{t},z^{t},y_{ag}^{t}\in y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$ for some $t\geq 0$ and prove it for $t+1$. Since $y_{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$ is a convex set and $y^{t+1}_{md}$ is a convex combination of $y^{t}_{ag}$ and $z^{t}$ we have $y^{t+1}_{md}\in y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. Next, the point $\frac{\alpha_{t}\lambda}{\lambda+\gamma_{t}}y^{t+1}_{md}+\frac{(1-\alpha_{t})\lambda+\gamma_{t}}{\lambda+\gamma_{t}}z^{t}$ also lies in $y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$ since it is convex combination of the points lying in this set. Due to (44), (45) and (46) we have that $\nabla\Psi_{k}(y_{md}^{t+1},\boldsymbol{\xi}^{t})=Ax(A^{\top}y_{md}^{t+1},\boldsymbol{\xi}^{t})+\lambda(y_{md}^{t+1}-y^{0})+\lambda\sum_{l=1}^{k}2^{l}(y_{md}^{t+1}-\hat{y}^{l})$. The first term lies in $\left(\text{Ker}(A^{\top})\right)^{\perp}$ since $\text{Im}(A)=\left(\text{Ker}(A^{\top})\right)^{\perp}$ and the second and the third terms also lie in $\left(\text{Ker}(A^{\top})\right)^{\perp}$ since $y_{md}^{t+1},y^{0},\hat{y}^{1},\ldots,\hat{y}^{k}\in y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. Putting all together we get $z^{t+1}\in y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. Finally, $y_{ag}^{t+1}$ lies in $y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$ as a convex combination of points from this set. ∎

We prove this result by induction. For $t=0$ the statement is trivial since $\hat{y}^{0}=y^{0}$. Next, assume that $\hat{y}^{0},\hat{y}^{1},\ldots,\hat{y}^{k}\in y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$ and prove that $\hat{y}^{k+1}\in y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. Our assumption implies that the assumptions from Theorem 5.4 and applying the result of the theorem we get that $y^{1}$ and $y^{2}$ from the method AC-SA2 applied to the $\psi_{k}$ also lie in $y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. That is, the output of AC-SA2 applied for $\psi_{k}$ lies in $y^{0}+\left(\text{Ker}(A^{\top})\right)^{\perp}$. ∎

Now we are ready to present our approach which was sketched in of constructing an accelerated method for the strongly convex dual problem using restarts of RRMA-AC-SA2. To explain the main idea we start with the simplest case: $\sigma_{\psi}^{2}=0$, $r=0$. It means that there is no stochasticity in the method and the bound (47) can be rewritten in the following form:

In the case when $\sigma_{\psi}^{2}\neq 0$ we need to modify this approach. The first ingredient to handle the stochasticity is large enough batch size for the $l$-th restart: $r_{l}$ should be $\Omega\left(\nicefrac{{\sigma_{\psi}^{2}}}{{(\bar{N}\|\nabla\psi(\bar{y}^{l-1})\|_{2}^{2})}}\right)$. However, in the stochastic case we do not have an access to the $\nabla\psi(\bar{y}^{l-1})$, so, such batch size is impractical. One possible way to fix this issue is to independently sample large enough number $\hat{r}_{l}\sim\nicefrac{{R_{y}^{2}}}{{\varepsilon^{2}}}$ of stochastic gradients additionally, which is the second ingredient of our approach, in order to get good enough approximation $\nabla\Psi(\bar{y}^{l-1},\boldsymbol{\xi}^{l-1},\hat{r}_{l})$ of $\nabla\psi(\bar{y}^{l-1})$ and use the norm of such an approximation which is close to the norm of the true gradient with big enough probability in order to estimate needed batch size $r^{l}$ for the optimization procedure. Using this, we can get the bound of the following form:

The third ingredient is the amplification trick: we run $p_{l}=\Omega(\ln(\nicefrac{{1}}{{\beta}}))$ independent trajectories of RRMA-AC-SA2, get points $\bar{y}^{l,1},\ldots,\bar{y}^{l,p_{l}}$ and choose such $\bar{y}^{l,p(l)}$ among of them that $\|\nabla\psi(\bar{y}^{l,p(l)})\|_{2}$ is close enough to $\min\limits_{p=1,\ldots,p_{l}}\|\nabla\psi(\bar{y}^{l,p})\|_{2}$ with high probability, i.e. $\|\nabla\psi(\bar{y}^{l,p(l)})\|_{2}^{2}\leq 2\min\limits_{p=1,\ldots,p_{l}}\|\nabla\psi(\bar{y}^{l,p})\|_{2}^{2}+\nicefrac{{\varepsilon^{2}}}{{8R_{y}^{2}}}$ with probability at least $1-\beta$ for fixed $\nabla\Psi(\bar{y}^{l-1},\boldsymbol{\xi}^{l-1},\hat{r}_{l})$. We achieve it due to additional sampling of $\bar{r}_{l}\sim\nicefrac{{R_{y}^{2}}}{{\varepsilon^{2}}}$ stochastic gradients at $\bar{y}^{l,p}$ for each trajectory and choosing such $p(l)$ corresponding to the smallest norm of the obtained batched stochastic gradient. By Markov’s inequality for all $p=1,\ldots,p_{l}$

That is, for $p_{l}=\log_{2}(\nicefrac{{1}}{{\beta}})$ we have that with probability at least $1-2\beta$

for fixed $\nabla\Psi(\bar{y}^{l-1},\boldsymbol{\xi}^{l-1},\hat{r}_{l})$ which means that

with probability at least $1-3\beta$. Therefore, after $l=\log_{2}(\nicefrac{{2R_{y}^{2}\|\nabla\psi(y^{0})\|_{2}^{2}}}{{\varepsilon^{2}}})$ of such restarts our method provide the point $\bar{y}^{l,p(l)}$ such that with probability at least $1-3l\beta$

The approach informally described above is stated as Algorithm 6.

Assume that $\psi$ is $\mu_{\psi}$-strongly convex and $L_{\psi}$-smooth. If Algorithm 6 is run with

for all $k=1,\ldots,l$ where $\bar{N}>1$ is such that $\frac{CL_{\psi}^{2}\ln^{4}\bar{N}}{\mu_{\psi}^{2}\bar{N}^{4}}\leq\frac{1}{32}$, $\beta\in(0,\nicefrac{{1}}{{3}})$ and $\varepsilon>0$, then with probability at least $1-3\beta$

and the total number of the oracle calls equals

Under assumptions of Theorem 5.6 we get that with probability at least $1-3\beta$

where $\beta\in(0,\nicefrac{{1}}{{3}})$ the total number of the oracle calls is defined in (51).

Inequalities (50) and $\mu_{\psi}\|y-y^{*}\|_{2}\leq\|\nabla\psi(y)\|_{2}$ which follows from $\mu_{\psi}$-strong convexity of $\psi$ imply that

Now we are ready to present convergence guarantees for the primal function and variables.

Let the assumptions of Theorem 5.6 hold. Assume that $f$ is $L_{f}$-Lipschitz continuous on $B_{R_{f}}(0)$ where

and $R_{x}=\|x(A^{\top}y^{*})\|_{2}$. Then, with probability at least $1-4\beta$

where $\beta\in(0,\nicefrac{{1}}{{4}})$, $\varepsilon\in(0,\mu_{\psi}R_{y}^{2})$ $x^{l}\stackrel{{\scriptstyle\text{def}}}{{=}}x(A^{\top}\bar{y}^{l,p(l)},\boldsymbol{\xi}^{l,p(l)},\bar{r}_{l})$ and to achieve it we need the total number of oracle calls equals

3 Direct Acceleration for Strongly Convex Dual Function

We consider first the following minimization problem:

We use Stochastic Similar Triangles Method which is stated in this section as Algorithm 7 to solve problem (55). To define the iterate $z^{k+1}$ we use the following sequence of functions:

Assume that Algorithm 7 is run to solve problem (55) with $\psi(y)$ being $\mu_{\psi}$-strongly convex and $L_{\psi}$-smooth. Then, for all $k\geq 0$ we have

for all $l=1,\ldots,N$, where $h,\delta,u$ and $c$ are some non-negative constants and $A_{k+1}=A_{k}+\alpha_{k+1}$, $\alpha_{k+1}\leq DA_{k}$ for some $D\geq 1$, $A_{0}=\alpha_{0}>0$. Assume that for each $k\geq 1$ vector $a^{k}$ is a function of $\eta^{0},\ldots,\eta^{k-1}$, $a^{0}$ is a deterministic vector, $u\geq 1$, sequence of random vectors $\{\eta^{k}\}_{k\geq 0}$ satisfy

hold for all $l=1,\ldots,N$ simultaneously, where $C_{1}$ is some positive constant, $g(N)=\frac{\ln\left(\frac{N}{\beta}\right)+\ln\ln\left(\frac{B}{b}\right)}{\left(1+\sqrt{3\ln\left(\frac{N}{\beta}\right)}\right)^{2}}$,

$b=2\sigma_{0}^{2}\alpha_{1}^{2}R_{0}^{2}$ and

Assume that the function $\psi$ is $\mu_{\psi}$-strongly convex and $L_{\psi}$-smooth,

i.e. $r_{k}\geq\frac{1}{C}\max\left\{1,\left(\frac{\mu_{\psi}}{L_{\psi}}\right)^{\nicefrac{{3}}{{2}}}\frac{N^{2}\sigma_{\psi}^{2}\left(1+\sqrt{3\ln\frac{N}{\beta}}\right)^{2}}{\varepsilon}\right\}$ with positive constants $C>0$, $\varepsilon>0$ and $N\geq 1$. If additionally $\delta\leq\frac{GR_{0}}{N\sqrt{A_{N}}}$ and $\varepsilon\leq\frac{HR_{0}^{2}}{A_{N}}$ where $R_{0}=\|y^{*}-y^{0}\|_{2}$ and Algorithm 7 is run for $N$ iterations, then with probability at least $1-3\beta$

and $C_{1}$ is some positive constant. In other words, to achieve $\|y^{N}-y^{*}\|_{2}^{2}\leq\varepsilon$ with probability at least $1-3\beta$ Algorithm 7 needs $N=\widetilde{O}\left(\sqrt{\frac{L_{\psi}}{\mu_{\psi}}}\right)$ iterations and $\widetilde{O}\left(\max\left\{\sqrt{\frac{L_{\psi}}{\mu_{\psi}}},\frac{\sigma_{\psi}^{2}}{\varepsilon}\right\}\right)$ oracle calls where $\widetilde{O}(\cdot)$ hides polylogarithmic factors depending on $L_{\psi},\mu_{\psi},R_{0},\varepsilon$ and $\beta$.

Next, we apply the SSTM_sc for the problem (18) when the objective of the primal problem (17) is $L$-smooth, $\mu$-strongly convex and $L_{f}$-Lipschitz continuous on some ball which will be specified next, i.e. we consider the same setup as in Section 5 but we additionally assume that the primal functional $f$ has $L$-Lipschitz continuous gradient. As in Section 5 we also consider the case when the gradient of the dual functional is known only through biased stochastic estimators, see (32)–(39) and the paragraphs containing these formulas.

This theorem makes it possible to apply the result from Theorem 5.11 for SSTM_sc which is run on the problem (18).

Under assumptions of Theorem 5.11 we get that after $N=\widetilde{O}\left(\sqrt{\frac{L_{\psi}}{\mu_{\psi}}}\ln\frac{1}{\varepsilon}\right)$ iterations of Algorithm 7 which is run on the problem (18) with probability at least $1-3\beta$

where $\beta\in\left(0,\nicefrac{{1}}{{3}}\right)$ and the total number of oracles calls equals

If additionally $\varepsilon\leq\mu_{\psi}R_{y}^{2}$, then with probability at least $1-3\beta$

Theorem 5.11 implies that with probability at least $1-3\beta$ we have

Using this and $L_{\psi}$-smoothness of $\psi$ we get that with probability $\geq 1-3\beta$

Since A\overset{\eqref{eq:A_k_lower_bound_str_cvx}}{\geq}\frac{1}{L_{\psi}}\left(1+\frac{1}{2}\sqrt{\frac{\mu_{\psi}}{L_{\psi}}}\right)^{2k}, it implies that after $N=\widetilde{O}\left(\sqrt{\frac{L_{\psi}}{\mu_{\psi}}}\ln\frac{1}{\varepsilon}\right)$ iterations of SSTM_sc we will get (65) with probability at least $1-3\beta$ and the number of oracle calls will be

Next, from $\mu_{\psi}$-strong convexity of $\psi(y)$ we have that with probability at least $1-3\beta$

and from this we obtain that with probability at least $1-3\beta$

Let the assumptions of Theorem 5.11 hold. Assume that $f$ is $L_{f}$-Lipschitz continuous on $B_{R_{f}}(0)$ where

$R_{x}=\|x(A^{\top}y^{*})\|_{2}$, $\varepsilon\leq\mu_{\psi}R_{y}^{2}$ and $\delta_{y}\leq\frac{G_{1}\varepsilon}{NR_{y}}$ for some positive constant $G_{1}$. Assume additionally that the last batch-size $r_{N}$ is slightly bigger than other batch-sizes, i.e.

Then, with probability at least $1-4\beta$

Applications to Decentralized Distributed Optimization

In this section we apply our results to the decentralized optimization problems. But let us consider first the centralized or parallel architecture. As we mentioned in the introduction, when the objective function is $L$-smooth one can compute batches in parallel in order to accelerate the work of the method and (14)-(16) imply that

number of workers in such a parallel scheme gives the method with working time proportional to the number of iterations defined in (14). However, number of workers defined in (73) could be too big in order to use such an approach in practice. But still computing the batches in parallel even with much smaller number of workers could reduce the working time of the method if the communication is fast enough and it follows from (16).

Besides the computation of batches in parallel for the general type of problem (1)+(2), parallel optimization is often applied to the finite-sum minimization problems (1)+(3) or (1)+(6) that we rewrite here in the following form:

We notice that in this section $m$ is a number of workers and $f_{k}(x)$ is known only for the $k$-th worker. Consider the situation when workers are connected in a network and one can construct a spanning tree for this network. Assume that the diameter of the obtained graph equals $d$, i.e. the height of the tree — maximal distance (in terms of connections) between the root and a leaf . If we run STM on such a spanning tree then we will get that the number of communication rounds will be $d$ times larger than number of iterations defined in (14).

For simplicity, we also call $W$ as a Laplacian matrix and it does not lead to misunderstanding since everywhere below we use $W$ instead of $\overline{W}$. The key observation here that computation of $Wx$ requires one round of communications when the $k$-th worker sends $x_{k}$ to all its neighbours and receives $x_{j}$ for all $j$ such that $(k,j)\in E$, i.e. $k$-th worker gets vectors from all its neighbours. Note, that $W$ is symmetric and positive semidefinite and, as a consequence, $\sqrt{W}$ exists. Moreover, we can replace $W$ by $\sqrt{W}$ in (77) and get the equivalent statement:

Using this we can rewrite the problem (74) in the following way:

with $\sigma_{f}^{2}=O\left(\nicefrac{{\sigma^{2}}}{{m}}\right)$.

Now it should become clear why in Section 4 we paid most of our attention on number of $A^{\top}A\mathbf{x}$ calculations. In this particular scenario $A^{\top}A\mathbf{x}=\sqrt{W}^{\top}\sqrt{W}x=Wx$ which can be computed via one round of communications of each node with its neighbours as it was mentioned earlier in this section. That is, for the primal approach we can simply use the results discussed in Section 4. For convenience, we summarize them in Tables 3 and 4 which are obtained via plugging the parameters that we obtained above in the bounds from Section 4. Note that the results presented in this match the lower bounds obtained in in terms of the number of communication rounds up to logarithmic factors and and there is a conjecture that these bounds are also optimal in terms of number of oracle calls per node for the class of methods that require optimal number of communication rounds. Recently, the very similar result about the optimal balance between number of oracle calls per node and number of communication round was proved for the case when the primal functional is convex and $L$-smooth and deterministic first-order oracle is available .

where $[\sqrt{W}\mathbf{x}]_{k}$ is the $k$-th $n$-dimensional block of $\sqrt{W}x$. Note that

Consider the stochastic function $f_{k}(x_{k},\xi_{k})$ which is defined implicitly as follows:

it is natural to define the stochastic gradient $\nabla\Phi(\mathbf{y},\xi)$ as follows:

with $\delta_{\Phi}=m\delta_{\varphi}$ and $\sigma_{\Phi}^{2}=O\left(m\sigma^{2}\right)$. Using this, we define the stochastic gradient of $\Psi(\mathbf{y})$ as $\nabla\Psi(\mathbf{y},\xi)\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{W}\nabla\Phi(\sqrt{W}\mathbf{y},\xi)=\sqrt{W}\mathbf{x}(\sqrt{W}\mathbf{y},\xi)$ and, as a consequence, we get

with $\delta_{\Psi}=\sqrt{\lambda_{\max}(W)}\delta_{\Phi}$ and $\sigma_{\Psi}=\sqrt{\lambda_{\max}(W)}\sigma_{\Phi}$.

Taking all of this into account we conclude that problem (82) is a special case of (18) with $A=\sqrt{W}$. To make the algorithms from Section 5 distributed we should change the variables in those methods via multiplying them by $\sqrt{W}$ from the left , e.g. for the iterates of SPDSTM we will get

which means that it is needed to multiply lines 4-6 of Algorithm 2 by $\sqrt{W}$ from the left. After such a change of variables all methods from Section 5 become suitable to run them in the distributed fashion. Besides that, it does not spoil the ability of recovering the primal variables since before the change of variables all of the methods mentioned in Section 5 used $\mathbf{x}(\sqrt{W}\mathbf{y})$ or $\mathbf{x}(\sqrt{W}\mathbf{y},\xi)$ where points $y$ were some dual iterates of those methods, so, after the change of variables we should use $\mathbf{x}(\mathbf{y})$ or $\mathbf{x}(\mathbf{y},\xi)$ respectively. Moreover, it is also possible to compute $\|\sqrt{W}x\|_{2}^{2}=\langle\mathbf{x},W\mathbf{x}\rangle$ in the distributed fashion using consensus type algorithms: one communication step is needed to compute $W\mathbf{x}$, then each worker computes $\langle x_{k},[W\mathbf{x}]_{k}\rangle$ locally and after that it is needed to run consensus algorithm. We summarize the results for this case in Tables 5 and 6. Note that the proposed bounds are optimal in terms of the number of communication rounds up to polylogarithmic factors . Note that the lower bounds from are presented for the convolution of two criteria: number of oracle calls per node and communication rounds. One can obtain lower bounds for the number of communication rounds itself using additional assumption that time needed for one communication is big enough and the term which corresponds to the number of oracle calls can be neglected. Regarding the number of oracle calls there is a conjecture that the bounds that we present in this paper are also optimal up to polylogarithmic factors for the class of methods that require optimal number of communication rounds.

Discussion

In this section we want to discuss some aspects of the proposed results that were not covered in the main part of this paper. First of all, we should say that in the smooth case for the primal approach our bounds for the number of communication steps coincides with the optimal bounds for the number of communication steps for parallel optimization if we substitute the diameter $d$ of the spanning tree in the bounds for parallel optimization by $\widetilde{O}(\sqrt{\chi(W)})$.

However, we want to discuss another interesting difference between parallel and decentralized optimization in terms of the complexity results which was noticed in . From the line of works it is known that for the problem (1)+(6) (here we use $m$ instead of $q$ and iterator $k$ instead of $i$ for consistency) with $L$-smooth and $\mu$-strongly convex $f_{k}$ for all $k=1,\ldots,m$ the optimal number of oracle calls, i.e. calculations of of the stochastic gradients of $f_{k}$ with $\sigma^{2}$-subgaussian variance is

The bad news is that (88) does not work with full parallelization trick and the best possible way to parallelize it is described in . However, standard accelerated scheme using mini-batched versions of the stochastic gradients without variance-reduction technique and incremental oracles which gives the bound

for the number of oracle calls and it admits full parallelization. It means that in the parallel optimization setup when we have computational network with $m$ nodes and the spanning tree for it with diameter $d$ the number of oracle calls per node is

However, for the decentralized setup the second row of Table 4 states that the number of communication rounds is the same as in (91) up to substitution of $d$ by $\sqrt{\chi(W)}$ and the number of oracle calls per node is

which has $m$ times bigger statistical term under the maximum than in (90). What is more, recently it was shown that there exists such a decentralized distributed method that requires

stochastic gradient oracle calls per node , but it is not optimal in terms of the number of communications. Moreover, there is a hypothesis that in the smooth case the bounds from Tables 3 and 4 (rows 2 and 3) are optimal in terms of the number of oracle calls per node for the class of methods that require optimal number of communication rounds up to polylogarithmic factors.

The same claim but for Table 5 was also presented in as a hypothesis and in this paper we propose the same hypothesis for the result stated Table 6 up to polylogarithmic and additionally we hypothesise that the noise level that we obtained is also unimprovable up to polylogarithmic factors.

As it was mentioned in Section 4, the recurrence technique that we use in Sections E and 5 can be very useful in the generalization of the results for STM from Section 4 for the case when instead of $\nabla f(x)$ only stochastic gradient $\nabla f(x,\xi)$ (see inequalities (10)-(11)) is available, $f$ is $L$-smooth and proximal step is computed in an inexact manner. It would be nice also to compare proposed methods for the case when $\delta$ with the results from . For the convex but non-strongly convex case one can also try to combine Nesterov’s smoothing technique with D-MASG from .

We emphasize that in our results we assume that each $f_{i}$ from (79) is $L$-smooth and $\mu$-strongly convex. When each $f_{i}$ is $L_{i}$-smooth and $\mu_{i}$-strongly convex, it means that in order to satisfy the assumption we use in our paper we need to choose $L=\max_{1\leq i\leq m}L_{i}$ and $\mu=\min_{1\leq i\leq m}\mu_{i}$. This choice can lead to a very slow rate in some situations, e.g. the worst-case $L$ can be $m$ times larger than $L$ for $f$ as for the case when $m=d$ and $f(x)=\nicefrac{{\|x\|_{2}^{2}}}{{2m}}=\nicefrac{{1}}{{m}}\sum_{i=1}^{m}f_{i}(x),$ $f_{i}(x)=\nicefrac{{x_{i}^{2}}}{{2}}$ where $L_{i}=1$ for all $i$ but $f$ is $\nicefrac{{1}}{{d}}$-smooth . It was shown that instead of worst-case $\mu$ and $L$ one can use $\bar{\mu}=\nicefrac{{1}}{{m}}\sum_{i=1}^{m}\mu_{i}$ and $\hat{L}$ to be some weighted average of $L_{i}$, but such techniques can spoil number of communication rounds needed to achieve desired accuracy.

It would be also interesting to generalize the proposed results for the case of more general stochastic gradients .

Application for Population Wasserstein Barycenter Calculation

In this section we consider the problem of calculation of population Wasserstein barycenter since this example hides different interesting details connected with the theory discussed in this paper. In our presentation of this example we rely mostly on the recent work .

Next, we consider the entropic OT problem (see )

For a given set of samples $q^{1},\ldots,q^{m}$ we introduce empirical barycenter as

We consider the problem (95) of finding population barycenter with some accuracy and discuss possible approaches to solve this problem in the following subsections.

However, before that, we need to mention some useful properties of ${\cal W}_{\mu}(p,q)$. First of all, one can write explicitly the dual function of $W_{\mu}(p,q)$ for a fixed $q\in S_{n}(1)$ (see ):

Using this representation one can deduce the following theorem.

We will also use another useful relation (see ):

where the gradient $\nabla{\cal W}_{\mu}(p,q)$ is taken w.r.t. the first argument.

2 SA Approach

Assume that one can obtain and use fresh samples $q^{1},q^{2},\ldots$ in online regime. This approach is called Stochastic Approximation (SA). It implies that at each iteration one can draw a fresh sample $q^{k}$ and compute the gradient w.r.t. $p$ of function ${\cal W}_{\mu}(p,q^{k})$ which is $\mu$-strongly convex and $M$-Lipschitz continuous with $M=\widetilde{O}(\sqrt{n}\|C\|_{\infty})$. Optimal methods for this case are based on iterations of the following form

for some $\delta\geq 0$ and for all $p,q\in S_{n}(1)$ after $N$ calls of this oracle produces such a point $p^{N}$ that with probability at least $1-\beta$ the following inequalities hold:

and, as a consequence of $\mu$-strong convexity of ${\cal W}_{\mu}(p,q)$ for all $q$,

with probability at least $1-\beta$, R-SGD requires

under additional assumption that $\delta=O(\mu\varepsilon^{2})$.

3 SAA Approach

Now let us assume that large enough collection of samples $q^{1},\ldots,q^{m}$ is available. Our goal is to find such $p\in S_{n}(1)$ that $\|\hat{p}-p_{\mu}^{*}\|_{2}\leq\varepsilon$ with high probability, i.e. $\varepsilon$-approximation of the population barycenter, via solving empirical barycenter problem (96). This approach is called Stochastic Average Approximation (SAA). Since ${\cal W}_{\mu}(p,q^{i})$ is $\mu$-strongly convex and $M$-Lipschitz in $p$ with $M=\widetilde{O}(\sqrt{n}\|C\|_{\infty})$ for all $i=1,\ldots,m$ we can conclude that with probability $\geq 1-\beta$

where we use that the diameter of $S_{n}(1)$ is $O(1)$. Moreover, in it was shown that one can guarantee that with probability $\geq 1-\beta$

Taking advantages of both inequalities we get that if

then with probability at least $1-\frac{\beta}{2}$

Assuming that we have such $\hat{p}\in S_{n}(1)$ that with probability at least $1-\frac{\beta}{2}$ the inequality

holds, we apply the union bound and get that with probability $\geq 1-\beta$

It remains to describe the approach that finds such $\hat{p}\in S_{n}(1)$ that satisfies (111) with probability at least $1-\beta$. Recall that in this subsection we consider the following problem

For each summand ${\cal W}_{\mu}(p,q^{i})$ in the sum above we have the explicit formula (98) for the dual function ${\cal W}_{q^{i},\mu}^{*}(\lambda)$. Note that one can compute the gradient of ${\cal W}_{q^{i},\mu}^{*}(\lambda)$ via $O(n^{2})$ arithmetical operations. What is more, ${\cal W}_{q^{i},\mu}^{*}(\lambda)$ has a finite-sum structure, so, one can sample $j$-th component of $q^{i}$ with probability $q_{j}^{i}$ and get stochastic gradient

which requires $O(n)$ arithmetical operations to be computed.

where $W=\overline{W}\otimes I_{n}$, this approach requires $\widetilde{O}\left(\sqrt{\frac{n\|C\|_{\infty}^{2}}{\mu\hat{\varepsilon}}\chi(W)}\right)$ communication rounds and $\widetilde{O}\left(n^{2.5}\sqrt{\frac{\|C\|_{\infty}^{2}}{\mu\hat{\varepsilon}}\chi(W)}\right)$ arithmetical operations per node to find gradients $\nabla{\cal W}_{q^{i},\mu}^{*}(\lambda)$. If instead of full gradients workers use stochastic gradients $\nabla{\cal W}_{q^{i},\mu}^{*}(\lambda,j)$ defined in (113) and these stochastic gradients have light-tailed distribution, i.e. satisfy the condition (86) with parameter $\sigma>0$, then to guarantee (114) with probability $\geq 1-\frac{\beta}{2}$ the aforementioned approach needs the same number of communications rounds and $\widetilde{O}\left(n\max\left\{\sqrt{\frac{n\|C\|_{\infty}^{2}}{\mu\hat{\varepsilon}}\chi(W)},\frac{m\sigma^{2}n\|C\|_{\infty}^{2}}{\hat{\varepsilon}^{2}}\chi(W)\right\}\right)$ arithmetical operations per node to find gradients $\nabla{\cal W}_{q^{i},\mu}^{*}(\lambda,j)$. Using $\mu$-strong convexity of ${\cal W}_{\mu}(p,q^{i})$ for all $i=1,\ldots,m$ and taking $\hat{\varepsilon}=\frac{\mu\varepsilon^{2}}{8}$ we get that our approach finds such a point $\hat{p}$ that satisfies (110) with probability at least $1-\frac{\beta}{2}$ using

arithmetical operations per node to find gradients in the deterministic case and

arithmetical operations per node to find stochastic gradients in the stochastic case. However, the state-of-the-art theory of learning states (see (108)) that $m$ should so large that in the stochastic case the second term in the bound for arithmetical operations typically dominates the first term and the dimensional dependence reduction from $n^{2.5}$ in the deterministic case to $n^{1.5}$ in the stochastic case is typically negligible in comparison with how much $\frac{m\sigma^{2}\sqrt{n}\|C\|_{\infty}^{2}}{\mu^{2}\varepsilon^{4}}\chi(W)$ is larger than $\frac{\|C\|_{\infty}}{\mu\varepsilon}\sqrt{\chi(W)}$. That is, our theory says that it is better to use full gradients in the particular example considered in this section (see also Section 7). Therefore, further in the section we will assume that $\sigma^{2}=0$, i.e. workers use full gradients of dual functions ${\cal W}_{q^{i},\mu}^{*}(\lambda)$.

However, bounds (115)-(116) were obtained under very restrictive at the first sight assumption that we have $m$ workers and each worker stores only one measure which is unrealistic. One can relax this assumption in the following way. Assume that we have $\hat{l}<m$ machines connected in a network with Laplacian matrix $\hat{W}$ and $j$-th machine stores $\hat{m}_{j}\geq 1$ measures for $j=1,\ldots,\hat{l}$ and $\sum_{j=1}^{\hat{l}}\hat{m}_{j}=m$. Next, for $j$-th machine we introduce $\hat{m}_{j}$ virtual workers also connected in some network that $j$-th machine can emulate along with communication between virtual workers and for every virtual worker we arrange one measure, e.g. it can be implemented as an array-like data structure with some formal rules for exchanging the data between cells that emulates communications. We also assume that inside the machine we can set the preferable network for the virtual nodes in such a way that each machine emulates communication between virtual nodes and computations inside them fast enough. Let us denote the Laplacian matrix of the obtained network of $m$ virtual nodes as $\overline{W}$. Then, our approach finds such a point $\hat{p}$ that satisfies (110) with probability at least $1-\frac{\beta}{2}$ using

time for arithmetical operations per machine to find gradients where $T_{\text{cm},j}$ is time needed for $j$-th machine to emulate communication between corresponding virtual nodes at each iteration and $T_{\text{cp},j}$ is time required by $j$-th machine to perform $1$ arithmetical operation for all corresponding virtual nodes in the gradients computation process at each iteration. For example, if we have only one machine and network of virtual nodes forms a complete graph than $\chi(W)=1$, but $T_{\text{cm},\max}$ and $T_{\text{cp},\max}$ can be large and to reduce the running time one should use more powerful machine. In contrast, if we have $m$ machines connected in a star-graph than $T_{\text{cm},\max}$ and $T_{\text{cp},\max}$ will be much smaller, but $\chi(W)$ will be of order $m$ which is large. Therefore, it is very important to choose balanced architecture of the network at least for virtual nodes per machine if it is possible. This question requires a separate thorough study and lies out of scope of this paper.

4 SA vs SAA comparison

Recall that in SA approach we assume that it is possible to sample new measures in online regime which means that the computational process is performed on one machine, whereas in SAA approach we assume that large enough collection of measures is distributed among the network of machines that form some computational network. In practice measures from $S_{n}(1)$ correspond to some images. As one can see from the complexity bounds, both SA and SAA approaches require large number of samples to learn the population barycenter defined in (95). If these samples are images, then they typically cannot be stored in RAM of one computer. Therefore, it is natural to use distributed systems to store the data.

Now let us compare complexity bounds for SA and SAA. We summarize them in Table 7.

When the communication is fast enough and $\mu$ is small we typically have that SAA approach significantly outperforms SA approach in terms of the complexity as well even for communication architectures with big $\chi(W)$. Therefore, for balanced architecture one can expect that SAA approach will outperform SA even more.

To conclude, we state that population barycenter computation is a natural example when it is typically much more preferable to use distributed algorithms with dual oracle instead of SA approach in terms of memory and complexity bounds.

Acknowledgments

We would like to thank F. Bach, P. Dvurechensky, M. Gürbüzbalaban, D. Kovalev, A. Nemirovski, A. Olshevsky, N. Srebro, A. Taylor and C. Uribe for useful discussions. The work of E. Gorbunov was supported by RFBR, project number 19-31-51001. The work of D. Dvinskikh was supported by Russian Science Foundation (project 18-71-10108). The work of A. Gasnikov was supported by RFBR, project number 19-31-51001.

Appendix A Basic Facts

In this section we enumerate for convenience basic facts that we use many times in our proofs.

Squared norm of the sum.

Appendix B Useful Facts about Duality

This section contains several useful results that we apply in our analysis.

Consider the function $f(x)$ defined on a closed convex set $Q\subseteq R^{n}$ and linear operator $A$ such that $\text{Ker}A\neq\{0\}$ and its dual function $\psi(y)$ defined as $\psi(y)=\max_{x\in Q}\left\{\langle y,Ax\rangle-f(x)\right\}$. Then

From Demyanov–Danskin theorem we have that $\nabla\psi(y)=Ax(A^{\top}y)$ which implies

Appendix C Auxiliary Results

In this section, we present the results from other papers that we rely on in our proofs.

where $\sigma_{k}^{2}$ belongs to the filtration $\sigma(\xi_{1},\ldots,\xi_{k-1})$ for all $k=1,\ldots,N$. Let $S_{N}=\sum\limits_{k=1}^{N}\xi_{k}$. Then there exists an absolute constant $C_{1}$ such that for any fixed $\beta>0$ and $B>b>0$ with probability at least $1-\beta$:

and let $S_{N}=\sum\limits_{k=1}^{N}\xi_{k}$. Assume that the sequence $\{\xi_{k}\}_{k=1}^{N}$ satisfy “light-tail” assumption:

where $\sigma_{1},\ldots,\sigma_{N}$ are some positive numbers. Then for all $\gamma\geq 0$

Appendix D Technical Results

For the sequence $\alpha_{k+1}\geq 0$ such that

Moreover, $A_{k}=\Omega\left(\frac{N^{2}}{L}\right)$.

We prove (125) by induction. For $k=0$ equation (124) gives us $\alpha_{1}=2L\alpha_{1}^{2}\Longleftrightarrow\alpha_{1}=\frac{1}{2L}$. Next we assume that (125) holds for all $k\leq l-1$ and prove it for $k=l$:

This quadratic inequality implies that $\alpha_{k+1}\leq\frac{1+\sqrt{4k^{2}+12k+1}}{4L}\leq\frac{1+\sqrt{(2k+3)^{2}}}{4L}\leq\frac{2k+4}{4L}=\frac{k+2}{2L}$.

Finally, the relation $A_{k}=\Omega\left(\frac{N^{2}}{L}\right)$ is proved in Lemma 1 from (see also ). ∎

For the sequence $\alpha_{k+1}\geq 0$ such that

If we solve quadratic equation $A_{k+1}(1+A_{k}\mu)=L\alpha_{k+1}^{2}$, $A_{k+1}=A_{k}+\alpha_{k+1}$ with respect to $\alpha_{k+1}$, we will get (127). Inequality (128) was established in Lemma 3 from and Lemma 4 from . It remains to prove (129). Since $\sqrt{a^{2}+b^{2}}\leq a+b$ for all $a,b\geq 0$ and $A_{k}\geq A_{0}=\frac{1}{L}$ we have

Let $A,B,D,r_{0},r_{1},\ldots,r_{N}$, where $N\geq 1$, be non-negative numbers such that

where $C$ is such positive number that $C^{2}\geq\max\{2A+2(B+D)C,1\}$, i.e. one can choose $C=\max\{B+D+\sqrt{(B+D)^{2}+2A},1\}$.

We prove (131) by induction. For $l=0$ the inequality $r_{l}\leq Cr_{0}$ trivially follows since $C\geq 1$. Next we assume that (131) holds for some $l<N$ and prove it for $l+1$:

Let $C,r_{0},r_{1},\ldots,r_{N}$, where $N\geq 1$, be non-negative numbers such that

and $C\in(0,\nicefrac{{1}}{{4}})$. Then for all $l=0,\ldots,N$ we have

We prove (133) by induction. For $l=0$ the inequality $r_{l}\leq 2r_{0}$ trivially follows. Next we assume that (133) holds for some $l\leq N-1$ and prove it for $l+1$. From (132), $C<\nicefrac{{1}}{{4}}$, $N\geq 1$ and $l\leq N-1$ we have

and $r_{0}\leq\frac{Cr_{0}}{\sqrt{A_{0}}}$ where $C$ is such positive number that

i.e. one can choose $C=\max\left\{\sqrt{A_{0}},\frac{3BD+\sqrt{9B^{2}D^{2}+4A}}{2}\right\}$.

since $C\geq\sqrt{A_{0}}$ and $C\geq\sqrt{A+BCD}\geq\sqrt{A+BD\sqrt{A_{0}}}$. Note that we also have $r_{0}\leq\frac{Cr_{0}}{\sqrt{A_{0}}}$. Next we assume that (135) holds for some $l\leq N-1$ and prove it for $l+1$:

That is, we proved the statement of the lemma for

w.r.t. $C$ one can show that the choice $C=\max\left\{\sqrt{A_{0}},\frac{3BD+\sqrt{9B^{2}D^{2}+4A}}{2}\right\}$ satisfies the assumption of the lemma on $C$.

Appendix E Similar Triangles Method with Inexact Proximal Step

In this section we focus on the composite optimization problem. i.e. problems of the type

where $f(x)$ is convex and $L$-smooth and $h(x)$ is convex and $L_{h}$-smooth. Before we present our method, let us introduce new notation.