Stochastic Optimization from Distributed, Streaming Data in Rate-limited Networks

TL;DR

This paper introduces distributed stochastic convex optimization algorithms designed for high-rate data streams over rate-limited networks, providing convergence guarantees and demonstrating their effectiveness through numerical experiments.

Contribution

It proposes two novel algorithms, D-SAMD and AD-SAMD, with convergence analysis and conditions for optimal performance in distributed streaming data settings.

Findings

Distributed algorithms achieve order-optimum convergence under certain network conditions.

Accelerated methods significantly improve convergence regimes compared to non-accelerated variants.

Algorithms perform effectively in numerical experiments simulating real-world network scenarios.

Abstract

Motivated by machine learning applications in networks of sensors, internet-of-things (IoT) devices, and autonomous agents, we propose techniques for distributed stochastic convex learning from high-rate data streams. The setup involves a network of nodes---each one of which has a stream of data arriving at a constant rate---that solve a stochastic convex optimization problem by collaborating with each other over rate-limited communication links. To this end, we present and analyze two algorithms---termed distributed stochastic approximation mirror descent (D-SAMD) and accelerated distributed stochastic approximation mirror descent (AD-SAMD)---that are based on two stochastic variants of mirror descent and in which nodes collaborate via approximate averaging of the local, noisy subgradients using distributed consensus. Our main contributions are (i) bounds on the convergence rates of D-SAMD and AD-SAMD in terms of the number of nodes, network topology, and ratio of the data streaming and communication rates, and (ii) sufficient conditions for order-optimum convergence of these algorithms. In particular, we show that for sufficiently well-connected networks, distributed learning schemes can obtain order-optimum convergence even if the communications rate is small. Further we find that the use of accelerated methods significantly enlarges the regime in which order-optimum convergence is achieved; this is in contrast to the centralized setting, where accelerated methods usually offer only a modest improvement. Finally, we demonstrate the effectiveness of the proposed algorithms using numerical experiments.

Figures

Equations

$$\underset{\mathbf{x}\in X}{\min }\psi (\mathbf{x})\triangleq \underset{\mathbf{x}\in X}{\min }{E}_{\xi }[\varphi (\mathbf{x},\xi )]$$

$$E[\psi (\mathbf{x}(T))-\psi ({\mathbf{x}}^{\ast })]\le O(1)\left[\frac{L}{T^2}+\frac{\mathcal{M}+\sigma }{\sqrt{T}}\right],$$

$$\psi (\mathbf{x})=f(\mathbf{x})+h(\mathbf{x}),$$

$${\mathbf{x}}^{\ast }\triangleq \arg \underset{\mathbf{x}\in X}{\min }\psi (\mathbf{x}),\text{and}{\psi }^{\ast }\triangleq \psi ({\mathbf{x}}^{\ast }).$$

$${\Vert \nabla f(\mathbf{x})-\nabla f(\mathbf{y})\Vert }_{\ast }\le L\Vert \mathbf{x}-\mathbf{x}\Vert ,\forall \mathbf{x},\mathbf{y}\in X,$$

$${\Vert \mathbf{g}\Vert }_{\ast }\triangleq \underset{\Vert \mathbf{x}\Vert \le 1}{\sup }⟨\mathbf{g},\mathbf{x}⟩.$$

$$\Vert h(\mathbf{x})-h(\mathbf{y})\Vert \le \mathcal{M}\Vert \mathbf{x}-\mathbf{y}\Vert ,\forall \mathbf{x},\mathbf{y}\in X.$$

$$\partial h(\mathbf{y})=\{\mathbf{g}:h(\mathbf{z})\ge h(\mathbf{y})+{\mathbf{g}}^T(\mathbf{z}-\mathbf{y}),\forall \mathbf{z}\in X\}.$$

$${\Vert \mathbf{g}\Vert }_{\ast }\le \mathcal{M},\forall \mathbf{g}\in \partial h(\mathbf{y}),\mathbf{y}\in X.$$

$${\Vert \partial \psi (\mathbf{x})-\partial \psi (\mathbf{y})\Vert }_{\ast }$$

$$\le {\Vert \nabla f(\mathbf{x})-\nabla f(\mathbf{y})\Vert }_{\ast }+{\Vert {\mathbf{g}}_{\mathbf{x}}-{\mathbf{g}}_{\mathbf{y}}\Vert }_{\ast }$$

$$\le L\Vert \mathbf{x}-\mathbf{y}\Vert +2\mathcal{M}.$$

$${\mathcal{N}}_i=\{j\in V:(i,j)\in E\},$$

$$E[G(\mathbf{x},\xi )]$$

$$E[{\Vert G(\mathbf{x},\xi )-\mathbf{g}(\mathbf{x})\Vert }_{\ast }^2]$$

$${\theta }_i(s)=\frac{1}{b}\sum _{t=(s-1)b+1}^{sb}G({\mathbf{x}}_i(s),{\xi }_i(t)).$$

$$E\left[{\Vert \frac{1}{k}\sum _{i=1}^k{\mathbf{z}}_i\Vert }_{\ast }^2\right]\le \frac{C_{\ast }{\sigma }^2}{k}.$$

$$r\le b\rho .$$

$$⟨\nabla \omega (\mathbf{x})-\nabla \omega (\mathbf{y}),\mathbf{x}-\mathbf{y}⟩\ge \alpha {\Vert \mathbf{x}-\mathbf{y}\Vert }^2,\forall \mathbf{x},\mathbf{y}\in X.$$

$$D_{\omega }\triangleq \sqrt{\underset{\mathbf{x}\in X}{\max }\omega (x)-\underset{\mathbf{x}\in X}{\min }\omega (x)},{\Omega }_{\omega }\triangleq \sqrt{\frac{2}{\alpha }{D}_{\omega }}.$$

$$V(\mathbf{x},\mathbf{z})=\omega (\mathbf{z})-(\omega (\mathbf{x})+⟨\nabla \omega (\mathbf{x}),\mathbf{z}-\mathbf{x}⟩).$$

$$P_{\mathbf{x}}(\mathbf{y})=\underset{\mathbf{z}\in X}{\mathrm{argmin}}⟨\mathbf{y},\mathbf{z}-\mathbf{x}⟩+V(\mathbf{x},\mathbf{z}).$$

$$\mathbf{x}(s+1)$$

$$\Vert P_{\mathbf{x}}(\mathbf{y})-P_{{\mathbf{x}}^{\prime }}({\mathbf{y}}^{\prime })\Vert \le \Vert \mathbf{x}-{\mathbf{x}}^{\prime }\Vert +\Vert \mathbf{y}-\mathbf{y}\prime \Vert .$$

$$\mathbf{H}(s)$$

$$\mathbf{G}(s)$$

$$\mathbf{Z}(s)$$

$$\overline{\mathbf{g}}(s)$$

$$\overline{\mathbf{z}}(s)$$

$$\mathbf{E}(s)$$

$$\tilde{\mathbf{Z}}(s)$$

$$\mathbf{H}(s)=\overline{\mathbf{G}}(s)+\mathbf{E}(s)+\overline{\mathbf{Z}}(s)+\tilde{\mathbf{Z}}(s).$$

$$E[{\Vert \overline{\mathbf{z}}(s)\Vert }_{\ast }^2]\le \frac{C_{\ast }{\sigma }^2}{mb}.$$

$${\Vert {\mathbf{e}}_i(s)\Vert }_{\ast }\le \underset{j,k}{\max }{m}^2\sqrt{C_{\ast }}{\lambda }_2^r{\Vert {\mathbf{g}}_j(s)-{\mathbf{g}}_k(s)\Vert }_{\ast }$$

Full text

I Introduction

Machine learning at its core involves solving stochastic optimization (SO) problems of the form

[TABLE]

to learn a “model” $\mathbf{x}\in X\subset\mathbb{R}^{n}$ that is then used for tasks such as dimensionality reduction, classification, clustering, regression, and/or prediction. A primary challenge of machine learning is to find a solution to the SO problem (1) without knowledge of the distribution $P(\xi)$. This involves finding an approximate solution to (1) using a sequence of $T$ training samples $\{\xi(t)\in\Upsilon\}_{t=1}^{T}$ drawn independently from the distribution $P(\xi)$, which is supported on a subset of $\Upsilon$. There are, in particular, two main categorizations of training data that, in turn, determine the types of methods that can be used to find approximate solutions to the SO problem. These are ($i$) batch training data and ($ii$) streaming training data.

In the case of batch training data, where all $T$ samples $\{\xi(t)\}$ are pre-stored and simultaneously available, a common strategy is sample average approximation (SAA) (also referred to as empirical risk minimization (ERM)), in which one minimizes the empirical average of the “risk” function $\phi(\cdot,\cdot)$ in lieu of the true expectation. In the case of streaming data, by contrast, the samples $\{\xi(t)\}$ arrive one-by-one, cannot be stored in memory for long, and should be processed as soon as possible. In this setting, stochastic approximation (SA) methods—the most well known of which is stochastic gradient descent (SGD)—are more common. Both SAA and SA have a long history in the literature; see [2] for a historical survey of SA methods, [3] for a comparative review of SAA and SA techniques, and [4] for a recent survey of SO techniques.

Among other trends, the rapid proliferation of sensing and wearable devices, the emergence of the internet-of-things (IoT), and the storage of data across geographically-distributed data centers have spurred a renewed interest in development and analysis of new methods for learning from fast-streaming and distributed data. The goal of this paper is to find a fast and efficient solution to the SO problem (1) in this setting of distributed, streaming data. In particular, we focus on geographically-distributed nodes that collaborate over rate-limited communication links (e.g., wireless links within an IoT infrastructure) and obtain independent streams of training data arriving at a constant rate.

The relationship between the rate at which communication takes place between nodes and the rate at which streaming data arrive at individual nodes plays a critical role in this setting. If, for example, data samples arrive much faster than nodes can communicate among themselves, it is difficult for the nodes to exchange enough information to enable an SA iteration on existing data in the network before new data arrives, thereby overwhelming the network. In order to address the challenge of distributed SO in the presence of a mismatch between the communications and streaming rates, we propose and analyze two distributed SA techniques, each based on distributed averaging consensus and stochastic mirror descent. In particular, we present bounds on the convergence rates of these techniques and derive conditions—involving the number of nodes, network topology, the streaming rate, and the communications rate—under which our solutions achieve order-optimum convergence speed.

I-A Relationship to Prior Work

SA methods date back to the seminal work of Robbins and Monro [5], and recent work shows that, for stochastic convex optimization, SA methods can outperform SAA methods [6, 7]. Lan [8] proposed accelerated stochastic mirror descent, which achieves the best possible convergence rate for general stochastic convex problems. This method, which makes use of noisy subgradients of $\psi(\cdot)$ computed using incoming training samples, satisfies

[TABLE]

where $\mathbf{x}^{*}$ denotes the minimizer of (1), $\sigma^{2}$ denotes variance of the subgradient noise, and $\mathcal{M}$ and $L$ denote the Lipschitz constants associated with the non-smooth (convex) component of $\psi$ and the gradient of the smooth (convex) component of $\psi$, respectively. Further assumptions such as smoothness and strong convexity of $\psi(\cdot)$ and/or presence of a structured regularizer term in $\psi(\cdot)$ can remove the dependence of the convergence rate on $\mathcal{M}$ and/or improve the convergence rate to $O(\sigma/T)$ [9, 6, 10, 11].

The problem of distributed SO goes back to the seminal work of Tsitsiklis et al. [12], which presents distributed first-order methods for SO and gives proofs of their asymptotic convergence. Myriad works since then have applied these ideas to other settings, each with different assumptions about the type of data, how the data are distributed across the network, and how distributed units process data and share information among themselves. In order to put our work in context, we review a representative sample of these works. A recent line of work was initiated by distributed gradient descent (DGD) [13], in which nodes descend using gradients of local data and collaborate via averaging consensus [14]. More recent works incorporate accelerated methods, time-varying or directed graphs, data structure, etc. [15, 16, 17, 18, 19]. These works tend not to address the SA problem directly; instead, they suppose a linearly separable function consistent with SAA using local, independent and identically distributed (i.i.d.) data. The works [20, 21, 22, 23] do consider SA directly, but suppose that nodes engage in a single round of message passing per stochastic subgradient sample.

We conclude by discussing two lines of works [24, 25, 26] that are most closely related to this work. In [24], nodes perform distributed SA by forming distributed mini-batch averages of stochastic gradients and using stochastic dual averaging. The main assumption in this work is that nodes can compute exact stochastic gradient averages (e.g., via AllReduce in parallel computing architectures). Under this assumption, it is shown in this work that there is an appropriate mini-batch size for which the nodes’ iterates converge at the optimum (centralized) rate. However, the need for exact averages in this work is not suited to rate-limited (e.g., wireless) networks, in which mimicking the AllReduce functionality can be costly and challenging.

The need for exact stochastic gradient averages in [24] has been relaxed recently in [26], in which nodes carry out distributed stochastic dual averaging by computing approximate mini-batch averages of dual variables via distributed consensus. In addition, and similar to our work, [26] allows for a mismatch between the communications rate and the data streaming rate. Nonetheless, there are four main distinctions between [26] and our work. First, we provide results for stochastic composite optimization, whereas [24, 26] suppose a differentiable objective. Second, we consider distributed mirror descent, which allows for a limited generalization to non-Euclidean settings. Third, we explicitly examine the impact of slow communications rate on performance, in particular highlighting the need for large mini-batches and their impact on convergence speed when the communications rate is slow. In [26], the optimum mini-batch size is first derived from [24], after which the communications rate needed to facilitate distributed consensus at the optimum mini-batch size is specified. While it appears to be possible to derive some of our results from a re-framing of the results of [26], it is crucial to highlight the trade-offs necessary under slow communications, which is not done in prior works. Finally, this work also presents a distributed accelerated mirror descent approach to distributed SA; a somewhat surprising outcome is that acceleration substantially improves the convergence rate in networks with slow communications.

I-B Our Contributions

In this paper, we present two strategies for distributed SA over networks with fast streaming data and slow communications links: distributed stochastic approximation mirror descent (D-SAMD) and accelerated distributed stochastic approximation mirror descent (AD-SAMD). In both cases, nodes first locally compute mini-batch stochastic subgradient averages to accommodate a fast streaming rate (or, equivalently, a slow communications rate), and then they collaboratively compute approximate network subgradient averages via distributed consensus. Finally, nodes individually employ mirror descent and accelerated mirror descent, respectively, on the approximate averaged subgradients for the next set of iterates.

Our main theoretical contribution is the derivation of upper bounds on the convergence rates of D-SAMD and AD-SAMD. These bounds involve a careful analysis of the impact of imperfect subgradient averaging on individual nodes’ iterates. In addition, we derive sufficient conditions for order-optimum convergence of D-SAMD and AD-SAMD in terms of the streaming and communications rates, the size and topology of the network, and the data statistics.

Two key findings of this paper are that distributed methods can achieve order-optimum convergence with small communication rates, as long as the number of nodes in the network does not grow too quickly as a function of the number of data samples each node processes, and that accelerated methods seem to offer order-optimum convergence in a larger regime than D-SAMD, thus potentially accommodating slower communications links relative to the streaming rate. By contrast, the convergence speeds of centralized stochastic mirror descent and accelerated stochastic mirror descent typically differ only in higher-order terms. We hasten to point out that we do not claim superior performance of D-SAMD and AD-SAMD versus other distributed methods. Instead, the larger goal is to establish the existence of methods for order-optimum stochastic learning in the fast-streaming, rate-limited regimes. D-SAMD and AD-SAMD should be best regarded as a proof of concept towards this end.

I-C Notation and Organization

We typically use boldfaced lowercase and boldfaced capital letters (e.g., $\mathbf{x}$ and $\mathbf{W}$) to denote (possibly random) vectors and matrices, respectively. Unless otherwise specified, all vectors are assumed to be column vectors. We use $(\cdot)^{T}$ to denote the transpose operation and $\mathbf{1}$ to denote the vector of all ones. Further, we denote the expectation operation by $E[\cdot]$ and the field of real numbers by $\mathbb{R}$. We use $\nabla$ to denote the gradient operator, while $\odot$ denotes the Hadamard product. Finally, given two functions $p(r)$ and $q(r)$, we write $p(r)=O(q(r))$ if there exists a constant $C$ such that $\forall r,p(r)\leq Cq(r)$, and we write $p(r)=\Omega(q(r))$ if $q(r)=O(p(r))$.

The rest of this paper is organized as follows. In Section II, we formalize the problem of distributed stochastic composite optimization. In Sections III and IV, we describe D-SAMD and AD-SAMD, respectively, and also derive performance guarantees for these two methods. We examine the empirical performance of the proposed methods via numerical experiments in Section V, and we conclude the paper in Section VI. Proofs are provided in the appendix.

II Problem Formulation

The objective of this paper is order-optimal, distributed minimization of the composite function

[TABLE]

where $\mathbf{x}\in X\subset\mathbb{R}^{n}$ and $X$ is convex and compact. The space $\mathbb{R}^{n}$ is endowed with an inner product $\langle\cdot,\cdot\rangle$ that need not be the usual one and a norm $\left\|{\cdot}\right\|$ that need not be the one induced by the inner product. In the following, the minimizer and the minimum value of $\psi$ are denoted as:

[TABLE]

We now make a few assumptions on the smooth ($f(\cdot)$) and non-smooth ($h(\cdot)$) components of $\psi$. The function $f:X\to\mathbb{R}$ is convex with Lipschitz continuous gradients, i.e.,

[TABLE]

where $\left\|{\cdot}\right\|_{*}$ is the dual norm associated with $\langle\cdot,\cdot\rangle$ and $\left\|{\cdot}\right\|$:

[TABLE]

The function $h:X\to\mathbb{R}$ is convex and Lipschitz continuous:

[TABLE]

Note that $h$ need not have gradients; however, since it is convex we can consider its subdifferential, denoted by $\partial h(\mathbf{y})$:

[TABLE]

An important fact that will be used in this paper is that the subgradient $\mathbf{g}\in\partial h$ of a Lipschitz-continuous convex function $h$ is bounded [27, Lemma 2.6]:

[TABLE]

Consequently, the gap between the subgradients of $\psi$ is bounded: $\forall\mathbf{x},\mathbf{y}\in X$ and $\mathbf{g}_{\mathbf{x}}\in\partial h(\mathbf{x})$, $\mathbf{g}_{\mathbf{y}}\in\partial h(\mathbf{y})$, we have

[TABLE]

II-A Distributed Stochastic Composite Optimization

Our focus in this paper is minimization of $\psi(\mathbf{x})$ over a network of $m$ nodes, represented by the undirected graph $G=(V,E)$. To this end, we suppose that nodes minimize $\psi$ collaboratively by exchanging subgradient information with their neighbors at each communications round. Specifically, each node $i\in V$ transmits a message at each communications round to each of its neighbors, defined as

[TABLE]

where we suppose that a node is in its own neighborhood, i.e., $i\in\mathcal{N}_{i}$. We assume that this message passing between nodes takes place without any error or distortion. Further, we constrain the messages between nodes to be members of the dual space of $X$ and to satisfy causality; i.e., messages transmitted by a node can depend only on its local data and previous messages received from its neighbors.

Next, in terms of data generation, we suppose that each node $i\in V$ queries a first-order stochastic “oracle” at a fixed rate—which may be different from the rate of message exchange—to obtain noisy estimates of the subgradient of $\psi$ at different query points in $X$. Formally, we use ‘$t$’ to index time according to data-acquisition rounds and define $\{\xi_{i}(t)\in\Upsilon\}_{t\geq 1}$ to be a sequence of independent (with respect to $i$ and $t$) and identically distributed (i.i.d.) random variables with unknown probability distribution $P(\xi)$. At each data-acquisition round $t$, node $i$ queries the oracle at search point $\mathbf{x}_{i}(s)$ to obtain a point $G(\mathbf{x}_{i}(s),\xi_{i}(t))$ that is a noisy version of the subgradient of $\psi$ at $\mathbf{x}_{i}(s)$. Here, we use ‘$s$’ to index time according to search-point update rounds, with possibly multiple data-acquisition rounds per search-point update. The reason for allowing the search-point update index $s$ to be different from the data-acquisition index $t$ is to accommodate the setting in which data (equivalently, subgradient estimates) arrive at a much faster rate than the rate at which nodes can communicate with each other; we will elaborate further on this in the next subsection.

Formally, $G(\mathbf{x},\xi)$ is a Borel function that satisfies the following properties:

[TABLE]

where the expectation is with respect to the distribution $P(\xi)$. We emphasize that this formulation is equivalent to that in which the objective function is $\psi(\mathbf{x})\triangleq E[\phi(\mathbf{x},\xi)]$, and where nodes in the network acquire data point $\{\xi_{i}(t)\}_{i\in V}$ at each data-acquisition round $t$ that are then used to compute the subgradients of $\phi(\mathbf{x},\xi_{i}(t))$, which—in turn—are noisy subgradients of $\psi(\mathbf{x})$.

II-B Mini-batching for Rate-Limited Networks

A common technique to reduce the variance of the (sub)gradient noise and/or reduce the computational burden in centralized SO is to average “batches” of oracle outputs into a single (sub)gradient estimate. This technique, which is referred to as mini-batching, is also used in this paper; however, its purpose in our distributed setting is to both reduce the subgradient noise variance and manage the potential mismatch between the communications rate and the data streaming rate. Before delving into the details of our mini-batch strategy, we present a simple model to parametrize the mismatch between the two rates. Specifically, let $\rho>0$ be the communications ratio, i.e. the fixed ratio between the rate of communications and the rate of data acquisition. That is, $\rho\geq 1$ implies nodes engage in $\rho$ rounds of message exchanges for every data-acquisition round. Similarly, $\rho<1$ means there is one communications round for every $1/\rho$ data-acquisition rounds. We ignore rounding issues for simplicity.

The mini-batching in our distributed problem proceeds as follows. Each mini-batch round spans $b\geq 1$ data-acquisition rounds and coincides with the search-point update round, i.e., each node $i$ updates its search point at the end of a mini-batch round. In each mini-batch round $s$, each node $i$ uses its current search point $\mathbf{x}_{i}(s)$ to compute an average of oracle outputs

[TABLE]

This is followed by each node computing a new search point $\mathbf{x}_{i}(s+1)$ using $\theta_{i}(s)$ and messages received from its neighbors.

In order to analyze the mini-batching distributed SA techniques proposed in this work, we need to generalize the usual averaging property of variances to non-Euclidean norms.

Lemma 1

Let $\mathbf{z}_{1},\dots,\mathbf{z}_{k}$ be i.i.d. random vectors in $\mathbb{R}^{n}$ with $E[\mathbf{z}_{i}]=0$ and $E[\left\|{\mathbf{z}_{i}}\right\|^{2}_{*}]\leq\sigma^{2}$. There exists a constant $C_{*}\geq 0$, which depends only on $\left\|{\cdot}\right\|$ and $\langle\cdot,\cdot\rangle$, such that

[TABLE]

Proof:

This follows directly from the property of norm equivalence in finite-dimensional spaces. ∎

In order to illustrate Lemma 1, notice that when $\left\|{\cdot}\right\|=\left\|{\cdot}\right\|_{1}$, i.e., the $\ell_{1}$ norm, and $\langle\cdot,\cdot\rangle$ is the standard inner product, the associated dual norm is the $\ell_{\infty}$ norm: $\left\|{\cdot}\right\|_{*}=\left\|{\cdot}\right\|_{\infty}$. Since $\left\|{\mathbf{x}}\right\|^{2}_{\infty}\leq\left\|{\mathbf{x}}\right\|^{2}_{2}\leq n\left\|{\mathbf{x}}\right\|_{\infty}^{2}$, we have $C_{*}=n$ in this case. Thus, depending on the norm in use, the extent to which averaging reduces subgradient noise variance may depend on the dimension of the optimization space.

In the following, we will use the notation $\mathbf{z}_{i}(s)\triangleq\theta_{i}(s)-\mathbf{g}(\mathbf{x}_{i}(s))$. Then, $E[\left\|{\mathbf{z}_{i}(s)}\right\|_{*}^{2}]\leq C_{*}\sigma^{2}/b$. We emphasize that the subgradient noise vectors $\mathbf{z}_{i}(s)$ depend on the search points $\mathbf{x}_{i}(s)$; we suppress this notation for brevity.

II-C Problem Statement

It is straightforward to see that mini-batching induces a performance trade-off: Averaging reduces subgradient noise and processing time, but it also reduces the rate of search-point updates (and hence slows down convergence). This trade-off depends on the relationship between the streaming and communications rates. In order to carry out distributed SA in an order-optimal manner, we will require that the nodes collaborate by carrying out $r\geq 1$ rounds of averaging consensus on their mini-batch averages $\theta_{i}(s)$ in each mini-batch round $s$ (see Section III for details). In order to complete the $r$ communication rounds in time for the next mini-batch round, we have the constraint

[TABLE]

If communications is faster, or if the mini-batch rounds are longer, nodes can fit in more rounds of information exchange between each mini-batch round or, equivalently, between each search-point update. But when the mismatch factor $\rho$ is small, the mini-batch size $b$ needed to enable sufficiently many consensus rounds may be so large that the reduction in subgradient noise is outstripped by the reduction in search-point updates and the resulting convergence speed is sub-optimum. In this context, our main goal is specification of sufficient conditions for $\rho$ such that the resulting convergence speeds of the proposed distributed SA techniques are optimum.

III Distributed Stochastic Approximation

Mirror Descent

In this section we present our first distributed SA algorithm, called distributed stochastic approximation mirror descent (D-SAMD). This algorithm is based upon stochastic approximated mirror descent, which is a generalized version of stochastic subgradient descent. Before presenting D-SAMD, we review a few concepts that underlie mirror descent.

III-A Stochastic Mirror Descent Preliminaries

Stochastic mirror descent, presented in [8], is a generalization of stochastic subgradient descent. This generalization is characterized by a distance-generating function $\omega:X\to\mathbb{R}$ that generalizes the Euclidean norm. The distance-generating function must be continuously differentiable and strongly convex with modulus $\alpha$, i.e.

[TABLE]

In the convergence analysis, we will require two measures of the “radius” of $X$ that will arise in the convergence analysis, defined as follows:

[TABLE]

The distance-generating function induces the prox function, or the Bregman divergence $V:X\times X\to\mathbb{R}_{+}$, which generalizes the Euclidean distance:

[TABLE]

The prox function $V(\mathbf{x},\cdot)$ inherits strong convexity from $\omega(\cdot)$, but it need not be symmetric or satisfy the triangle inequality. We define the prox mapping $P_{\mathbf{x}}:\mathbb{R}^{n}\to X$ as

[TABLE]

The prox mapping generalizes the usual subgradient descent step, in which one minimizes the local linearization of the objective function regularized by the Euclidean distance of the step taken. In (centralized) stochastic mirror descent, one computes iterates of the form

[TABLE]

where $\mathbf{g}(s)$ is a stochastic subgradient of $\psi(\mathbf{x}(s))$, and $\gamma_{s}$ is a step size. These iterates have the same form as (stochastic) subgradient descent; indeed, choosing $\omega(\mathbf{x})=\frac{1}{2}\left\|{\mathbf{x}}\right\|^{2}_{2}$ as well as $\langle\cdot,\cdot\rangle$ and $\left\|{\cdot}\right\|$ to be the usual ones results in subgradient descent iterations.

One can speed up convergence by choosing $\omega(\cdot)$ to match the structure of $X$ and $\psi$. For example, if the optimization space $X$ is the unit simplex over $\mathbb{R}^{n}$, one can choose $\omega(\mathbf{x})=\sum_{i}x_{i}\log(x_{i})$ and $\left\|{\cdot}\right\|$ to be the $\ell_{1}$ norm. This leads to $V(\mathbf{x},\mathbf{z})=D(\mathbf{z}||\mathbf{x}))$, where $D(\cdot||\cdot)$ denotes the Kullback-Leibler (K-L) divergence between $\mathbf{z}$ and $\mathbf{x}$. This choice speeds up convergence on the order of $O\left(\sqrt{n/\log(n)}\right)$ over using the Euclidean norm throughout. Along a similar vein, when $\psi$ includes an $\ell_{1}$ regularizer to promote a sparse minimizer, one can speed up convergence by choosing $\omega(\cdot)$ to be a $p$-norm with $p=\log(n)/(\log(n)-1)$.

In order to guarantee convergence for D-SAMD, we need to restrict further the distance-generating function $\omega(\cdot)$. In particular, we require that the resulting prox mapping be 1-Lipschitz continuous in $\mathbf{x},\mathbf{y}$ pairs, i.e., $\forall\ \mathbf{x},\mathbf{x}^{\prime},\mathbf{y},\mathbf{y}^{\prime}\in\mathbb{R}^{n}$,

[TABLE]

This condition is in addition to the conditions one usually places on the Bregman divergence for stochastic optimization; we will use it to guarantee that imprecise gradient averages make a bounded perturbation in the iterates of stochastic mirror descent. The condition holds whenever the prox mapping is the projection of a 1-Lipschitz function of $\mathbf{x},\mathbf{y}$ onto $X$. For example, it is easy to verify that this condition holds in the Euclidean setting. One can also show that when the distance-generating function $\omega(\mathbf{x})$ is an $\ell_{p}$ norm for $p>1$, the resulting prox mapping is 1-Lipschitz continuous in $\mathbf{x}$ and $\mathbf{y}$ as required.

However, not all Bregman divergences satisfy this condition. One can show that the K-L divergence results in a prox mapping that is not Lipschitz. Consequently, while we present our results in terms of general prox functions, we emphasize that the results do not apply in all cases.111We note further that it is possible to relax the constraint that the best Lipschitz constant be no larger than unity. This worsens the scaling laws—in particular, the required communications ratio $\rho$ grows in $T$ rather than decreases—and we omit this case for brevity’s sake. One can think of the results primarily in the setting of Euclidean (accelerated) stochastic subgradient descent—for which case they are guaranteed to hold—with the understanding that one can check on a case-by-case basis to see if they hold for a particular non-Euclidean setting.

III-B Description of D-SAMD

Here we present in detail D-SAMD, which generalizes stochastic mirror descent to the setting of distributed, streaming data. In D-SAMD, nodes carry out iterations similar to stochastic mirror descent as presented in [8], but instead of using local stochastic subgradients associated with the local search points, they carry out approximate consensus to estimate the average of stochastic subgradients across the network. This reduces the subgradient noise at each node and speeds up convergence.

Let $\mathbf{W}$ be a symmetric, doubly-stochastic matrix consistent with the network graph $G$, i.e., $[\mathbf{W}]_{ij}\triangleq w_{ij}=0$ if $(i,j)\notin E$. Further suppose that $\mathbf{W}-\mathbf{1}\mathbf{1}^{T}/n$ has spectral radius strictly less than one, i.e. the second-largest eigenvalue magnitude is strictly less than one. This condition is guaranteed by choosing the diagonal elements of $\mathbf{W}$ to be strictly greater than zero.

Next, we focus on the case of constant step size $\gamma$ and set it as $0<\gamma\leq\alpha/(2L)$.222It is shown in [8] that a constant step size is sufficient for order-optimal performance, so we adopt such a rule here. For simplicity, we suppose that there is a predetermined number of data-acquisition rounds $T$, which leads to $S=T/b$ mini-batch rounds. We detail the steps of D-SAMD in Algorithm 1. Further, in Figure 1 we illustrate the data acquisition round, mini-batch round, communication round, and search point update counters and their role in the D-SAMD algorithm.

In D-SAMD, each node $i$ initializes its iterate at the minimizer of $\omega(\cdot)$, which is guaranteed to be unique due to strong convexity. At each mini-batch round $s$, each node $i$ obtains its mini-batched subgradient and nodes engage in $r$ rounds of averaging consensus to produce the (approximate) average subgradients $\mathbf{h}_{i}^{r}(s)$. Then, each node $i$ takes a mirror prox step, using $\mathbf{h}^{r}_{i}(s)$ instead of its own mini-batched estimate. Finally, each node keeps a running average of its iterates, which is well-known to speed up convergence [28].

III-C Convergence Analysis

The convergence rate of D-SAMD depends on the bias and variance of the approximate subgradient averages $\mathbf{h}_{i}^{r}(s)$. In principle, averaging subgradients together reduces the noise variance and speeds up convergence. However, because averaging consensus using only $r$ communications rounds results in approximate averages, each node takes a slightly different mirror prox step and therefore ends up with a different iterate. At each mini-batch round $s$, nodes then compute subgradients at different search points, leading to bias in the averages $\mathbf{h}^{r}_{i}(s)$. This bias accumulates at a rate that depends on the subgradient noise variance, the topology of the network, and the number of consensus rounds per mini-batch round.

Therefore, the first step in bounding the convergence speed of D-SAMD is to bound the bias and the variance of the subgradient estimates $\mathbf{h}^{r}_{i}(s)$, which we do in the following lemma.

Lemma 2

Let $0\leq\lambda_{2}<1$ denote the magnitude of the second-largest (ordered by magnitude) eigenvalue of $\mathbf{W}$. Define the matrices

[TABLE]

recalling that the subgradient noise $\mathbf{z}_{i}(s)$ is defined with respect to the mini-batched subgradient $\theta_{i}(s)$. Also define

[TABLE]

where the matrices $\overline{\mathbf{G}}(s),\overline{\mathbf{Z}}(s)\in\mathbb{R}^{n\times n}$ have identical columns. Finally, define the matrices

[TABLE]

of average consensus error on the subgradients and subgradient noise, respectively. Then, the following facts are true. First, one can write $\mathbf{H}(s)$ as

[TABLE]

Second, the columns of $\overline{\mathbf{Z}}(s)$ satisfy

[TABLE]

Finally, the $i$th columns of $\mathbf{E}(s)$ and $\tilde{\mathbf{Z}}(s)$, denoted by $\mathbf{e}_{i}(s)$ and $\tilde{\mathbf{z}}_{i}(s)$, respectively, satisfy

[TABLE]

and

[TABLE]

where we have used $\mathbf{g}_{j}(s)$ as a shorthand for $\mathbf{g}(\mathbf{x}_{j}(s))$.

The next step in the convergence analysis is to bound the distance between iterates at different nodes. As long as iterates are not too far apart, the subgradients computed at different nodes have sufficiently similar means that averaging them together reduces the overall subgradient noise.

Lemma 3

Let $a_{s}\triangleq\max_{i,j}\left\|{\mathbf{x}_{i}(s)-\mathbf{x}_{j}(s)}\right\|$. The moments of $a_{s}$ follow:

[TABLE]

Now, we bound D-SAMD’s expected gap to optimality.

Theorem 1

For D-SAMD, the expected gap to optimality at each node $i$ satisfies

[TABLE]

where

[TABLE]

and

[TABLE]

quantify the moments of the effective subgradient noise.

The convergence rate proven in Theorem 1 is akin to that provided in [8], with $\Delta_{s}^{2}$ taking the role of the subgradient noise variance. A crucial difference is the presence of the final term involving $\Xi_{s}$. In [8], this term vanishes because the intrinsic subgradient noise has zero mean. However, the equivalent gradient error in D-SAMD does not have zero mean in general. As nodes’ iterates diverge, their subgradients differ, and the nonlinear mapping between iterates and subgradients results in noise with nonzero mean.

The critical question is how fast communication needs to be for order-optimum convergence speed, i.e., the convergence speed that one would obtain if nodes had access to other nodes’ subgradient estimates at each round. After $S$ mini-batch rounds, the network has processed $mT$ data samples. Centralized mirror descent, with access to all $mT$ data samples in sequence, achieves the convergence rate [8]

[TABLE]

The final term dominates the error as a function of $m$ and $T$ if $\sigma^{2}>0$. In the following corollary we derive conditions under which the convergence rate of D-SAMD matches this term.

Corollary 1

The optimality gap for D-SAMD satisfies

[TABLE]

provided the mini-batch size $b$, the communications ratio $\rho$, the number of users $m$, and the Lipschitz constant $\mathcal{M}$ satisfy

[TABLE]

III-D Discussion

Corollary 1 gives new insights into influences of the communications and streaming rates, network topology, and mini-batch size on the convergence rate of distributed stochastic learning. In [24], a mini-batch size of $b=O(T^{1/2})$ is prescribed—which is sufficient whenever gradient averages are perfect—and in [26] the number of imperfect consensus rounds needed to facilitate the mini-batch size $b$ prescribed in [24] is derived. By contrast, we derive a mini-batch condition sufficient to drive the effective noise variance to $O(\sigma^{2}/(mT))$ while taking into consideration the impact of imperfect subgradient averaging. This condition depends not only on $T$ but also on $m$, $\rho$, $\lambda_{2}$, and $\sigma^{2}$—indeed, for all else constant, the optimum mini-batch size is merely $\Omega(\log(T))$. Then, the condition on $\rho$ essentially ensures that $b=O(T^{1/2})$ as specified in [24].

We note that Corollary 1 imposes a strict requirement on $\mathcal{M}$, the Lipschitz constant of the non-smooth part of $\psi$. Essentially the non-smooth part must vanish as $m$, $T$, or $\sigma^{2}$ becomes large. This is because the contribution of $h(\mathbf{x})$ to the convergence rate depends only on the number of iterations taken, not on the noise variance. Reducing the effective subgradient noise via mini-batching has no impact on this contribution, so we require the Lipschitz constant $\mathcal{M}$ to be small to compensate.

Finally, we note that Corollary 1 dictates the relationship between the size of the network and the number of data samples obtained at each node. Leaving the terms besides $m$ and $T$ constant, Corollary 1 requires $T=\Omega(m)$, i.e. the number of nodes in the network should scale no faster than the number of data samples processed per node. This is a relatively mild condition for big data applications; many applications involve data streams that are large relative to the size of the network. Furthermore, ignoring the $\log(mT)$ term and assuming $\lambda_{2}$ and $\sigma$ to be fixed, Corollary 1 indicates that a communication ratio of \rho=\Omega\big{(}\sqrt{m/T}\big{)} is sufficient for order optimality; i.e., nodes need to communicate at least \Omega\big{(}\sqrt{m/T}\big{)} times per data sample. This means that if $T$ scales faster than $\Omega(m)$ then the required communications ratio approaches zero in this case as $m,T\to\infty$. In particular, fast stochastic learning is possible in expander graphs, for which the spectral gap $1-\lambda_{2}$ is bounded away from zero, even in communication rate-limited scenarios. For graph families that are poor expanders, however, the required communications ratio depends on the scaling of $\lambda_{2}$ as a function of $m$.

IV Accelerated Distributed Stochastic Approximation Mirror Descent

In this section, we present accelerated distributed stochastic approximation mirror descent (AD-SAMD), which distributes the accelerated stochastic approximation mirror descent proposed in [8]. The centralized version of accelerated mirror descent achieves the optimum convergence rate of

[TABLE]

Consequently, we will see that the convergence rate of AD-SAMD has $1/S^{2}$ as its first term. This faster convergence in $S$ allows for more aggressive mini-batching, and the resulting conditions for order-optimal convergence are less stringent.

IV-A Description of AD-SAMD

The setting for AD-SAMD is the same as in Section III. We again suppose a distance function $\omega:X\to\mathbb{R}$, its associated prox function/Bregman divergence $V:X\times X\to\mathbb{R}$, and the resulting (Lipschitz) prox mapping $P_{x}:\mathbb{R}^{n}\to X$.

As in Section III-B, we suppose a mixing matrix $\mathbf{W}\in\mathbb{R}^{m\times m}$ that is symmetric, doubly stochastic, consistent with $G$, and has nonzero spectral gap. The main distinction between accelerated and standard mirror descent is the way one averages iterates. Rather than simply average the sequence of iterates, one maintains several distinct sequences of iterates, carefully averaging them along the way. This involves two sequences of step sizes $\beta_{s}\in[1,\infty)$ and $\gamma_{s}\in\mathbb{R}$, which are not held constant. Again we suppose that the number of mini-batch rounds $S=T/b$ is predetermined. We detail the steps of AD-SAMD in Algorithm 2.

The sequences of iterates $\mathbf{x}_{i}(s)$, $\mathbf{x}_{i}^{\mathrm{md}}(s)$, and $\mathbf{x}^{\mathrm{ag}}_{i}(s)$ are interrelated in complicated ways; we refer the reader to [8] for an intuitive explanation of these iterations.

IV-B Convergence Analysis

As with D-SAMD, the convergence analysis relies on bounds on the bias and variance of the averaged subgradients. To this end, we note first that Lemma 2 also holds for AD-SAMD, where $\mathbf{H}(s)$ has columns corresponding to noisy subgradients evaluated at $\mathbf{x}_{i}^{\mathrm{md}}(s)$. Next, we bound the distance between iterates at different nodes. This analysis is somewhat more complicated due to the relationships between the three iterate sequences.

Lemma 4

Let

[TABLE]

Then, the moments of $a_{s}$, $b_{s}$, and $c_{s}$ satisfy:

[TABLE]

Now, we bound the expected gap to optimality of the AD-SAMD iterates.

Theorem 2

For AD-SAMD, there exist step size sequences $\beta_{s}$ and $\gamma_{s}$ such that the expected gap to optimality satisfies

[TABLE]

where

[TABLE]

and

[TABLE]

As with D-SAMD, we study the conditions under which AD-SAMD achieves order-optimum convergence speed. The centralized version of accelerated mirror descent, after processing the $mT$ data samples that the network sees after $S$ mini-batch rounds, achieves the convergence rate

[TABLE]

This is the optimum convergence rate under any circumstance. In the following corollary, we derive the conditions under which the convergence rate matches the second term, which usually dominates the error when $\sigma^{2}>0$.

Corollary 2

The optimality gap satisfies

[TABLE]

provided

[TABLE]

IV-C Discussion

The crucial difference between the two schemes is that AD-SAMD has a convergence rate of $1/S^{2}$ in the absence of noise and non-smoothness. This faster term, which is often negligible in centralized mirror descent, means that AD-SAMD tolerates more aggressive mini-batching without impact on the order of the convergence rate. As a result, while the condition on the mini-batch size $b$ is the same in terms of $\rho$, the condition on $\rho$ is relaxed by $1/4$ in the exponents of $m$ and $T$. This is because the condition $b=O(T^{1/2})$, which holds for standard stochastic SO methods, is relaxed to $b=O(T^{3/4})$ for accelerated stochastic mirror descent.

Similar to Corollary 1, Corollary 2 prescribes a relationship between $m$ and $T$, but the relationship for AD-SAMD is $T~{}=~{}\Omega(m^{1/3})$, holding all but $m,T$ constant. This again is due to the relaxed mini-batch condition $b=O(T^{3/4})$ for accelerated stochastic mirror descent. Furthermore, ignoring the $\log$ term, Corollary 2 indicates that a communications ratio $\rho=\Omega\left(\frac{m^{1/4}}{T^{3/4}}\right)$ is needed for well-connected graphs such as expander graphs. In this case, as long as $T$ grows faster than the cube root of $m$, order-optimum convergence rates can be obtained even for small communications ratio. Thus, the use of accelerated methods increases the domain in which order optimum rate-limited learning is guaranteed.

V Numerical Example: Logistic Regression

To demonstrate the scaling laws predicted by Corollaries 1 and 2 and to investigate the empirical performance of D-SAMD and AD-SAMD, we consider supervised learning via binary logistic regression. Specifically, we assume each node observes a stream of pairs $\xi_{i}(t)=(y(t),l(t))$ of data points $y_{i}(t)\in\mathbb{R}^{d}$ and their labels $l_{i}(t)\in\{0,1\}$, from which it learns a classifier with the log-likelihood function

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{d}$ and $x_{0}\in\mathbb{R}$ are regression coefficients.

The SO task is to learn the optimum regression coefficients $\mathbf{x},x_{0}$. In terms of the framework of this paper, $\Upsilon=(\mathbb{R}^{d}\times\{0,1\})$, and $X=\mathbb{R}^{d+1}$ (i.e., $n=d+1$). We use the Euclidean norm, inner product, and distance-generating function to compute the prox mapping. The convex objective function is the negative of the log-likelihood function, averaged over the unknown distribution of the data, i.e.,

[TABLE]

Minimizing $\psi$ is equivalent to performing maximum likelihood estimation of the regression coefficients [29].

We examine performance on synthetic data so that there exists a “ground truth” distribution with which to compute $\psi(\mathbf{x})$ and evaluate empirical performance. We suppose that the data follow a Gaussian distribution. For $l_{i}(t)~{}\in~{}\{0,1\}$, we let $\mathbf{y}_{i}(t)~{}\sim~{}\mathcal{N}(\mu_{l_{i}(t)},\sigma_{r}^{2}\mathbf{I})$, where $\mu_{l(t)}$ is one of two mean vectors, and $\sigma_{r}^{2}>0$ is the noise variance.333The variance $\sigma_{r}^{2}$ is distinct from the resulting gradient noise variance $\sigma^{2}$. For this experiment, we draw the elements $\mu_{0}$ and $\mu_{1}$ randomly from the standard normal distribution, let $d=20$, and choose $\sigma_{r}^{2}=2$. We consider several network topologies, as detailed in the next subsections.

We compare the performance of D-SAMD and AD-SAMD against several other schemes. As a best-case scenario, we consider centralized mirror descent, meaning that at each data-acquisition round $t$ all $m$ data samples and their associated gradients are available at a single machine, which carries out stochastic mirror descent and accelerated stochastic mirror descent. These algorithms naturally have the best average performance. As a baseline, we consider local (accelerated) stochastic mirror descent, in which nodes simply perform mirror descent on their own data streams without collaboration. This scheme does benefit from an insensitivity to the communications ratio $\rho$, and no mini-batching is required, but it represents a minimum standard for performance in the sense that it does not require collaboration among nodes.

Finally, we consider a communications-constrained adaptation of distributed gradient descent (DGD), introduced in [13], where local subgradient updates are followed by a single round of consensus averaging on the search points $\mathbf{x}_{i}(s)$. DGD implicitly supposes that $\rho=1$. To handle the $\rho<1$ case, we consider two adaptations: naive DGD, in which data samples that arrive between communications rounds are simply discarded, and mini-batched DGD, in which nodes compute local mini-batches of size $b=1/\rho$, take gradient updates with the local mini-batch, and carry out a consensus round. While it is not designed for the communications rate-limited scenario, DGD has good performance in general, so it represents a natural alternative against which to compare the performance of D-SAMD and AD-SAMD.

V-A Fully Connected Graphs

First, we consider the simple case of a fully connected graph, in which $E$ is the set of all possible edges, and the obvious mixing matrix is $\mathbf{W}=\mathbf{1}\mathbf{1}^{T}/n$, which has $\lambda_{2}=0$. This represents a best-case scenario in which to validate the theoretical claims made above. We choose $\rho=1/2$ to examine the regime of low communications ratio, and we let $m$ and $T$ grow according to two regimes: $T=m$, and $T=\sqrt{m}$, which are the regimes in which D-SAMD and AD-SAMD are predicted to give order-optimum performance, respectively. The constraint on mini-batch size per Corollaries 1 and 2 is trivial, so we take $b=2$ to ensure that nodes can average each mini-batch gradient via (perfect) consensus. We select the following step-size parameter $\gamma$: $0.5$ and $2$ for (local and centralized) stochastic mirror descent (MD) and accelerated stochastic mirror descent (A-MD), respectively; 5 for both variants of DGD; and $5$ and $20$ for D-SAMD and AD-SAMD, respectively.444While the accelerated variant of stochastic mirror descent makes use of two sequences of step sizes, $\beta_{s}$ and $\gamma_{s}$, these two sequences can be expressed as a function of a single parameter $\gamma$; see, e.g., the proof of Theorem 2. These values were selected via trial-and-error to give good performance for all algorithms; future work involves the use of adaptive step size rules such as AdaGrad and ADAM [30, 31].

In Figure 2(a) and Figure 2(b), we plot the performance averaged over 1200 and 2400 independent instances of the problem, respectively. We also plot the order-wise theoretical performance $1/\sqrt{mT}$, which has a constant slope on the log-log axis. As expected, the distributed methods significantly outperform the local methods. The performance of the distributed methods is on par with the asymptotic theoretical predictions, as seen by the slope of the performance curves, with the possible exception of D-SAMD for $T=m$. However, we observe that D-SAMD performance is at least as good as predicted by theory for $T=\sqrt{m}$, a regime in which optimality is not guaranteed for D-SAMD. This suggests the possibility that the requirement that $T=\Omega(m)$ for D-SAMD is an artifact of the analysis, at least for this problem.

V-B Expander Graphs

For a more realistic setting, we consider expander graphs, which are families of graphs that have spectral gap $1-\lambda_{2}$ bounded away from zero as $m\to\infty$. In particular, we use 6-regular graphs, i.e., regular graphs in which each node has six neighbors, drawn uniformly from the ensemble of such graphs. Because the spectral gap is bounded away from zero for expander graphs, one can more easily examine whether performance of D-SAMD and AD-SAMD agrees with the ideal scaling laws discussed in Corollaries 1 and 2. At the same time, because D-SAMD and AD-SAMD make use of imperfect averaging, expander graphs also allow us to examine non-asymptotic behavior of the two schemes. Per Corollaries 1 and 2, we choose $b=\frac{1}{10}\frac{\log(mT)}{\rho\log(1/\lambda_{2})}$. While this choice is guaranteed to be sufficient for optimum asymptotic performance, we chose the multiplicative constant $1/10$ via trial-and-error to give good non-asymptotic performance.

In Figure 3 we plot the performance averaged over 600 problem instances. We again take $\rho=1/2$, and consider the regimes $T=m$ and $T=\sqrt{m}$. The step sizes are the same as in the previous subsection except that $\gamma=2.5$ for D-SAMD when $T=\sqrt{m}$, $\gamma=28$ for AD-SAMD when $T=m$, and $\gamma=8$ for AD-SAMD when $T=\sqrt{m}$. Again, we see that AD-SAMD and D-SAMD outperform local methods, while their performance is roughly in line with asymptotic theoretical predictions. The performance of DGD, on the other hand, depends on the regime: For $T=m$, it appears to have order-optimum performance, whereas for $T=\sqrt{m}$ it has suboptimum performance on par with local methods. The reason for the dependency on regime is not immediately clear and suggests the need for further study into DGD-style methods in the case of rate-limited networks.

V-C Erdős-Renyi Graphs

Finally, we consider Erdős-Renyi graphs, in which a random fraction (in this case $0.1$) of possible edges are chosen. These graphs are not expanders, and their spectral gaps are not bounded. Therefore, order-optimum performance is not easy to guarantee, since the conditions on the rate and the size of the network depend on $\lambda_{2}$, which is not guaranteed to be well behaved. We again take $\rho=1/2$, consider the regimes $T=m$ and $T=\sqrt{m}$, and again we choose $b=\frac{1}{10}\frac{\log(mT)}{\rho\log(1/\lambda_{2})}$. The step sizes are chosen to be the same as for expander graphs in both regimes.

Once again, we observe a clear distinction in performance between local and distributed methods; in particular, all distributed methods (including DGD) appear to show near-optimum performance in both regimes. However, as expected the performance is somewhat more volatile than in the case of expander graphs, especially for the case of $T=m$, and it is possible that the trends seen in these plots will change as $T$ and $m$ increase.

VI Conclusion

We have presented two distributed schemes, D-SAMD and AD-SAMD, for convex stochastic optimization over networks of nodes that collaborate via rate-limited links. Further, we have derived sufficient conditions for the order-optimum convergence of D-SAMD and AD-SAMD, showing that accelerated mirror descent provides a foundation for distributed SO that better tolerates slow communications links. These results characterize relationships between network communications speed and the convergence speed of stochastic optimization.

A limitation of this work is that we are restricted to settings in which the prox mapping is Lipschitz continuous, which excludes important Bregman divergences such as the Kullbeck-Liebler divergence. Further, the conditions for optimum convergence restrict the Lipschitz constant of non-smooth component of the objective function to be small. Future work in this direction includes study of the limits on convergence speed for more general divergences and composite objective functions.

-A Proof of Lemma 2

The first claim follows immediately from the definitions of the constituent matrices. The second claim follows from Lemma 1. To establish the final claim, we first bound the norm of the columns of $\mathbf{E}(s)$:

[TABLE]

where the first equality follows from $\overline{\mathbf{G}}(s)=\mathbf{G}(s)\mathbf{1}\mathbf{1}^{T}/n$ by definition, and the second equality follows from the fact that the dominant eigenspace of a column-stochastic matrix $\mathbf{W}$ is the subspace of constant vectors, thus the rows of $\bar{\mathbf{G}}(s)$ lie in the null space of $\mathbf{W}^{r}-\mathbf{1}\mathbf{1}^{T}/n$. We bound the norm of the columns of $\mathbf{E}(s)$ via the Frobenius norm of the entire matrix:

[TABLE]

where $C_{1}$ bounds the gap between the Euclidean and dual norms, $C_{2}$ bounds the gap in the opposite direction, and thus $C_{*}\leq C_{1}C_{2}$. Then, we bound the variance of the subgradient noise. Similar to the case of $\mathbf{E}(s)$, we can rewrite

[TABLE]

We bound the expected squared norm of each column $\tilde{\mathbf{z}}_{i}(s)$ in terms of the expected Frobenius norm:

[TABLE]

-B Proof of Lemma 3

By the Lipschitz condition on $P_{\mathbf{x}}(\mathbf{y})$,

[TABLE]

where the final inequality is due to the constraint $\gamma\leq\alpha/(2L)$. Applying Lemma 2, we obtain

[TABLE]

which, applying (22), yields

[TABLE]

Applying (10), we get

[TABLE]

where $i^{*}\in V$ is the index that maximizes $\left\|{\tilde{\mathbf{z}}_{i}}\right\|_{*}$.

The final expression (32) gives the recurrence relationship that governs the divergence of the iterates, and it takes the form of a standard first-order difference equation. Recalling the base case $a_{1}=0$, we obtain the solution

[TABLE]

To prove the lemma, we bound the moments of $a_{s}$ using the solution (-B). First, we bound $E[a_{s}^{2}]$:

[TABLE]

which we can rewrite via a change of variables as

[TABLE]

Then, we apply the Cauchy-Schwartz inequality and (23) to obtain

[TABLE]

Finally, we apply the geometric sum identity and simplify:

[TABLE]

This establishes the result for $E[a_{s}^{2}]$, and the bound on $E[a_{s}]$ follows a fortiori from Jensen’s inequality.

-C Proof of Theorem 1

The proof involves an analysis of mirror descent, with the subgradient error quantified by Lemmas 2 and 3. As necessary, we will cite results from the proof of [8, Theorem 1] rather than reproduce its analysis.

Define

[TABLE]

Applying Lemmas 2 and 3, we bound $E[\left\|{\boldsymbol{\eta}_{i}(\tau)}\right\|_{*}]$:

[TABLE]

Next, we bound the expected value of $\zeta_{i}(\tau)$:

[TABLE]

where the first inequality is due to the definition of the dual norm, and the third equality follows because $\bar{\mathbf{z}}(\tau)$ and $\tilde{\mathbf{z}}_{i}(\tau)$ have zero mean and are independent of $\mathbf{x}^{*}-\mathbf{x}_{i}(\tau)$.

Next, we bound $E[\left\|{\boldsymbol{\delta}_{i}(\tau)}\right\|^{2}_{*}]$. By Lemma 2,

[TABLE]

which we can bound via repeated uses of the triangle inequality:

[TABLE]

Next, we apply Lemma 2:

[TABLE]

Finally, we apply Lemma 3:

[TABLE]

Applying the proof of Theorem 1 of [8] to the D-SAMD iterates, we have that

[TABLE]

We take the expectation to obtain

[TABLE]

Because $\{\Xi_{\tau}\}$ and $\{\Delta_{\tau}^{2}\}$ are increasing sequences, we can bound the preceding by

[TABLE]

Minimizing this bound over $\gamma$ in the interval $(0,\alpha/(2L)]$, we obtain the following optimum value:

[TABLE]

which results in the bound

[TABLE]

-D Proof of Corollary 1

First, we make a few simplifying approximations to $\Xi_{\tau}$:

[TABLE]

where the first inequality is trivial, the next two statements are due to the binomial theorem and the exponential upper bound on the binomial coefficient, and the final inequality is true if we constrain $S\alpha m^{2}\sqrt{C_{*}}\lambda_{2}^{r}\leq 1$, which is consistent with the order-wise constraints listed in the statement of this corollary. Recalling that $r\leq b\rho$, we obtain

[TABLE]

Now, we find the optimum scaling law on the mini-batch size $b$. To ensure $S\alpha m^{2}\sqrt{C_{*}}\lambda_{2}^{r}\leq 1$ for all $s$, we need

[TABLE]

In order to get optimum scaling we need $\Xi_{S}~{}=~{}O(\sigma/\sqrt{mT}+\mathcal{M})$, which yields a slightly stronger necessary condition:

[TABLE]

where the first term is necessary because $b$ cannot be smaller than unity.

Similar approximations establish that

[TABLE]

Thus, (36) ensures $\Delta_{S}^{2}=O(\sigma^{2}/(mb)+\mathcal{M}^{2})$. The resulting gap to optimality scales as

[TABLE]

Substituting (36) yields

[TABLE]

In order to achieve order optimality, we need the first term to be $O((\mathcal{M}+\sigma)/\sqrt{mT})$, which requires

[TABLE]

Finally, we also need the second and fourth terms to be $O((\mathcal{M}~{}+~{}\sigma)/\sqrt{mT})$, which requires

[TABLE]

This establishes the result.

-E Proof of Lemma 4

The sequences $a_{s}$, $b_{s}$, $c_{s}$ are interdependent. Rather than solving the complex recursion directly, we bound the dominant sequence. Define $d_{s}\triangleq\max\{a_{s},b_{s},c_{s}\}$. First, we bound $a_{s}$:

[TABLE]

which, via the triangle inequality, is bounded by

[TABLE]

Because the prox-mapping is Lipschitz,

[TABLE]

Applying the first part of Lemma 2 yields

[TABLE]

and applying the second part and the triangle inequality yields

[TABLE]

We apply (10) and collect terms to obtain

[TABLE]

Next, we bound $b_{s}$ via similar steps:

[TABLE]

Finally, we bound $c_{s}$:

[TABLE]

We bound this sequence in terms of $a_{s}$ and $b_{s}$:

[TABLE]

Comparing the three bounds, we observe that

[TABLE]

Solving the recursion yields

[TABLE]

Following the same argument as in Lemma 3, we obtain

[TABLE]

The bound on $E[d_{s}]$ follows a fortiori.

-F Proof of Theorem 2

Similar to the proof of Theorem 1, the proof involves an analysis of accelerated mirror descent, and as necessary we will cite results from [8, Theorem 2]. Define

[TABLE]

Applying (10) and Lemmas 2 and 4, we bound $E[\left\|{\boldsymbol{\eta}_{i}(\tau)}\right\|]$:

[TABLE]

By the definition of the dual norm:

[TABLE]

We bound $E[\left\|{\boldsymbol{\delta}_{i}(\tau)}\right\|^{2}]$ via the triangle inequality:

[TABLE]

which by (10) and Lemma 2 is bounded by

[TABLE]

Applying Lemma 4, we obtain

[TABLE]

From the proof of [8, Theorem 2], we observe that

[TABLE]

Taking the expectation yields

[TABLE]

Letting $\beta_{\tau}=\frac{\tau+1}{2}$ and $\gamma_{\tau}=\frac{\tau+1}{2}\gamma$, observing that $\Delta_{\tau}$ and $\Xi_{\tau}$ are increasing in $\tau$, and simplifying, we obtain

[TABLE]

Solving for the optimum $\gamma$ in the interval $0\leq\gamma\leq\alpha/(2L)$, we obtain:

[TABLE]

This gives the bound

[TABLE]

which simplifies to the desired bound

[TABLE]

-G Proof of Corollary 2

First, we bound $\Xi_{\tau}$ and $\Delta_{\tau}^{2}$ order-wise. Using steps similar to those in the proof of Corollary 1, one can show that

[TABLE]

and

[TABLE]

In order to achieve order-optimum convergence, we need $\Delta_{S}^{2}~{}=~{}O(\sigma^{2}/(mb)+\mathcal{M}^{2})$. Combining the above equation with the constraint $r\leq b\rho$, we obtain the condition

[TABLE]

This condition also ensures that $\Xi_{S}=O(\sigma/\sqrt{mT}+\mathcal{M})$, as is necessary for optimum convergence speed. This leads to a convergence speed of

[TABLE]

In order to make the first term $O((\mathcal{M}+\sigma)/\sqrt{mT})$, we need

[TABLE]

To make the second and fourth terms $O((\mathcal{M}+\sigma)/\sqrt{mT})$, we need, as before,

[TABLE]