# Communication-Efficient Variance-Reduced Decentralized Stochastic
Optimization over Time-Varying Directed Graphs

Yiyue Chen, Abolfazl Hashemi, Haris Vikalo Yiyue Chen and Haris Vikalo are
with the Department of Electrical and Computer Engineering, University of
Texas at Austin, Austin, TX 78712 USA. Abolfazl Hashemi is with the School of
Electrical and Computer Engineering, Purdue University, West Lafayette, IN
47907, USA. A preliminary version of this article is presented at the 2021
International Conference on Acoustics, Speech, and Signal Processing (ICASSP)
[1].

###### Abstract

We consider the problem of decentralized optimization over time-varying
directed networks. The network nodes can access only their local objectives,
and aim to collaboratively minimize a global function by exchanging messages
with their neighbors. Leveraging sparsification, gradient tracking and
variance-reduction, we propose a novel communication-efficient decentralized
optimization scheme that is suitable for resource-constrained time-varying
directed networks. We prove that in the case of smooth and strongly-convex
objective functions, the proposed scheme achieves an accelerated linear
convergence rate. To our knowledge, this is the first decentralized
optimization framework for time-varying directed networks that achieves such a
convergence rate and applies to settings requiring sparsified communication.
Experimental results on both synthetic and real datasets verify the
theoretical results and demonstrate efficacy of the proposed scheme.

##  1 Introduction

Decentralized optimization problems are encountered in a number of settings in
control, signal processing, and machine learning [2, 3, 4]. Formally, the goal
of a decentralized optimization task is to minimize global objective in the
form of a finite sum

$\min_{\mathbf{x}\in\mathcal{X}}\left[f(\mathbf{x}):=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\mathbf{x})\right],$
(1)

where
$f_{i}(\mathbf{x})=\frac{1}{m_{i}}\sum_{j=1}^{m_{i}}f_{i,j}(\mathbf{x}):\mathbb{R}^{d}\to\mathbb{R}$
for $i\in[n]:=\left\\{1,...,n\right\\}$ denotes the local objective function
that averages loss of $m_{i}$ data points at node $i$, and $\mathcal{X}$
denotes a convex compact constraint set. The $n$ nodes of the network exchange
messages to collaboratively solve (1). Since the communication links between
nodes in real-world networks are often uni-directional and dynamic (i.e.,
time-varying), we model the network by a sequence of directed graphs
${\mathcal{G}}(t)=(|n|,{\mathcal{E}}(t))$, where the existence of an edge
$e_{i,j}\in{\mathcal{E}}(t)$ implies that node $i$ can send messages to node
$j$ at time step $t$.

As the networks and datasets keep increasing in size, computational complexity
and communication cost of decentralized optimization start presenting major
challenges. To reduce the complexity of computing gradients, decentralized
stochastic methods that allow each agent to perform gradient estimation by
processing small subset of local data are preferred [5]. However, such
techniques exhibit low convergence rates and suffer from high variance of the
local stochastic gradients. Remedies for these impediments in decentralized
optimization over networks and in stochastic optimization in centralized
settings include gradient tracking [6, 7] and variance reduction [8, 9],
respectively; however, no such remedies have been developed for decentralized
optimization over time-varying directed networks. On another note,
communication constraints often exacerbate large-scale decentralized
optimization problems where the size of the network or the dimension of local
model parameters may be on the order of millions. This motivates design of
communication-efficient algorithms that compress messages exchanged between
network nodes yet preserve fast convergence. In this paper, we study the
general setting where a network is time-varying and directed, and present, to
the best of our knowledge, the first variance-reduced communication-sparsified
algorithm for decentralized convex optimization over such networks. Moreover,
we theoretically establish that when the local objective functions are smooth
and strongly-convex, the proposed algorithm achieves accelerated linear
convergence rate. Note that while our focus is on communication-constrained
settings, the proposed scheme readily applies to decentralized optimization
problems where time-varying directed graphs operate under no communication
constraints; in fact, to our knowledge, this is the first variance-reduced
stochastic algorithm for such problems.

###  1.1 Related work

The first work on decentralized optimization over networks dates back to the
1980s [10]. A number of studies that followed in subsequent years was focused
on the decentralized average consensus problem, where the network nodes work
collaboratively to compute the average value of local vectors. Convergence
conditions for achieving consensus over directed and undirected time-varying
graphs were established in [3, 2, 11, 12, 13]. The first communication-
efficient algorithm that achieves linear convergence over time-invariant
(static) undirected graphs was proposed in [14].

The consensus problem can be viewed as a stepping stone towards more general
decentralized optimization problems, where the nodes in a network aim to
collaboratively minimize the sum of local objective functions. A number of
solutions to this problem has been proposed for the setting where the network
is undirected, including the well-known distributed (sub)gradient descent
algorithm (DGD) [4, 15], distributed alternating direction method of
multipliers (D-ADMM) [16], and decentralized dual averaging methods [17, 18,
19]. Recently, [20, 14] proposed a novel communication-efficient algorithm for
decentralized convex optimization problems; the provably convergent algorithm
relies on a message-passing scheme with memory and biased compression.

A key technical property required to ensure convergence of decentralized
convex optimization algorithms over undirected networks is that the so-called
mixing matrix characterizing the network connectivity is doubly-stochastic.
However, in directed networks characterized by communication link asymmetry,
doubly-stochastic mixing matrices are atypical. This motivated algorithms that
rely on auxiliary variables to cancel out the imbalance in asymmetric directed
networks in order to achieve convergence. For instance, the subgradient-push
algorithm [21, 22] works with column-stochastic mixing matrices and introduces
local normalization scalars to ensure converge. The directed distributed
gradient descent (D-DGD) algorithm [23], on the other hand, keeps track of the
variations of local models by utilizing auxiliary variables of the same
dimension as the local model parameters. For convex objective functions, both
algorithms achieve ${\mathcal{O}}(\frac{\mathrm{ln}T}{\sqrt{T}})$ convergence
rate. When the objectives are strongly-convex with Lipshitz gradients, and
assuming availability of only the stochastic gradient terms, the stochastic
gradient-push algorithm proposed in [24] achieves
${\mathcal{O}}(\frac{\mathrm{ln}T}{T})$ convergence rate. A common feature of
these algorithms is their reliance upon diminishing stepsizes to achieve
convergence to the optimal solution; in comparison, using a fixed stepsize can
accelerate the gradient search but cannot guarantee the exact convergence,
only the convergence to a neighborhood of the optimal solution. The implied
exactness-speed dilemma can be overcome using schemes that deploy gradient
tracking (see, e.g., [6, 7, 25]). These schemes utilize fixed step sizes to
achieve linear convergence rate when the objective functions are both smooth
and strongly-convex. Among them, the Push-DIGing algorithm [7] follows the
same basic ideas of the subgradient-push algorithm, while TV-AB [25] relies on
column- and row-stochastic matrices to update model parameters and gradient
terms, respectively.

The aforementioned linearly convergent methods rely on full gradient, i.e.,
each node is assumed to use all of its data to compute the local gradient.
However, if the number of data points stored at each node is large, full
gradient computation becomes infeasible. To this end, stochastic gradient
descent was adapted to decentralized settings, but the resulting computational
savings come at the cost of sublinear convergence rate [26, 24]. To accelerate
traditional stochastic gradient methods in centralized settings, variance-
reduction algorithms such as SAG [9] and SVRG [8] have been proposed; these
schemes enable linear convergence when the objective functions are smooth and
strongly-convex. In decentralized settings, GT-SVRG [26] and Network-SVRG [27]
leverage variance-reduction techniques to achieve linear convergence rate.
However, these algorithms are restricted to static and undirected networks,
and a narrow class of directed networks where the mixing matrices can be
rendered doubly stochastic.111To have doubly stochastic mixing matrices,
directed graphs require weight balance, i.e., at each node of a graph, the sum
of the weights from in-coming edges should be equal to that of the weights
from out-coming edges [28]. The existing algorithms for decentralized
optimization over directed networks, such as Push-SAGA [29], are restricted to
static networks.

In recent years, decentralized learning tasks have experienced rapid growth in
the amount of data and the dimension of the optimization problems, which may
lead to practically infeasible demands for communication between network
nodes. To this end, various communication compression schemes have been
proposed; among them, the most frequently encountered are quantization and
sparsification. Quantization schemes limit the number of bits encoding the
messages, while the sparsification schemes select a subset of features and
represent messages in lower dimensional space. For instance, [30, 31, 32, 20,
33, 34] propose algorithms for distributed training of (non)convex machine
learning models in static master-worker settings (i.e., star graph topology)
using quantized/compressed information, while [35, 14, 36, 37] develop
communication-efficient algorithms for decentralized optimization over static
and undirected networks. However, directed networks in general, and time-
varying ones in particular, have received considerably less attention.
Decentralized optimization over such networks faces technical challenges of
developing an algorithmic framework conducive to theoretical analysis and
establishing convergence guarantees, which is further exacerbated when the
communication is compressed. Early steps in this direction were made in [38]
by building upon the subgradient-push algorithm to develop a quantized
communication framework for decentralized optimization over a static network.

Our proposed algorithm utilizes gradient tracking and variance reduction to
achieve fast convergence, and relies on stochastic gradients to solve the
decentralized convex optimization task at feasible computational cost.
Preliminary results of this work, focused on a significantly slower full
gradient framework (${\mathcal{O}}(\frac{1}{\epsilon^{2}})$ vs.
${\mathcal{O}}(\ln\frac{1}{\epsilon})$), were presented in [1]. In Table 1, we
briefly summarize and contrast several algorithms for decentralized
optimization over directed graphs.

###  1.2 Notation

We use lowercase bold letters to represent vectors and uppercase letters to
represent matrices. $[A]_{ij}$ denotes the $(i,j)$ entry of matrix $A$, while
$\|\cdot\|$ denotes the standard Euclidean norm. The spectral radius of a
matrix $A$ is denoted by $\rho(A)$. The weighted infinity norm of
${\mathbf{x}}$ given a positive vector $\mathbf{w}$ is
$\|{\mathbf{x}}\|^{\mathbf{w}}_{\infty}=\max_{i}|x_{i}|/w_{i}$ and the induced
matrix norm is $\||\cdot|\|_{\infty}^{\mathbf{w}}$. Finally, $I$ denotes the
identity matrix whose dimension is inferred from the context.

Table 1: The settings and convergence rates of algorithms for decentralized optimization over directed graphs.  Algorithm | Convergence | Digraph | Gradient | Convex Objective Setting | Compress
---|---|---|---|---|---
Subgradient-push [24] | ${\mathcal{O}}(\frac{1}{\epsilon})$ | Time-varying | Stochastic | Strong convexity | ✗
Push-DIGing [7] | ${\mathcal{O}}(\ln\frac{1}{\epsilon})$ | Time-varying | Full | Strong convexity and smoothness | ✗
TV-AB [25] | ${\mathcal{O}}(\ln\frac{1}{\epsilon})$ | Time-varying | Full | Strong convexity and smoothness | ✗
Quantized Push-sum [38] | ${\mathcal{O}}(\frac{1}{\epsilon^{2}})$ | Static | Full | - | ✓
This work | ${\mathcal{O}}(\ln\frac{1}{\epsilon})$ | Time-varying | Stochastic | Strong convexity and smoothness | ✓

##  2 Preliminaries

###  2.1 Problem Formulation

For convenience, we remind the reader of the problem formulation (1): In a
network of $n$ agents, where each node maintains a local model consisting of
$d$ parameters, the agents collaboratively solve the decentralized convex
optimization

$\min_{\mathrm{\mathbf{x}\in\mathbb{R}^{d}}}\left[f(\mathbf{x}):=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\mathbf{x})\right],$
(2)

where
$f_{i}(\mathbf{x})=\frac{1}{m_{i}}\sum_{j=1}^{m_{i}}f_{i,j}(\mathbf{x}):\mathbb{R}^{d}\to\mathbb{R}$
for $i\in[n]:=\left\\{1,...,n\right\\}$ denotes the local objective function
at node $i$. Each component of $f_{i}$ is assumed to be smooth and not
accessible to nodes other than the $i^{th}$ one, and the global objective $f$
is assumed to be strongly-convex. We further assume existence of a unique
optimal solution ${\mathbf{x}}^{*}\in\mathbb{R}^{d}$ and that each node can
communicate to its neighbors; the nodes identify ${\mathbf{x}}^{*}$ by
exchanging messages over a time-varying directed network. The network’s
connectivity properties are elaborated upon in Section 3.

###  2.2 Communication-Efficient Methods

In practice, bandwidth limitations may restrict the amount of data that the
network nodes can communicate to each other; this is typical of high-
dimensional scenarios where the dimension $d$ of local parameters
${\mathbf{x}}_{i}$ is exceedingly large. To handle communication constraints,
network nodes may employ sparsification to reduce the size of the messages.
Typically, there are two approaches to sparsification: (i) each node selects
and communicates $d_{q}$ out of $d$ components of a $d$-dimensional message;
or (ii) each component of a $d$-dimensional message is selected to be
communicated independently with probability $d_{q}/d$. Note that the former
imposes a hard constraint on the number of communicated entries while the
latter results in $d_{q}$ communicated entries on expectation; both select a
specific entry with probability $d_{q}/d$. Throughout this paper, we focus on
and study the first approach.

Let $Q:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ denote the sparsification
operator; we allow for biased $Q$ with variance proportional to the argument
norm, i.e., $\mathbb{E}[Q({\mathbf{x}})]\neq{\mathbf{x}}$ and
$\mathbb{E}[\|Q({\mathbf{x}})-{\mathbf{x}}\|^{2}]\propto\|{\mathbf{x}}\|^{2}$.
This stands in contrast to typical compression operators which aim to achieve
no bias and have bounded variance (see, e.g., [30]). More recent works [20,
14, 36, 38] do consider biased compression operators but only for time-
invariant communication networks – a setting that is more restrictive than the
one considered in this paper.

##  3 Algorithm development

In this section, we first introduce a novel average consensus algorithm, an
intermediate step towards the main (optimization) framework; then we present
the optimization algorithm consisting of the consensus and the gradient
components.

###  3.1 Foundations: Decentralized average consensus

We start by specifying a procedure for decentralized average consensus, an
intermediate step towards decentralized optimization and, ultimately, an
integral part thereof. The decentralized average consensus problem is
formulated as the computation
$\bar{{\mathbf{x}}}=\frac{1}{n}\sum_{i=1}^{n}{\mathbf{x}}_{i}$, where
${\mathbf{x}}_{i}\in\mathcal{R}^{d}$ is the parameter vector at node $i$.
Following the idea of [12], for each node of the network we define the so-
called surplus vector, i.e., an auxiliary variable
$\mathbf{y}_{i}\in\mathbb{R}^{d}$ which tracks local state vector variations
over consecutive time steps; as shown later in this section, one can use the
surplus vector to help provide guarantees of convergence to the optimal
solution of the decentralized problem. The surplus vector is exchanged along
with the state vector, i.e., at time $t$, node $i$ sends both
$\mathbf{y}_{i}^{t}$ and the state vector $\mathbf{x}_{i}^{t}$ to its out-
neighbors. For the sake of having compact notation, let us introduce
$\mathbf{z}_{i}^{t}\in\mathbb{R}^{d}$,

$\mathbf{z}_{i}^{t}=\begin{cases}\mathbf{x}_{i}^{t},&i\in\left\\{1,...,n\right\\}\\\
\mathbf{y}_{i-n}^{t},&i\in\left\\{n+1,...,2n\right\\},\end{cases}$ (3)

to represent messages node $i$ communicates to its neighbors in the network at
time $t$.

We assume that the time-varying graph is $\mathcal{B}$-jointly connected,
i.e., that there exists a window size $\mathcal{B}\geq 1$ such that the
aggregate graph $\bigcup_{l=t}^{t+\mathcal{B}-1}{\mathcal{G}}_{l}$ is strongly
connected for all $t=k\mathcal{B}$, $k\in\mathcal{N}$. Note that if
$\mathcal{B}=1$, each instance of the graph is strongly connected. This is a
more general assumption than the often used $\mathcal{B}$-bounded strong-
connectivity (see, e.g. [25]) which requires strong connectivity of the union
graph $\bigcup_{l=t}^{t+\mathcal{B}-1}{\mathcal{G}}_{l}$ for all $t\geq 0$.

For $\mathcal{B}$-jointly connected graphs, the product of mixing matrices of
graph instances over $\mathcal{B}$ consecutive time steps has a non-zero
spectral gap. To formalize this statement, let us construct two weight
matrices that reflect the network topology; in particular, let $W_{in}^{t}$
(row-stochastic) and $W_{out}^{t}$ (column-stochastic) denote the in-neighbor
and out-neighbor weight matrices at time $t$, respectively. It holds that
$[W_{in}^{t}]_{ij}>0$ if and only if $j\in\mathcal{N}_{{in},i}^{t}$ and
$[W_{out}^{t}]_{ij}>0$ if and only if $i\in\mathcal{N}_{out,j}^{t}$, where
$\mathcal{N}_{{in},i}^{t}$ denotes the set of nodes that may send information
to node $i$ (including $i$) whereas $\mathcal{N}_{out,j}^{t}$ denotes the set
of nodes that may receive information from node $j$ (including $j$) at time
$t$. We assume $W_{in}^{t}$ and $W_{out}^{t}$ are given and that both
$\mathcal{N}_{{in},i}^{t}$ and $\mathcal{N}_{out,i}^{t}$ are known to node
$i$. A common policy for designing $W^{t}_{in}$ and $W^{t}_{out}$ is to assign

$[W^{t}_{in}]_{ij}=1/|\mathcal{N}^{t}_{in,i}|,\qquad[W^{t}_{out}]_{ij}=1/|\mathcal{N}^{t}_{out,j}|.$
(4)

Recall that we are interested in sparsifying messages exchanged between nodes
of a network; clearly, the sparsification should impact the structure of a
mixing matrix. Indeed, if one attempts sparsifying messages used by existing
methods, e.g. [21, 12, 13, 22], without any modifications of the mixing
matrices therein, non-vanishing error terms induced by the compression
operator will prevent those methods from converging. We note that the impact
of sparsification on the components of a message vector is similar to the
impact of link failures, and may therefore be captured by the structure of the
weight matrices. To elaborate on this, observe that the vector-valued problem
at time $t$ can essentially be decomposed to $d$ individual scalar-valued
tasks with weight matrices $\\{W_{in,m}^{t}\\}_{m=1}^{d}$ and
$\\{W_{out,m}^{t}\\}_{m=1}^{d}$. For the sparsified components of the message
vector, i.e., those that are set to zero and not communicated, the
corresponding entries in the weight matrices can be replaced by zeros; on the
other hand, the entries in the weight matrices corresponding to the
communicated components of the message vector remain unchanged, leading to the
violation of the stochasticity of the weight matrices. To address this, we re-
normalize the weight matrices $\\{W_{in,m}^{t}\\}_{m=1}^{d}$ and
$\\{W_{out,m}^{t}\\}_{m=1}^{d}$, thus ensuring their row and column
stochasticity. Note that the re-normalization of the $i\textsuperscript{th}$
row of $\\{W_{in,m}^{t}\\}_{m=1}^{d}$ ($i\textsuperscript{th}$ column of
$\\{W_{out,m}^{t}\\}_{m=1}^{d}$) is performed by the $i^{th}$ network node.

To specify the normalization rule, we first need to define the sparsification
operation. Sparsification of $\mathbf{x}_{i}^{t}$ (and, consequently,
$\mathbf{y}_{i}^{t}$) is done via the compression operator $Q(\cdot)$ applied
to $\mathbf{z}_{i}^{t}$; we denote the result by $Q(\mathbf{z}_{i}^{t})$. Let
$[Q({{\mathbf{z}}}_{i}^{t})]_{m}$ denote the $m\textsuperscript{th}$ component
of $Q(\mathbf{z}_{i}^{t})$. Let $\\{A_{m}^{t}\\}_{m=1}^{d}$ and
$\\{B_{m}^{t}\\}_{m=1}^{d}$ be the weight matrices obtained after normalizing
$\\{W_{in,m}^{t}\\}_{m=1}^{d}$ and $\\{W_{out,m}^{t}\\}_{m=1}^{d}$,
respectively. To formalize the normalization procedure, we introduce the
weight matrix

$[A^{t}_{m}]_{ij}=\begin{cases}\frac{[W^{t}_{in,m}]_{ij}}{\sum_{j\in\mathcal{S}_{m}^{t}(i,j)}[W^{t}_{in,m}]_{ij}}&\text{if
}j\in\mathcal{S}_{m}^{t}(i,j)\\\ 0&\mathrm{otherwise},\end{cases}$ (5)

where
$\mathcal{S}_{m}^{t}(i,j):=\\{j|j\in\mathcal{N}^{t}_{in,i},[Q({{\mathbf{z}}}_{j}^{t})]_{m}\neq
0\\}\cup\\{i\\}$. Likewise, the weight matrix $B^{t}_{m}$ is defined as

$[B^{t}_{m}]_{ij}=\begin{cases}\frac{[W^{t}_{out,m}]_{ij}}{\sum_{i\in\mathcal{T}_{m}^{t}(i,j)}[W^{t}_{out,m}]_{ij}}&\text{if
}i\in\mathcal{T}_{m}^{t}(i,j)\\\ 0&\mathrm{otherwise},\end{cases}$ (6)

where
$\mathcal{T}_{m}^{t}(i,j):=\\{i|i\in\mathcal{N}^{t}_{out,j},[Q({{\mathbf{z}}}_{i}^{t})]_{m}\neq
0\\}\cup\\{j\\}$.

We can now define the mixing matrix of a directed network with sparsified
messages.

######  Definition 1.

The $m\textsuperscript{th}$ mixing matrix at time $t$ of a time-varying
directed network with sparsified messages,
$\bar{M}_{m}^{t}\in\mathbb{R}^{2n\times 2n}$, is a matrix whose columns sum up
to $1$ and whose eigenvalues satisfy
$1=|\lambda_{1}(\bar{M}_{m}^{t})|=|\lambda_{2}(\bar{M}_{m}^{t})|\geq|\lambda_{3}(\bar{M}_{m}^{t})|\geq\cdots|\lambda_{2n}(\bar{M}_{m}^{t})|$,
constructed from the current network topology as

$\bar{M}_{m}^{t}=\left[\begin{matrix}A_{m}^{t}&\mathbf{0}\\\
I-A_{m}^{t}&B_{m}^{t}\\\ \end{matrix}\right],$ (7)

where $A_{m}^{t}$ and $B_{m}^{t}$ denote the $m\textsuperscript{th}$
normalized in-neighbor and out-neighbor weight matrices at time $t$,
respectively.

Given $\mathbf{z}_{i}^{t}$ and $\bar{M}_{m}^{t}$ in (3) and (7), respectively,
we can formulate a compact recursive update rule for $\mathbf{z}_{i}^{t}$ as

$\displaystyle z_{im}^{t+1}$
$\displaystyle=\sum_{j=1}^{2n}[\bar{M}^{t}_{m}]_{ij}[Q({\mathbf{z}}_{j}^{t})]_{m}+\mathbbm{1}_{\left\\{t\
\text{mod}\
\mathcal{B}=\mathcal{B}-1\right\\}}\gamma[F]_{ij}z_{jm}^{\mathcal{B}\lfloor
t/\mathcal{B}\rfloor},$ (8)

where $F=\left[\begin{matrix}\mathbf{0}&I\\\ \mathbf{0}&-I\end{matrix}\right]$
and $m$ denotes the coordinate index.

As seen in (8), vectors ${\mathbf{z}}_{i}^{t}$ (which contain
${\mathbf{x}}_{i}^{t}$, the quantities to be averaged) are updated in a
straightforward manner via sparsification and multiplication with the mixing
matrix at all times $t$ except those that satisfy

$t\mod\mathcal{B}=\mathcal{B}-1.$ (9)

In particular, when (9) holds, vectors ${\mathbf{z}}_{i}^{\mathcal{B}\lfloor
t/\mathcal{B}\rfloor}$, stored at time $\mathcal{B}\lfloor
t/\mathcal{B}\rfloor$, are also used to update ${\mathbf{z}}_{i}^{t}$. The
usage of the stored vectors is motivated by the observation that
$\bar{M}_{m}^{t}$ may have zero spectral gap, which is undesirable since for
such mixing matrices the convergence of the consensus algorithms will not be
guaranteed. However, for a judiciously chosen perturbation parameter
$\epsilon$, which determines to what extent
$\sum_{j=1}^{2n}[F]_{ij}z_{jm}^{\mathcal{B}\lfloor t/\mathcal{B}\rfloor}$
affects the update, we can ensure a nonzero spectral gap of the product of
$\mathcal{B}$ consecutive mixing matrices starting from $t=k\mathcal{B}$. The
described communication-sparsified average consensus procedure over directed
graphs, referred to for convenience as Di-CS-AC, is formalized as Algorithm 1.

Algorithm 1 Directed Communication-Sparsified Average Consensus (Di-CS-AC)

1: Input: $T$, $\mathbf{x}^{0}$, $\mathbf{y}^{0}=\mathbf{0}$, $\gamma$

2: set $\mathbf{z}^{0}=[\mathbf{x}^{0};\mathbf{y}^{0}]$,
$\tilde{w}^{0}={\mathbf{z}}^{0}$ and

3: for each $s\in[0,1,...,S]$ do

4: generate non-negative matrices $\\{W_{in,m}^{t}\\}_{m=1}^{d}$ and
$\\{W_{out,m}^{t}\\}_{m=1}^{d}$

5: for each $m\in[1,...,d]$ do

6: construct a row-stochastic $A^{t}_{m}$ and a column-stochastic $B^{t}_{m}$
according to (5) and (6)

7: construct $\bar{M}^{t}_{m}$ according to (7)

8: for each $i\in[1,...,2n]$ do

9: update $z_{im}^{t+1}$ according to (8)

10: end for

11: end for

12: end for

###  3.2 Decentralized gradient component

Going beyond the simple consensus problem and towards solving optimization
(2), we re-define the recursive update rule for $\mathbf{z}_{i}^{t}$ as

$\displaystyle z_{im}^{t+1}$
$\displaystyle=\sum_{j=1}^{2n}[\bar{M}^{t}_{m}]_{ij}[Q({\mathbf{z}}_{j}^{t})]_{m}+\mathbbm{1}_{\left\\{t\
\text{mod}\
\mathcal{B}=\mathcal{B}-1\right\\}}\gamma[F]_{ij}z_{jm}^{\mathcal{B}\lfloor
t/\mathcal{B}\rfloor}-\mathbbm{1}_{\left\\{t\ \text{mod}\
\mathcal{B}=\mathcal{B}-1\right\\}}\alpha g_{im}^{\mathcal{B}\lfloor
t/\mathcal{B}\rfloor},$ (10)

where $F$ and $m$ denote the same objects as in (8), and $g_{im}$ combines
global gradient tracking with local stochastic reduction to achieve
accelerated convergence (elaborated upon shortly). Note that (10) implies the
following element-wise update rules for state and surplus vectors,
respectively:

$\displaystyle x_{im}^{t+1}$
$\displaystyle=\sum_{j=1}^{n}[A_{m}^{t}]_{ij}[Q({\mathbf{x}}_{j}^{t})]_{m}+\mathbbm{1}_{\left\\{t\
\text{mod}\ \mathcal{B}=\mathcal{B}-1\right\\}}\gamma
y_{im}^{\mathcal{B}\lfloor t/\mathcal{B}\rfloor}-\mathbbm{1}_{\left\\{t\
\text{mod}\ \mathcal{B}=\mathcal{B}-1\right\\}}\alpha
g_{im}^{\mathcal{B}\lfloor t/\mathcal{B}\rfloor},$ (11) $\displaystyle
y_{im}^{t+1}$
$\displaystyle=\sum_{j=1}^{n}[B_{m}^{t}]_{ij}[Q({\mathbf{y}}_{j}^{t})]_{m}-(x_{im}^{t+1}-x_{im}^{t}).$
(12)

Paralleling the basic consensus task discussed in the previous section,
vectors ${\mathbf{z}}_{i}^{t}$ (containing state vectors to be averaged) are
updated via sparsification and multiplication with the mixing matrix at all
times $t$ except those that satisfy (9). When (9) does hold, vectors
${\mathbf{z}}_{i}^{\mathcal{B}\lfloor t/\mathcal{B}\rfloor}$, stored at times
$\mathcal{B}\lfloor t/\mathcal{B}\rfloor$, are also used to update
${\mathbf{z}}_{i}^{t}$; the motivation and reasoning for such special
treatment are as same as in the consensus algorithm.222Note that $F$ has all-
zero matrices for its $(1,1)$ and $(2,1)$ blocks and thus we only need to
store ${\mathbf{z}}_{i}^{\mathcal{B}\lfloor t/\mathcal{B}\rfloor}$
(equivalently, ${\mathbf{y}}_{i-n}^{\mathcal{B}\lfloor
t/\mathcal{B}\rfloor}$), where $n+1\leq i\leq 2n$.

In the proposed algorithm, updates of the gradient term ${\mathbf{g}}_{i}^{t}$
combine global gradient tracking with local stochastic variance reduction. In
particular, the updates of ${\mathbf{g}}_{i}^{t}$ mix gradient messages while
keeping track of the changes in gradient estimates ${\mathbf{v}}_{i}^{t}$;
this guides ${\mathbf{g}}_{i}^{t}$ towards the gradient of the global
objective, ultimately ensuring convergence to the optimal solution
$\mathbf{x}^{*}$ (i.e., global gradient tracking helps avoid the pitfall of
non-vanishing local gradients which would otherwise lead the search only to a
neighborhood of $\mathbf{x}^{*}$). The $m\textsuperscript{th}$ entry of
${\mathbf{g}}_{i}^{t}$, $g_{im}^{t}$, is updated as

$\displaystyle g_{im}^{\mathcal{B}(\lfloor
t/\mathcal{B}\rfloor)}=\begin{cases}\sum_{j=1}^{n}[B_{m}(k\mathcal{B}-1:(k-1)\mathcal{B})]_{ij}g_{jm}^{(k-1)\mathcal{B}}+v_{im}^{\mathcal{B}(\lfloor
t/\mathcal{B}\rfloor)}-v_{im}^{\mathcal{B}(\lfloor
t/\mathcal{B}\rfloor-1)}&i\leq n\\\ 0&\mathrm{otherwise},\end{cases}$ (13)

where $k=\lfloor t/\mathcal{B}\rfloor$. The gradient estimate
${\mathbf{v}}_{i}^{t}$ in (13) is updated via the stochastic variance-
reduction method [8]. Specifically,

$\mathbf{v}_{i}^{\mathcal{B}(\lfloor t/\mathcal{B}\rfloor+1)}=\nabla
f_{i,l_{i}}({\mathbf{z}}_{i}^{\mathcal{B}\lfloor t/\mathcal{B}\rfloor})-\nabla
f_{i,l_{i}}(\tilde{w}_{i})+\tilde{\mu}_{i},\quad\forall i\in[n].$ (14)

One can interpret this update as being executed in a double loop fashion: when
a local full gradient at node $i$, $\tilde{\mu}_{i}$, is computed (in what can
be considered an outer loop), it is retained in the subsequent $T$ iterations
(the inner loop). In each iteration of the inner loop, if the time step
satisfies (9), node $i$ uniformly at random selects a local sample, $l_{i}$,
for the calculation of two stochastic gradient estimates – an estimate of the
current state, $\nabla f_{i,l_{i}}({\mathbf{z}}_{i}^{\mathcal{B}\lfloor
t/\mathcal{B}\rfloor})$, and an estimate of the state from the last outer
loop, $\nabla f_{i,l_{i}}(\tilde{w}_{i})$ – the terms needed to perform update
of $\mathbf{v}_{i}^{t}$. By computing a full gradient periodically in the
outer loop and estimating gradient stochastically in the inner loop, the
described procedure trades computational cost for convergence speed,
ultimately achieving linear convergence at fewer gradient computations per
sample than the full gradient techniques.

The described procedure is formalized as Algorithm 2.

######  Remark 1.

We highlight a few important observations regarding Algorithm 2.

  1. (a)

When there are no communication constraint and each agent in the network can
send full information to out-neighboring agents, Algorithm 2 reduces to a
novel stochastic variance-reduced scheme for decentralized convex optimization
over such networks.

  2. (b)

For $\mathcal{B}=1$, the problem reduces to decentralized optimization over
networks that are strongly connected at all time steps, a typical connectivity
assumption for many decentralized optimization algorithms [23, 14].

  3. (c)

Algorithm 2 requires each node in the network to store local vectors of size
$4d$, including the current state vector, current and past surplus vector, and
local gradient vector. While the current state vector and current surplus
vector may be communicated to the neighboring nodes, past surplus vectors are
only used locally to add local perturbations at the time steps satisfying (9).

  4. (d)

The columns of $\bar{M}_{m}^{t}$ sum up to one. However, $\bar{M}_{m}^{t}$ is
not column-stochastic as it has negative entries, which stands in contrast to
the stochasticity property of the mixing matrices appearing in the average
consensus algorithms [39, 14].

Algorithm 2 Directed Communication-Sparsified Stochastic Variance-Reduced
Gradient Descent (Di-CS-SVRG)

1: Input: $T$, $\mathbf{x}^{0}$, $\mathbf{y}^{0}=\mathbf{0}$, $\alpha$,
$\gamma$

2: set $\mathbf{z}^{0}=[\mathbf{x}^{0};\mathbf{y}^{0}]$,
$\tilde{w}^{0}={\mathbf{z}}^{0}$ and
${\mathbf{g}}_{i}^{0}={\mathbf{v}}_{i}^{0}=\nabla{\mathbf{f}}_{i}({\mathbf{x}}_{i}^{0})\quad\forall
i\in[n]$

3: for each $s\in[0,1,...,S]$ do

4: $\tilde{w}=\tilde{w}^{s}$

5: $\tilde{\mu}_{i}=\nabla
f_{i}(\tilde{w})=\frac{1}{m_{i}}\sum_{j=1}^{m_{i}}\nabla f_{i,j}(\tilde{w})$

6: for each $t\in[sT+1,...,(s+1)T-1]$ do

7: generate non-negative matrices $\\{W_{in,m}^{t}\\}_{m=1}^{d}$ and
$\\{W_{out,m}^{t}\\}_{m=1}^{d}$

8: for each $m\in[1,...,d]$ do

9: construct a row-stochastic $A^{t}_{m}$ and a column-stochastic $B^{t}_{m}$
according to (5) and (6)

10: construct $\bar{M}^{t}_{m}$ according to (7)

11: for each $i\in[1,...,2n]$ do

12: update $z_{im}^{t+1}$ according to (10)

13: end for

14: if $t\mod\mathcal{B}=\mathcal{B}-1$ then

15: for each $i\in[1,...,n]$ do

16: select $l_{i}$ uniformly randomly from $[m_{i}]$:

17: update $\mathbf{v}_{i}^{\mathcal{B}(\lfloor t/\mathcal{B}\rfloor+1)}$
according to (14)

18: update $g_{im}^{\mathcal{B}(\lfloor t/\mathcal{B}\rfloor+1)}$ according to
(13)

19: end for

20: end if

21: end for

22: end for

23: $\tilde{w}^{s+1}={\mathbf{z}}^{(s+1)T}$

24: end for

##  4 Convergence Analysis

For convenience, let us denote the product of a sequence of mixing matrices
from time step $s$ to $T$ as

$\bar{M}_{m}(T:s)=\bar{M}_{m}^{T}\bar{M}_{m}^{T-1}\cdots\bar{M}_{m}^{s}.$ (15)

To further simplify notation, we also introduce

$M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})=\bar{M}_{m}((k+1)\mathcal{B}-1:k\mathcal{B})+\gamma
F,$ (16)

and

$\displaystyle
M_{m}(t:k_{1}\mathcal{B})=\bar{M}_{m}(t:k_{2}\mathcal{B})M_{m}(k_{2}\mathcal{B}-1:(k_{2}-1)\mathcal{B})\cdots
M_{m}((k_{1}+1)\mathcal{B}-1:k_{1}\mathcal{B}),$ (17)

where $k_{2}\mathcal{B}\leq t\leq(k_{2}+1)\mathcal{B}-1$ and
$k_{1},k_{2}\in\mathcal{N},k_{1}\leq k_{2}$. Note that
$M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$ is formed by adding a perturbation
matrix $\gamma F$ to the product
$\bar{M}_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$. Finally, we also introduce
shorthand notation for the product of the weight matrices $B_{m}$ from time
$s$ to $T$,

$B_{m}(T:s)=B_{m}^{T}B_{m}^{T-1}\cdots B_{m}^{s}.$ (18)

Our analysis relies on several standard assumptions about the graph and
network connectivity matrices as well as the characteristics of the local and
global objective functions. These are given next.

######  Assumption 1.

Suppose the following conditions hold:

  1. (a)

The product of consecutive mixing matrices
$M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$ in (16) has a non-zero spectral gap
for all $k\geq 0$, $1\leq m\leq d$, and all $0<\gamma<\gamma_{0}$ for some
$0<\gamma_{0}<1$.

  2. (b)

The collection of all possible mixing matrices $\\{\bar{M}_{m}^{t}\\}$ is a
finite set.

  3. (c)

Each component of the local objective function $f_{i,j}$ is $L$-smooth and the
global objective $f$ is $\mu$-strongly-convex333This implies that
$\forall{\mathbf{x}}_{1},{\mathbf{x}}_{2}\in{\mathbb{R}}^{d}$ there exists
$L>0$ such that $\|\nabla f_{i,j}({\mathbf{x}}_{1})-\nabla
f_{i,j}({\mathbf{x}}_{2})\|\leq L\|{\mathbf{x}}_{1}-{\mathbf{x}}_{2}\|$.
Furthermore, $\forall{\mathbf{x}}_{1},{\mathbf{x}}_{2}\in{\mathbb{R}}^{d}$
there exists $\mu>0$ such that $f({\mathbf{x}}_{2})\geq
f({\mathbf{x}}_{1})+\langle\nabla
f({\mathbf{x}}_{1}),{\mathbf{x}}_{2}-{\mathbf{x}}_{1}\rangle+\frac{\mu}{2}\|{\mathbf{x}}_{2}-{\mathbf{x}}_{1}\|^{2}$.
.

######  Remark 2.

Assumption 1(a) is readily satisfied for a variety of graph structures such as
the $\mathcal{B}$-strongly connected directed graph introduced in [7], i.e.,
the setting where the union of graphs over $B$ consecutive instances starting
from $k\mathcal{B}$ forms a strongly connected graph for any non-negative
integer $k$.444There are two versions of the definition of
$\mathcal{B}$-strongly connected directed graphs, the difference being the
window starting time. As noted in Section II, we consider the definition where
the window may start at any time $t=k\mathcal{B}$; this differs from the (more
demanding in regards to connectivity) definition in [25] where the starting
time is an arbitrary non-negative integer. Furthermore, one can readily verify
that Assumption 1(b) holds for the weight matrices defined in (4).

Before stating the main theorem, we provide the following lemma which, under
Assumption 1, establishes the consensus contraction of the product of mixing
matrices and the product of normalized weight matrices.

######  Lemma 1.

Suppose Assumptions 1(a) and 1(b) hold. Let
$\sigma=\max(|\lambda_{M,2}|,|\lambda_{B,2}|)$ denote the larger of the second
largest eigenvalues of $M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$ and
$B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$. Then,

$\displaystyle\begin{split}\|M_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{z}}-\bar{{\mathbf{z}}}\|\leq\sigma\|{\mathbf{z}}-\bar{{\mathbf{z}}}\|,\
\forall{\mathbf{z}}\in{\mathbb{R}}^{2n}\\\ \mbox{and}\\\
\|B_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{y}}-\bar{{\mathbf{y}}}\|\leq\sigma\|{\mathbf{y}}-\bar{{\mathbf{y}}}\|,\
\forall{\mathbf{y}}\in{\mathbb{R}}^{n},\\\ \end{split}$ (19)

where
$\bar{{\mathbf{z}}}=[\frac{1}{n}\sum_{i=1}^{2n}z_{i},\cdots,\frac{1}{n}\sum_{i=1}^{2n}z_{i}]^{T}$
and
$\bar{{\mathbf{y}}}=[\frac{1}{n}\sum_{i=1}^{n}y_{i},\cdots,\frac{1}{n}\sum_{i=1}^{n}y_{i}]^{T}$.

###### Proof.

The proof of the lemma is in the supplementary document. ∎

Our main result, stated in Theorem 1 below, establishes that Algorithm 2
provides linear convergence of local parameters to their average values, which
itself converges linearly to the optimal solution of (1).

######  Theorem 1.

Suppose Assumption 1 holds. Denote the condition number of $f$ by
$\tilde{Q}=\frac{L}{\mu}$. If the step size $\alpha$ is chosen according to

$\alpha=\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L},$ (20)

the iterates of Algorithm 2 satisfy

$\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\|\bar{{\mathbf{z}}}^{ST}-{\mathbf{z}}_{i}^{ST}\|^{2}+\mathbb{E}\|\bar{{\mathbf{z}}}^{ST}-{\mathbf{x}}^{\ast}\|^{2}\leq
2\lambda^{S}\times\Bigg{(}\frac{1}{n}\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{0}-{\mathbf{z}}_{i}^{0}\|^{2}+\|\bar{{\mathbf{z}}}^{0}-{\mathbf{x}}^{\ast}\|^{2}\\\
+\frac{(1-\sigma^{2})^{2}}{1457nL^{2}}\sum_{i=1}^{n}\sum_{m=1}^{d}\mathbb{E}|g_{im}^{0}-\bar{g}_{m}^{0}|^{2}\Bigg{)},$
(21)

where

$\lambda=8\tilde{Q}^{2}\exp{(-\frac{(1-\sigma^{2})^{2}T}{748\tilde{Q}^{2}})}+0.66,$
(22)

$\bar{{\mathbf{z}}}^{t}=\frac{1}{n}\sum_{i=1}^{2n}{\mathbf{z}}_{i}^{t}$, and
$\bar{{\mathbf{g}}}^{t}=\frac{1}{n}\sum_{i=1}^{n}{\mathbf{g}}_{i}^{t}$.

######  Corollary 1.1.

Instate the notation and hypotheses of Theorem 1. If, in addition, the inner-
loop duration $T$ is chosen as

$T=\mathcal{B}\lceil\frac{1496\tilde{Q}^{2}}{(1-\sigma^{2})^{2}\mathcal{B}}\ln(200\tilde{Q}^{2})\rceil,$
(23)

the proposed algorithm achieves a linear convergence rate. Furthermore, to
reach an $\epsilon$-accurate solution, Algorithm 2 takes at most
${\mathcal{O}}(\frac{\tilde{Q}^{2}\mathcal{B}\ln\tilde{Q}}{(1-\sigma^{2})^{2}}\ln
1/\epsilon)$ communication rounds and performs
${\mathcal{O}}((\frac{\tilde{Q}^{2}\ln\tilde{Q}}{(1-\sigma^{2})^{2}}+\max_{i}\left\\{m_{i}\right\\})\ln
1/\epsilon)$ stochastic gradient computations.

###### Proof.

It is straightforward to verify that for the stated value of $T$, $\lambda\leq
0.7<1$ and thus the algorithm converges linearly. ∎

Note that due to the gradient tracking step in Algorithm 2, in particular when
constructing the linear system of inequalities (which includes the gradient
tracking error), the rate of convergence is dependent upon $\tilde{Q}^{2}$
(through the factor in the coefficient matrix).

######  Remark 3.

Clearly, the communication cost of Algorithm 2 depends on the level of
sparsification, i.e., the value of parameter $d_{q}$. Intuitively, if agents
communicate fewer entries in each round, the communication cost per round
decreases but the algorithm may take more rounds to reach the same accuracy.
Therefore, the total communication cost until reaching a pre-specified
$\epsilon$-accuracy, found as the product of the number of communication
rounds and the cost per round, is of interest. Let $q$ denote the fraction of
entries being communicated per iteration; smaller $q$ implies more aggressive
sparsification. This compression level parameter, $q$, impacts $\sigma$ in
Theorem 1; in particular, for a fixed network connectivity parameter
$\mathcal{B}$, smaller $q$ leads to sparser mixing matrices and, consequently,
greater $\sigma$. Note that large $\mathcal{B}$ may be caused by sparsity of
the instances of a time-varying network, thus leading to large values of
$\sigma$.

######  Remark 4.

It is worthwhile discussing and comparing the constants in Corollary 1.1 to
those in the original SVRG [8] (centralized optimization) and GT-SVRG [26]
(decentralized optimization over undirected graphs). For SVRG, this constant
is $O(\frac{1}{\mu\alpha(1-2L\alpha)T}+\frac{2L\alpha}{1-2L\alpha})$, where
$\alpha$ denotes the step size, $L$ is the smoothness parameter, $\mu$ is the
strong convexity parameter and $T$ denotes the inner loop duration [8]. For
both GT-SVRG and our proposed algorithm, the inner loop duration is
$T=O(\frac{\tilde{Q}^{2}\log\tilde{Q}}{(1-\sigma)^{2}})$ and the linear
convergence constant
$O(\tilde{Q}^{2}\exp{(-\frac{(1-\sigma^{2})^{2}T}{\tilde{Q}^{2}})})$, where
$\tilde{Q}$ is the condition number and $\sigma$ is specified by the network
topology and the applied compression, i.e., sparsification of the communicated
quantities.

###  4.1 Proof of Theorem 1

In this section, we prove Theorem 1 by analyzing various error terms that
collectively impact the convergence rate of Algorithm 2. The main technical
challenge in deriving the linear convergence result is the analysis of the
vanishing errors formally introduced in the next paragraph: the consensus
error, the optimality error and the gradient tracking error. Note that unlike
in undirected graphs, the mixing matrices of directed graphs are not
necessarily doubly stochastic; as a result, decentralized optimization schemes
may produce state vectors that converge to a weighted average, rather than a
consensus average. Furthermore, recall that in order to accelerate the
convergence, we deploy two techniques: global gradient tracking and local
variance reduction. Both the gradient tracking technique, which relies on the
communication of gradient information over the network, and the variance
reduction trick increase the difficulty of analyzing the vanishing gradient
tracking error. The combination of these issues renders the analysis of the
aforementioned errors challenging.

Specifically, the convergence rate depends on: (i) the expected consensus
error, i.e., the expected squared difference between local vectors and the
average vectors at time $(k+1)\mathcal{B}$,
$\mathbb{E}[|z_{im}^{(k+1)\mathcal{B}}-\bar{z}_{m}^{(k+1)\mathcal{B}}|^{2}]$;
(ii) the expected optimality error, i.e., the expected squared difference
between the average vectors and the optimal vector,
$\mathbb{E}[\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$;
and (iii) the expected gradient tracking error,
$\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}$.
Hence, it is critical to determine evolution of these sequences. Note that
compared to the gradient-tracking based work, e.g. [7, 26], analyzing the
proposed scheme is more involved due to its reliance upon a combination of
variance reduction techniques and a communication-sparsified consensus;
showing that the novel scheme achieves linear convergence on general directed
time-varying graphs despite sparsified communication calls for a careful
examination of the error terms in a manner distinct from the analysis found in
prior work.

Dynamics of the aforementioned errors are clearly interconnected.
Consequently, our analysis relies on deriving recursive bounds for the errors
in terms of the linear combinations of their past values. The results are
formally stated in Lemma 2, Lemma 3, and Lemma 4. Proofs of these lemmas are
provided in the supplementary document.

######  Lemma 2.

Suppose Assumption 1 holds. Then $\forall i\leq n$, $k\geq 0$ and $0<m\leq d$,
updates generated by Algorithm 2 satisfy

$\displaystyle\mathbb{E}[|z_{im}^{(k+1)\mathcal{B}}-\bar{z}_{m}^{(k+1)\mathcal{B}}|^{2}]$
$\displaystyle\leq\frac{1+\sigma^{2}}{2}\mathbb{E}[|z_{im}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}|^{2}]+\frac{2\alpha^{2}}{1-\sigma^{2}}\mathbb{E}[|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}|^{2}].$
(24)

Having established in Lemma 2 a recursive bound on the expected consensus
error, we proceed by stating in Lemmas 3 and 4 recursive bounds on the
expected optimality and gradient tracking errors, respectively. First, let us
introduce (for $k\geq 0$)

$\displaystyle\bar{\tau}^{k\mathcal{B}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\tau_{i}^{k\mathcal{B}},$ (25)
$\displaystyle\tau_{i}^{(k+1)\mathcal{B}}$
$\displaystyle=\begin{cases}{\mathbf{x}}_{i}^{(k+1)\mathcal{B}}&\mathrm{if}\quad(k+1)\mathcal{B}\mod
T=0\\\ \tilde{w}_{i}&\mathrm{otherwise}.\end{cases}$

######  Lemma 3.

Suppose Assumption 1 holds and let $0<\alpha<\frac{\mu}{8L^{2}}$. Then for all
$k>0$ it holds that

$\displaystyle\begin{split}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]&\leq\frac{2L^{2}\alpha}{\mu}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]\\\
&\quad+(1-\frac{\mu\alpha}{2})\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].\end{split}$
(26)

######  Lemma 4.

Suppose Assumption 1 holds. Then

$\displaystyle\frac{1}{L^{2}}\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}]$
$\displaystyle\leq\frac{120}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]$
(27)
$\displaystyle\quad+\frac{89}{1-\sigma^{2}}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$
$\displaystyle\quad+\frac{3+\sigma^{2}}{4}\mathbb{E}[\frac{\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}|^{2}}{L^{2}}]$
$\displaystyle\quad+\frac{38}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]$
$\displaystyle\quad+\frac{38}{1-\sigma^{2}}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].$

We proceed by defining a system of linear inequalities involving the three
previously discussed error terms; study of the conditions for the geometric
convergence of the powers of the resultant matrix in the system of linear
inequalities leads to the linear convergence result in Theorem 1. To this end,
we first state Proposition 1 whose proof follows by combining and re-
organizing the inequalities in Lemmas 2-4 in a matrix form.

######  Proposition 1.

Suppose Assumption 1 holds. Define

${\mathbf{u}}^{k\mathcal{B}}=\begin{bmatrix}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]\\\
\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
\mathbb{E}[\frac{\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}|^{2}}{L^{2}}]\end{bmatrix},$
(28)

$\tilde{{\mathbf{u}}}^{k\mathcal{B}}=\begin{bmatrix}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]\\\
\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
\mathbf{0}\end{bmatrix},$ (29)

$J_{\alpha}=\begin{bmatrix}\frac{1+\sigma^{2}}{2}&0&\frac{2\alpha^{2}L^{2}}{1-\sigma^{2}}\\\
\frac{2L^{2}\alpha}{\mu}&1-\frac{\mu\alpha}{2}&0\\\
\frac{120}{1-\sigma^{2}}&\frac{89}{1-\sigma^{2}}&\frac{3+\sigma^{2}}{4}\end{bmatrix},$
(30)

$H_{\alpha}=\begin{bmatrix}0&0&0\\\
\frac{4L^{2}\alpha^{2}}{n}&\frac{4L^{2}\alpha^{2}}{n}&0\\\
\frac{38}{1-\sigma^{2}}&\frac{38}{1-\sigma^{2}}&0\end{bmatrix}.$ (31)

If $0\leq\alpha\leq\frac{\mu(1-\sigma^{2})}{14\sqrt{2}L^{2}}$, then for any
$k\geq 0$ it holds that

${\mathbf{u}}^{(k+1)\mathcal{B}}\leq
J_{\alpha}{\mathbf{u}}^{k\mathcal{B}}+H_{\alpha}\tilde{{\mathbf{u}}}^{k\mathcal{B}}.$
(32)

As a direct consequence of Proposition 1, for the iterations of the inner loop
of Algorithm 2, for all $k\in[s\lfloor T/\mathcal{B}\rfloor,(s+1)\lfloor
T/\mathcal{B}\rfloor-1]$ it holds that

${\mathbf{u}}^{(k+1)\mathcal{B}}\leq
J_{\alpha}{\mathbf{u}}^{k\mathcal{B}}+H_{\alpha}{\mathbf{u}}^{sT},$ (33)

while for the outer loop of Algorithm 2 it holds for all $\forall s\geq 0$,

${\mathbf{u}}^{(s+1)T}\leq(J_{\alpha}^{T}+\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}){\mathbf{u}}^{sT}.$
(34)

Now, to guarantee linear decay of the outer loop sequence, we restrict the
range of the inner loop duration $T$ and the step size $\alpha$ according to

$\rho(J_{\alpha}^{T}+\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha})<1,$ (35)

where $\rho(\cdot)$ denotes the spectral radius of its argument.

In Lemma 5 below, we establish the range of $\alpha$ such that the weighted
matrix norms of $J_{\alpha}^{T}$ and
$\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}$ are small, thereby ensuring the
geometric convergence of the powers of these matrices to $\mathbf{0}$.

######  Lemma 5.

Suppose Assumption 1 holds and let
$0<\alpha\leq\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L}$, where
$\tilde{Q}=\frac{L}{\mu}$. Then

$\rho(J_{\alpha})<\||J_{\alpha}|\|^{\mathbf{\delta}}_{\infty}<1-\frac{\mu\alpha}{4}$
(36)

and

$\||\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}|\|^{{\mathbf{q}}}_{\infty}\leq\||(I-J_{\alpha})^{-1}H_{\alpha}|\|^{{\mathbf{q}}}_{\infty}<0.66,$
(37)

where
$\mathbf{\delta}=\begin{bmatrix}1,8\tilde{Q}^{2},\frac{6656\tilde{Q}^{2}}{(1-\sigma^{2})^{2}}\end{bmatrix}$
and ${\mathbf{q}}=[1,1,\frac{1457}{(1-\sigma^{2})^{2}}]$.

Essentially, Lemma 5 establishes the range of the stepsize $\alpha$ such that
the matrices involved in the system of linear inequalities (34) have small
norm. Upon setting $\alpha=\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L}$ (i.e.,
assigning $\alpha$ the largest feasible value), we proceed to determine the
number of iterations in the outer loop such that the powers of matrix
$J_{\alpha}^{T}+\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}$ in (34), and hence
the components of ${\mathbf{u}}$ (i.e. the error terms), converge to zero at a
geometric rate.

To this end, note that since $J_{\alpha}$ is non-negative,
$\sum_{l=0}^{T-1}J_{\alpha}^{l}\leq\sum_{l=0}^{\infty}J_{\alpha}^{l}=(I-J_{\alpha})^{-1}$
. Hence, for all $s\geq 0$ it holds

${\mathbf{u}}^{(s+1)T}\leq(J_{\alpha}^{T}+(I-J_{\alpha})^{-1}H_{\alpha}){\mathbf{u}}^{sT}.$
(38)

Since $\alpha=\frac{(1-\sigma^{2})^{2}}{187QL}$, we may write

$\displaystyle\begin{split}\|{\mathbf{u}}^{(s+1)T}\|_{\infty}^{{\mathbf{q}}}&\leq\||J_{\alpha}^{T}+(I-J_{\alpha})^{-1}H_{\alpha}|\|_{\infty}^{{\mathbf{q}}}\|{\mathbf{u}}^{sT}\|_{\infty}^{{\mathbf{q}}}\\\
&\leq(\||J_{\alpha}^{T}|\|_{\infty}^{{\mathbf{q}}}+0.66)\|{\mathbf{u}}^{sT}\|_{\infty}^{{\mathbf{q}}}\\\
&\leq(8\tilde{Q}^{2}(\||J_{\alpha}^{T}|\|_{\infty}^{\delta})^{T}+0.66)\|{\mathbf{u}}^{sT}\|_{\infty}^{{\mathbf{q}}}\\\
&\leq(8\tilde{Q}^{2}\exp{(-\frac{(1-\sigma^{2})^{2}T}{748\tilde{Q}^{2}})}+0.66)\|{\mathbf{u}}^{sT}\|_{\infty}^{{\mathbf{q}}}\\\
&:=\lambda\|{\mathbf{u}}^{sT}\|_{\infty}^{{\mathbf{q}}}.\end{split}$ (39)

The result in (21) follows simply by noting the definitions of
${\mathbf{u}}^{sT}$ and the $\|\cdot\|^{\mathbf{q}}_{\infty}$ norm. Therefore,
the proof of Theorem 1 is complete.

##  5 Experiments

In this section, we report results of benchmarking the proposed Algorithm 2;
for convenience, we refer to Algorithm 2 as Di-CS-SVRG (Directed
Communication-Sparsified Stochastic Variance-Reduced Gradient Descent). The
results for the proposed Algorithm 1 are presented in the supplementary
material.

We start with a network consisting of $10$ nodes with randomly generated time-
varying connections while ensuring strong connectivity at each time step. The
construction begins with the Erdős–Rényi model [40] where each edge exists
independently with probability $0.9$; then, $2$ directed edges are dropped
from each strongly connected graph, leading to directed graphs. Building upon
this basic structure, we can design networks with different connectivity
profiles. Recall that the window size parameter $\mathcal{B}$, introduced in
Assumption 1(a), may imply that the union graph over $\mathcal{B}$ consecutive
instances, starting from any instance that is a multiple of $\mathcal{B}$,
forms an almost-surely strongly connected Erdős–Rényi graph. When
$\mathcal{B}=1$, the network is strongly connected at each time step. The
parameter $q$, the fraction of entries being communicated to neighboring
nodes, characterizes the level of message sparsification; $q=1$ implies
communication without compression, while $q=0$ indicates there is no
communication in the network.

###  5.1 Decentralized Optimization Problem

We test the proposed Di-CS-SVRG on two tasks, linear and logistic regression,
and benchmark it against several baseline algorithms. In particular, we
compare Di-CS-SVRG with Decentralized Stochastic Gradient Descent (De-Stoc,
stochastic variant), Push-DIGing (Push-DIG-Full) [7], Push-DIGing Stochastic
(Push-DIG-Stoc, stochastic variant) and the TV-AB algorithm (AB Algorithm)
[25]. In addition, we include comparisons to the full gradient schemes we
considered in our preliminary work [1], Decentralized Full Gradient Descent
(De-Full) and its communication-sparsified variant Sparsified De-Full (S-De-
Full).

####  5.1.1 Decentralized Linear Regression

We consider the setting where $n$ nodes collaboratively solve the optimization
problem

$\min_{{\mathbf{x}}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\|\mathbf{y}_{i}-D_{i}\mathbf{x}\|^{2}\right\\},$
(40)

where for the $i^{th}$ node $D_{i}\in\mathbb{R}^{200\times 64}$ denotes the
matrix with $200$ local samples of size $d=64$, and
$\mathbf{y}_{i}\in\mathbb{R}^{200}$ denotes the corresponding measurement
vector. The true value of $\mathbf{x}$, $\mathbf{x}^{*}$, is generated from a
normal distribution, and the samples are generated synthetically. The
measurements are generated as $\mathbf{y}_{i}=M_{i}\mathbf{x}^{*}+\eta_{i}$,
where the entries of $M_{i}$ are generated randomly from the standard normal
distribution and then $M_{i}$ is normalized such that its rows sum up to one.
The local noise vector $\eta_{i}$ is drawn at random from a zero-mean Gaussian
distribution with variance $0.01$. All algorithms are initialized using
randomly generated local vectors $\mathbf{x}_{i}^{0}$, and utilize constant
step size $\alpha_{t}=0.002$ except the AB algorithm for which we followed the
recommendation in [25] and set $\alpha_{t}=0.0025$.

Performance of the algorithms is characterized using three metrics: residual
over iterations, residual over average gradient computation and residual over
communication cost, where the residual is computed as
$\frac{\|\mathbf{x}^{t}-\mathbf{x}^{*}\|}{\|\mathbf{x}^{0}-\mathbf{x}^{*}\|}$.
The results are shown in Fig. 1. As seen in Fig. 1(a), Di-CS-SVRG (i.e., our
Algorithm 2) with $q=1$ converges at a linear rate and, while being a
stochastic gradient algorithm, reaches the same residual floor as the full
gradient method Push-DIG-Full and the AB algorithm. Di-CS-SVRG converges much
faster than the two baseline algorithms employing SGD, Push-DIG-Stoc and De-
Stoc. Fig. 1(b) shows that Di-CS-SVRG with varied compression levels $q$
converges to the same residual floor, and that (as expected) larger $q$ leads
to faster convergence. Moreover, the figure shows that for a fixed $q$, Di-CS-
SVRG achieves faster convergence than the benchmark algorithms.

Fig. 1(c) compares different algorithms in terms of the number of gradients
computed per sample, demonstrating computation efficiency of Di-CS-SVRG.
Finally, Fig. 1(d) shows the communication cost for varied $q$, computed as
the total number of the (state, surplus and gradient) vector entries
communicated across the network. As seen in the figure, to achieve a pre-
specified level of the residual, Di-CS-SVRG with $q=0.05$ incurs smaller
communication cost than any other considered algorithm.

(a) Residual vs. iterations: full communication schemes.

(b) Residual vs. iterations: compressed communication schemes.

(c) Residual vs. gradient computations per sample.

(d) Residual vs. communication cost.

Figure 1: Linear regression, $\mathcal{B}=5$. (a) The residual achieved by
full communication schemes and the residual of Di-CS-SVRG (Algorithm 2) with
$q=1$ vs. iterations. (b) The residual achieved by Di-CS-SVRG with different
compression levels: $q=1$, $q=0.08$, and $q=0.05$ vs. iterations. (c) The
cumulative number of gradient computations needed to reach given level of the
residual; showing both the compressed as well as full communication schemes.
(d) The cumulative communication cost to reach given level of the residual;
showing both the compressed as well as full communication schemes.

(a) Correct rate vs. iterations: full communication schemes.

(b) Correct rate vs. iterations: compressed schemes.

(c) Correct rate vs. gradient computations per sample.

(d) Correct rate vs. communication cost.

Figure 2: Logistic regression, $\mathcal{B}=1$. (a) The correct classification
rate achieved by full communication schemes and the correct classification
rate of Di-CS-SVRG vs. iterations. (b) The correct classification rate for Di-
CS-SVRG with varied compression levels: $q=1$, $q=0.12$, $q=0.08$ vs.
iterations. (c) The cumulative number of gradient computations required to
reach given levels of correct classification rate. (d) The cumulative
communication cost required to reach given level of the correct classification
rate.

####  5.1.2 Decentralized Logistic Regression

To perform benchmarking on a logistic regression task, we solve a multi-class
classification problem on the Stackoverflow dataset [41].

$\min_{{\mathbf{x}}}\left\\{\frac{\mu}{2}\|\mathbf{x}\|^{2}+\sum_{i=1}^{n}\sum_{j=1}^{N}\mathrm{ln}(1+\mathrm{exp}(-(\mathbf{m}_{ij}^{T}\mathbf{x}){\mathbf{y}}_{ij}))\right\\},$
(41)

where the training samples
$(\mathbf{m}_{ij},{\mathbf{y}}_{ij})\in\mathbb{R}^{400+5}$, $\mathbf{m}_{ij}$
represents a vectorized text feature and ${\mathbf{y}}_{ij}$ represents the
corresponding tag vector. We compare the performance of Di-CS-SVRG with the
same benchmarking algorithms as in the linear regression problem, and use the
same initialization setup. The logistic regression experiment is run with the
stepsize $\alpha_{t}=0.01$; the regularization parameter is set to
$\mu=10^{-5}$.

Performance of the algorithms on the logistic regression task is characterized
by the classification correct rate. In particular, we evaluate the following
three metrics: the correct rate vs. iterations, the correct rate vs. average
gradient computation, and the correct rate vs. communication cost; they are
all shown in Fig. 2. As seen in Fig. 2(a), Di-CS-SVRG converges and reaches
the same residual floor as the full gradient method Push-DIG-Full and the AB
Algorithm. Di-CS-SVRG converges much faster than then the two algorithms that
rely on SGD, Push-DIG-Stoc and De-Stoc. Fig. 2(b) shows that for varied
compression levels $q$, Di-CS-SVRG converges to the same residual floor. As
expected, larger $q$ leads to faster convergence. For a fixed $q$, Di-CS-SVRG
converges faster than the benchmark algorithm.

Fig. 2(c) reports the average gradient computation, i.e., the number of
gradients computed per sample. As can be seen, Di-CS-SVRG with different
compression levels uses fewer gradient computation than the full gradient
schemes (Push-DIG-Full, the AB algorithm and De-Full) to reach $90\%$ correct
classification rate.

Fig. 2(d) shows the communication cost, defined as the total number of the
(state, surplus and gradient) vector entries exchanged across the network, for
various values of $q$. Among the considered schemes, Di-CS-SVRG with $q=0.08$
reaches the pre-specified residual level with smaller communication cost than
any stochastic scheme (Push-DIG-Stoc, De-Stoc as well as Di-CS-SVRG with other
$q$’s). Even though Push-DIG-Stoc reaches various accuracies slightly faster
than Di-CS-SVRG, it requires considerably higher amount of gradient
computation to do so (see Fig. 2).

###  5.2 Results on different network topologies

To further test Di-CS-SVRG in different settings, we apply it to decentralized
optimization over networks with varied connectivity and sizes.

####  5.2.1 Varied network connectivity

We consider the linear regression problem and vary the values of the joint
connectivity parameter $\mathcal{B}$. Fig. 3(a) shows the resulting residuals;
as seen there, larger $\mathcal{B}$, implying the network takes longer time
before the union of its instances forms a strongly connected graph, leads to
slower convergence.

(a)

(b)

Figure 3: Varied network connectivity and size.

####  5.2.2 Varied network size

We now consider the logistic regression problem over networks of varied sizes.
In particular, we fix the total number of data points and vary the number of
nodes in the network. As the network grows, i.e., the number of nodes in the
network becomes larger, each agent has fewer locally available data points.
Fig. 3(b) shows the correct rate for compression levels $q=1$, $q=0.12$ and
$q=0.08$ as $n$ grows from $10$ to $30$. For a pre-specified sparsification
level, larger networks, in which each agent has fewer data points to train its
local model, requires more communication rounds and therefore takes longer to
converge.

##  6 Conclusion

In this paper we studied decentralized convex optimization problems over time-
varying directed networks and proposed a stochastic variance-reduced algorithm
for solving them. The algorithm sparsifies messages exchanged between network
nodes thus enabling collaboration in resource-constrained settings. We proved
that the proposed algorithm, Di-CS-SVRG, enjoys linear convergence rate, and
demonstrated its efficacy through simulation studies. As part of the future
work, it is of interest to extend this work to decentralized non-convex
optimization problems.

## References

  * [1] Yiyue Chen, Abolfazl Hashemi and Haris Vikalo  “Decentralized Optimization on Time-Varying Directed Graphs Under Communication Constraints”  In _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2021, pp. 3670–3674  IEEE
  * [2] Wei Ren and Randal W Beard  “Consensus seeking in multiagent systems under dynamically changing interaction topologies”  In _IEEE Transactions on automatic control_ 50.5 IEEE, 2005, pp. 655–661
  * [3] Ali Jadbabaie, Jie Lin and A Stephen Morse  “Coordination of groups of mobile autonomous agents using nearest neighbor rules”  In _IEEE Transactions on automatic control_ 48.6 IEEE, 2003, pp. 988–1001
  * [4] Angelia Nedic and Asuman Ozdaglar  “Distributed subgradient methods for multi-agent optimization”  In _IEEE Transactions on Automatic Control_ 54.1 IEEE, 2009, pp. 48–61
  * [5] S Sundhar Ram, Angelia Nedić and Venugopal V Veeravalli  “Distributed stochastic subgradient projection algorithms for convex optimization”  In _Journal of optimization theory and applications_ 147.3 Springer, 2010, pp. 516–545
  * [6] Guannan Qu and Na Li  “Harnessing smoothness to accelerate distributed optimization”  In _IEEE Transactions on Control of Network Systems_ 5.3 IEEE, 2017, pp. 1245–1260
  * [7] Angelia Nedic, Alex Olshevsky and Wei Shi  “Achieving geometric convergence for distributed optimization over time-varying graphs”  In _SIAM Journal on Optimization_ 27.4 SIAM, 2017, pp. 2597–2633
  * [8] Rie Johnson and Tong Zhang  “Accelerating stochastic gradient descent using predictive variance reduction”  In _Advances in neural information processing systems_ , 2013, pp. 315–323
  * [9] Mark Schmidt, Nicolas Le Roux and Francis Bach  “Minimizing finite sums with the stochastic average gradient”  In _Mathematical Programming_ 162.1-2 Springer, 2017, pp. 83–112
  * [10] John Nikolas Tsitsiklis  “Problems in decentralized decision making and computation.”, 1984
  * [11] Wei Ren, Randal W Beard and Ella M Atkins  “Information consensus in multivehicle cooperative control”  In _IEEE Control systems magazine_ 27.2 IEEE, 2007, pp. 71–82
  * [12] Kai Cai and Hideaki Ishii  “Average consensus on general strongly connected digraphs”  In _Automatica_ 48.11 Elsevier, 2012, pp. 2750–2761
  * [13] Kai Cai and Hideaki Ishii  “Average consensus on arbitrary strongly connected digraphs with time-varying topologies”  In _IEEE Transactions on Automatic Control_ 59.4 IEEE, 2014, pp. 1066–1071
  * [14] Anastasia Koloskova, Sebastian Stich and Martin Jaggi  “Decentralized Stochastic Optimization and Gossip Algorithms with Compressed Communication”  In _International Conference on Machine Learning_ , 2019, pp. 3478–3487
  * [15] Björn Johansson, Maben Rabi and Mikael Johansson  “A randomized incremental subgradient method for distributed optimization in networked systems”  In _SIAM Journal on Optimization_ 20.3 SIAM, 2010, pp. 1157–1170
  * [16] Ermin Wei and Asuman Ozdaglar  “Distributed alternating direction method of multipliers”  In _2012 IEEE 51st IEEE Conference on Decision and Control (CDC)_ , 2012, pp. 5445–5450  IEEE
  * [17] John C Duchi, Alekh Agarwal and Martin J Wainwright  “Dual averaging for distributed optimization: Convergence analysis and network scaling”  In _IEEE Transactions on Automatic control_ 57.3 IEEE, 2011, pp. 592–606
  * [18] Angelia Nedić, Soomin Lee and Maxim Raginsky  “Decentralized online optimization with global objectives and local communication”  In _2015 American Control Conference (ACC)_ , 2015, pp. 4497–4503  IEEE
  * [19] Lie He, An Bian and Martin Jaggi  “Cola: Decentralized linear learning”  In _Advances in Neural Information Processing Systems_ , 2018, pp. 4536–4546
  * [20] Sebastian U Stich, Jean-Baptiste Cordonnier and Martin Jaggi  “Sparsified SGD with memory”  In _Advances in Neural Information Processing Systems_ , 2018, pp. 4447–4458
  * [21] David Kempe, Alin Dobra and Johannes Gehrke  “Gossip-based computation of aggregate information”  In _44th Annual IEEE Symposium on Foundations of Computer Science, 2003. Proceedings._ , 2003, pp. 482–491  IEEE
  * [22] Angelia Nedić and Alex Olshevsky  “Distributed optimization over time-varying directed graphs”  In _IEEE Transactions on Automatic Control_ 60.3 IEEE, 2014, pp. 601–615
  * [23] Chenguang Xi, Qiong Wu and Usman A Khan  “On the distributed optimization over directed networks”  In _Neurocomputing_ 267 Elsevier, 2017, pp. 508–515
  * [24] Angelia Nedić and Alex Olshevsky  “Stochastic gradient-push for strongly convex functions on time-varying directed graphs”  In _IEEE Transactions on Automatic Control_ 61.12 IEEE, 2016, pp. 3936–3947
  * [25] Fakhteh Saadatniaki, Ran Xin and Usman A Khan  “Optimization over time-varying directed graphs with row and column-stochastic matrices”  In _arXiv preprint arXiv:1810.07393_ , 2018
  * [26] Ran Xin, Usman A Khan and Soummya Kar  “Variance-reduced decentralized stochastic optimization with accelerated convergence”  In _arXiv preprint arXiv:1912.04230_ , 2019
  * [27] Boyue Li, Shicong Cen, Yuxin Chen and Yuejie Chi  “Communication-efficient distributed optimization in networks with gradient tracking and variance reduction”  In _International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 1662–1672  PMLR
  * [28] Bahman Gharesifard and Jorge Cortés  “When does a digraph admit a doubly stochastic adjacency matrix?”  In _Proceedings of the 2010 American Control Conference_ , 2010, pp. 2440–2445  IEEE
  * [29] Muhammad I Qureshi, Ran Xin, Soummya Kar and Usman A Khan  “Push-SAGA: A decentralized stochastic algorithm with variance reduction over directed graphs”  In _IEEE Control Systems Letters_ IEEE, 2021
  * [30] Hanlin Tang et al.  “Communication compression for decentralized training”  In _Advances in Neural Information Processing Systems_ , 2018, pp. 7652–7662
  * [31] Wei Wen et al.  “Terngrad: Ternary gradients to reduce communication in distributed deep learning”  In _Advances in neural information processing systems_ , 2017, pp. 1509–1519
  * [32] Hantian Zhang et al.  “Zipml: Training linear models with end-to-end low precision, and a little bit of deep learning”  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , 2017, pp. 4035–4043  JMLR. org
  * [33] Rudrajit Das, Abolfazl Hashemi, Sujay Sanghavi and Inderjit S Dhillon  “Improved Convergence Rates for Non-Convex Federated Learning with Compression”  In _arXiv preprint arXiv:2012.04061_ , 2020
  * [34] Amirhossein Reisizadeh et al.  “Robust and communication-efficient collaborative learning”  In _Advances in Neural Information Processing Systems_ , 2019, pp. 8386–8397
  * [35] Zebang Shen et al.  “Towards More Efficient Stochastic Decentralized Learning: Faster Convergence and Sparse Communication”  In _International Conference on Machine Learning_ , 2018, pp. 4624–4633
  * [36] Anastasia Koloskova, Tao Lin, Sebastian U Stich and Martin Jaggi  “Decentralized deep learning with arbitrary communication compression”  In _arXiv preprint arXiv:1907.09356_ , 2019
  * [37] Abolfazl Hashemi et al.  “On the Benefits of Multiple Gossip Steps in Communication-Constrained Decentralized Optimization”  In _arXiv preprint arXiv:2012.04061_ , 2020
  * [38] Hossein Taheri, Aryan Mokhtari, Hamed Hassani and Ramtin Pedarsani  “Quantized Decentralized Stochastic Learning over Directed Graphs”  In _International Conference on Machine Learning (ICML)_ , 2020
  * [39] Lin Xiao and Stephen Boyd  “Fast linear iterations for distributed averaging”  In _Systems & Control Letters_ 53.1 Elsevier, 2004, pp. 65–78
  * [40] Paul Erdös and Alfréd Rényi  “On random graphs”  In _Publicationes mathematicae_ 6.26, 1959, pp. 290–297
  * [41] “The Stack Overflow Data”  URL: https://www.kaggle.com/stackoverflow/stackoverflow
  * [42] Ran Xin, Anit Kumar Sahu, Usman A Khan and Soummya Kar  “Distributed stochastic optimization with gradient tracking over strongly-connected networks”  In _arXiv preprint arXiv:1903.07266_ , 2019

Supplementary Material

In this document, we include detailed proofs of the auxiliary lemmas stated in
Section 1 of the main manuscript, and present experimental results for the
consensus algorithm, Algorithm 1.

##  Appendix A Analysis

We start by proving auxiliary lemmas utilized in the proof of Theorem 1.

######  Lemma $1$.

Suppose Assumptions 1 (a) and (b) hold. Let
$\sigma=\max(|\lambda_{M,2}|,|\lambda_{B,2}|)$ denote the larger of the second
largest eigenvalues of $M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$ and
$B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$. Then,

$\displaystyle\begin{split}\|M_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{z}}-\bar{{\mathbf{z}}}\|\leq\sigma\|{\mathbf{z}}-\bar{{\mathbf{z}}}\|,\
\forall{\mathbf{z}}\in{\mathbb{R}}^{2n}\\\ \mbox{and}\\\
\|B_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{y}}-\bar{{\mathbf{y}}}\|\leq\sigma\|{\mathbf{y}}-\bar{{\mathbf{y}}}\|,\
\forall{\mathbf{y}}\in{\mathbb{R}}^{n},\\\ \end{split}$ (42)

where
$\bar{{\mathbf{z}}}=[\frac{1}{n}\sum_{i=1}^{2n}z_{i},\cdots,\frac{1}{n}\sum_{i=1}^{2n}z_{i}]^{T}$
and
$\bar{{\mathbf{y}}}=[\frac{1}{n}\sum_{i=1}^{n}y_{i},\cdots,\frac{1}{n}\sum_{i=1}^{n}y_{i}]^{T}$.

###### Proof.

To prove Lemma 1, we first need to establish the following.

######  Lemma $1.1$.

Assume that $M_{m}((s+1)\mathcal{B}-1:s\mathcal{B})$ has non-zero spectral gap
for each $m$. Then the following statements hold:

  1. (a)

The sequence of matrix products $M_{m}((s+1)\mathcal{B}-1:s\mathcal{B})$
converges to the limit matrix

$\mathrm{lim}_{t\to\infty}(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))^{t}=\left[\begin{matrix}\frac{\mathbf{1}_{n}\mathbf{1}^{T}_{n}}{n}&\frac{\mathbf{1}_{n}\mathbf{1}^{T}_{n}}{n}\\\
\mathbf{0}&\mathbf{0}\\\ \end{matrix}\right].$ (43)

  2. (b)

Let
$1=|\lambda_{1}(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))|>|\lambda_{2}(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))|\geq\cdots\geq|\lambda_{2n}(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))|$
be the eigenvalues of $M_{m}((s+1)\mathcal{B}-1:s\mathcal{B})$, and let
$\sigma_{m}=|\lambda_{2}(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))|$; then there
exists $\Gamma_{m}^{\prime}>0$ such that

$\displaystyle\|(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))^{t}-\mathcal{I}\|_{\infty}\leq\Gamma_{m}^{\prime}\sigma_{m}^{t},$
(44)

where $\mathcal{I}:=\frac{1}{n}[\mathbf{1}^{T}\
\mathbf{0}^{T}]^{T}[\mathbf{1}^{T}\ \mathbf{1}^{T}]$.

###### Proof.

For each $m$, $M_{m}((s+1)\mathcal{B}-1:s\mathcal{B})$ has column sum equal to
$1$. According to Assumption 1, definition of the mixing matrix and the
construction of the product, $M_{m}((s+1)\mathcal{B}-1:s\mathcal{B})$ has a
simple eigenvalue $1$ with the corresponding left eigenvector
$[\mathbf{1}^{T}\ \mathbf{1}^{T}]$ and right eigenvector $[\mathbf{1}^{T}\
\mathbf{0}^{T}]^{T}$. Following Jordan matrix decomposition for the simple
eigenvalue, there exist some $P,Q\in\mathcal{R}^{(2n-1)\times(2n-1)}$ such
that

$\displaystyle(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))^{t}$
$\displaystyle=\mathcal{I}^{t}+PJ_{m}^{t}Q=\mathcal{I}+PJ_{m}^{t}Q.$ (45)

Let $\gamma_{m}$ be the second largest eigenvalue magnitude of
$M_{m}((s+1)\mathcal{B}-1:s\mathcal{B})$; then, $\gamma_{m}$ is also the
spectral norm of $J_{m}$. The proof of part (a) follows by noting that
$\mathrm{lim}_{t\to\infty}J_{m}^{t}=\mathbf{0}.$ Since $\|P\|$, $\|Q\|$ and
$\|J_{m}\|$ are finite, there exists some $\Gamma_{m}^{\prime}>0$ such that

$\displaystyle\quad\|(M_{m}((s+1)\mathcal{B}-1:s\mathcal{B}))^{t}-\mathcal{I}\|_{\infty}\leq\|PJ_{m}^{t}Q\|_{\infty}\leq\Gamma_{m}^{\prime}\sigma_{m}^{t}$
(46)

which completes the proof of part (b). ∎

Then let $\sigma^{\prime}=\max_{m}\sigma_{m}$, where $\sigma_{m}$ is as
defined in Lemma Lemma $1.1$ and, by mathematical induction, for each $m$ it
holds that

$\rho(M_{m}(T\mathcal{B}-1:0)-\frac{1}{n}[\mathbf{1}^{T}\
\mathbf{0}^{T}]^{T}[\mathbf{1}^{T}\ \mathbf{1}^{T}])\leq\sigma^{\prime T}.$
(47)

Referring to the fact that $M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$ and
$B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$ both have column sums equal to $1$
and defining $\sigma$ as stated in the lemma, we can conclude the proof of
Lemma 1.

∎

The following lemma, restated for convenience, establishes an upper bound on
the consensus error.

######  Lemma $2$.

Suppose Assumption 1 holds. Then, $\forall i\leq n$, $k\geq 0$, and $0<m\leq
d$, the updates generated by Algorithm 1 satisfy

$\displaystyle\mathbb{E}[|z_{im}^{(k+1)\mathcal{B}}-\bar{z}_{m}^{(k+1)\mathcal{B}}|^{2}]$
$\displaystyle\leq\frac{1+\sigma^{2}}{2}\mathbb{E}[|z_{im}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}|^{2}]$
(48)
$\displaystyle\quad+\frac{2\alpha^{2}}{1-\sigma^{2}}\mathbb{E}[|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\|^{2}].$

###### Proof.

By constructing normalized weight matrices and relying on the definition of
$M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})$, we can simplify the update as

$\displaystyle z_{im}^{(k+1)\mathcal{B}}$
$\displaystyle=\sum_{j=1}^{2n}[M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{ij}[Q(z_{j}^{k\mathcal{B}})]_{m}-\alpha
g_{im}^{k\mathcal{B}}$ (49)
$\displaystyle=\sum_{j=1}^{2n}[M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{ij}z_{jm}^{k\mathcal{B}}-\alpha
g_{im}^{k\mathcal{B}}.$

Since for any $m$ we have
$\bar{g}_{m}^{k\mathcal{B}}=\frac{1}{n}\sum_{j=1}^{n}g_{jm}^{k\mathcal{B}}$,
it holds that

$\displaystyle|z_{im}^{(k+1)\mathcal{B}}-\bar{z}_{m}^{(k+1)\mathcal{B}}|^{2}=|\sum_{j=1}^{2n}[M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{ij}z_{jm}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}-\alpha(g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}})|^{2}.$
(50)

By Young’s inequality,

$\|{\mathbf{a}}+{\mathbf{b}}\|^{2}\leq(1+\eta)\|{\mathbf{a}}\|^{2}+(1+\frac{1}{\eta})\|{\mathbf{b}}\|^{2},\quad\forall\eta>0,{\mathbf{a}},{\mathbf{b}}.$

Using Lemma 1, $\forall i\leq n$ and $0<m\leq d$ it holds that

$|\sum_{j=1}^{2n}[M_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{ij}z_{jm}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}|\leq\sigma|z_{im}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}|.$
(51)

Then for all $k\geq 0$,

$\displaystyle|z_{im}^{(k+1)\mathcal{B}}-\bar{z}_{m}^{(k+1)\mathcal{B}}|^{2}$
$\displaystyle\leq(1+\eta)\sigma^{2}|z_{im}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}|^{2}+(1+\frac{1}{\eta})\alpha^{2}|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}|^{2}.$
(52)

Setting $\eta=\frac{1-\sigma^{2}}{2\sigma^{2}}$ completes the proof.

∎

Next, we provide proofs of two lemmas stating upper bounds on the optimality
gap and the gradient tracking error. For convenience, we first introduce

$\bar{\tau}^{k\mathcal{B}}=\frac{1}{n}\sum_{i=1}^{n}\tau_{i}^{k\mathcal{B}}$
(53)

and

$\tau_{i}^{(k+1)\mathcal{B}}=\begin{cases}{\mathbf{x}}_{i}^{(k+1)\mathcal{B}}&\mathrm{if}\quad(k+1)\mathcal{B}\mod
T=0\\\ \tilde{w}_{i}&\mathrm{otherwise}.\end{cases}$ (54)

######  Lemma $3$.

Suppose Assumption 1 holds and let $0<\alpha<\frac{\mu}{8L^{2}}$. Then for all
$k>0$ it holds that

$\displaystyle\begin{split}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]&\leq\frac{2L^{2}\alpha}{\mu}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+(1-\frac{\mu\alpha}{2})\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].\end{split}$
(55)

######  Lemma $4$.

Suppose Assumption 1 holds. Then,

$\displaystyle\frac{1}{L^{2}}\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}]$
$\displaystyle\leq\frac{120}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+\frac{89}{1-\sigma^{2}}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$
(56)
$\displaystyle\quad+\frac{3+\sigma^{2}}{4}\mathbb{E}[\frac{\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}|^{2}}{L^{2}}]$
$\displaystyle\quad+\frac{38}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+\frac{38}{1-\sigma^{2}}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].$

Proving Lemmas 3 and 4 requires a series of auxiliary lemmas, Lemma 3.1-3.4.
We start with Lemma 3.1, which states an upper bound on
$\mathbb{E}[\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$.

######  Lemma $3.1$.

Suppose Assumption 1 holds. Let $0<\alpha<\frac{1}{L}$, where $L$ is the
smoothness parameter. For all $k>0$, it holds that

$\displaystyle\mathbb{E}[\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$
$\displaystyle\leq\frac{L^{2}\alpha}{n\mu}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+(1-\mu\alpha)\mathbb{E}[\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]+\frac{\alpha^{2}}{n^{2}}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla\mathbf{f}({\mathbf{x}}^{k\mathcal{B}})\|^{2}],$
(57)

where $\mu$ is the strong convexity parameter and
$\nabla\mathbf{f}({\mathbf{x}}^{k\mathcal{B}})=[\nabla
f_{1}({\mathbf{x}}_{1}^{k\mathcal{B}});\cdots;\nabla
f_{n}({\mathbf{x}}_{n}^{k\mathcal{B}})]$.

###### Proof.

By definition, $\bar{z}_{m}^{t}=\frac{1}{n}\sum_{i=1}^{2n}z_{im}^{t}$. Let us
denote
$\nabla\bar{f}({\mathbf{x}}^{t})=\nabla\bar{f}({\mathbf{z}}^{t})=\frac{1}{n}\sum_{i=1}^{n}\nabla
f_{i}({\mathbf{z}}_{i}^{t})$. By induction,

$\bar{g}_{m}^{k\mathcal{B}}=\bar{v}_{m}^{k\mathcal{B}}.$ (58)

Next, we have that for any $m$

$\bar{z}_{m}^{(k+1)\mathcal{B}}=\bar{z}_{m}^{k\mathcal{B}}-\alpha\bar{g}_{m}^{k\mathcal{B}}=\bar{z}_{m}^{k\mathcal{B}}-\alpha\bar{v}_{m}^{k\mathcal{B}},$
(59)

which implies that

$\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}=\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\bar{{\mathbf{g}}}^{k\mathcal{B}}=\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\bar{{\mathbf{v}}}^{k\mathcal{B}}.$
(60)

Note that the randomness in Algorithm 1 originates from a set of independent
random variables $\\{\omega_{i}^{t}\\}_{i\in[2n]}^{t\geq 0}$. We rely on the
$\sigma$-algebra ${\mathcal{F}}^{k\mathcal{B}}$ to characterize the history of
the dynamical system generated by $\\{\omega_{i}^{t}\\}_{i\in[2n]}^{t\leq
k\mathcal{B}-1}$,

$\displaystyle\begin{split}&\mathbb{E}[\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]=\mathbb{E}[\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\bar{{\mathbf{v}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&=\mathbb{E}[\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-{\mathbf{x}}^{*}+\alpha(\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}})\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&=\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-{\mathbf{x}}^{*}\|^{2}+\alpha^{2}\mathbb{E}[\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&\quad+2\alpha\langle\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\nabla
f(\bar{z}^{k\mathcal{B}})-{\mathbf{x}}^{*},\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\rangle\\\
&=\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-{\mathbf{x}}^{*}\|^{2}+\alpha^{2}\mathbb{E}[\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&\quad+2\alpha\langle\bar{{\mathbf{z}}}^{k\mathcal{B}}-\alpha\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-{\mathbf{x}}^{*},\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{\mathcal{B}})\rangle.\end{split}$
(61)

We then proceed by considering $\mathbb{E}[\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]$,

$\displaystyle\begin{split}\mathbb{E}[\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]&=\mathbb{E}[\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})+\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&=\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})\|^{2}+\mathbb{E}[\|\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}],\end{split}$
(62)

where the fact that
$\mathbb{E}[\bar{{\mathbf{v}}}^{t}|{\mathcal{F}}^{t}]=\nabla\bar{f}({\mathbf{x}}^{t})$
is used. Furthermore, note that

$\displaystyle\mathbb{E}[\|\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})-\bar{{\mathbf{v}}}^{k\mathcal{B}}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}[\|\sum_{i=1}^{n}({\mathbf{v}}_{i}^{k\mathcal{B}}-\nabla
f_{i}({\mathbf{z}}_{i}^{k\mathcal{B}}))\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]$
(63)
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla\mathbf{f}({\mathbf{x}}^{k\mathcal{B}})\|^{2}|{\mathcal{F}}^{k\mathcal{B}}],$

since $\left\\{{\mathbf{v}}_{i}^{t}\right\\}_{i=1}^{n}$ are independent given
${\mathcal{F}}^{t}$ and $\mathbb{E}[\sum_{i\neq
j}\langle{\mathbf{v}}_{i}^{t}-\nabla
f_{i}({\mathbf{z}}_{i}^{t}),{\mathbf{v}}_{j}^{t}-\nabla
f_{j}({\mathbf{z}}_{j}^{t})\rangle|{\mathcal{F}}^{t}]=0$. Recall the strong
convexity of the objective, i.e., we have that if $0<\alpha\leq\frac{1}{L}$,
$\forall{\mathbf{x}}$

$\|{\mathbf{x}}-\alpha\nabla{\mathbf{f}}({\mathbf{x}})-{\mathbf{x}}^{*}\|\leq(1-\mu\alpha)\|{\mathbf{x}}-{\mathbf{x}}^{*}\|.$
(64)

It follows that

$\displaystyle\begin{split}\mathbb{E}[\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]&\leq(1-\mu\alpha)^{2}\|\bar{z}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}+\alpha^{2}\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})\|^{2}\\\
&\quad+2\alpha(1-\mu\alpha)\|\bar{z}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})\|\\\
&\quad+\frac{\alpha^{2}}{n^{2}}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}|{\mathcal{F}}^{k\mathcal{B}}].\end{split}$
(65)

Using Young’s inequality, we readily obtain that

$\displaystyle
2\alpha\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})\|\leq\mu\alpha\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}+\frac{\alpha}{\mu}\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})\|^{2}.$
(66)

On the other hand, assuming convexity and smoothness, we have that $\forall
k\geq 0$,

$\displaystyle\|\nabla
f(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla\bar{f}({\mathbf{x}}^{k\mathcal{B}})\|$
$\displaystyle=\|\sum_{i=1}^{n}\frac{\nabla
f_{i}(\bar{{\mathbf{z}}}^{k\mathcal{B}})-\nabla
f_{i}({\mathbf{z}}_{i}^{k\mathcal{B}})}{n}\|$ (67) $\displaystyle\leq
L\sum_{i=1}^{n}\frac{\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|}{n}$
(68) $\displaystyle\leq
L\sqrt{\sum_{i=1}^{n}\frac{\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}}{n}}$
(69)
$\displaystyle=\frac{L}{\sqrt{n}}\sqrt{\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}}.$
(70)

Then by taking the total expectation,

$\displaystyle\mathbb{E}[\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$
$\displaystyle\leq\frac{L^{2}\alpha}{n\mu}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+(1-\mu\alpha)\mathbb{E}[\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]+\frac{\alpha^{2}}{n^{2}}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla\mathbf{f}({\mathbf{x}}^{k\mathcal{B}})\|^{2}].$

∎

The following lemma helps establish an upper bound on the expected gradient
tracking error.

######  Lemma $3.2$.

Suppose the objective function $f$ is $\mu$-strongly-convex and that each
component of the local objective function $f_{i,j}$ is $L$-smooth. If
$0<\alpha<\frac{1}{4\sqrt{2}L}$,

$\displaystyle\begin{split}\mathbb{E}[\|{\mathbf{g}}_{\sim
n}^{(k+1)\mathcal{B}}-\mathbf{1}_{n}\bar{{\mathbf{g}}}^{(k+1)\mathcal{B}}\|^{2}]&\leq\frac{33L^{2}}{1-\sigma^{2}}E[\|{\mathbf{z}}^{k\mathcal{B}}-\mathbf{1}_{2n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}\|^{2}]+\frac{2L^{2}}{1-\sigma^{2}}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+(\frac{1+\sigma^{2}}{2}+\frac{32\alpha^{2}L^{2}}{1-\sigma^{2}})\mathbb{E}[\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}\bar{{\mathbf{g}}}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{5}{1-\sigma^{2}}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}]\\\
&\quad+\frac{4}{1-\sigma^{2}}\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})\|^{2}],\end{split}$
(71)

where ${\mathbf{g}}_{\sim
n}^{t}=[{\mathbf{g}}_{1}^{t};\cdots;{\mathbf{g}}_{n}^{t}]\in{\mathbb{R}}^{n\times
d}.$

###### Proof.

For all $i\leq n$ and $0<m\leq d$, it holds that

$\displaystyle|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}$
$\displaystyle=|\sum_{j=1}^{n}[B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{ij}g_{jm}^{k\mathcal{B}}+v_{im}^{(k+1)\mathcal{B}}-v_{im}^{k\mathcal{B}}$
(72)
$\displaystyle\quad-\frac{1}{n}\sum_{l=1}^{n}(\sum_{j=1}^{n}[B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{lj}g_{jm}^{k\mathcal{B}}+v_{lm}^{(k+1)\mathcal{B}}-v_{lm}^{k\mathcal{B}})|^{2}.$

Denoting
${\mathbf{g}}_{:m}^{(k+1)\mathcal{B}}=[g_{1m}^{(k+1)\mathcal{B}},\cdots,g_{nm}^{(k+1)\mathcal{B}}]^{T}$
and
${\mathbf{v}}_{:m}^{(k+1)\mathcal{B}}=[v_{1m}^{(k+1)\mathcal{B}},\cdots,v_{nm}^{(k+1)\mathcal{B}}]$,

$\displaystyle\begin{split}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}&=\|{\mathbf{g}}_{:m}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}\mathbf{1}_{n}\|^{2}\\\
&=\sum_{i=1}^{n}|\sum_{j=1}^{n}[B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{ij}g_{jm}^{k\mathcal{B}}+v_{im}^{(k+1)\mathcal{B}}-v_{im}^{k\mathcal{B}}-\\\
&\quad\frac{1}{n}\sum_{l=1}^{n}(\sum_{j=1}^{n}[B_{m}((k+1)\mathcal{B}-1:k\mathcal{B})]_{lj}g_{jm}^{k\mathcal{B}}+v_{lm}^{(k+1)\mathcal{B}}-v_{lm}^{k\mathcal{B}})|^{2}\\\
&=\|B_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{g}}_{:m}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\mathbf{1}_{n}+({\mathbf{v}}_{:m}^{(k+1)\mathcal{B}}-\bar{v}_{m}^{(k+1)\mathcal{B}}\mathbf{1}_{n})-({\mathbf{v}}_{:m}^{k\mathcal{B}}-\bar{v}_{m}^{k\mathcal{B}}\mathbf{1}_{n})\|^{2}.\end{split}$
(73)

Once again applying Young’s inequality yields

$\displaystyle\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}$
$\displaystyle\leq(1+\frac{1-\sigma^{2}}{2\sigma^{2}})\|B_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{g}}_{:m}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\mathbf{1}_{n}\|^{2}$
(74)
$\displaystyle\quad+(1+\frac{2\sigma^{2}}{1-\sigma^{2}})\|({\mathbf{v}}_{:m}^{(k+1)\mathcal{B}}-\bar{v}_{m}^{(k+1)\mathcal{B}}\mathbf{1}_{n})-({\mathbf{v}}_{:m}^{k\mathcal{B}}-\bar{v}_{m}^{k\mathcal{B}}\mathbf{1}_{n})\|^{2}.$

Summing up the above objects over all $m$, we obtain

$\displaystyle\begin{split}\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}&\leq\sum_{m=1}^{d}(1+\frac{1-\sigma^{2}}{2\sigma^{2}})\|B_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{g}}_{:m}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\mathbf{1}_{n}\|^{2}\\\
&\quad+(1+\frac{2\sigma^{2}}{1-\sigma^{2}})\|({\mathbf{v}}_{:m}^{(k+1)\mathcal{B}}-\bar{v}_{m}^{(k+1)\mathcal{B}}\mathbf{1}_{n})-({\mathbf{v}}_{:m}^{k\mathcal{B}}-\bar{v}_{m}^{k\mathcal{B}}\mathbf{1}_{n})\|^{2}\\\
&\leq\sum_{m=1}^{d}\frac{1+\sigma^{2}}{2}\|{\mathbf{g}}_{:m}^{k\mathcal{B}}-\bar{g}_{m}^{t}\mathbf{1}_{n}\|^{2}\\\
&\quad+\frac{2}{1-\sigma^{2}}\|({\mathbf{v}}_{:m}^{(k+1)\mathcal{B}}-\bar{v}_{m}^{(k+1)\mathcal{B}}\mathbf{1}_{n})-({\mathbf{v}}_{:m}^{k\mathcal{B}}-\bar{v}_{m}^{k\mathcal{B}}\mathbf{1}_{n})\|^{2}\\\
&\leq\frac{1+\sigma^{2}}{2}\sum_{m=1}^{d}\|{\mathbf{g}}_{:m}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\mathbf{1}_{n}\|^{2}+\frac{2}{1-\sigma^{2}}\|{\mathbf{v}}^{(k+1)\mathcal{B}}-{\mathbf{v}}^{k\mathcal{B}}\|^{2}\end{split}$
(75)

Taking the total expectation yields

$\displaystyle\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}]$
$\displaystyle\leq\frac{1+\sigma^{2}}{2}\mathbb{E}[\sum_{m=1}^{d}\|{\mathbf{g}}_{:m}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\mathbf{1}_{n}\|^{2}]+\frac{2}{1-\sigma^{2}}\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-{\mathbf{v}}^{k\mathcal{B}}\|^{2}].$
(76)

Next, we derive an upper bound on
$\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-{\mathbf{v}}^{k\mathcal{B}}\|^{2}]$
as

$\displaystyle\begin{split}\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-{\mathbf{v}}^{k\mathcal{B}}\|^{2}]&\leq
2\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-{\mathbf{v}}^{k\mathcal{B}}-(\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}}))\|^{2}]\\\
&\quad+2\mathbb{E}[\|\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}]\\\
&\leq
2\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})\|^{2}]\\\
&\quad+2\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}]+2L^{2}\mathbb{E}[\|{\mathbf{x}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{k\mathcal{B}}\|^{2}]\\\
&\leq
2\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})\|^{2}]\\\
&\quad+2\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}]+2L^{2}\mathbb{E}[\|{\mathbf{z}}^{(k+1)\mathcal{B}}-{\mathbf{z}}^{k\mathcal{B}}\|^{2}],\end{split}$
(77)

where $\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})=[\nabla
f_{1}({\mathbf{x}}_{1}^{(k+1)\mathcal{B}});\cdots;\nabla
f_{n}({\mathbf{x}}_{n}^{(k+1)\mathcal{B}})]$. To proceed, let us derive an
upper bound on
$E[\|{\mathbf{z}}^{(k+1)\mathcal{B}}-{\mathbf{z}}^{k\mathcal{B}}\|^{2}]$.
First, consider each column of ${\mathbf{z}}^{(k+1)\mathcal{B}}$ and
${\mathbf{z}}^{k\mathcal{B}}$ (i.e., ${\mathbf{z}}_{:m}^{(k+1)\mathcal{B}}$
and ${\mathbf{z}}_{:m}^{k\mathcal{B}}$) separately and observe that

$\displaystyle\begin{split}\|{\mathbf{z}}_{:m}^{(k+1)\mathcal{B}}-{\mathbf{z}}_{:m}^{k\mathcal{B}}\|^{2}&=\|M_{m}((k+1)\mathcal{B}-1:k\mathcal{B}){\mathbf{z}}_{:m}^{k\mathcal{B}}-\alpha{\mathbf{g}}_{:m}^{k\mathcal{B}}-{\mathbf{z}}_{:m}^{k\mathcal{B}}\|^{2}\leq
8\|{\mathbf{z}}_{:m}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}\mathbf{1}_{2n}\|^{2}+2\alpha^{2}\|{\mathbf{g}}_{:m}^{k\mathcal{B}}\|^{2}.\end{split}$
(78)

Then

$\|{\mathbf{z}}^{(k+1)\mathcal{B}}-{\mathbf{z}}^{k\mathcal{B}}\|^{2}\leq
8\sum_{m=1}^{d}\|{\mathbf{z}}_{:m}^{k\mathcal{B}}-\bar{z}_{m}^{k\mathcal{B}}\mathbf{1}_{2n}\|^{2}+2\alpha^{2}\|{\mathbf{g}}^{k\mathcal{B}}\|^{2}.$
(79)

To derive an upper bound on $\|{\mathbf{g}}_{\sim n}^{k\mathcal{B}}\|$, let
$\bar{{\mathbf{z}}}^{k\mathcal{B}}=[\bar{z}_{1}^{k\mathcal{B}}\mathbf{1}_{2n},\cdots,\bar{z}_{d}^{k\mathcal{B}}\mathbf{1}_{2n}]$
and note that

$\displaystyle\begin{split}\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}\|&=\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}+\mathbf{1}_{n}(\bar{{\mathbf{v}}}^{k\mathcal{B}})^{\prime}-\mathbf{1}_{n}(\nabla\bar{{\mathbf{f}}}({\mathbf{x}}^{k\mathcal{B}}))^{\prime}+\mathbf{1}_{n}(\nabla\bar{{\mathbf{f}}}({\mathbf{x}}^{k\mathcal{B}}))^{\prime}-\mathbf{1}_{n}(\nabla\bar{{\mathbf{f}}}({\mathbf{x}}^{*}))^{\prime}\|\\\
&\leq\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}\|+\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|+L\|{\mathbf{x}}^{k\mathcal{B}}-\mathbf{1}_{n}({\mathbf{x}}^{*})^{\prime}\|\\\
&\leq\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}\|+\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|+L\|{\mathbf{x}}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}+\mathbf{1}_{n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}-\mathbf{1}_{n}({\mathbf{x}}^{*})^{\prime}\|\\\
&\leq\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}\|+\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|+L\|{\mathbf{z}}^{k\mathcal{B}}-\mathbf{1}_{2n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}\|+\sqrt{2n}L\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|.\end{split}$
(80)

Squaring both sides of the above inequality yields

$\displaystyle\|{\mathbf{g}}_{\sim n}^{k\mathcal{B}}\|^{2}$ $\displaystyle\leq
4L^{2}\|{\mathbf{z}}^{k\mathcal{B}}-\mathbf{1}_{2n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}\|^{2}+8nL^{2}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}+4\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}\|^{2}+4\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}.$
(81)

Imposing $0<\alpha<\frac{1}{4\sqrt{2}L}$,

$\displaystyle\mathbb{E}[\|{\mathbf{z}}^{(k+1)\mathcal{B}}-{\mathbf{z}}^{k\mathcal{B}}\|^{2}]$
$\displaystyle\leq
8.25\mathbb{E}[\|{\mathbf{z}}^{k\mathcal{B}}-\mathbf{1}_{2n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}\|^{2}]+0.5\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$
(82) $\displaystyle\quad+8\alpha^{2}\mathbb{E}[\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}\|^{2}]+8\alpha^{2}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}].$

Then

$\displaystyle\begin{split}\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-{\mathbf{v}}^{k\mathcal{B}}\|^{2}]&\leq
16.5L^{2}\mathbb{E}[\|{\mathbf{z}}^{k\mathcal{B}}-\mathbf{1}_{2n}(\bar{{\mathbf{z}}}^{k\mathcal{B}})^{\prime}\|^{2}]+L^{2}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+16\alpha^{2}L^{2}\mathbb{E}[\|{\mathbf{g}}_{\sim
n}^{k\mathcal{B}}-\mathbf{1}_{n}(\bar{{\mathbf{g}}}^{k\mathcal{B}})^{\prime}\|^{2}]+2.5\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{k\mathcal{B}})\|^{2}]\\\
&\quad+2\mathbb{E}[\|{\mathbf{v}}^{(k+1)\mathcal{B}}-\nabla{\mathbf{f}}({\mathbf{x}}^{(k+1)\mathcal{B}})\|^{2}].\end{split}$
(83)

The proof is completed by combining (76) with (83).

∎

The gradient estimate error, $\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla
f({\mathbf{x}}^{k\mathcal{B}})\|^{2}]$, appearing on the right-hand side of
the inequalities in Lemmas 3.1 and 3.2, is analyzed in the following lemma.

######  Lemma $3.3$.

Suppose the objective function $f$ is $\mu$-strongly-convex, and let $\tau$,
$\bar{\tau}$ be defined as above. Then $\forall k\geq 0$,

$\displaystyle\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla
f({\mathbf{x}}^{k\mathcal{B}})\|^{2}]$ $\displaystyle\leq
4L^{2}\sum_{i=1}^{n}\mathbb{E}[\|{\mathbf{x}}_{i}^{k\mathcal{B}}-\bar{{\mathbf{z}}}^{k\mathcal{B}}\|^{2}]+4L^{2}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]$
(84)
$\displaystyle\quad+4L^{2}\sum_{i=1}^{n}\mathbb{E}[\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]+4L^{2}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].$

###### Proof.

For all $i\leq n$, it holds that

$\displaystyle\begin{split}\mathbb{E}[\|{\mathbf{v}}_{i}^{k\mathcal{B}}-\nabla
f_{i}({\mathbf{x}}_{i}^{k\mathcal{B}})\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]&=\mathbb{E}[\|\nabla
f_{i,l_{i}^{k\mathcal{B}}}({\mathbf{x}}_{i}^{k\mathcal{B}})-f_{i,l_{i}^{k\mathcal{B}}}(\tau_{i}^{k\mathcal{B}})-(\nabla
f_{i}({\mathbf{x}}_{i}^{k\mathcal{B}})-\nabla
f_{i}(\tau_{i}^{k\mathcal{B}}))\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&\leq\mathbb{E}[\|\nabla
f_{i,l_{i}^{k\mathcal{B}}}({\mathbf{x}}_{i}^{k\mathcal{B}})-f_{i,l_{i}^{k\mathcal{B}}}(\tau_{i}^{k\mathcal{B}})\|^{2}|{\mathcal{F}}^{k\mathcal{B}}]\\\
&=\frac{1}{m_{i}}\sum_{j=1}^{m_{i}}\|\nabla
f_{i,j}({\mathbf{x}}_{i}^{k\mathcal{B}})-\nabla
f_{i,j}({\mathbf{x}}^{*})+(\nabla f_{i,j}({\mathbf{x}}^{*})-\nabla
f_{i,j}(\tau_{i}^{k\mathcal{B}}))\|^{2}\\\ &\leq
2L^{2}\|{\mathbf{x}}_{i}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}+2L^{2}\|\tau_{i}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}\\\
&\leq
4L^{2}\|{\mathbf{x}}_{i}^{k\mathcal{B}}-\bar{{\mathbf{z}}}^{k\mathcal{B}}\|^{2}+4L^{2}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}+4L^{2}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}+4L^{2}\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}.\end{split}$
(85)

The proof of the lemma is completed by summing over $i$ from $1$ to $n$ and
taking the total expectation.

∎

Combining the results of Lemma 2, 3.1 and 3.3, we obtain the following result.

######  Lemma $3.4$.

Suppose the objective function $f$ is $\mu$-strongly-convex. If
$0<\alpha\leq\frac{1}{8L}$, then for all $k\geq 0$ it holds

$\displaystyle\begin{split}\mathbb{E}[\|{\mathbf{v}}_{i}^{(k+1)\mathcal{B}}-\nabla
f_{i}({\mathbf{x}}_{i}^{(k+1)\mathcal{B}})\|^{2}]&\leq
16.75L^{2}\mathbb{E}[\|{\mathbf{x}}_{i}^{(k+1)\mathcal{B}}-\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}\|^{2}]\\\
&\quad+16L^{2}\alpha^{2}\mathbb{E}[\|{\mathbf{g}}_{i}^{k\mathcal{B}}-\mathbf{1}_{n}\bar{g}^{k\mathcal{B}}\|^{2}]+16.5L^{2}\mathbb{E}[\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+4.5L^{2}\mathbb{E}[\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]+4.5L^{2}\mathbb{E}[\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].\end{split}$
(86)

###### Proof.

The proof is completed by combining Lemma 1, 2 and 3.3. ∎

We can now present an argument proving Lemmas 3 and 4 in the main paper. In
particular, combining Lemma 3.1

$\displaystyle\mathbb{E}[\|\bar{{\mathbf{z}}}^{(t+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\leq\frac{L^{2}\alpha}{n\mu}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+(1-\mu\alpha)\mathbb{E}[\|\bar{z}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]+\frac{\alpha^{2}}{n^{2}}\mathbb{E}[\|{\mathbf{v}}^{k\mathcal{B}}-\nabla
f({\mathbf{x}}^{k\mathcal{B}})\|^{2}]$

and the results in Lemma 3.2, we obtain

$\displaystyle\begin{split}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]&\leq
L^{2}\alpha(\frac{1}{\mu}+\frac{4\alpha}{n})\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+(1-\mu\alpha+\frac{4L^{2}\alpha^{2}}{n})\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].\end{split}$
(87)

By letting $0<\alpha\leq\frac{\mu}{8L^{2}}$, we can further bound

$\displaystyle\begin{split}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{(k+1)\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]&\leq\frac{2L^{2}\alpha}{\mu}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]+(1-\frac{\mu\alpha}{2})\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]+\frac{4L^{2}\alpha^{2}}{n}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}],\end{split}$

which completes the proof of Lemma 3. Moreover,

$\displaystyle\begin{split}\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}]&\leq\frac{120L^{2}}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{89L^{2}}{1-\sigma^{2}}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+(\frac{1+\sigma^{2}}{2}+\frac{96L^{2}\alpha^{2}}{1-\sigma^{2}})\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}|^{2}]\\\
&\quad+\frac{38L^{2}}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{38L^{2}}{1-\sigma^{2}}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].\end{split}$
(88)

For $0<\alpha\leq\frac{1-\sigma^{2}}{14\sqrt{2}L}$, we have
$\frac{1+\sigma^{2}}{2}+\frac{98L^{2}\alpha^{2}}{1-\sigma^{2}}\leq\frac{3+\sigma^{2}}{4}$;
this helps complete the proof of Lemma 4,

$\displaystyle\begin{split}\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}|g_{im}^{(k+1)\mathcal{B}}-\bar{g}_{m}^{(k+1)\mathcal{B}}|^{2}]&\leq\frac{120L^{2}}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{89L^{2}}{1-\sigma^{2}}\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
&\quad+\frac{3+\sigma^{2}}{4}\mathbb{E}[\sum_{m=1}^{d}\sum_{i=1}^{n}\|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{38L^{2}}{1-\sigma^{2}}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]\\\
&\quad+\frac{38L^{2}}{1-\sigma^{2}}\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}].\end{split}$

Using the inequalities shown in Lemma 2, 3 and 4, we can construct a dynamic
system and continue the proof of linear convergence of Algorithm 1. To this
end, we first define

${\mathbf{u}}^{k\mathcal{B}}=\begin{bmatrix}\mathbb{E}[\sum_{i=1}^{n}\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{z}}_{i}^{k\mathcal{B}}\|^{2}]\\\
\mathbb{E}[n\|\bar{{\mathbf{z}}}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
\mathbb{E}[\frac{\sum_{m=1}^{d}\sum_{i=1}^{n}\|g_{im}^{k\mathcal{B}}-\bar{g}_{m}^{k\mathcal{B}}\|^{2}}{L^{2}}]\end{bmatrix}$
(89)

$\tilde{{\mathbf{u}}}^{k\mathcal{B}}=\begin{bmatrix}\mathbb{E}[\sum_{i=1}^{n}\|\tau_{i}^{k\mathcal{B}}-\bar{\tau}^{k\mathcal{B}}\|^{2}]\\\
\mathbb{E}[n\|\bar{\tau}^{k\mathcal{B}}-{\mathbf{x}}^{*}\|^{2}]\\\
\mathbf{0}\end{bmatrix}$ (90)

$J_{\alpha}=\begin{bmatrix}\frac{1+\sigma^{2}}{2}&0&\frac{2\alpha^{2}L^{2}}{1-\sigma^{2}}\\\
\frac{2L^{2}\alpha}{\mu}&1-\frac{\mu\alpha}{2}&0\\\
\frac{120}{1-\sigma^{2}}&\frac{89}{1-\sigma^{2}}&\frac{3+\sigma^{2}}{4}\end{bmatrix}$
(91)

$H_{\alpha}=\begin{bmatrix}0&0&0\\\
\frac{4L^{2}\alpha^{2}}{n}&\frac{4L^{2}\alpha^{2}}{n}&0\\\
\frac{38}{1-\sigma^{2}}&\frac{38}{1-\sigma^{2}}&0\end{bmatrix}$ (92)

and then formally state the dynamic system in Proposition 1.

######  Proposition 2.

Suppose Assumption 1 holds, the objective function $f$ is $\mu$-strongly-
convex and each component of the local objective function $f_{i,j}$ is
$L$-smooth. If $0\leq\alpha\leq\frac{\mu(1-\sigma^{2})}{14\sqrt{2}L^{2}}$,
then for any $k\geq 0$

${\mathbf{u}}^{(k+1)\mathcal{B}}\leq
J_{\alpha}{\mathbf{u}}^{k\mathcal{B}}+H_{\alpha}\tilde{{\mathbf{u}}}^{k\mathcal{B}}.$
(93)

It follows that for the inner loop, for all $k\in[sT,(s+1)T-1]$

${\mathbf{u}}^{(k+1)\mathcal{B}}\leq
J_{\alpha}{\mathbf{u}}^{k\mathcal{B}}+H_{\alpha}{\mathbf{u}}^{sT}.$ (94)

For the outer loop, for all $s\geq 0$, it holds

${\mathbf{u}}^{(s+1)T}\leq(J_{\alpha}^{T}+\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}){\mathbf{u}}^{sT}.$
(95)

To guarantee linear decay of the outer loop sequence and ultimately show
linear convergence of Algorithm 1, we require that the spectral radius of
$J_{\alpha}^{T}+\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}$ is small. In Lemma
5, we compute the range of the step size, $\alpha$, such that the weighted
matrix norms of both $J_{\alpha}^{T}$ and
$\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}$ are small.

######  Lemma $5$.

Suppose Assumption 1 holds and assume that
$0<\alpha\leq\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L}$, where
$\tilde{Q}=\frac{L}{\mu}$. Then,

$\rho(J_{\alpha})<\||J_{\alpha}|\|^{\mathbf{\delta}}_{\infty}<1-\frac{\mu\alpha}{4},$
(96)

and

$\||\sum_{l=0}^{T-1}J_{\alpha}^{l}H_{\alpha}|\|^{{\mathbf{q}}}_{\infty}\leq\||(I-J_{\alpha})^{-1}H_{\alpha}|\|^{{\mathbf{q}}}_{\infty}<0.66,$
(97)

where
$\mathbf{\delta}=\begin{bmatrix}1,8\tilde{Q}^{2},\frac{6656\tilde{Q}^{2}}{(1-\sigma^{2})^{2}}\end{bmatrix}$
and ${\mathbf{q}}=[1,1,\frac{1457}{(1-\sigma^{2})^{2}}]$.

###### Proof.

Following Lemma 10 in [42], consider a matrix $A\in{\mathbb{R}}^{d\times d}$
and a positive vector ${\mathbf{x}}\in{\mathbb{R}}^{d}$, and note that if
$A{\mathbf{x}}\leq\beta{\mathbf{x}}$ for $\beta>0$ we have that
$\rho(A)\leq\||A|\|^{{\mathbf{x}}}_{\infty}\leq\beta$. Using this lemma we
solve for a range of $\alpha$ and a positive vector
$\mathbf{\delta}\in{\mathbb{R}}^{3}$ such that

$J_{\alpha}\mathbf{\delta}\leq(1-\frac{\mu\alpha}{4})\mathbf{\delta},$ (98)

which is equivalent to the element-wise inequalities

$\displaystyle\frac{1+\sigma^{2}}{2}\delta_{1}+\frac{2\alpha^{2}L^{2}}{1-\sigma^{2}}\delta_{3}\leq(1-\frac{\mu\alpha}{4})\delta_{1}$
(99)
$\displaystyle\frac{2L^{2}\alpha}{\mu}\delta_{1}+(1-\frac{\mu\alpha}{2})\delta_{2}\leq(1-\frac{\mu\alpha}{4})\delta_{2}$
$\displaystyle\frac{120}{1-\sigma^{2}}\delta_{1}+\frac{89}{1-\sigma^{2}}\delta_{2}+\frac{3+\sigma^{2}}{4}\delta_{3}\leq(1-\frac{\mu\alpha}{4})\delta_{3}.$

To solve for a meaningful $\mathbf{\delta}$, we set $\delta_{1}=1$ and
$\delta_{2}=8\tilde{Q}^{2}$; then
$\mathbf{\delta}=\begin{bmatrix}1,8\tilde{Q}^{2},\frac{6656\tilde{Q}^{2}}{(1-\sigma^{2})^{2}}\end{bmatrix}$
and $0<\alpha\leq\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L}$ are sufficient to
satisfy the first inequality. Since $J_{\alpha}$ is non-negative,
$\sum_{l=0}^{T-1}J_{\alpha}^{l}\leq\sum_{l=0}^{\infty}J_{\alpha}^{l}=(I_{3}-J_{\alpha})^{-1}$;
this yields

${\mathbf{u}}^{(s+1)T}\leq(J_{\alpha}^{T}+(I_{3}-J_{\alpha})^{-1}H_{\alpha}){\mathbf{u}}^{sT},$
(100)

where

$I_{3}-J_{\alpha}=\begin{bmatrix}\frac{1-\sigma^{2}}{2}&0&-\frac{2\alpha^{2}L^{2}}{1-\sigma^{2}}\\\
-\frac{2L^{2}\alpha}{\mu}&\frac{\mu\alpha}{2}&0\\\
-\frac{120}{1-\sigma^{2}}&-\frac{89}{1-\sigma^{2}}&\frac{1-\sigma^{2}}{4}\end{bmatrix}$
(101)

and its determinant is

$\det(I_{3}-J_{\alpha})=\frac{(1-\sigma^{2})^{2}\mu\alpha}{16}-\frac{356L^{4}\alpha^{3}}{\mu(1-\sigma^{2})^{2}}-\frac{120\alpha^{3}\mu
L^{2}}{(1-\sigma^{2})^{2}}.$ (102)

When $0<\alpha\leq\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L}$,
$\det(I_{3}-J_{\alpha})\geq\frac{(1-\sigma^{2})^{2}\mu\alpha}{32}$, and

$\displaystyle[\mathrm{adj}(I_{3}-J_{\alpha})]_{1,2}=\frac{178L^{2}\alpha^{2}}{(1-\sigma^{2})^{2}},$
$\displaystyle[\mathrm{adj}(I_{3}-J_{\alpha})]_{1,3}=\frac{\mu
L^{2}\alpha^{3}}{(1-\sigma^{2})},$
$\displaystyle[\mathrm{adj}(I_{3}-J_{\alpha})]_{2,2}\leq\frac{(1-\sigma^{2})^{2}}{8},$
$\displaystyle[\mathrm{adj}(I_{3}-J_{\alpha})]_{2,3}=\frac{4L^{4}\alpha^{3}}{\mu(1-\sigma^{2})},$
$\displaystyle[\mathrm{adj}(I_{3}-J_{\alpha})]_{3,2}=44.5,$
$\displaystyle[\mathrm{adj}(I_{3}-J_{\alpha})]_{3,3}=\frac{\mu\alpha(1-\sigma^{2})}{4}.$

Next, we derive a matrix upper-bounding (element-wise)
$(I_{3}-J_{\alpha})^{-1}H_{\alpha}=\frac{\mathrm{adj}(I_{3}-J_{\alpha})}{\det(I_{3}-J_{\alpha})}H_{\alpha}$.
For $0\leq\alpha\leq\frac{(1-\sigma^{2})^{2}}{187\tilde{Q}L}$,

$(I_{3}-J_{\alpha})^{-1}H_{\alpha}\leq\begin{bmatrix}0.039&0.039&0\\\
0.23&0.23&0\\\
\frac{335}{(1-\sigma^{2})^{2}}&\frac{335}{(1-\sigma^{2})^{2}}&0\end{bmatrix}.$
(103)

If ${\mathbf{q}}=[1,1,\frac{1457}{(1-\sigma^{2})^{2}}]$, we have

$((I_{3}-J_{\alpha})^{-1}H_{\alpha}){\mathbf{q}}\leq 0.66{\mathbf{q}}.$ (104)

Finally, invoking the definition of the weighted matrix norm, we complete the
proof of the second inequality in the lemma. ∎

This completes the presentation of auxiliary results that support the proof of
Theorem 1 in the main paper.

##  Appendix B Experimental Results for the Decentralized Average Consensus
Problem

(a) Consensus residual: $\mathcal{B}=1$

(b) Consensus residual: $\mathcal{B}=10$

Figure 4: Average consensus on a jointly connected network with
$\mathcal{B}=1,10,\epsilon=0.05$. In each of the subplots, we show the
performance of Di-CS-AC, i.e., Algorithm 1, for $6$ different sparsification
levels and compare it to $2$ benchmark quantization algorithms, Q-Push-sum and
Q-Push-Gossip. The quantization level is chosen such that the number of
communicated bits for the benchmark algorithms is equal to that of Di-CS-AC
when $q=0.078$.

In this section, we refer to the proposed Algorithm 1 in the main paper as Di-
CS-AC (Directed Communication-Sparsifying Average Consensus) and present its
performance.

We consider an average consensus problem where the dimension of a local
parameter vector at each node is $d=64$. The initial state
$\mathbf{x}_{i}^{0}$ is randomly generated from the normal distribution; the
goal of the network is to reach the average consensus vector, i.e., compute
$\bar{\mathbf{x}}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i}^{0}$. The network
setups are exactly the same as the decentralized optimization experiments in
Section V.

For benchmarking purposes, we consider two quantized versions of the push-sum
algorithm: (i) Q-Push-sum, obtained by applying simple quantization to the
push-sum scheme [21, 22], and (ii) Q-Push-Gossip, a quantized push-sum for
gossip algorithm recently proposed in [38]. The former was originally
developed for unconstrainted communication settings, while the latter
originally targeted static networks; in the absence of prior work on
communication-constrained consensus over time-varying directed networks, we
adopt these two as the benchmarking schemes.

We compare the performance of different algorithms by computing the residual
value
$\frac{\|\mathbf{x}^{t}-\bar{\mathbf{x}}\|}{\|\mathbf{x}^{0}-\bar{\mathbf{x}}\|}$;
the results are shown in Fig. 4. As the figures demonstrate, at the considered
levels of sparsification $q$ and values of the connectivity parameter
$\mathcal{B}$, Di-CS-AC converges to the same limit as the full communication
schemes. The convergence rate is linear in the number of iterations $t$ but
smaller compression level and larger connectivity period slow the convergence
down. In Fig. 4 (a) and (b), the two benchmarking quantization algorithms
cannot reach the desired consensus accuracy in the time-varying directed
network while the proposed Di-CS-AC achieves considerably smaller consensus
error.

In regards to communication cost, to reach the same residual threshold (i.e.,
$10^{-12}$), our consensus algorithm, Di-CS-AC, with $q=0.078$, incurs only
$7.8\%$ per-iteration cost as compared to the non-compression consensus
algorithm (with $q=1$) and needs fewer than double the number of iterations
(see Fig. 4 (a) and (b)). This demonstrates communication efficiency of Di-CS-
AC.

##  Appendix C Experimental results for larger size of network

In this section, we consider a network with $100$ nodes and repeat the
logistic regression experiment. We still use the stackoverflow dataset
introduced in Section V of the main paper. The network is constructed
following the same steps introduced in Section V of the main paper with
parameter $B=1$. In Fig. 5, we show the correct rate for three different
sparsification levels, i.e., $q=0.25,0.5$ and $1$.

Figure 5: Di-CS-SVRG performance on a time-varying directed network with $100$
nodes.

As the plot shows, schemes with different sparsification levels have different
convergence rates – more aggressive sparsification leads to slower
convergence.