# Zero-Order One-Point Estimate with Distributed Stochastic Gradient-Tracking
Technique

Elissa Mhanna<EMAIL_ADDRESS>
Mohamad Assaad<EMAIL_ADDRESS>
Laboratoire des Signaux & Systèmes
CentraleSupélec
3 rue Joliot Curie
91190 Gif-sur-Yvette, France

###### Abstract

In this work, we consider a distributed multi-agent stochastic optimization
problem, where each agent holds a local objective function that is smooth and
convex, and that is subject to a stochastic process. The goal is for all
agents to collaborate to find a common solution that optimizes the sum of
these local functions. With the practical assumption that agents can only
obtain noisy numerical function queries at exactly one point at a time, we
extend the distributed stochastic gradient-tracking method to the bandit
setting where we don’t have an estimate of the gradient, and we introduce a
zero-order (ZO) one-point estimate (1P-DSGT). We analyze the convergence of
this novel technique for smooth and convex objectives using stochastic
approximation tools, and we prove that it converges almost surely to the
optimum. We then study the convergence rate for when the objectives are
additionally strongly convex. We obtain a rate of $O(\frac{1}{\sqrt{k}})$
after a sufficient number of iterations $k>K_{2}$ which is usually optimal for
techniques utilizing one-point estimators. We also provide a regret bound of
$O(\sqrt{k})$, which is exceptionally good compared to the aforementioned
techniques. We further illustrate the usefulness of the proposed technique
using numerical experiments.

Keywords: Gradient-free optimization, stochastic framework, distributed
algorithms, consensus, gradient tracking, convergence analysis

##  1 Introduction

Gradient-free optimization is an old topic in the research community; however,
there has been an increased interest recently, especially in machine learning
applications, where optimization problems are typically solved with gradient
descent algorithms. Successful applications of gradient-free methods in
machine learning include competing with an adversary in bandit problems
(Flaxman et al., 2004; Agarwal et al., 2010), generating adversarial attacks
for deep learning models (Chen et al., 2019; Liu et al., 2019) and
reinforcement learning (Vemula et al., 2019). Gradient-free optimization aims
to solve optimization problems with only functional (zero-order) information
rather than first-order (FO) gradient information. These techniques are
essential in settings where explicit gradient computation may be impractical,
expensive, or impossible. Instances of such settings include high data
dimensionality, time or resource straining function differentiation, or the
cost function not having a closed-form. Zero-order information-based methods
include direct search methods (Golovin et al., 2019), 1-point methods (Li and
Assaad, 2021; Flaxman et al., 2004; Bach and Perchet, 2016; Vemula et al.,
2019) where the function is evaluated at a single point with some
randomization to estimate the gradient, 2- or more point methods (Duchi et
al., 2015; Nesterov and Spokoiny, 2017; Gorbunov et al., 2018; Bach and
Perchet, 2016; Hajinezhad et al., 2019; Kumar Sahu et al., 2018; Agarwal et
al., 2010; Chen et al., 2019; Liu et al., 2019; Vemula et al., 2019), where
functional difference at various points is employed for estimation, and other
methods such as sign information of gradient estimates (Liu et al., 2019).

Another area of great interest is distributed multi-agent optimization, where
agents try to cooperatively solve a problem with information exchange only
limited to immediate neighbors in the network. Distributed computing and data
storing are particularly essential in fields such as vehicular communications
and coordination, data processing and distributed control in sensor networks
(Shi and Eryilmaz, 2020), big-data analytics (Daneshmand et al., 2015), and
federated learning (McMahan et al., 2017). More specifically, one direction of
research integrates (sub)gradient-based methods with a consensus/averaging
strategy; the local agent incorporates one or multiple consensus steps
alongside evaluating the local gradient during optimization. Hence, these
algorithms can tackle a fundamental challenge: overcoming differences between
agents’ local data distributions.

###  1.1 Problem Description

Consider a set of agents $\mathcal{N}=\\{1,2,\ldots,n\\}$ connected by a
communication network. Each agent $i$ is associated with a local objective
function $f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$. The global goal of the
agents is to collaboratively locate the decision variable $x\in\mathbb{R}^{d}$
that solves the stochastic optimization problem:

$min_{x\in\mathbb{R}^{d}}\mathcal{F}(x)=\frac{1}{n}\sum_{i=1}^{n}F_{i}(x)$ (1)

where

$F_{i}(x)=\mathbb{E}_{S}f_{i}(x,S),$

with $S\in\mathcal{S}$ denoting an i.i.d. ergodic stochastic process
describing uncertainties in the communication system.

We assume that at each time step, agent $i$ can only query the function values
of $f_{i}$ at exactly one point, and can only communicate with its neighbors.
Further, we assume that the function queries are noisy
$\tilde{f}_{i}=f_{i}+\zeta_{i}$ with $\zeta_{i}$ some additive noise. Agent
$i$ must then employ this query to estimate the gradient of the form
$g_{i}(x,S_{i})$.

One efficient algorithm with a straightforward averaging scheme to solve this
problem is the gradient-tracking (GT) technique, which has proven to achieve
rates competing with its centralized counterparts. For example, the acquired
error bound under a distributed stochastic variant was found to decrease with
the network size n (Pu and Nedić, 2018). In most work, this technique proved
to converge linearly to the optimal solution with constant step size (Qu and
Li, 2018; Nedić et al., 2017; Pu, 2020), which is also a unique attribute
among other distributed stochastic gradient algorithms. It has been extended
to time-varying (undirected or directed) graphs (Nedić et al., 2017), a
gossip-like method which is efficient in communication (Pu and Nedić, 2018),
and nonconvex settings (Tang et al., 2021; Lorenzo and Scutari, 2016; Jiang et
al., 2022; Lu et al., 2019). All references mentioned above consider the case
where an accurate gradient computation or a non-biased gradient estimation
with bounded variance (BV) is available.

###  1.2 Function Classes

Consider the following five classes of functions:

  * •

The convex class $\mathcal{C}_{cvx}$ containing all functions
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that are convex.

  * •

The strongly convex class $\mathcal{C}_{sc}$ containing all functions
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that are continuously differentiable
and admit a constant $\mu_{f}$ such that

$\langle\nabla f(x)-\nabla f(y),x-y\rangle\|\geq\mu_{f}\|x-y\|^{2},\;\;\forall
x,y\in\mathbb{R}^{d}.$

  * •

The Lipschitz continuous class $\mathcal{C}_{lip}$ containing all functions
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that admit a constant $L_{f}$ such
that

$|f(x)-f(y)|\leq L_{f}\|x-y\|,\;\;\forall x,y\in\mathbb{R}^{d}.$

  * •

The smooth class $\mathcal{C}_{smo}$ containing all functions
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that are continuously differentiable
and admit a constant $G_{f}$ such that

$\|\nabla f(x)-\nabla f(y)\|\leq G_{f}\|x-y\|,\;\;\forall
x,y\in\mathbb{R}^{d}.$

  * •

The gradient dominated class $\mathcal{C}_{gd}$ containing all functions
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that are differentiable, have a global
minimizer $x^{*}$, and admit a constant $\nu_{f}$ such that

$2\nu_{f}(f(x)-f(x^{*}))\leq\|\nabla f(x)\|^{2},\;\;\forall
x\in\mathbb{R}^{d}.$

|  GRADIENT ESTIMATE | OP |  FUNCTION CLASS |  CONS- ENSUS |  REGRET BOUND |  CONVERGENCE RATE
---|---|---|---|---|---|---
ZO | One-point | Cent. | $\mathcal{C}_{cvx}\bigcap\mathcal{C}_{lip}$ | - | $O(k^{\frac{3}{4}})$ |  $O(\frac{1}{\sqrt[4]{k}})$ Flaxman et al. (2004)
Dist. | $\mathcal{C}_{sc}\bigcap\mathcal{C}_{lip}\bigcap\mathcal{C}_{smo}$ | None | - |  $O(\frac{1}{\sqrt{k}})$ Li and Assaad (2021)
Dist. | $\mathcal{C}_{sc}\bigcap\mathcal{C}_{lip}\bigcap\mathcal{C}_{smo}$ | GT | $O(\sqrt{k})$ |  $O(\frac{1}{\sqrt{k}})$ 1P-DSGT
Two-point | Cent. | $\mathcal{C}_{cvx}\bigcap\mathcal{C}_{lip}$ | - | $O(\sqrt{k})$ |  $O(\frac{1}{\sqrt{k}})$ Agarwal et al. (2010)
Cent. | $\mathcal{C}_{sc}\bigcap\mathcal{C}_{lip}$ | - | $O(\log k)$ |  $O(\frac{\log k}{k})$Agarwal et al. (2010)
Dist. | $\mathcal{C}_{lip}\bigcap\mathcal{C}_{smo}$ | None | - |  $O(\frac{1}{\sqrt{k}}\log k)$ Tang et al. (2021)
Dist. | $\mathcal{C}_{smo}\bigcap\mathcal{C}_{gd}$ | None | - |  $O(\frac{1}{k})$ Tang et al. (2021)
(d+1)-point | Cent. | $\mathcal{C}_{sc}\bigcap\mathcal{C}_{lip}\bigcap\mathcal{C}_{smo}$ | - | $O(\log k)$ |  $O(\frac{\log k}{k})$ Agarwal et al. (2010)
2d-point | Dist. | $\mathcal{C}_{smo}$ | GT | - |  $O(\frac{1}{k})$ Tang et al. (2021)
Dist. | $\mathcal{C}_{smo}\bigcap\mathcal{C}_{gd}$ | GT | - |  $O(\lambda^{k})$ Tang et al. (2021)
Kernel-based | Cent. | $\mathcal{C}_{cvx}\bigcap\mathcal{C}_{lip}$ | - | $O(\sqrt{k})$ |  $O(\frac{1}{\sqrt{k}})$ Bubeck et al. (2021)
FO | Unbiased/BV | Dist. | $\mathcal{C}_{sc}\bigcap\mathcal{C}_{smo}$ | GT | - |  $O(\lambda^{k})$ Pu and Nedić (2018);
|  |  |  | Xin et al. (2019); Pu (2020)
Dist. | $\mathcal{C}_{smo}$ | GT | - |  $O(\frac{1}{\sqrt{k}})$ Lu et al. (2019)
Table 1: Convergence rates for various algorithms related to our work,
classified according to the nature of the gradient estimate, whether the
optimization problem (OP) is centralized or distributed, the assumptions on
the objective function, whether consensus is aimed at or not, and the achieved
regret bound and convergence rate

###  1.3 Related Work

An adversarial convex bandit problem in a zero-order framework is studied in
Flaxman et al. (2004); Agarwal et al. (2010); Bubeck et al. (2021). In Flaxman
et al. (2004), a one point estimator is proposed and a regret bound of
$O(k^{\frac{3}{4}})$ is achieved for Lipschitz continuous functions. In
Agarwal et al. (2010), they extend the problem to the multi-point setting,
where the player queries each loss function at many points; hence, their
estimates are unbiased with respect to the true gradient and have bounded
variance. When the number of points is two, they prove regret bounds of
$\tilde{O}(\sqrt{k})$ with high probability and of $O(\log(k))$ in expectation
for strongly convex loss functions. When the number is $d+1$ point, they prove
regret bounds of $O(\sqrt{k})$ and of $O(\log(k))$ with strong convexity. As
mentioned in Agarwal et al. (2010) and the references therein, these bounds
achieved are even optimal for the full-information setting. In Bubeck et al.
(2021), a kernelized loss estimator is proposed where a generalization of
Bernoulli convolutions is adopted, and an annealing schedule for exponential
weights is used to control the estimator’s variance in a focus region for
dimensions higher than 1. Bubeck et al. (2021) achieves a regret bound of
$O(\sqrt{k})$. However, all these references solve the bandit problem in a
centralized framework or the case where a single agent has all the data and is
the only one responsible for the decision-making.

Li and Assaad (2021) consider a distributed stochastic optimization problem in
networks under the assumption of smoothness and concavity. They solve it using
a stochastic gradient descent method, where the exchange between the nodes is
limited to a partial trade of the observation of their local utilities, which
each agent then uses to estimate the global utility. With this estimated
global utility, each agent can construct a one-point gradient estimator based
on a stochastic perturbation. The convergence rate is proved to scale as
$O(\frac{1}{\sqrt{k}})$, under the further assumption of a strongly concave
objective function. We must remark that in Li and Assaad (2021), there is a
distributed solution, the optimization variable is a scalar, not a vector, and
each agent can only update its own scalar. Hence, there is no consensus to be
aimed at, as in our case. In the context of the derivative-free stochastic
_centralized_ optimization, all Duchi et al. (2015); Jamieson et al. (2012);
Shamir (2013) have derived several lower bounds to confirm that the
convergence rate cannot be better than $O(\frac{1}{\sqrt{k}})$ after k
iterations for strongly convex and smooth objective functions.

All Qu and Li (2018); Lorenzo and Scutari (2016); Nedić et al. (2017); Shi et
al. (2015); Li et al. (2022); Jiang et al. (2022) present a distributed
gradient-tracking method that employs local auxiliary variables to track the
average of all agents’ gradients, considering the availability of accurate
gradient information. In both Li et al. (2022); Jiang et al. (2022), each
local objective function is an average of finite instantaneous functions.
Thus, they incorporate the gradient-tracking algorithm with stochastic
averaging gradient technology (Li et al., 2022) (smooth convex optimization)
or with variance reduction techniques (Jiang et al., 2022) (smooth nonconvex
optimization). At each iteration, they randomly select only one gradient of an
instantaneous function to approximate the local batch gradient. In Li et al.
(2022), this is an unbiased estimate of the local gradient, whereas, in Jiang
et al. (2022), it is biased. Nevertheless, both references assume access to an
exact gradient oracle.

In Tang et al. (2021), they develop two algorithms for nonconvex multi-agent
optimization. One is a gradient-tracking algorithm based on a noise-free
2d-point estimator of the gradient, achieving a rate of $O(\frac{1}{k})$ with
smoothness assumptions and a linear rate for an extra $\nu$-gradient dominated
objective assumption. The other is based on a 2-point unbiased estimator
without global gradient tracking and achieves a rate of
$O(\frac{1}{\sqrt{k}}\log k)$ under Lipschitz continuity and smoothness
conditions and $O(\frac{1}{k})$ under an extra gradient dominated function
structure. Both these gradient estimators are unbiased with respect to the
exact gradient and even have a vanishing variance.

All Pu and Nedić (2018); Xin et al. (2019); Pu (2020); Lu et al. (2019) assume
access to local stochastic first-order oracles where the gradient is unbiased
and with a bounded variance. In the first three, they additionally assume
smooth and strongly-convex local objectives, and they all accomplish a linear
convergence rate under a constant step size. Pu and Nedić (2018) proposes a
distributed stochastic gradient-tracking method (DSGT) and a gossip-like
stochastic gradient-tracking method (GSGT) where at each round, each agent
wakes up with a certain probability. Further, in Pu and Nedić (2018), when the
step-size is diminishing, the convergence rate is that of $O(\frac{1}{k})$.
Xin et al. (2019) employs a gradient-tracking algorithm for agents
communicating over a strongly-connected graph. Pu (2020) introduces a robust
gradient-tracking method (R-Push-Pull) in the context of noisy information
exchange between agents and with a directed network topology. In Lu et al.
(2019), they propose a gradient-tracking based nonconvex stochastic
decentralized (GNSD) algorithm for nonconvex optimization problems in machine
learning, and they fulfill a convergence rate of $O(\frac{1}{\sqrt{k}})$ under
constant step size.

###  1.4 Contributions

In this paper, we consider smooth and convex local objectives, and we extend
the gradient-tracking algorithm to the case where we do not have an estimation
of the gradient. Under the realistic assumption that the agent only has access
to a single noisy function value at each time, without necessarily knowing the
form of this function, we propose a one-point estimator in a stochastic
framework. Naturally, one-point estimators are biased with respect to the true
gradient and suffer from high variance (Liu et al., 2020) hence they do not
match the assumptions for convergence presented in Tang et al. (2021); Pu and
Nedić (2018); Xin et al. (2019); Pu (2020); Lu et al. (2019). However, in this
work, we analyze and indeed prove the convergence of the algorithm with a
biased estimate. We also consider that a stochastic process influences the
objective function from one iteration to the other, which is not the case in
the aforementioned references. We then study the convergence rate and prove
that a rate of $O(\frac{1}{\sqrt{k}})$ after a sufficient number of iterations
$k>K_{2}$ is attainable for smooth and strongly convex objectives. This rate
satisfies the lower bounds achieved by its centralized counterparts in the
same derivative-free setting (Duchi et al., 2015; Jamieson et al., 2012;
Shamir, 2013). Finally, we show that a regret bound of $O(\sqrt{k})$ is
achieved for this algorithm.

###  1.5 Notation

In all that follows, vectors are column-shaped unless defined otherwise and
$\mathbf{1}$ denotes the vector of all entries equal to $1$. For two vectors
$a$, $b$ of the same dimension, $\langle a,b\rangle$ is the inner product. For
two matrices $A$, $B\in\mathbb{R}^{n\times d}$, we define

$\langle A,B\rangle=\sum_{i=1}^{n}\langle A_{i},B_{i}\rangle$

where $A_{i}$ (respectively, $B_{i}$) represents the $i$-th row of $A$
(respectively, $B$). $\|.\|$ denotes the $2$-norm for vectors and the
Frobenius norm for matrices.

We assume that each agent $i$ maintains a local copy $x_{i}\in\mathbb{R}^{d}$
of the decision variable and another auxiliary variable
$y_{i}\in\mathbb{R}^{d}$ and each agent’s local function is subject to the
stochastic variable $S_{i}\in\mathbb{R}^{m}$. At iteration $k$, the respective
values are denoted as $x_{i,k}$, $y_{i,k}$, and $S_{i,k}$. Bold notations
denote the concatenated version of the variables, i.e.,

$\begin{split}\mathbf{x}:=[x_{1},x_{2},\ldots,x_{n}]^{T},\;\mathbf{y}:=[y_{1},y_{2},\ldots,y_{n}]^{T}\in\mathbb{R}^{n\times
d},\;\text{and}\;\mathbf{S}:=[S_{1},S_{2},\ldots,S_{n}]^{T}\in\mathbb{R}^{n\times
m}.\end{split}$

We then define the means of the previous two variables as
$\bar{x}:=\frac{1}{n}\mathbf{1}^{T}\mathbf{x}$ and
$\bar{y}:=\frac{1}{n}\mathbf{1}^{T}\mathbf{y}\in\mathbb{R}^{1\times d}$.

We define the gradient of $F_{i}$ at the local variable $\nabla
F_{i}(x_{i})\in\mathbb{R}^{d}$ and its Hessian matrix
$\nabla^{2}F_{i}(x_{i})\in\mathbb{R}^{d\times d}$ and we let

$\nabla F(\mathbf{x}):=[\nabla F_{1}(x_{1}),\nabla F_{2}(x_{2}),\ldots,\nabla
F_{n}(x_{n})]^{T}\in\mathbb{R}^{n\times d}$

and

$\begin{split}\mathbf{g}:=g(\mathbf{x},\mathbf{S}):=[g_{1}(x_{1},S_{1}),g_{2}(x_{2},S_{2}),\ldots,g_{n}(x_{n},S_{n})]^{T}\in\mathbb{R}^{n\times
d}.\end{split}$

We define its mean
$\bar{g}:=\frac{1}{n}\mathbf{1}^{T}\mathbf{g}\in\mathbb{R}^{1\times d}$ and we
denote each agent’s gradient estimate at time $k$ by
$g_{i,k}=g_{i}(x_{i,k},S_{i,k})$.

###  1.6 Basic Assumptions

In this subsection, we introduce the fundamental assumptions that ensure the
performance of the 1P-DSGT algorithm.

###### Assumption 1

(on the graph) The topology of the network is represented by the graph
$\mathcal{G}=(\mathcal{N},\mathcal{E})$ where the edges in
$\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ represent communication
links. The graph $\mathcal{G}$ is undirected, i.e., $(i,j)\in\mathcal{E}$ iff
$(j,i)\in\mathcal{E}$, and connected (there exists a path of links between any
two agents).

$W=[w_{ij}]\in\mathbb{R}^{n\times n}$ denotes the agents’ coupling matrix,
where agents $i$ and $j$ are connected iff $w_{ij}=w_{ji}>0$
($w_{ij}=w_{ji}=0$ otherwise). $W$ is a nonnegative matrix and doubly
stochastic, i.e., $W\mathbf{1}=\mathbf{1}$ and
$\mathbf{1}^{T}W=\mathbf{1}^{T}$. All diagonal elements $w_{ii}$ are strictly
positive.

###### Assumption 2

(on the objective function) We assume the existence and the continuity of both
$\nabla F_{i}(x)$ and $\nabla^{2}F_{i}(x)$. Let $x^{*}\in\mathbb{R}^{d}$
denote the solution of the problem (1), then $\nabla F_{i}(x^{*})=0$ and
$\det(\nabla^{2}F_{i}(x^{*}))>0$, $\forall i\in\mathcal{N}$. To insure the
existence of $x^{*}$, we let the objective function be strictly convex, i.e.,

$\langle x-x^{*},\nabla\mathcal{F}(x)\rangle\geq 0,\forall
x\in\mathbb{R}^{d}.$ (2)

We further assume the boundedness of the local Hessian where there exists a
constant $\alpha_{1}\in\mathbb{R}^{+}$ such that

$\|\nabla^{2}F_{i}(x)\|_{2}\leq\alpha_{1},\;\forall i\in\mathcal{N},$

where here it suffices to use the Euclidean norm for matrices (keeping in mind
for a matrix $A$, $\|A\|_{2}\leq\|A\|_{F}$).

###### Assumption 3

(on the local functions) All local functions $x\longmapsto f_{i}(x,S)$ are
Lipschitz continuous with Lipschitz constant $L_{S}$,

$\|f_{i}(x,S)-f_{i}(x^{\prime},S)\|\leq L_{S}\|x-x^{\prime}\|,\;\forall
i\in\mathcal{N}.$

In addition, we assume $\mathbb{E}_{S}f_{i}(x,S)<\infty$, $\forall
i\in\mathcal{N}$, to guarantee the boundedness of the objective
$\mathcal{F}(x)$.

###### Assumption 4

(on the additive noise) $\zeta_{i,k}$ is a zero-mean uncorrelated noise with
bounded variance, where $E(\zeta_{i,k})=0$,
$E(\zeta_{i,k}^{2})=\alpha_{4}<\infty$, $\forall i\in\mathcal{N}$, and
$E(\zeta_{i,k}\zeta_{j,k})=0$ if $i\neq j$.

###### Lemma 1

(Qu and Li, 2018) Let $\rho_{w}$ be the spectral norm
$W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}$. When Assumption 1 is satisfied, we
have the following inequality

$\|W\omega-\mathbf{1}\bar{\omega}\|\leq\rho_{w}\|\omega-\mathbf{1}\bar{\omega}\|,\;\forall\omega\in\mathbb{R}^{n\times
d}\text{ and }\bar{\omega}=\frac{1}{n}\mathbf{1}^{T}\omega,$

and $\rho_{w}<1$.

###### Lemma 2

Define $h(\mathbf{x}):=\frac{1}{n}\mathbf{1}^{T}\nabla
F(\mathbf{x})\in\mathbb{R}^{1\times d}$. Due to the boundedness of the second
derivative in Assumption 2, the objective function is thus $L$-smooth and we
have

$\|\nabla\mathcal{F}(\bar{x}_{k})-h(\mathbf{x}_{k})\|\leq\frac{L}{\sqrt{n}}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|.$

##  2 Distributed Stochastic Gradient-Tracking Method

We propose to employ a zero-order one-point estimate of the gradient subject
to the stochastic process $S$ and an additive noise $\zeta$ while a stochastic
perturbation and a step size are introduced, and we assume that each agent can
perform this estimation at each iteration. To elaborate, let $g_{i,k}$ denote
the aforementioned gradient estimate for agent $i$ at time $k$, then we define
it as

$\begin{split}g_{i,k}&=\Phi_{i,k}\tilde{f}_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})\\\
&=\Phi_{i,k}(f_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})+\zeta_{i,k}),\end{split}$
(3)

where $\gamma_{k}>0$ is a vanishing step size and
$\Phi_{i,k}\in\mathbb{R}^{d}$ is a perturbation randomly and independently
generated by each agent $i$. $g_{i,k}$ is in fact a biased estimation of the
gradient $\nabla F_{i}(x_{i,k})$ and the algorithm can converge under the
condition that all parameters are properly chosen. For clarification on the
form of this bias and more on the properties of this estimate, refer to
Appendix A.

###  2.1 The 1P-DSGT Algorithm

The following distributed stochastic gradient-tracking method is considered in
this part making use of the gradient estimate presented in (3).

Every agent $i$ initializes its variables with an arbitrary value $x_{i,0}$
and $y_{i,0}=g_{i,0}$. Then, at each time $k\in\mathbb{N}$, agent $i$ updates
its variables independently according to the following 3 steps:

$\begin{split}&x_{i,k+1}=\sum_{j=1}^{n}w_{ij}(x_{j,k}-\alpha_{k}y_{j,k})\\\
&\text{perform the action: }x_{i,k+1}+\gamma_{k+1}\Phi_{i,k+1}\\\
&y_{i,k+1}=\sum_{j=1}^{n}w_{ij}y_{j,k}+g_{i,k+1}-g_{i,k}\end{split}$ (4)

where $\alpha_{k}>0$ is a vanishing step size. Algorithm (4) can then be
written in the following compact matrix form for clarity of analysis:

$\begin{split}&\mathbf{x}_{k+1}=W(\mathbf{x}_{k}-\alpha_{k}\mathbf{y}_{k})\\\
&\text{perform the action:
}\mathbf{x}_{k+1}+\gamma_{k+1}\mathbf{\Phi}_{k+1}\\\
&\mathbf{y}_{k+1}=W\mathbf{y}_{k}+\mathbf{g}_{k+1}-\mathbf{g}_{k}\end{split}$
(5)

where $\mathbf{\Phi}_{k}\in\mathbb{R}^{n\times d}$ is defined as
$\mathbf{\Phi}_{k}=[\Phi_{1,k},\Phi_{2,k},\ldots,\Phi_{n,k}]^{T}$.

As is evident from the update of the variables, the exchange between agents is
limited to neighboring nodes, and it encompasses the decision variable
$\mathbf{x}_{k+1}$ and the auxiliary variable $\mathbf{y}_{k+1}$.

By construction of Algorithm (5), we note that the mean of the auxiliary
variable $\mathbf{y}_{k}$ is equal to that of the gradient estimate
$\mathbf{g}_{k}$ at every iteration $k$ since $\mathbf{y}_{0}=\mathbf{g}_{0}$,
and by recursion, we obtain
$\bar{y}_{k}=\frac{1}{n}\mathbf{1}^{T}\mathbf{g}_{k}=\bar{g}_{k}$.

###### Assumption 5

(on the step-sizes) Both $\alpha_{k}$ and $\gamma_{k}$ vanish to $0$ as
$k\rightarrow\infty$, and satisty the the following sums

$\sum_{k=1}^{\infty}\alpha_{k}\gamma_{k}=\infty,\;\text{and}\;\sum_{k=1}^{\infty}\alpha_{k}^{2}<\infty.$

###### Assumption 6

(on the random perturbation) Let
$\Phi_{i,k}=(\phi_{i,k}^{1},\phi_{i,k}^{2},\ldots,\phi_{i,k}^{d})^{T}$.

Each agent $i$ chooses its $\Phi_{i,k}$ vector independently from other agents
$j\neq i$. In addition, the elements of $\Phi_{i,k}$ are assumed i.i.d with
$\mathbb{E}(\phi_{i,k}^{d_{1}}\phi_{i,k}^{d_{2}})=0$ for $d_{1}\neq d_{2}$ and
there exists $\alpha_{2}>0$ such that
$\mathbb{E}(\phi_{i,k}^{d_{j}})^{2}=\alpha_{2}$, $\forall{d_{j}}$, $\forall
i$. We further assume that there exists a constant $\alpha_{3}>0$ where
$\|\Phi_{i,k}\|\leq\alpha_{3}$, $\forall i.$

###### Example 1

One example is to take $\alpha_{k}=\alpha_{0}(k+1)^{-\upsilon_{1}}$ and
$\gamma_{k}=\gamma_{0}(k+1)^{-\upsilon_{2}}$ with the constants $\alpha_{0}$,
$\gamma_{0}$, $\upsilon_{1}$, $\upsilon_{2}$ $\in\mathbb{R}^{+}$. As
$\sum_{k=1}^{\infty}\alpha_{k}\gamma_{k}$ diverges for
$\upsilon_{1}+\upsilon_{2}\leq 1$ and $\sum_{k=1}^{\infty}\alpha_{k}^{2}$
converges for $\upsilon_{1}>0.5$, we can find pairs of $\upsilon_{1}$ and
$\upsilon_{2}$ so that Assumption 5 is satisfied.

To achieve the conditions in Assumption 6, we can choose the probability
distribution of $\phi_{i,k}^{d_{j}}$ to be the symmetrical Bernoulli
distribution where
$\phi_{i,k}^{d_{j}}\in\\{-\frac{1}{\sqrt{d}},\frac{1}{\sqrt{d}}\\}$ with
$\mathbb{P}(\phi_{i,k}^{d_{j}}=-\frac{1}{\sqrt{d}})=\mathbb{P}(\phi_{i,k}^{d_{j}}=\frac{1}{\sqrt{d}})=0.5$,
$\forall d_{j}$, $\forall i$.

###  2.2 Convergence Results

In this part, we analyze the asymptotic behavior of Algorithm (5). We start
the analysis by defining $\mathcal{H}_{k}$ as the history sequence
$\\{x_{0},y_{0},S_{0},\ldots,x_{k-1},y_{k-1},S_{k-1},x_{k}\\}$ and denoting by
$\mathbb{E}[.|\mathcal{H}_{k}]$ as the conditional expectation given
$\mathcal{H}_{k}$.

We define $\tilde{g}_{k}$ to be the expected value of $\bar{g}_{k}$ with
respect to all the stochastic terms $S,\Phi,\zeta$ given $\mathcal{H}_{k}$,
i.e.,

$\tilde{g}_{k}=\mathbb{E}_{S,\Phi,\zeta}[\bar{g}_{k}|\mathcal{H}_{k}]$

In what follows, we use
$\tilde{g}_{k}=\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]$ for shorthand
notation.

We define the error $e_{k}$ to be the difference between the value of a single
realization of $\bar{g}_{k}$ and its conditional expectation $\tilde{g}_{k}$,
i.e.,

$e_{k}=\bar{g}_{k}-\tilde{g}_{k},$

where $e_{k}$ can be seen as a stochastic noise. The following lemma
describing the vanishing of the stochastic noise is essential for our main
result.

###### Lemma 3

If all Assumptions $2-6$ hold and $\|\bar{x}_{k}\|<\infty$ almost surely, then
for any constant $\nu>0$, we have

$\lim_{K\rightarrow\infty}\mathbb{P}(\sup_{K^{\prime}\geq
K}\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|\geq\nu)=0,\;\forall\nu>0.$

Proof: See Appendix B.

For any integer $k\geq 0$, we define the divergence, or the error between the
average action taken by the agents $\bar{x}_{k}$ and the optimal solution
$x^{*}$ as

$d_{k}=\|\bar{x}_{k}-x^{*}\|^{2}.$ (6)

The following theorem describes the main convergence result.

###### Theorem 4

If all Assumptions $1-6$ hold and $\|\mathbf{x}_{k}\|<\infty$ almost surely,
then as $k\rightarrow\infty$, $d_{k}\rightarrow 0$, $\bar{x}_{k}\rightarrow
x^{*}$, and $x_{i,k}\rightarrow\bar{x}_{k}$ for all $i\in\mathcal{N}$ almost
surely by applying the Algorithm.

Proof: See Appendix C.

###  2.3 Convergence Rate

This part deals with how fast the expected divergence vanishes to find the
proposed algorithm’s expected convergence rate. To do so, we define the
expected divergence to be $D_{k}=\mathbb{E}[\|\bar{x}_{k}-x^{*}\|^{2}]$. The
goal is to bound this divergence from above by sequences whose convergence
rate is known. The analysis is highly associated with the parameters
$\alpha_{k}$ and $\gamma_{k}$ that play a significant role in determining this
upper bound. Hence, in what follows, the analysis starts with a general form
of $\alpha_{k}$ and $\gamma_{k}$, then a particular case is considered.

####  2.3.1 General Form of $\alpha_{k}$ and $\gamma_{k}$

We start by considering an additional assumption on the objective function for
what follows.

###### Assumption 7

Let $\mathcal{F}(x)$ be strongly convex, then there exists $\lambda>0$ such
that

$\langle\nabla\mathcal{F}(x),x-x^{*}\rangle\geq\lambda\|x-x^{*}\|^{2}.$

Our main result regarding the convergence rate is summarized in the following
theorem.

###### Theorem 5

Let $R=\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}$,
$\delta_{k}=(\frac{1+\rho_{w}^{2}}{2})^{k}$, and
$\beta_{k}=\sum_{j=1}^{k}\delta_{j}\alpha_{k-j}^{2}$. Let $\bar{M}$ denote the
upper bound of $\mathbb{E}[\|\bar{g}_{k}\|^{2}]$ and $G$ that of
$\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|$ and define
$\bar{G}=\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}G$ (Refer to Appendix A, Appendix
C.1, and Appendix C.2 for proof of boundedness).

Next, we define the constants $A=\lambda\alpha_{2}$,
$B=\alpha_{1}\alpha_{3}^{3}$, $C=\alpha_{2}\frac{L^{2}}{\lambda n}$, and

$K_{0}=\arg\min_{A\alpha_{k}\gamma_{k}<1}k.$

We finally define the following parameters:

$\begin{matrix}\kappa_{k}=\frac{1-(\frac{\gamma_{k+1}}{\gamma_{k}})^{2}}{\alpha_{k}\gamma_{k}},&\sigma_{1}=\underset{k\geq
K_{0}}{\max}\;\kappa_{k},&\sigma_{2}=\underset{k\geq
K_{0}}{\max}\;\frac{\delta_{k}}{\gamma_{k}^{2}},&\sigma_{3}=\underset{k\geq
K_{0}}{\max}\;\frac{\beta_{k}}{\gamma_{k}^{2}},&\sigma_{4}=\underset{k\geq
K_{0}}{\max}\;\frac{\alpha_{k}}{\gamma_{k}^{3}},\\\
\tau_{k}=\frac{1-\frac{\alpha_{k+1}\gamma_{k+1}^{-1}}{\alpha_{k}\gamma_{k}^{-1}}}{\alpha_{k}\gamma_{k}},&\sigma_{5}=\underset{k\geq
K_{0}}{\max}\;\tau_{k},&\sigma_{6}=\underset{k\geq
K_{0}}{\max}\;\sqrt{\frac{\gamma_{k}^{3}}{\alpha_{k}}},&\sigma_{7}=\underset{k\geq
K_{0}}{\max}\;\frac{\gamma_{k}}{\alpha_{k}}\delta_{k},&\sigma_{8}=\underset{k\geq
K_{0}}{\max}\;\frac{\gamma_{k}}{\alpha_{k}}\beta_{k}.\end{matrix}$ (7)

If $\kappa_{k}<A$ for any $k\geq K_{0}$, then

$D_{k}\leq\varsigma_{1}^{2}\gamma_{k}^{2},\;\forall k\geq K_{0},$ (8)

with

$\varsigma_{1}\geq\max\Bigg{\\{}\frac{\sqrt{D_{K_{0}}}}{\gamma_{K_{0}}},\frac{B}{2(A-\sigma_{1})}+\sqrt{\Big{(}\frac{B}{2(A-\sigma_{1})}\Big{)}^{2}+\frac{CR\sigma_{2}+C\bar{G}\sigma_{3}+\bar{M}\sigma_{4}}{(A-\sigma_{1})}}\Bigg{\\}}.$
(9)

If $\tau_{k}<A$ for any $k\geq K_{0}$, then

$D_{k}\leq\varsigma_{2}^{2}\frac{\alpha_{k}}{\gamma_{k}},\;\forall k\geq
K_{0},$ (10)

with

$\varsigma_{2}\geq\max\Bigg{\\{}\sqrt{\frac{D_{K_{0}}\gamma_{K_{0}}}{\alpha_{K_{0}}}},\frac{B\sigma_{6}}{2(A-\sigma_{5})}+\sqrt{\bigg{(}\frac{B\sigma_{6}}{2(A-\sigma_{5})}\bigg{)}^{2}+\frac{C(R\sigma_{7}+\bar{G}\sigma_{8})+\bar{M}}{(A-\sigma_{5})}}\Bigg{\\}}.$
(11)

Proof: See Appendix D.1.

####  2.3.2 A Special Case of $\alpha_{k}$ and $\gamma_{k}$

We now consider the special case mentioned in Example 1:

$\alpha_{k}=\alpha_{0}(k+1)^{-\upsilon_{1}}\text{ and
}\gamma_{k}=\gamma_{0}(k+1)^{-\upsilon_{2}},$ (12)

where $0.5<\upsilon_{1}<1$ and $0<\upsilon_{2}\leq 1-\upsilon_{1}$.

Before stating the main result, we consider the following lemma.

###### Lemma 6

(Study of $\beta_{k}$) Let $K_{1}$ be such that

$K_{1}=\arg\min_{\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{k}\leq\alpha_{k}^{2}}k.$
(13)

Then, the convergence rate of $\beta_{k}$ is at least that of order
$\frac{1}{k^{3\upsilon_{1}-1}}$.

Proof: See Appendix D.4

We then let $K_{2}$ be such that $K_{2}\geq\max\\{K_{0},K_{1}\\}$ and state
our next theorem.

###### Theorem 7

Let $\alpha_{k}$ and $\gamma_{k}$ have the forms given in (12), if
$\alpha_{0}\gamma_{0}\geq\max\\{2\upsilon_{2},\upsilon_{1}-\upsilon_{2}\\}/A$,
then we can say that there exists $\Upsilon<\infty$, where

$D_{k}\leq\Upsilon(k+1)^{-\min\\{2\upsilon_{2},\upsilon_{1}-\upsilon_{2}\\}},\;\forall
k\geq K_{2}.$

Proof: See Appendix D.5.

The parameters clearly affect the upper bound of the convergence rate or rate
of expected divergence decay in Theorem 7. As it is evident that

$\max\\{2\upsilon_{2},\upsilon_{1}-\upsilon_{2}\\}\leq 0.5,$

the best choice is when equality holds for $\upsilon_{1}=0.75$ and
$\upsilon_{2}=0.25$. With the sufficient condition on the parameters in
Theorem 7, we can finally state that our algorithm converges with a rate of
$O(\frac{1}{\sqrt{k}})$ after a sufficient number of iterations $k>K_{2}$.

###  2.4 Regret Bound

To further examine the performance of our algorithm, we present the following
theorem on the achieved regret bound.

###### Theorem 8

Let $\alpha_{k}$ and $\gamma_{k}$ have the forms of (12), with
$\upsilon_{1}=0.75$ and $\upsilon_{2}=0.25$. Then, the regret bound is given
by

$\mathbb{E}_{\mathcal{H}_{k}}\bigg{[}\sum_{k=1}^{K}\mathcal{F}(\bar{x}_{k})-\mathcal{F}(x^{*})\bigg{]}\leq\Upsilon
L(\sqrt{K+1}-1).$

Proof: See Appendix E.

##  3 Numerical Results

In this section, we give numerical examples to illustrate the performance of
the 1P-DSGT algorithm. We compare it with a general DSGT algorithm based on an
unbiased estimator with bounded variance. For this unbiased estimator, we
calculate the exact gradient and add white noise to it. The network topology
is a connected Erdős-Rényi random graph with a probability of $0.3$.

We consider a logistic classification problem to classify $m$ images of the
two digits, labeled as $y_{ij}=+1$ or $-1$ from the MNIST data set (LeCun and
Cortes, 2005). Each image, $X_{ij}$, is a $785$-dimensional vector and is
compressed using a lossy autoencoder to become $10$-dimensional, i.e., $d=10$.
The total images are split equally among the agents such that each agent has
$m_{i}=\frac{m}{n}$ images and no access to other ones for privacy
constraints. However, the goal is still to make use of all images and to solve
collaboratively

$min_{\theta\in\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{m}\sum_{j=1}^{m_{i}}\mathbb{E}_{u\sim\mathcal{N}(1,\sigma_{u})}\ln(1+\exp(-u_{ij}y_{ij}.X_{ij}^{T}\theta))+c\|\theta\|^{2},$

while reaching consensus on the decision variable $\theta\in\mathbb{R}^{d}$.
We note here that $u$ models some perturbation on the local querying of every
example to add to the randomization of the communication process.

For the first example, we consider classifying the digits $2$ and $9$ where
$m=11907$ images. There are $n=21$ agents in the network and thus each has a
local batch of $m_{i}=567$ images. We take $\sigma_{u}=0.01$ and let
$\alpha_{k}=0.2(k+1)^{-0.75}$, $\gamma_{k}=1.3(k+1)^{-0.25}$, and $\Phi_{k}\in
1.5\times\\{-\frac{1}{\sqrt{d}},\frac{1}{\sqrt{d}}\\}^{d}$ with equal
probability. Also, every function query is subject to a white noise generated
by the standard normal distribution. For the general DSGT algorithm, we let
the step size to $\alpha_{k}=4(k+1)^{-0.75}$, and we don’t consider the
perturbation on the objective function nor the noise on the objective
function, only the noise on the exact gradient. We let $c=0.1$, and the
initialization be the same for both algorithms, with $\theta_{i,0}$ uniformly
chosen from $[-0.5,0.5]^{d}$, $\forall i\in\mathcal{N}$, per instance. We
finally average the simulations over $50$ instances.

The expected evolution of the loss objective function is present in Figure 1.
Since the DSGT algorithm with vanishing step size converges at a rate of
$O(\frac{1}{k})$ (Pu and Nedić, 2018), the result we obtain is actually
anticipated. Next, at every iteration, we measure the accuracy of the
classification against an independent test set of $2041$ images using the
updated mean vector $\bar{\theta}_{k}=\frac{1}{n}\sum_{i=1}^{n}\theta_{i,k}$
of the local decision variables and we present the results in Figure 2.
1P-DSGT achieves an accuracy of $95.3\%$ whereas DSGT achieves $97.34\%$ at
the final iteration, which are both quite good results and an especially
important one for 1P-DSGT considering the difference in the convergence rate.

Figure 1: Expected loss function evolution of the algorithms 1P-DSGT vs. DSGT
considering images of the digits $2$ and $9$.

Figure 2: Expected accuracy evolution of the algorithms 1P-DSGT vs. DSGT
considering images of the digits $2$ and $9$.

We display the same result for the last two Figures, 3 and 4; only the second
Figure is shown for later iterations to zoom in on the curves. The curves are
those of the evolution of the consensus error, or
$\sum_{i=1}^{n}\|\theta_{i,k}-\theta_{k}\|^{2}$ which is the error between the
local decision variables and their average. For both algorithms, the error
decreases quite fast. Nonetheless, a swift consensus is reached by 1P-DSGT
compared to that of DSGT, which is unpredictable considering all the
randomizations. This is an outstanding result for parallel training and
computation in networked resource nodes, as consensus no longer raises an
issue.

Figure 3: Expected consensus error evolution of the algorithms 1P-DSGT vs.
DSGT considering images of the digits $2$ and $9$.

Figure 4: Expected consensus error evolution of the algorithms 1P-DSGT vs.
DSGT considering images of the digits $2$ and $9$.

The second example is the classification of the images of the digits $3$ and
$7$. We have $m=12396$ and we consider $n=6$ agents, thus $m_{i}=2066$,
$\forall i$. Here, we take $\alpha_{k}=0.15(k+1)^{-0.75}$, and everything else
remains exactly the same as in the previous example. The results are shown in
Figures 5 through 8. We note that 1P-DSGT achieves $92.46\%$ accuracy, and
DSGT achieves $97.44\%$ at the last iteration.

Figure 5: Expected loss function evolution of the algorithms 1P-DSGT vs. DSGT
considering images of the digits $3$ and $7$.

Figure 6: Expected accuracy evolution of the algorithms 1P-DSGT vs. DSGT
considering images of the digits $3$ and $7$.

Figure 7: Expected consensus error evolution of the algorithms 1P-DSGT vs.
DSGT considering images of the digits $3$ and $7$.

Figure 8: Expected consensus error evolution of the algorithms 1P-DSGT vs.
DSGT considering images of the digits $3$ and $7$.

One ending remark we can bring attention to is that the 1P-DSGT algorithm is
clearly more robust against the distortion brought by the compression of the
images. While the DSGT algorithm introduces noticeable fluctuations to all
results, 1P-DSGT demonstrates a steady evolution towards the optimum.

##  4 Conclusion

In this work, we extended the gradient-tracking algorithm to present a
practical solution to a relevant problem with realistic assumptions. A
distributed stochastic gradient-tracking algorithm was studied and proved to
converge with a biased and high variance one-point gradient estimate and a
stochastic perturbation on the objective function. The convergence rate was
proved to be $O(\frac{1}{\sqrt{k}})$ and the regret bound that of
$O(\sqrt{k})$. A numerical application confirmed the success and efficiency of
the algorithm. For further work, we wish to relax some of the assumptions.

##  Appendix Appendix A Estimated Gradient

In this section, we derive the bias of the gradient estimate with respect to
the real gradient of the local objective function. Let

$\breve{g}_{i,k}=\mathbb{E}_{S,\Phi,\zeta}[g_{i,k}|\mathcal{H}_{k}].$

Thus,

$\begin{split}\breve{g}_{i,k}&=\mathbb{E}_{S,\Phi,\zeta}[\Phi_{i,k}(f_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})+\zeta_{i,k})|\mathcal{H}_{k}]\\\
&=\mathbb{E}_{S,\Phi}[\Phi_{i,k}f_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})|\mathcal{H}_{k}]\\\
&=\mathbb{E}_{\Phi}[\Phi_{i,k}F_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k})|\mathcal{H}_{k}].\end{split}$

By Taylor’s theorem and the mean-valued theorem, there exists
$\tilde{x}_{i,k}$ located between $x_{i,k}$ and $x_{i,k}+\gamma_{k}\Phi_{i,k}$
where

$F_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k})=F_{i}(x_{i,k})+\gamma_{k}\langle\Phi_{i,k},\nabla
F_{i}(x_{i,k})\rangle+\frac{\gamma_{k}^{2}}{2}\langle\Phi_{i,k},\nabla^{2}F_{i}(\tilde{x}_{i,k})\Phi_{i,k}\rangle,$

substituting in the previous definition,

$\begin{split}\breve{g}_{i,k}&=F_{i}(x_{i,k})\mathbb{E}_{\Phi}[\Phi_{i,k}]+\gamma_{k}\mathbb{E}_{\Phi}[\Phi_{i,k}\Phi_{i,k}^{T}]\nabla
F_{i}(x_{i,k})+\frac{\gamma_{k}^{2}}{2}\mathbb{E}_{\Phi}[\Phi_{i,k}\Phi_{i,k}^{T}\nabla^{2}F_{i}(\tilde{x}_{i,k})\Phi_{i,k}|\mathcal{H}_{k}]\\\
&=\alpha_{2}\gamma_{k}[\nabla F_{i}(x_{i,k})+b_{i,k}].\end{split}$

Thus, the estimation bias has the form

$\begin{split}b_{i,k}&=\frac{\breve{g}_{i,k}}{\alpha_{2}\gamma_{k}}-\nabla
F_{i}(x_{i,k})\\\
&=\frac{\gamma_{k}}{2\alpha_{2}}\mathbb{E}_{\Phi}[\Phi_{i,k}\Phi_{i,k}^{T}\nabla^{2}F_{i}(\tilde{x}_{i,k})\Phi_{i,k}|\mathcal{H}_{k}].\end{split}$

Let Assumptions 2 and 6 hold. Then, we can bound the bias as

$\begin{split}\|b_{i,k}\|&\leq\frac{\gamma_{k}}{2\alpha_{2}}\mathbb{E}_{\Phi}[\|\Phi_{i,k}\|_{2}\|\Phi_{i,k}^{T}\|_{2}\|\nabla^{2}F_{i}(\tilde{x}_{i,k})\|_{2}\|\Phi_{i,k}\|_{2}|\mathcal{H}_{k}]\\\
&\leq\gamma_{k}\frac{\alpha_{3}^{3}\alpha_{1}}{2\alpha_{2}}.\end{split}$

We can see $\|b_{i,k}\|\rightarrow 0$ as $k\rightarrow\infty$ since
$\gamma_{k}$ is vanishing. We remark that

$\begin{split}\tilde{g}_{k}&=\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\\\
&=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[g_{i,k}|\mathcal{H}_{k}]\\\
&=\frac{1}{n}\sum_{i=1}^{n}\alpha_{2}\gamma_{k}[\nabla
F_{i}(x_{i,k})+b_{i,k}]\\\
&=\alpha_{2}\gamma_{k}[h(\mathbf{x}_{k})+\bar{b}_{k}]\end{split}$ (14)

is also a biased estimator of $h(\mathbf{x}_{k})$ with

$\begin{split}\|\bar{b}_{k}\|&=\|\frac{1}{n}\sum_{i=1}^{n}b_{i,k}\|\\\
&\leq\frac{1}{n}\sum_{i=1}^{n}\|b_{i,k}\|\\\
&\leq\frac{1}{n}\sum_{i=1}^{n}\gamma_{k}\frac{\alpha_{3}^{3}\alpha_{1}}{2\alpha_{2}}\\\
&=\gamma_{k}\frac{\alpha_{3}^{3}\alpha_{1}}{2\alpha_{2}}.\end{split}$ (15)

###### Lemma 9

Let all Assumptions $2-6$ hold and $\|\mathbf{x}_{k}\|<\infty$ almost surely,
then there exists a bounded constant $M>0$, such that
$E[\|\mathbf{g}_{k}\|^{2}]<M$ almost surely.

Proof  $\forall i\in\mathcal{N}$, we have

$\begin{split}\mathbb{E}[\|g_{i,k}\|^{2}|\mathcal{H}_{k}]&=\mathbb{E}[\|\Phi_{i,k}(f_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})+\zeta_{i,k})\|^{2}|\mathcal{H}_{k}]\\\
&=\mathbb{E}[\|\Phi_{i,k}\|^{2}\|f_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})+\zeta_{i,k}\|^{2}|\mathcal{H}_{k}]\\\
&\overset{(a)}{\leq}\alpha_{3}^{2}\mathbb{E}[(f_{i}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})+\zeta_{i,k})^{2}|\mathcal{H}_{k}]\\\
&\overset{(b)}{=}\alpha_{3}^{2}\mathbb{E}[f_{i}^{2}(x_{i,k}+\gamma_{k}\Phi_{i,k},S_{i,k})|\mathcal{H}_{k}]+\alpha_{3}^{2}\alpha_{4}\\\
&\overset{(c)}{\leq}\alpha_{3}^{2}\mathbb{E}[(\|f_{i}(0,S_{i,k})\|+L_{S_{i,k}}\|x_{i,k}+\gamma_{k}\Phi_{i,k}\|)^{2}|\mathcal{H}_{k}]+\alpha_{3}^{2}\alpha_{4}\\\
&\overset{(d)}{\leq}2\alpha_{3}^{2}\mathbb{E}[\mu_{S_{i,k}}^{2}+L_{S_{i,k}}^{2}(\|x_{i,k}\|+\gamma_{k}\alpha_{3})^{2}|\mathcal{H}_{k}]+\alpha_{3}^{2}\alpha_{4}\\\
&\overset{(e)}{=}2\alpha_{3}^{2}(\mu+L^{\prime}(\|x_{i,k}\|+\gamma_{k}\alpha_{3})^{2})+\alpha_{3}^{2}\alpha_{4}\\\
&<\infty,\end{split}$

where $(a)$ is due to Assumption 6, $(b)$ Assumption 4, and $(c)$ Assumption
3. We denote $\|f_{i}(0,S_{i,k})\|=\mu_{S_{i,k}}$ in $(d)$ and the inequality
is due to $\frac{x+y}{2}\leq\sqrt{\frac{x^{2}+y^{2}}{2}}$, $\forall
x,y\in\mathbb{R}$. In $(e)$, $\mu=\mathbb{E}[\mu_{S_{i,k}}^{2}]$ and
$L^{\prime}=\mathbb{E}[L_{S_{i,k}}^{2}]$.

##  Appendix Appendix B Stochastic Noise

To prove Lemma 3, we begin by demonstrating that the sequence
$\\{\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\\}_{K^{\prime}\geq K}$ is a
martingale. Since $\bar{g}_{k}$ and $\bar{g}_{k^{\prime}}$ are independent if
$k\neq k^{\prime}$ and

$\begin{split}\mathbb{E}[e_{k}]=&\mathbb{E}[\bar{g}_{k}-\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]]\\\
=&\mathbb{E}_{\mathcal{H}_{k}}[\mathbb{E}[\bar{g}_{k}-\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]|\mathcal{H}_{k}]]\\\
=&0\end{split}$

by the law of total expectation, the sequence is a martingale. Therefore, for
any constant $\nu>0$, we can state

$\begin{split}\mathbb{P}(\sup_{K^{\prime}\geq
K}\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|\geq\nu)&\overset{(a)}{\leq}\mathbb{E}(\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|^{2})\\\
&=\frac{1}{\nu^{2}}\mathbb{E}(\sum_{k=K}^{K^{\prime}}\sum_{k^{\prime}=K}^{K^{\prime}}\alpha_{k}\alpha_{k^{\prime}}\langle
e_{k},e_{k^{\prime}}\rangle)\\\
&\overset{(b)}{=}\frac{1}{\nu^{2}}\mathbb{E}(\sum_{k=K}^{K^{\prime}}\|\alpha_{k}e_{k}\|^{2})\\\
&\leq\frac{1}{\nu^{2}}\sum_{k=K}^{\infty}\mathbb{E}(\alpha_{k}^{2}\|\bar{g}_{k}-\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\|^{2})\\\
&=\frac{1}{\nu^{2}}\sum_{k=K}^{\infty}\alpha_{k}^{2}\mathbb{E}(\|\bar{g}_{k}\|^{2})-\mathbb{E}_{\mathcal{H}_{k}}(\|\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\|^{2})\\\
&\leq\frac{1}{\nu^{2}}\sum_{k=K}^{\infty}\alpha_{k}^{2}\mathbb{E}(\|\bar{g}_{k}\|^{2})\\\
&\overset{(c)}{\leq}\frac{M}{\nu^{2}}\sum_{k=K}^{\infty}\alpha_{k}^{2},\\\
\end{split}$

where $(a)$ is due to Doob’s martingale inequality (Doob, 1953), $(b)$ is
since $\mathbb{E}[\langle e_{k},e_{k^{\prime}}\rangle]=0$ for any $k\neq
k^{\prime}$, and $(c)$ is by Lemma 9.

Since $M$ is a bounded constant and
$\lim_{K\rightarrow\infty}\sum_{k=K}^{\infty}\alpha_{k}^{2}=0$ by Assumption
5, we get
$\lim_{K\rightarrow\infty}\frac{M}{\nu^{2}}\sum_{k=K}^{\infty}\alpha_{k}^{2}=0$
for any bounded constant $\nu$. Hence, the probability that
$\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|\geq\nu$ also vanishes as
$K\rightarrow\infty$, which concludes the proof.

##  Appendix Appendix C Proof of Convergence

We start by stating the following lemma that will be useful for the proof of
convergence.

###### Lemma 10

If all Assumptions $1-6$ hold and $\|\mathbf{x}_{k}\|<\infty$ almost surely,
then
$\lim_{k\rightarrow\infty}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}=0$. In
fact, we have
$\sum_{k=0}^{\infty}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}<\infty$ as
well as

$\sum_{k=0}^{\infty}\gamma_{k}\alpha_{k}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|<\infty$

almost surely.

Proof: See Appendix C.2.

###  Appendix C.1 Proof of Theorem 4

The goal is to write the divergence in terms of its previous term and to prove
that it’s finally vanishing. We know that
$\bar{x}_{k+1}=\bar{x}_{k}-\alpha_{k}\bar{g}_{k}$. With this equation, the
divergence at time $k+1$ can be written as

$\begin{split}d_{k+1}&=\|\bar{x}_{k+1}-x^{*}\|^{2}\\\
&=\|\bar{x}_{k}-\alpha_{k}\bar{g}_{k}-x^{*}\|\\\
&=\|\bar{x}_{k}-x^{*}\|^{2}-2\alpha_{k}\langle\bar{x}_{k}-x^{*},\bar{g}_{k}-\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]+\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle+\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}\\\
&=\|\bar{x}_{k}-x^{*}\|^{2}-2\alpha_{k}\langle\bar{x}_{k}-x^{*},\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle-2\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle+\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}\\\
&=\|\bar{x}_{k}-x^{*}\|^{2}-2\alpha_{k}\langle\bar{x}_{k}-x^{*},\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle-2\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle+\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}\\\
&\overset{(a)}{=}d_{k}-2\alpha_{2}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},h(\mathbf{x}_{k})+\bar{b}_{k}\rangle-2\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle+\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}\\\
&=d_{k}-2\alpha_{2}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle+2\alpha_{2}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})-h(\mathbf{x}_{k})\rangle\\\
&\hskip
14.22636pt-2\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle+\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}\\\
&\overset{(b)}{\leq}d_{k}-2\alpha_{2}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle+\frac{2\alpha_{2}L\gamma_{k}\alpha_{k}}{\sqrt{n}}\|\bar{x}_{k}-x^{*}\|\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|\\\
&\hskip
14.22636pt-2\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle+\alpha_{k}^{2}\|\bar{g}_{k}\|^{2},\end{split}$

where $(a)$ is due to (14) and $(b)$ is due to Lemma 2. By recursion, we have

$\begin{split}d_{K+1}&\leq
d_{0}-2\alpha_{2}\sum_{k=0}^{K}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle+\frac{2\alpha_{2}L}{\sqrt{n}}\sum_{k=0}^{K}\gamma_{k}\alpha_{k}\|\bar{x}_{k}-x^{*}\|\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|\\\
&\hskip
14.22636pt-2\sum_{k=0}^{K}\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle+\sum_{k=0}^{K}\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}.\\\
\end{split}$ (16)

By Lemma 3, we have
$\lim_{K\rightarrow\infty}\|\sum_{k=0}^{K}\alpha_{k}e_{k}\|<\infty$ almost
surely. Since
$\|\bar{x}_{k}-x^{*}\|=\|\frac{1}{n}\sum_{i=1}^{n}(x_{i,k}-x^{*})\|\leq\frac{1}{n}\sum_{i=1}^{n}\|x_{i,k}-x^{*}\|<\infty$
almost surely, hence

$\lim_{K\rightarrow\infty}\|\sum_{k=0}^{K}\alpha_{k}\langle\bar{x}_{k}-x^{*},e_{k}\rangle\|<\infty.$
(17)

From Lemma 9, we have

$\|\bar{g}_{k}\|^{2}=\|\frac{1}{n}\sum_{i=1}^{n}g_{i,k}\|^{2}\leq\frac{n}{n^{2}}\sum_{i=1}^{n}\|g_{i,k}\|^{2}=\frac{1}{n}\sum_{i=1}^{n}\|g_{i,k}\|^{2}<\infty,\;\text{almost
surely.}$ (18)

Then, by Assumption 5,

$\lim_{K\rightarrow\infty}\sum_{k=0}^{K}\alpha_{k}^{2}\|\bar{g}_{k}\|^{2}<\infty.$
(19)

As stated in Lemma 10, we have
$\sum_{k=0}^{\infty}\gamma_{k}\alpha_{k}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|<\infty$,
adding to $\|\bar{x}_{k}-x^{*}\|<\infty$ almost surely, then

$\lim_{K\rightarrow\infty}\sum_{k=0}^{K}\gamma_{k}\alpha_{k}\|\bar{x}_{k}-x^{*}\|\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|<\infty.$
(20)

From the above equations (16)-(20), we conclude that there exists $D$ such
that $d_{K+1}\leq D+z_{K}$, with

$z_{K}=-2\alpha_{2}\sum_{k=0}^{K}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle.$
(21)

Carefully examine (15), we can say that there exists $K_{b}$, such that for
$k\geq K_{b}$, $\bar{b}_{k}$ becomes very small where we can find an arbitrary
small positive value $\epsilon_{b}$ such that,

$\|\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\|\geq(1-\epsilon_{b})\|\nabla\mathcal{F}(\bar{x}_{k})\|,\;\forall
k\geq K_{b},$ (22)

leading to
$-\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle\leq
0$, due to the convexity of $\mathcal{F}$ in (2).

Consequently, for any big $K$, $0\leq d_{K+1}<\infty$ and the limit
$\lim_{K\rightarrow\infty}d_{K+1}=\bar{d}$ exists.

Thus, there are 2 cases: $\bar{d}>0$ or $\bar{d}=0$. Assume hypothesis H1)
$\bar{d}>0$ to be valid, i.e., $\bar{x}_{k}$ does not converge to $x^{*}$,
then $\forall\epsilon_{h}>0$, $\exists K_{h}$ such that

$-\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})\rangle<-\epsilon_{h},\;\forall
k\geq K_{h}.$

From (21) and (22), we get that $\forall k\geq K_{m}=\max\\{K_{b},K_{h}\\}$,

$-\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle<-\epsilon_{h}(1-\epsilon_{b}),$

implying

$\begin{split}\lim_{K\rightarrow\infty}-\sum_{k=K_{m}}^{K}\gamma_{k}\alpha_{k}\langle\bar{x}_{k}-x^{*},\nabla\mathcal{F}(\bar{x}_{k})+\bar{b}_{k}\rangle<-\epsilon_{h}(1-\epsilon_{b})\lim_{K\rightarrow\infty}\sum_{k=K_{m}}^{K}\gamma_{k}\alpha_{k}<-\infty\end{split}$

since $\sum\gamma_{k}\alpha_{k}$ diverges by Assumption 5. As a result, we get
$z_{K}<-\infty$ and $d_{K+1}<-\infty$. However, by definition in (6),
$d_{K+1}>0$. Accordingly, the hypothesis H1 cannot be true and the case
$\bar{d}=0$ is the valid one. We conclude that
$\lim_{k\rightarrow\infty}d_{k}=0$,
$\lim_{k\rightarrow\infty}\nabla\mathcal{F}(\bar{x}_{k})=0$, and
$\lim_{k\rightarrow\infty}\bar{x}_{k}=x^{*}$ almost surely.

###  Appendix C.2 Proof of Lemma 10

We start by replacing the variables with their algorithmic updates.

$\begin{split}\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|^{2}&=\|W\mathbf{x}_{k}-\alpha_{k}W\mathbf{y}_{k}-\mathbf{1}\bar{x}_{k}+\alpha_{k}\mathbf{1}\bar{y}_{k}\|^{2}\\\
&=\|W\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}-2\alpha_{k}\langle
W\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k},W\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\rangle+\alpha_{k}^{2}\|W\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}\\\
&\overset{(a)}{\leq}\|W\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}+\alpha_{k}[\frac{1-\rho_{w}^{2}}{2\rho_{w}^{2}\alpha_{k}}\|W\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}+\frac{2\rho_{w}^{2}\alpha_{k}}{1-\rho_{w}^{2}}\|W\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}]\\\
&\hskip
14.22636pt+\alpha_{k}^{2}\|W\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}\\\
&\overset{(b)}{\leq}\rho_{w}^{2}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}+\rho_{w}^{2}\alpha_{k}[\frac{1-\rho_{w}^{2}}{2\rho_{w}^{2}\alpha_{k}}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}+\frac{2\rho_{w}^{2}\alpha_{k}}{1-\rho_{w}^{2}}\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}]\\\
&\hskip
14.22636pt+\rho_{w}^{2}\alpha_{k}^{2}\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}\\\
&=\frac{1+\rho_{w}^{2}}{2}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}+\alpha_{k}^{2}\frac{(1+\rho_{w}^{2})\rho_{w}^{2}}{1-\rho_{w}^{2}}\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}\end{split}$
(23)

where $(a)$ is by $-2\epsilon\times\frac{1}{\epsilon}\langle
a,b\rangle=-2\langle\epsilon
a,\frac{1}{\epsilon}b\rangle\leq\epsilon^{2}\|a\|^{2}+\frac{1}{\epsilon^{2}}\|b\|^{2}$
and $(b)$ is by Lemma 1. By induction, we have

$\begin{split}\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|^{2}&\leq(\frac{1+\rho_{w}^{2}}{2})^{k+1}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}+\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}\sum_{j=0}^{k}(\frac{1+\rho_{w}^{2}}{2})^{j+1}\alpha_{k-j}^{2}\|\mathbf{y}_{k-j}-\mathbf{1}\bar{y}_{k-j}\|^{2}.\end{split}$
(24)

Since $\sqrt{a+b}<\sqrt{a}+\sqrt{b}$,

$\begin{split}\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|&\leq(\frac{1+\rho_{w}^{2}}{2})^{\frac{k+1}{2}}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|+\sqrt{\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}}\sum_{j=0}^{k}(\frac{1+\rho_{w}^{2}}{2})^{\frac{j+1}{2}}\alpha_{k-j}\|\mathbf{y}_{k-j}-\mathbf{1}\bar{y}_{k-j}\|.\end{split}$
(25)

By repeatedly replacing the variables with the algorithm’s iterations, we see
that

$\begin{split}\mathbf{y}_{k+1}&=W\mathbf{y}_{k}+\mathbf{g}_{k+1}-\mathbf{g}_{k}\\\
&=W(W\mathbf{y}_{k-1}+\mathbf{g}_{k}-\mathbf{g}_{k-1})+\mathbf{g}_{k+1}-\mathbf{g}_{k}\\\
&=W^{2}\mathbf{y}_{k-1}-W\mathbf{g}_{k-1}+(W-I)\mathbf{g}_{k}+\mathbf{g}_{k+1}\\\
&=W^{2}(W\mathbf{y}_{k-2}+\mathbf{g}_{k-1}-\mathbf{g}_{k-2})-W\mathbf{g}_{k-1}+(W-I)\mathbf{g}_{k}+\mathbf{g}_{k+1}\\\
&=W^{3}\mathbf{y}_{k-2}-W^{2}\mathbf{g}_{k-2}+W(W-I)\mathbf{g}_{k-1}+(W-I)\mathbf{g}_{k}+\mathbf{g}_{k+1}\\\
&=W^{3}(W\mathbf{y}_{k-3}+\mathbf{g}_{k-2}-\mathbf{g}_{k-3})-W^{2}\mathbf{g}_{k-2}+W(W-I)\mathbf{g}_{k-1}+(W-I)\mathbf{g}_{k}+\mathbf{g}_{k+1}\\\
&=W^{4}\mathbf{y}_{k-3}-W^{3}\mathbf{g}_{k-3}+W^{2}(W-I)\mathbf{g}_{k-2}+W(W-I)\mathbf{g}_{k-1}+(W-I)\mathbf{g}_{k}+\mathbf{g}_{k+1}\\\
&=\ldots\\\
&=W^{k+1}\mathbf{y}_{0}-W^{k}\mathbf{g}_{0}+\sum_{j=0}^{k-1}W^{j}(W-I)\mathbf{g}_{k-j}+\mathbf{g}_{k+1}\\\
&=W^{k}(W-I)\mathbf{g}_{0}+\sum_{j=0}^{k-1}W^{j}(W-I)\mathbf{g}_{k-j}+\mathbf{g}_{k+1}\\\
&=\sum_{j=0}^{k}W^{j}(W-I)\mathbf{g}_{k-j}+\mathbf{g}_{k+1},\end{split}$
$\begin{split}\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}&=\sum_{j=0}^{k-1}W^{j}(W-I)\mathbf{g}_{k-1-j}+\mathbf{g}_{k}-\sum_{j=0}^{k-1}\frac{1}{n}\mathbf{1}\mathbf{1}^{T}W^{j}(W-I)\mathbf{g}_{k-1-j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{g}_{k}\\\
&=\sum_{j=0}^{k-1}W^{j}(W-I)\mathbf{g}_{k-1-j}+\mathbf{g}_{k}-\sum_{j=0}^{k-1}\frac{1}{n}\mathbf{1}\mathbf{1}^{T}(W-I)\mathbf{g}_{k-1-j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{g}_{k}\\\
&=\sum_{j=0}^{k-1}(W^{j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})(W-I)\mathbf{g}_{k-1-j}+\mathbf{g}_{k}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{g}_{k}\\\
&=\sum_{j=0}^{k-1}(W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})^{j}(W-I)\mathbf{g}_{k-1-j}+\mathbf{g}_{k}-\mathbf{1}\bar{g}_{k},\\\
\end{split}$

where the last equality can be proven by recursion and the fact that the
matrix $W$ is doubly stochastic by Assumption 1:

$(W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})^{j+1}=(W^{j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})(W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})=W^{j+1}-\frac{1}{n}W^{j}\mathbf{1}\mathbf{1}^{T}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}W+\frac{1}{n}\mathbf{1}\mathbf{1}^{T}=W^{j+1}-\frac{2}{n}\mathbf{1}\mathbf{1}^{T}+\frac{1}{n}\mathbf{1}\mathbf{1}^{T}=W^{j+1}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}$.

Thus,

$\begin{split}\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|&\leq\sum_{j=0}^{k-1}\|(W-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})^{j}(W-I)\mathbf{g}_{k-1-j}\|+\|\mathbf{g}_{k}-\mathbf{1}\bar{g}_{k}\|\\\
&\leq\sum_{j=0}^{k-1}\rho_{w}^{j}\|(W-I)\mathbf{g}_{k-1-j}\|+\|\mathbf{g}_{k}-\mathbf{1}\bar{g}_{k}\|.\\\
\end{split}$

From Lemma 9, we have $\|\mathbf{g}_{k}\|^{2}<\infty$ almost surely.

$\begin{split}\|\mathbf{g}_{k}-\mathbf{1}\bar{g}_{k}\|^{2}&=\sum_{i=1}^{n}\|g_{i,k}-\frac{1}{n}\sum_{j=1}^{n}g_{j,k}\|^{2}\\\
&=\sum_{i=1}^{n}\bigg{(}\|g_{i,k}\|^{2}-2\langle
g_{i,k},\frac{1}{n}\sum_{j=1}^{n}g_{j,k}\rangle+\|\bar{g}_{k}\|^{2}\bigg{)}\\\
&=\|\mathbf{g}_{k}\|^{2}-2n\|\bar{g}_{k}\|^{2}+n\|\bar{g}_{k}\|^{2}\\\
&=\|\mathbf{g}_{k}\|^{2}-n\|\bar{g}_{k}\|^{2}\\\
&\leq\|\mathbf{g}_{k}\|^{2}\\\ &\leq M^{\prime 2}<\infty.\\\ \end{split}$

Inserting in the previous inequality, we get

$\begin{split}\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|&\leq\frac{M^{\prime}}{1-\rho_{w}}\|(W-I)\|+M^{\prime}\\\
&=G<\infty,\end{split}$ (26)

where we have a geometric sum as $\rho_{w}<1$.

  1. 1.

Proving
$\lim_{k\rightarrow\infty}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}=0$

Reconsider (23),

$\begin{split}\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|^{2}&\leq\frac{1+\rho_{w}^{2}}{2}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}+\alpha_{k}^{2}\frac{(1+\rho_{w}^{2})\rho_{w}^{2}}{1-\rho_{w}^{2}}\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|^{2}\\\
\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}&\leq\frac{1+\rho_{w}^{2}}{2}\|\mathbf{x}_{k-1}-\mathbf{1}\bar{x}_{k-1}\|^{2}+\alpha_{k-1}^{2}\frac{(1+\rho_{w}^{2})\rho_{w}^{2}}{1-\rho_{w}^{2}}\|\mathbf{y}_{k-1}-\mathbf{1}\bar{y}_{k-1}\|^{2}\\\
&\ldots\\\
\|\mathbf{x}_{1}-\mathbf{1}\bar{x}_{1}\|^{2}&\leq\frac{1+\rho_{w}^{2}}{2}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}+\alpha_{0}^{2}\frac{(1+\rho_{w}^{2})\rho_{w}^{2}}{1-\rho_{w}^{2}}\|\mathbf{y}_{0}-\mathbf{1}\bar{y}_{0}\|^{2}.\\\
\end{split}$ (27)

Adding all inequalities in (27), we obtain

$\begin{split}\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|^{2}&\leq-\frac{1-\rho_{w}^{2}}{2}\sum_{i=1}^{k}\|\mathbf{x}_{i}-\mathbf{1}\bar{x}_{i}\|^{2}+\frac{1+\rho_{w}^{2}}{2}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}\\\
&\hskip
14.22636pt+\frac{(1+\rho_{w}^{2})\rho_{w}^{2}}{1-\rho_{w}^{2}}\sum_{i=0}^{k}\alpha_{i}^{2}\|\mathbf{y}_{i}-\mathbf{1}\bar{y}_{i}\|^{2}\\\
&\overset{(a)}{\leq}-\frac{1-\rho_{w}^{2}}{2}\sum_{i=1}^{k}\|\mathbf{x}_{i}-\mathbf{1}\bar{x}_{i}\|^{2}+\frac{1+\rho_{w}^{2}}{2}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}\\\
&\hskip
14.22636pt+G^{2}\frac{(1+\rho_{w}^{2})\rho_{w}^{2}}{1-\rho_{w}^{2}}\sum_{i=0}^{k}\alpha_{i}^{2},\end{split}$

with $(a)$ being due to (26). Let $k\rightarrow\infty$, then the second and
third terms are bounded due to Assumption 5. There are then 2 cases:
$\sum\|\mathbf{x}_{i}-\mathbf{1}\bar{x}_{i}\|^{2}$ either diverges or
converges. Assume the validity of the hypothesis H2)
$\sum\|\mathbf{x}_{i}-\mathbf{1}\bar{x}_{i}\|^{2}$ diverges, i.e.,
$\sum_{i=1}^{\infty}\|\mathbf{x}_{i}-\mathbf{1}\bar{x}_{i}\|^{2}\rightarrow\infty$.
This leads to

$\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|^{2}<-\infty,$

as $-\frac{1-\rho_{w}^{2}}{2}<0$. However,
$\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|^{2}$ should be positive. Thus,
hypothesis H2 cannot be true and
$\sum\|\mathbf{x}_{i}-\mathbf{1}\bar{x}_{i}\|^{2}$ converges. Hence,
$\lim_{k\rightarrow\infty}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}=0$
almost surely.

  2. 2.

Proving
$\sum_{k=0}^{\infty}\gamma_{k}\alpha_{k}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|<\infty$

Going back to (25),

$\begin{split}\|\mathbf{x}_{k+1}-\mathbf{1}\bar{x}_{k+1}\|&\leq(\frac{1+\rho_{w}^{2}}{2})^{\frac{k+1}{2}}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|+G\sqrt{\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}}\sum_{j=0}^{k}(\frac{1+\rho_{w}^{2}}{2})^{\frac{j+1}{2}}\alpha_{k-j},\\\
\end{split}$

then substituting into the sum
$\sum_{k=0}^{\infty}\gamma_{k}\alpha_{k}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|$,

$\begin{split}&\sum_{k=1}^{\infty}\gamma_{k}\alpha_{k}((\frac{1+\rho_{w}^{2}}{2})^{\frac{k}{2}}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|+G\sqrt{\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}}\sum_{j=0}^{k-1}(\frac{1+\rho_{w}^{2}}{2})^{\frac{j+1}{2}}\alpha_{k-1-j})\\\
\leq&\gamma_{0}\alpha_{0}\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|\frac{\sqrt{1+\rho_{w}^{2}}}{\sqrt{2}-\sqrt{1+\rho_{w}^{2}}}+G\sqrt{\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}}\sum_{k=1}^{\infty}\gamma_{k}\alpha_{k}\sum_{j=0}^{k-1}(\frac{1+\rho_{w}^{2}}{2})^{\frac{j+1}{2}}\alpha_{k-1-j},\\\
\end{split}$

where the inequality is due to the fact that $\gamma_{k}$ and $\alpha_{k}$ are
both decreasing step-sizes and we have a geometric sum of ratio
$\sqrt{\frac{1+\rho^{2}}{2}}<1$. We then study the sums in the second term,

$\begin{split}\sum_{k=1}^{\infty}\gamma_{k}\alpha_{k}\sum_{j=0}^{k-1}(\frac{1+\rho_{w}^{2}}{2})^{\frac{j+1}{2}}\alpha_{k-1-j}\leq&\sum_{k=1}^{\infty}\gamma_{k}\sum_{j=0}^{k-1}(\frac{1+\rho_{w}^{2}}{2})^{\frac{j+1}{2}}\alpha_{k-1-j}^{2}\\\
=&\sum_{k=1}^{\infty}\gamma_{k}\sum_{j=1}^{k}(\frac{1+\rho_{w}^{2}}{2})^{\frac{k-j+1}{2}}\alpha_{j-1}^{2}\\\
=&\sum_{j=1}^{\infty}\alpha_{j-1}^{2}\sum_{k=j}^{\infty}\gamma_{k}(\frac{1+\rho_{w}^{2}}{2})^{\frac{k-j+1}{2}}\\\
\leq&\gamma_{0}\sum_{j=1}^{\infty}\alpha_{j-1}^{2}\sum_{k=j}^{\infty}(\frac{1+\rho_{w}^{2}}{2})^{\frac{k-j+1}{2}}\\\
=&\gamma_{0}\frac{\sqrt{1+\rho_{w}^{2}}}{\sqrt{2}-\sqrt{1+\rho_{w}^{2}}}\sum_{j=1}^{\infty}\alpha_{j-1}^{2}\\\
<&\infty\end{split}$

as $\sum\alpha_{k}^{2}$ converges by Assumption 5.

Finally,
$\sum_{k=0}^{\infty}\gamma_{k}\alpha_{k}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|<\infty$.

##  Appendix Appendix D Convergence Rate

As a reminder of the definitions of the parameters in Theorem 5, consider
again $R=\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}$,
$\delta_{k}=(\frac{1+\rho_{w}^{2}}{2})^{k}$, and
$\beta_{k}=\sum_{j=1}^{k}\delta_{j}\alpha_{k-j}^{2}$. From Lemma 9 and (18),
we let $\bar{M}$ denote the upper bound of $\mathbb{E}[\|\bar{g}_{k}\|^{2}]$,
and from (26), $G$ that of $\|\mathbf{y}_{k}-\mathbf{1}\bar{y}_{k}\|$. Let
again $\bar{G}=\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}G$.

Our primary result, stated in the following Lemma, is based on finding a
relation between two successive iterations of the expected divergence.

###### Lemma 11

Let $A=\lambda\alpha_{2}$, $B=\alpha_{1}\alpha_{3}^{3}$, and
$C=\alpha_{2}\frac{L^{2}}{\lambda n}$. Whenever Assumption 7 holds, for $k>1$,
we get

$D_{k+1}\leq(1-A\alpha_{k}\gamma_{k})D_{k}+B\alpha_{k}\gamma_{k}^{2}\sqrt{D_{k}}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}.$
(28)

Proof: See Appendix D.2.

Next, we let

$K_{0}=\arg\min_{A\alpha_{k}\gamma_{k}<1}k.$

For the ensuing part, the purpose is to locate a vanishing upper bound of
$D_{k}$, making use of the inequality (28). The idea is to propose a
decreasing sequence $U_{k+1}\leq U_{k}$ and suppose that $D_{k}\leq U_{k}$,
$\forall k\geq K_{0}$, and then verify that $D_{k+1}\leq U_{k+1}$ by
induction. The choice of $U_{k}$ is the most difficult component as one has to
keep in mind the general forms of $\alpha_{k}$ and $\gamma_{k}$ in (28) and
what kind of decisions to take regarding these forms, alongside knowing the
exact rate of $\beta_{k}$. An essential property of $U_{k}$ is presented in
the subsequent lemma.

###### Lemma 12

If a decreasing sequence $U_{k+1}\leq U_{k}$ for $k\geq K_{0}$ exists such
that $D_{k+1}\leq U_{k+1}$ can be deduced from $D_{k}\leq U_{k}$ and (28),
then

$U_{k}\geq\Bigg{(}\frac{B}{2A}\gamma_{k}+\sqrt{\bigg{(}\frac{B}{2A}\bigg{)}^{2}\gamma_{k}^{2}+\frac{C}{A}[\delta_{k}R+\beta_{k}\bar{G}]+\frac{\bar{M}}{A}\frac{\alpha_{k}}{\gamma_{k}}}\Bigg{)}^{2}.$
(29)

An important remark is that the lower bound of $U_{k}$ in (29) is vanishing as
$\gamma_{k},\delta_{k},\beta_{k}$, and $\frac{\alpha_{k}}{\gamma_{k}}$ are all
vanishing. This lower bound provides an insight on the convergence rate of
$D_{k}$ as it cannot be better than that of
$\gamma_{k}^{2},\delta_{k},\beta_{k}$, and $\frac{\alpha_{k}}{\gamma_{k}}$.

Proof: See Appendix D.3.

The previous Lemma allows us to move forward in confirming the existence of
the constants $\varsigma_{1}$ and $\varsigma_{2}$ that permit
$D_{k}\leq\varsigma_{1}^{2}\gamma_{k}^{2}$ and
$D_{k}\leq\varsigma_{2}^{2}\frac{\alpha_{k}}{\gamma_{k}}$ in Theorem 5,
respectively.

###  Appendix D.1 Proof of Theorem 5

  1. 1.

Proof of (8)

By definition of $\varsigma_{1}$,
$D_{K_{0}}\leq\varsigma_{1}^{2}\gamma_{K_{0}}^{2}$. The next step is to make
sure that $D_{k+1}\leq U_{k+1}$ can be obtained from $D_{k}\leq U_{k}$,
$\forall k\geq K_{0}$. Take $U_{k}=\varsigma_{1}^{2}\gamma_{k}^{2}$, let
$D_{k}\leq U_{k}$ hold, and substitute in (28),

$D_{k+1}\leq(1-A\alpha_{k}\gamma_{k})\gamma_{k}^{2}\varsigma_{1}^{2}+B\alpha_{k}\gamma_{k}^{3}\varsigma_{1}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}.$

We solve $D_{k+1}\leq U_{k+1}$ for $\varsigma_{1}\in\mathbb{R}^{+}$

$(1-A\alpha_{k}\gamma_{k})\gamma_{k}^{2}\varsigma_{1}^{2}+B\alpha_{k}\gamma_{k}^{3}\varsigma_{1}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}\leq
U_{k+1}=\varsigma_{1}^{2}\gamma_{k+1}^{2}.$

Then, by considering
$\kappa_{k}=\frac{1-(\frac{\gamma_{k+1}}{\gamma_{k}})^{2}}{\alpha_{k}\gamma_{k}}>0$
as given in (7),

$(\kappa_{k}-A)\varsigma_{1}^{2}+B\varsigma_{1}+CR\delta_{k}\gamma_{k}^{-2}+C\bar{G}\beta_{k}\gamma_{k}^{-2}+\bar{M}\alpha_{k}\gamma_{k}^{-3}\leq
0,$

and assuming $\kappa_{k}-A<0$, we find a constant $\bar{\varsigma_{1}}$ such
that

$\varsigma_{1}\geq\bar{\varsigma_{1}}=\frac{B}{2(A-\kappa_{k})}+\sqrt{\Bigg{(}\frac{B}{2(A-\kappa_{k})}\Bigg{)}^{2}+\frac{CR\delta_{k}\gamma_{k}^{-2}+C\bar{G}\beta_{k}\gamma_{k}^{-2}+\bar{M}\alpha_{k}\gamma_{k}^{-3}}{(A-\kappa_{k})}},$

keeping in mind that $B$ and
$CR\delta_{k}\gamma_{k}^{-2}+C\bar{G}\beta_{k}\gamma_{k}^{-2}+\bar{M}\alpha_{k}\gamma_{k}^{-3}$
are positive by definition. Examine the parameters $\sigma_{1}$, $\sigma_{2}$,
$\sigma_{3}$ and $\sigma_{4}$ as they are introduced in (7), then

$\bar{\varsigma_{1}}\leq\frac{B}{2(A-\sigma_{1})}+\sqrt{\Bigg{(}\frac{B}{2(A-\sigma_{1})}\Bigg{)}^{2}+\frac{CR\sigma_{2}+C\bar{G}\sigma_{3}+\bar{M}\sigma_{4}}{(A-\sigma_{1})}}.$

We conclude that $D_{k}\leq\varsigma_{1}^{2}\gamma_{k}^{2}$ where
$\varsigma_{1}$ satisfies the definition (9).

  2. 2.

Proof of (10)

$D_{K_{0}}\leq\varsigma_{2}^{2}\frac{\gamma_{K_{0}}}{\alpha_{K_{0}}}$ by
definition of $\varsigma_{2}$. $\forall k\geq K_{0}$, let
$D_{k}\leq\varsigma_{2}^{2}\frac{\alpha_{k}}{\gamma_{k}}$, then

$D_{k+1}\leq(1-A\alpha_{k}\gamma_{k})\varsigma_{2}^{2}\frac{\alpha_{k}}{\gamma_{k}}+B\alpha_{k}\gamma_{k}^{2}\varsigma_{2}\sqrt{\frac{\alpha_{k}}{\gamma_{k}}}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}.$

Solving $D_{k+1}\leq\varsigma_{2}^{2}\frac{\alpha_{k+1}}{\gamma_{k+1}}$ for
$\varsigma_{2}\in\mathbb{R}^{+}$,

$(1-A\alpha_{k}\gamma_{k})\varsigma_{2}^{2}\frac{\alpha_{k}}{\gamma_{k}}+B\alpha_{k}\gamma_{k}^{2}\varsigma_{2}\sqrt{\frac{\alpha_{k}}{\gamma_{k}}}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}\leq\varsigma_{2}^{2}\frac{\alpha_{k+1}}{\gamma_{k+1}}.$

Take
$\tau_{k}=\frac{\frac{\alpha_{k}}{\gamma_{k}}-\frac{\alpha_{k+1}}{\gamma_{k+1}}}{\alpha_{k}^{2}}>0$
as given in (7), then

$(\tau_{k}-A)\varsigma_{2}^{2}+B\gamma_{k}^{\frac{3}{2}}\alpha_{k}^{-\frac{1}{2}}\varsigma_{2}+C\gamma_{k}\alpha_{k}^{-1}[\delta_{k}R+\beta_{k}\bar{G}]+\bar{M}\leq
0.$

If
$\frac{\alpha_{k}}{\gamma_{k}}-\frac{\alpha_{k+1}}{\gamma_{k+1}}<A\alpha_{k}^{2}$,
then $\exists$ $\bar{\varsigma_{2}}$ such that

$\varsigma_{2}\geq\bar{\varsigma_{2}}=\frac{B\gamma_{k}^{\frac{3}{2}}\alpha_{k}^{-\frac{1}{2}}}{2(A-\tau_{k})}+\sqrt{\bigg{(}\frac{B\gamma_{k}^{\frac{3}{2}}\alpha_{k}^{-\frac{1}{2}}}{2(A-\tau_{k})}\bigg{)}^{2}+\frac{C\gamma_{k}\alpha_{k}^{-1}[\delta_{k}R+\beta_{k}\bar{G}]+\bar{M}}{(A-\tau_{k})}}.$

Examine $\sigma_{5},\sigma_{6},\sigma_{7}$ and $\sigma_{8}$ that are defined
in (7), we can say

$\bar{\varsigma_{2}}\leq\frac{B\sigma_{6}}{2(A-\sigma_{5})}+\sqrt{\bigg{(}\frac{B\sigma_{6}}{2(A-\sigma_{5})}\bigg{)}^{2}+\frac{C(R\sigma_{7}+\bar{G}\sigma_{8})+\bar{M}}{(A-\sigma_{5})}}.$

We conclude that $D_{k}\leq\varsigma_{2}^{2}\frac{\alpha_{k}}{\gamma_{k}}$
with $\varsigma_{2}$ satisfying (11).

###  Appendix D.2 Proof of Lemma 11

We start by expressing the expected divergence in terms of its previous
iteration.

$\begin{split}D_{k+1}=&\mathbb{E}[\|\bar{x}_{k+1}-x^{*}\|^{2}]\\\
=&\mathbb{E}[\|\bar{x}_{k}-\alpha_{k}\bar{g}_{k}-x^{*}\|^{2}]\\\
=&D_{k}+\alpha_{k}^{2}\mathbb{E}[\|\bar{g}_{k}\|^{2}]-2\alpha_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\bar{g}_{k}\rangle]\\\
\overset{(a)}{=}&D_{k}+\alpha_{k}^{2}\mathbb{E}[\|\bar{g}_{k}\|^{2}]-2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},h(\mathbf{x}_{k})+\bar{b}_{k}\rangle]\\\
=&D_{k}+\alpha_{k}^{2}\mathbb{E}[\|\bar{g}_{k}\|^{2}]-2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathcal{F}(\bar{x}_{k})\rangle]+2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathcal{F}(\bar{x}_{k})-h(\mathbf{x}_{k})\rangle]\\\
&-2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\bar{b}_{k}\rangle]\\\
\end{split}$ (30)

where $(a)$ is due to both $\mathbb{E}[e_{k}|\mathcal{H}_{k}]=0$ and (14):

$\begin{split}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\bar{g}_{k}\rangle]&=\mathbb{E}[\langle\bar{x}_{k}-x^{*},\bar{g}_{k}-\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]+\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle]\\\
&=\mathbb{E}[\langle\bar{x}_{k}-x^{*},e_{k}\rangle]+\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle]\\\
&=\mathbb{E}_{\mathcal{H}_{k}}[\mathbb{E}[\langle\bar{x}_{k}-x^{*},e_{k}\rangle|\mathcal{H}_{k}]]+\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle]\\\
&=0+\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathbb{E}[\bar{g}_{k}|\mathcal{H}_{k}]\rangle].\end{split}$

From Lemma 9 and (18), we have $\mathbb{E}[\|\bar{g}_{k}\|^{2}]<\bar{M}$
almost surely, with $\bar{M}=\frac{1}{n}M$ a bounded constant.

By the strong convexity in Assumption 7, we have

$\begin{split}-2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathcal{F}(\bar{x}_{k})\rangle]&\leq-2\lambda\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\|\bar{x}_{k}-x^{*}\|^{2}]\\\
&=-2\lambda\alpha_{2}\alpha_{k}\gamma_{k}D_{k}.\end{split}$ (31)

Next, from Lemma 2, we have

$\begin{split}2\alpha_{2}\alpha_{k}\gamma_{k}\langle\bar{x}_{k}-x^{*},\mathcal{F}(\bar{x}_{k})-h(\mathbf{x}_{k})\rangle&\leq
2\alpha_{2}\alpha_{k}\gamma_{k}\frac{L}{\sqrt{n}}\|\bar{x}_{k}-x^{*}\|\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|\\\
&\overset{(a)}{\leq}\lambda\alpha_{2}\alpha_{k}\gamma_{k}\|\bar{x}_{k}-x^{*}\|^{2}+\alpha_{2}\alpha_{k}\gamma_{k}\frac{L^{2}}{\lambda
n}\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2},\end{split}$

where $(a)$ is due to $2\sqrt{\epsilon}\times\frac{1}{\sqrt{\epsilon}}\langle
a,b\rangle=2\langle\sqrt{\epsilon}a,\frac{1}{\sqrt{\epsilon}}b\rangle\leq\epsilon\|a\|^{2}+\frac{1}{\epsilon}\|b\|^{2}$.
From (24) and (26), we get

$\|\mathbf{x}_{k}-\mathbf{1}\bar{x}_{k}\|^{2}\leq\delta_{k}R+\beta_{k}\bar{G}$

where $R=\|\mathbf{x}_{0}-\mathbf{1}\bar{x}_{0}\|^{2}$,
$\delta_{k}=(\frac{1+\rho_{w}^{2}}{2})^{k}$,
$\beta_{k}=\sum_{j=1}^{k}\delta_{j}\alpha_{k-j}^{2}$, and
$\bar{G}=\frac{2\rho_{w}^{2}}{1-\rho_{w}^{2}}G$. Hence,

$2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\mathcal{F}(\bar{x}_{k})-h(\mathbf{x}_{k})\rangle]\leq\lambda\alpha_{2}\alpha_{k}\gamma_{k}D_{k}+\alpha_{2}\alpha_{k}\gamma_{k}\frac{L^{2}}{\lambda
n}[\delta_{k}R+\beta_{k}\bar{G}].$ (32)

From (15),

$\begin{split}-2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\langle\bar{x}_{k}-x^{*},\bar{b}_{k}\rangle]&\leq
2\alpha_{2}\alpha_{k}\gamma_{k}\mathbb{E}[\|\bar{x}_{k}-x^{*}\|\|\bar{b}_{k}\|]\\\
&\leq\alpha_{1}\alpha_{3}^{3}\alpha_{k}\gamma_{k}^{2}\sqrt{D_{k}}.\end{split}$
(33)

Finally, by combining (30), (31), (32), and (33), we get (28).

###  Appendix D.3 Proof of Lemma 12

Since $1-A\alpha_{k}\gamma_{k}>0$ when $k\geq K_{0}$, we may substitute
$D_{k}\leq U_{k}$ in (28),

$D_{k+1}\leq(1-A\alpha_{k}\gamma_{k})U_{k}+B\alpha_{k}\gamma_{k}^{2}\sqrt{U_{k}}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}.$

Testing $D_{k+1}\leq U_{k+1}$ in the previous inequality, we get

$(1-A\alpha_{k}\gamma_{k})U_{k}+B\alpha_{k}\gamma_{k}^{2}\sqrt{U_{k}}+C\alpha_{k}\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]+\alpha_{k}^{2}\bar{M}\leq
U_{k+1}\leq U_{k}$
$A\gamma_{k}U_{k}-B\gamma_{k}^{2}\sqrt{U_{k}}-C\gamma_{k}[\delta_{k}R+\beta_{k}\bar{G}]-\alpha_{k}\bar{M}\geq
0.$ (34)

Since $\sqrt{U_{k}}\geq 0$, the solution of (34) is given by (29), which
concludes the proof.

###  Appendix D.4 Proof of Lemma 6

We consider the following series, as it affects the maximum possible
convergence rate which appears in (29) and all inequalities related to
$D_{k}$,

$\beta_{k}=\sum_{j=1}^{k}\delta_{j}\alpha_{k-j}^{2}=\sum_{j=1}^{k}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j}\alpha_{k-j}^{2}.$

We then try to find $\beta_{k+1}$ in terms of $\beta_{k}$,

$\begin{split}\beta_{k+1}=&\sum_{j=1}^{k+1}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j}\alpha_{k+1-j}^{2}\\\
=&\sum_{j=2}^{k+1}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j}\alpha_{k+1-j}^{2}+\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}\alpha_{k}^{2}\\\
=&\sum_{j=1}^{k}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j+1}\alpha_{k-j}^{2}+\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}\alpha_{k}^{2}\\\
=&\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}(\beta_{k}+\alpha_{k}^{2}).\end{split}$

We know that this series has a Q-linear to a Q-sublinear convergence rate as

$\begin{split}\frac{\beta_{k+1}}{\beta_{k}}&=\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}\Bigg{(}1+\frac{\alpha_{k}^{2}}{\sum_{j=1}^{k}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{j}\alpha_{k-j}^{2}}\Bigg{)}\\\
&=\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}+\frac{\alpha_{k}^{2}}{\sum_{j=1}^{k}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{j-1}\alpha_{k-j}^{2}}\\\
&=\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}+\frac{\alpha_{k}^{2}}{\sum_{j=0}^{k-1}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{j}\alpha_{k-1-j}^{2}}\\\
&=\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}+\frac{1}{\sum_{j=0}^{k-1}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{j}\frac{\alpha_{k-1-j}^{2}}{\alpha_{k}^{2}}}\\\
&\overset{(a)}{\leq}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}+\frac{1}{\sum_{j=0}^{k-1}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{j}}\\\
&\overset{(b)}{=}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}+\frac{1-\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}}{1-\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{k}}\\\
\end{split}$

where $(a)$ is since $\frac{\alpha_{k-1-j}^{2}}{\alpha_{k}^{2}}\geq 1$ for
every $j\in\\{0,\ldots,k-1\\}$ and $(b)$ is due to the geometric sum
$\sum_{j=0}^{k-1}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{j}=\frac{1-\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{k}}{1-\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}}$.
Hence,

$\lim_{k\rightarrow\infty}\frac{\beta_{k+1}}{\beta_{k}}\leq\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}+1-\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}=1$

since $\lim_{k\rightarrow\infty}\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{k}=0$.
Next, to get an idea of how $\beta_{k}$ converges in terms of $k$, we consider
$K_{1}$ that is defined in (13). We know that $K_{1}$ exists and is finite as
$\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{k}$ decreases much faster than
$\alpha_{k}^{2}$. Taking $\alpha_{k}$ as that in (12), we know that
$0.5<\upsilon_{1}<1$ and we find the condition on $\upsilon_{1}$ such that

$\begin{split}\Bigg{(}\frac{1+\rho_{w}^{2}}{2}\Bigg{)}^{k}\leq&\alpha_{k}^{2}\\\
k\log\Bigg{(}\frac{1+\rho_{w}^{2}}{2}\Bigg{)}\leq&-2\upsilon_{1}\log(k+1)\\\
\upsilon_{1}\leq&\frac{1}{2}\log\bigg{(}\frac{2}{1+\rho_{w}^{2}}\bigg{)}\frac{k}{\log(k+1)},\\\
\end{split}$

which is feasible as $\frac{1}{2}\log\big{(}\frac{2}{1+\rho_{w}^{2}}\big{)}>0$
and $\frac{k}{\log(k+1)}$ grows very large for $k\geq K_{1}$, a simple
condition is that
$1<\frac{1}{2}\log(\frac{2}{1+\rho_{w}^{2}})\frac{k}{\log(k+1)}$ which gives

$K_{1}=\arg\min_{2\log^{-1}(\frac{2}{1+\rho_{w}^{2}})<k\log^{-1}(k+1)}k.$

Thus, we can write $\beta_{k}$ as

$\begin{split}\beta_{k}&=\sum_{j=1}^{K_{1}-1}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j}\alpha_{k-j}^{2}+\sum_{j=K_{1}}^{k}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j}\alpha_{k-j}^{2}\\\
&\overset{(a)}{\leq}\alpha_{k-K_{1}}^{2}\sum_{j=1}^{K_{1}-1}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{j}+\sum_{j=K_{1}}^{k}\alpha_{j}^{2}\alpha_{k-j}^{2}\\\
&\overset{(b)}{=}\alpha_{k-K_{1}}^{2}\bigg{(}\frac{1+\rho_{w}^{2}}{1-\rho_{w}^{2}}\bigg{)}\bigg{(}1-\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{K_{1}-1}\bigg{)}+\sum_{j=K_{1}}^{k}\alpha_{j}^{2}\alpha_{k-j}^{2}\\\
&=\beta^{\prime}\alpha_{k-K_{1}}^{2}+\sum_{j=K_{1}}^{k}\alpha_{j}^{2}\alpha_{k-j}^{2},\\\
\end{split}$ (35)

where in $(a)$ we used the definition of $K_{1}$ and the fact that
$\alpha_{k}^{2}$ is a decreasing step-size, and in $(b)$ the sum of a
geometric series. We know that
$\beta^{\prime}=\big{(}\frac{1+\rho_{w}^{2}}{1-\rho_{w}^{2}}\big{)}\big{(}1-\big{(}\frac{1+\rho_{w}^{2}}{2}\big{)}^{K_{1}-1}\big{)}$
is finite, then the first sum decreases as
$\alpha_{k-K_{1}}^{2}=\frac{\alpha_{0}^{2}}{(k-K_{1}+1)^{2\upsilon_{1}}}\sim\frac{1}{k^{2\upsilon_{1}}}$
with $1<2\upsilon_{1}<2$, which leaves us with the second sum,

$\begin{split}\sum_{j=K_{1}}^{k}\alpha_{j}^{2}\alpha_{k-j}^{2}&=\alpha_{0}^{2}\sum_{j=K_{1}}^{k}\frac{1}{(j+1)^{2\upsilon_{1}}}\frac{1}{(k-j+1)^{2\upsilon_{1}}}\\\
&=\alpha_{0}^{2}\sum_{j=K_{1}+1}^{k+1}\frac{1}{j^{2\upsilon_{1}}}\frac{1}{(k-j+2)^{2\upsilon_{1}}}\\\
&=\alpha_{0}^{2}\sum_{j=K_{1}+1}^{k+1}\frac{1}{j^{2\upsilon_{1}}}\frac{1}{(k+2)^{2\upsilon_{1}}(1-\frac{j}{k+2})^{2\upsilon_{1}}}\\\
&\overset{(a)}{=}\alpha_{0}^{2}\frac{1}{(k+2)^{4\upsilon_{1}}}\sum_{u\in\mathcal{U}_{k}}\Big{(}\frac{1}{u(1-u)}\Big{)}^{2\upsilon_{1}}\\\
&=\alpha_{0}^{2}\frac{|\mathcal{U}_{k}|}{(k+2)^{4\upsilon_{1}}}\sum_{u\in\mathcal{U}_{k}}\Big{(}\frac{1}{u^{2}(1-u)^{2}}\Big{)}^{\upsilon_{1}}\times\frac{1}{|\mathcal{U}_{k}|}\\\
&\overset{(b)}{\leq}\alpha_{0}^{2}\frac{|\mathcal{U}_{k}|}{(k+2)^{4\upsilon_{1}}}\bigg{(}\sum_{u\in\mathcal{U}_{k}}\frac{1}{u^{2}(1-u)^{2}}\times\frac{1}{|\mathcal{U}_{k}|}\bigg{)}^{\upsilon_{1}}\\\
&\overset{(c)}{=}\alpha_{0}^{2}\frac{|\mathcal{U}_{k}|^{1-\upsilon_{1}}}{(k+2)^{4\upsilon_{1}}}\bigg{(}\sum_{u\in\mathcal{U}_{k}}\frac{1}{u^{2}}+\frac{1}{(1-u)^{2}}+\frac{2}{u}+\frac{2}{(1-u)}\bigg{)}^{\upsilon_{1}}\\\
&=\alpha_{0}^{2}\frac{(k-K_{1}+1)^{1-\upsilon_{1}}}{(k+2)^{4\upsilon_{1}}}\bigg{(}\sum_{u\in\mathcal{U}_{k}}\frac{1}{u^{2}}+\frac{2}{u}+\sum_{u\in\mathcal{U}_{k}\setminus\\{\frac{k+1}{k+2}\\}}\frac{1}{(1-u)^{2}}+\frac{2}{(1-u)}\\\
&\hskip 253.22934pt+2(k+2)+(k+2)^{2}\bigg{)}^{\upsilon_{1}}\\\
&\overset{(d)}{\leq}\alpha_{0}^{2}\frac{(k-K_{1}+1)^{1-\upsilon_{1}}}{(k+2)^{4\upsilon_{1}}}\bigg{(}(k+2)\int_{\frac{K_{1}}{k+2}}^{\frac{k+1}{k+2}}\Big{(}\frac{1}{u^{2}}+\frac{2}{u}\Big{)}du\\\
&\hskip
59.75095pt+(k+2)\int_{\frac{K_{1}+1}{k+2}}^{\frac{k+1}{k+2}}\Big{(}\frac{1}{(1-u)^{2}}+\frac{2}{(1-u)}\Big{)}du+2(k+2)+(k+2)^{2}\bigg{)}^{\upsilon_{1}}\\\
&=\alpha_{0}^{2}\frac{(k-K_{1}+1)^{1-\upsilon_{1}}}{(k+2)^{3\upsilon_{1}}}\bigg{(}\frac{(k+2)(k-K_{1}+1)}{K_{1}(k+1)}\\\
&\hskip
85.35826pt+2\ln\frac{e(k+1)(k-K_{1}+1)}{K_{1}}+\frac{(k+2)(k-K_{1})}{k-K_{1}+1}+(k+2)\bigg{)}^{\upsilon_{1}}\\\
&\leq\alpha_{0}^{2}\frac{(k-K_{1}+1)^{1-\upsilon_{1}}}{(k+2)^{3\upsilon_{1}}}\bigg{(}\frac{(k+2)(k+1)}{K_{1}(k+1)}\\\
&\hskip
85.35826pt+2\ln\frac{e(k+1)(k-K_{1}+1)}{K_{1}}+\frac{(k+2)(k-K_{1})}{k-K_{1}}+(k+2)\bigg{)}^{\upsilon_{1}}\\\
&=\alpha_{0}^{2}\frac{(k-K_{1}+1)^{1-\upsilon_{1}}}{(k+2)^{3\upsilon_{1}}}\bigg{(}2\ln\frac{e(k+1)(k-K_{1}+1)}{K_{1}}+(k+2)\Big{(}2+\frac{1}{K_{1}}\Big{)}\bigg{)}^{\upsilon_{1}}\\\
&\leq\alpha_{0}^{2}\frac{(k-K_{1}+1)^{1-\upsilon_{1}}}{(k+2)^{3\upsilon_{1}}}\bigg{(}4\ln\frac{e(k+2)}{K_{1}}+(k+2)\Big{(}2+\frac{1}{K_{1}}\Big{)}\bigg{)}^{\upsilon_{1}}\\\
\end{split}$ (36)
$\begin{split}&\leq\alpha_{0}^{2}\frac{(k+2)^{1-\upsilon_{1}}}{(k+2)^{3\upsilon_{1}}}(k+2)^{\upsilon_{1}}\bigg{(}\frac{4}{(k+2)}\ln\frac{e}{K_{1}}+6+\frac{1}{K_{1}}\bigg{)}^{\upsilon_{1}}\\\
&\leq\beta^{\prime\prime}\frac{1}{(k+2)^{3\upsilon_{1}-1}},\\\ \end{split}$

where in $(a)$ we changed the summation variable to $u=\frac{j}{k+2}$ and
$\mathcal{U}_{k}=\\{\frac{K_{1}+1}{k+2},\frac{K_{1}+2}{k+2},\ldots,\frac{k+1}{k+2}\\}$,
in $(b)$ we used Jensen’s inequality
$\mathbb{E}[\varphi(x)]\leq\varphi(\mathbb{E}[x])$ for the concave function
$\varphi(x)=x^{\upsilon_{1}}$ with $0.5<\upsilon_{1}<1$, and $(c)$ is by
partial fraction decomposition. In $(d)$, we interpret the sum over
$\mathcal{U}_{k}$ as a Riemann sum in which the function $\frac{1}{u}$ is
evaluated at the right endpoint of the interval
$[\frac{K_{1}+i}{k+2},\frac{K_{1}+i+1}{k+2}]$, for $i=0,1,\ldots,k-K_{1}$.
Since the function $\frac{1}{u}$ is monotonically decreasing, it is in fact a
_lower_ Riemann sum and therefore bounded from above by the integral
$\int_{\frac{K_{1}}{k+2}}^{\frac{k+1}{k+2}}\frac{1}{u}du$. Analogously, we
estimate the sum of $\frac{1}{u^{2}}$ over $\mathcal{U}_{k}$. Mutatis
mutandis, the estimate for the monotonically increasing functions
$\frac{1}{1-u}$ and $\frac{1}{(1-u)^{2}}$ follows.

We have also let $\beta^{\prime\prime}=\max_{k\geq
K_{1}}\bigg{(}\frac{4}{(k+2)}\ln\frac{e}{K_{1}}+6+\frac{1}{K_{1}}\bigg{)}^{\upsilon_{1}}=\bigg{(}\frac{4}{(K_{1}+2)}\ln\frac{e}{K_{1}}+6+\frac{1}{K_{1}}\bigg{)}^{\upsilon_{1}}$.
Hence,
$\beta^{\prime\prime}\frac{1}{(k+2)^{3\upsilon_{1}-1}}\sim\frac{1}{k^{3\upsilon_{1}-1}}$
with $0.5<3\upsilon_{1}-1<2$.

Since $3\upsilon_{1}-1<2\upsilon_{1}$ for $0.5<\upsilon_{1}<1$, then
$\frac{1}{k^{2\upsilon_{1}}}<\frac{1}{k^{3\upsilon_{1}-1}}$ for all $k\geq
K_{1}$.

From (35) and (36),

$\begin{split}\beta_{k}&\leq\beta^{\prime}\frac{1}{(k-K_{1}+1)^{2\upsilon_{1}}}+\beta^{\prime\prime}\frac{1}{(k+2)^{3\upsilon_{1}-1}}\\\
&\leq\beta^{\prime}\frac{1}{(k-K_{1}+1)^{3\upsilon_{1}-1}}+\beta^{\prime\prime}\frac{1}{(k-K_{1}+1)^{3\upsilon_{1}-1}}\\\
&=(\beta^{\prime}+\beta^{\prime\prime})\frac{1}{(k-K_{1}+1)^{3\upsilon_{1}-1}},\end{split}$

and we deduce that $\beta_{k}$ has a convergence rate of at least
$\frac{1}{k^{3\upsilon_{1}-1}}$, which concludes the proof.

###  Appendix D.5 Proof of Theorem 7

Theorem 5 indicates that the convergence rate is a function of $\upsilon_{1}$
and $\upsilon_{2}$, as $\gamma_{k}^{2}\propto(k+1)^{-2\upsilon_{2}}$ and
$\frac{\alpha_{k}}{\gamma_{k}}\propto(k+1)^{-(-\upsilon_{1}-\upsilon_{2})}$.
Nonetheless, we must still verify the validity of the assumptions presented in
the theorem, meaning:

  * •

Are $\sigma_{1}<A$ and $\sigma_{5}<A$ fulfilled?

  * •

Are $\varsigma_{1}$ and $\varsigma_{2}$ bounded?

We must remark that in what follows, the analysis is done for $k\geq K_{2}$.

Let $\alpha_{k}$ and $\gamma_{k}$ have the forms given in (12).

  1. 1.

Verifying that $\sigma_{1}<A$ and $\sigma_{5}<A$

The idea is to find a bound on $\alpha_{0}$ and $\gamma_{0}$ to guarantee
$\sigma_{1}<A$ and $\sigma_{5}<A$. We start by bounding $\sigma_{1}$ and
$\sigma_{5}$ from above, i.e.,

$\begin{split}\sigma_{1}=\underset{k\geq
K_{2}}{\max}\;\frac{1-(\frac{\gamma_{k+1}}{\gamma_{k}})^{2}}{\alpha_{k}\gamma_{k}}=\underset{k\geq
K_{2}}{\max}\;\frac{1-(1+\frac{1}{k+1})^{-2\upsilon_{2}}}{\alpha_{0}\gamma_{0}(k+1)^{-\upsilon_{1}-\upsilon_{2}}}\end{split}$

and

$\sigma_{5}=\underset{k\geq
K_{2}}{\max}\;\frac{1-\frac{\alpha_{k+1}\gamma_{k+1}^{-1}}{\alpha_{k}\gamma_{k}^{-1}}}{\alpha_{k}\gamma_{k}}=\underset{k\geq
K_{2}}{\max}\;\frac{1-(1+\frac{1}{k+1})^{-(\upsilon_{1}-\upsilon_{2})}}{\alpha_{0}\gamma_{0}(k+1)^{-\upsilon_{1}-\upsilon_{2}}}.$

To do so, we define a function $q(x)=x^{-a}(1-(1-x)^{-b})$ with
$a,b,x\in(0,1]$. Since $x^{-a}\leq x^{-1}$, we have $q(x)\leq
x^{-1}(1-(1-x)^{-b})=r(x)$. To further bound $q(x)$, We study the derivative
of $r(x)$ as it’s simpler to do so,

$r^{\prime}(x)=x^{-2}\bigg{(}((b+1)x+1)(1+x)^{-b-1}-1\bigg{)}=x^{-2}s(x).$

Hence the sign of $r^{\prime}(x)$ is that of $s(x)$. We again calculate the
derivative of $s(x)$ to find its sign,

$s^{\prime}(x)=-b(b+1)x(1+x)^{-b-2}\leq 0$

since $b>0$ and $x>0$. Then, $s(x)$ is a decreasing function of $x$ over
$(0,1]$. We remark that $\lim_{x\rightarrow 0}s(x)=0$, meaning $s(x)<0$ and
$r^{\prime}(x)<0$, $\forall x\in(0,1]$. Finally,

$r(x)<\lim_{x\rightarrow 0}r(x)=\frac{1-(1+x)^{-b}}{x}=b,$

and $q(x)\leq r(x)<b$, noting that $\lim_{x\rightarrow 0}q(x)=b$ for $a=1$. We
conclude that $\sigma_{1}<\frac{2\upsilon_{2}}{\alpha_{0}\gamma_{0}}$ and
$\sigma_{1}<\frac{\upsilon_{1}-\upsilon_{2}}{\alpha_{0}\gamma_{0}}$. For
$\sigma_{1}<A$ and $\sigma_{5}<A$ to be valid, we must have

$\alpha_{0}\gamma_{0}\geq\max\\{2\upsilon_{2},\upsilon_{1}-\upsilon_{2}\\}/A.$
(37)

  2. 2.

Verifying that $\varsigma_{1}$ and $\varsigma_{2}$ are bounded

The goal is to verify that the constant term in the convergence rate is
bounded. Thus, we must check that the lower bounds given in (9) and (11) are
indeed finite. We begin by analyzing $\sigma_{2}$ and $\sigma_{7}$, i.e.,

$\sigma_{2}=\underset{k\geq
K_{2}}{\max}\;\frac{\delta_{k}}{\gamma_{k}^{2}}=\underset{k\geq
K_{2}}{\max}\;\gamma_{0}^{-2}(k+1)^{2\upsilon_{2}}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{k}$

and

$\sigma_{7}=\underset{k\geq
K_{2}}{\max}\;\frac{\gamma_{k}}{\alpha_{k}}\delta_{k}=\underset{k\geq
K_{2}}{\max}\;\gamma_{0}\alpha_{0}^{-1}(k+1)^{\upsilon_{1}-\upsilon_{2}}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{k}.$

To prove that $\sigma_{2}$ and $\sigma_{7}$ are finite, we define the function
$p(k)=a(k+1)^{b}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{k}$, with $a>0$ and
$0<b<1$. To find the maximum, we find $k$ such that $p^{\prime}(k)=0$, meaning

$a(k+1)^{b}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{k}\bigg{(}\frac{b}{k+1}+\ln\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}\bigg{)}=0.$

Hence,

$k=b\ln^{-1}\bigg{(}\frac{2}{1+\rho_{w}^{2}}\bigg{)}-1.$

     * •

From here, we can say if
$b\ln^{-1}\bigg{(}\frac{2}{1+\rho_{w}^{2}}\bigg{)}-1\geq K_{2}$, then

$\underset{k\geq
K_{2}}{\max}\;p(k)=\frac{ab^{b}}{\big{(}\frac{2}{1+\rho_{w}^{2}}\big{)}\ln^{b}\big{(}\frac{2}{1+\rho_{w}^{2}}\big{)}}e^{-b}<\infty.$

     * •

Otherwise, if $b\ln^{-1}\bigg{(}\frac{2}{1+\rho_{w}^{2}}\bigg{)}-1<K_{2}\leq
k$, then $\frac{b}{k+1}+\ln\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}<0$ which
gives $p^{\prime}<0$ for $k\geq K_{2}$, meaning $p$ is strictly decreasing and

$\underset{k\geq
K_{2}}{\max}\;p(k)=p(K_{2})=a(K_{2}+1)^{b}\bigg{(}\frac{1+\rho_{w}^{2}}{2}\bigg{)}^{K_{2}}<\infty.$

We conclude that $\sigma_{2}<\infty$ and $\sigma_{7}<\infty$.

Next, we study the finiteness of $\sigma_{3}$ and $\sigma_{8}$. From Lemma 6,
we can write

$\begin{split}\sigma_{3}=\underset{k\geq
K_{2}}{\max}\;\frac{\beta_{k}}{\gamma_{k}^{2}}&\leq\underset{k\geq
K_{2}}{\max}\;(\beta^{\prime}+\beta^{\prime\prime})\frac{\alpha_{k-K_{1}}^{2}}{\gamma_{k}^{2}}\\\
&=\underset{k\geq
K_{2}}{\max}\;(\beta^{\prime}+\beta^{\prime\prime})\alpha_{0}^{2}\gamma_{0}^{-2}\frac{(k+1)^{2\upsilon_{2}}}{(k-K_{1}+1)^{3\upsilon_{1}-1}}\end{split}$

and

$\begin{split}\sigma_{8}=\underset{k\geq
K_{2}}{\max}\;\frac{\gamma_{k}}{\alpha_{k}}\beta_{k}&\leq\underset{k\geq
K_{2}}{\max}\;(\beta^{\prime}+\beta^{\prime\prime})\frac{\gamma_{k}}{\alpha_{k}}=\frac{\gamma_{k}\alpha_{k-K_{1}}^{2}}{\alpha_{k}}\\\
&=\underset{k\geq
K_{2}}{\max}\;(\beta^{\prime}+\beta^{\prime\prime})\gamma_{0}\alpha_{0}\frac{(k+1)^{\upsilon_{1}-\upsilon_{2}}}{(k-K_{1}+1)^{3\upsilon_{1}-1}}.\end{split}$

We then define a function
$q(k)=a\frac{(k+1)^{b}}{(k-K_{1}+1)^{3\upsilon_{1}-1}}$ for $k\geq K_{2}$,
with $a>0$ and $0<b<1$, and we study its derivative

$q^{\prime}(k)=a\frac{(b-3\upsilon_{1}+1)k-bK_{1}+b-3\upsilon_{1}+1}{(k+1)^{1-b}(k-K_{1}+1)^{3\upsilon_{1}}}.$

We know that $q^{\prime}<0$, and thus $q$ is strictly decreasing for $k\geq
K_{2}$, when $b-3\upsilon_{1}+1\leq 0$.

Hence,

$\begin{split}\sigma_{3}=\Bigg{\\{}\begin{matrix}\alpha_{0}^{2}\gamma_{0}^{-2}\frac{(K_{2}+1)^{2\upsilon_{2}}}{(K_{2}-K_{1}+1)^{3\upsilon_{1}-1}},&\text{if}\;2\upsilon_{2}-3\upsilon_{1}+1\leq
0,\\\
\infty,&\text{if}\;2\upsilon_{2}-3\upsilon_{1}+1>0,\end{matrix}\end{split}$
(38)

and

$\begin{split}\sigma_{8}=(\beta^{\prime}+\beta^{\prime\prime})\gamma_{0}\alpha_{0}\frac{(K_{2}+1)^{\upsilon_{1}-\upsilon_{2}}}{(K_{2}-K_{1}+1)^{3\upsilon_{1}-1}}<\infty\end{split}$

since $b-3\upsilon_{1}+1=-\upsilon_{2}-2\upsilon_{1}+1<0$ always holds.

We end with the analysis of $\sigma_{4}$ and $\sigma_{6}$, i.e.,

$\sigma_{4}=\alpha_{0}\gamma_{0}^{-3}\underset{k\geq
K_{2}}{\max}\;(1+k)^{-(\upsilon_{1}-3\upsilon_{2})}=\Bigg{\\{}\begin{matrix}\alpha_{0}\gamma_{0}^{-3}(1+K_{2})^{-(\upsilon_{1}-3\upsilon_{2})},&\text{if}\;\upsilon_{1}\geq
3\upsilon_{2},\\\ \infty,&\text{if}\;\upsilon_{1}<3\upsilon_{2},\end{matrix}$

and

$\sigma_{6}=\alpha_{0}^{-\frac{1}{2}}\gamma_{0}^{\frac{3}{2}}\underset{k\geq
K_{2}}{\max}\;(1+k)^{\frac{\upsilon_{1}-3\upsilon_{2}}{2}}=\Bigg{\\{}\begin{matrix}\alpha_{0}^{-\frac{1}{2}}\gamma_{0}^{\frac{3}{2}}(1+K_{2})^{\frac{\upsilon_{1}-3\upsilon_{2}}{2}},&\text{if}\;\upsilon_{1}\leq
3\upsilon_{2},\\\ \infty,&\text{if}\;\upsilon_{1}>3\upsilon_{2}.\end{matrix}$

There are clearly 3 cases:

     * •

$\upsilon_{1}>3\upsilon_{2}$

Thus, $\sigma_{4}$ is bounded.

Since now
$2\upsilon_{2}-3\upsilon_{1}+1<\frac{2}{3}\upsilon_{1}-3\upsilon_{1}+1=-\frac{7}{3}\upsilon_{1}+1<0$
always holds, then $\sigma_{3}$ (38) and $\varsigma_{1}$ (by definition) are
also bounded provided that $\alpha_{0}\gamma_{0}\geq\frac{2\upsilon_{2}}{A}$
in (37).

However, $\varsigma_{2}\rightarrow\infty$ since $\sigma_{6}\rightarrow\infty$
resulting in a loose upper bound in (10).

To that end, we can write $D_{k}\leq\Upsilon_{1}(1+k)^{-2\upsilon_{2}}$ with
$\Upsilon_{1}$ a bounded constant.

     * •

$\upsilon_{1}<3\upsilon_{2}$

Similarly, $\sigma_{6}$ is bounded while $\sigma_{4}\rightarrow\infty$. Then,
$\exists$ $\Upsilon_{2}<\infty$, where
$D_{k}\leq\Upsilon_{2}(1+k)^{-(\upsilon_{1}-\upsilon_{2})}$ provided that
$\alpha_{0}\gamma_{0}\geq\frac{\upsilon_{1}-\upsilon_{2}}{A}$.

     * •

$\upsilon_{1}=3\upsilon_{2}$

Both $\sigma_{4}$ and $\sigma_{6}$ are bounded allowing both previous
inequalities corresponding to $D_{k}$ to be valid.

By this analysis, we conclude the proof of Theorem 7.

##  Appendix Appendix E Regret Analysis

Since, by Lemma 2, the objective function is $L$-smooth, we can write

$\mathcal{F}(y)\leq\mathcal{F}(x)+\langle\nabla\mathcal{F}(x),y-x\rangle+\frac{L}{2}\|y-x\|^{2},\;\forall
x,y\in\mathbb{R}^{d}.$ (39)

To find the regret bound, consider

$\begin{split}\mathbb{E}_{\mathcal{H}_{k}}\bigg{[}\sum_{k=1}^{K}\mathcal{F}(\bar{x}_{k})-\mathcal{F}(x^{*})\bigg{]}&\overset{(a)}{\leq}\frac{L}{2}\sum_{k=1}^{K}\mathbb{E}_{\mathcal{H}_{k}}\bigg{[}\|\bar{x}_{k}-x^{*}\|^{2}\bigg{]}\\\
&=\frac{L}{2}\sum_{k=1}^{K}D_{k}\\\
&\overset{(b)}{\leq}\Upsilon\frac{L}{2}\sum_{k=1}^{K}\frac{1}{\sqrt{k+1}}\\\
&\overset{(c)}{\leq}\Upsilon\frac{L}{2}\int_{0}^{K}\frac{1}{\sqrt{u+1}}du\\\
&=\Upsilon L(\sqrt{K+1}-1)\\\ \end{split}$

where $(a)$ is due to $\nabla\mathcal{F}(x^{*})=0$, by definition of $x^{*}$,
in (39), $(b)$ is by Theorem 7, and in $(c)$ we perform the same
interpretation of the monotonically decreasing function $\frac{1}{\sqrt{u+1}}$
as in (36).

## References

  * Agarwal et al. (2010) A. Agarwal, O. Dekel, and L. Xiao.  Optimal algorithms for online convex optimization with multi-point bandit feedback.  In _COLT_ , 2010.
  * Bach and Perchet (2016) F. Bach and V. Perchet.  Highly-smooth zero-th order online optimization vianney perchet, 2016\.  URL https://arxiv.org/abs/1605.08165.
  * Bubeck et al. (2021) S. Bubeck, R. Eldan, and Y. T. Lee.  Kernel-based methods for bandit convex optimization.  _J. ACM_ , 68(4), jun 2021.  ISSN 0004-5411.  doi: 10.1145/3453721.  URL https://doi.org/10.1145/3453721.
  * Chen et al. (2019) X. Chen, S. Liu, K. Xu, X. Li, X. Lin, M. Hong, and D. Cox.  Zo-adamm: Zeroth-order adaptive momentum method for black-box optimization.  In _NeurIPS_ , 2019.
  * Daneshmand et al. (2015) A. Daneshmand, F. Facchinei, V. Kungurtsev, and G. Scutari.  Hybrid random/deterministic parallel algorithms for convex and nonconvex big data optimization.  _IEEE Transactions on Signal Processing_ , 63(15):3914–3929, 2015.  doi: 10.1109/TSP.2015.2436357.
  * Doob (1953) J. L. Doob.  Stochastic processes.  1953\.
  * Duchi et al. (2015) J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono.  Optimal rates for zero-order convex optimization: The power of two function evaluations.  _IEEE Transactions on Information Theory_ , 61(5):2788–2806, 2015.  doi: 10.1109/TIT.2015.2409256.
  * Flaxman et al. (2004) A. Flaxman, A. T. Kalai, and H. B. McMahan.  Online convex optimization in the bandit setting: gradient descent without a gradient.  _CoRR_ , cs.LG/0408007, 2004.  URL http://arxiv.org/abs/cs.LG/0408007.
  * Golovin et al. (2019) D. Golovin, J. Karro, G. Kochanski, C. Lee, X. Song, and Q. R. Zhang.  Gradientless descent: High-dimensional zeroth-order optimization.  _CoRR_ , abs/1911.06317, 2019.  URL http://arxiv.org/abs/1911.06317.
  * Gorbunov et al. (2018) E. Gorbunov, P. Dvurechensky, and A. Gasnikov.  An accelerated method for derivative-free smooth stochastic convex optimization, 2018.  URL https://arxiv.org/abs/1802.09022.
  * Hajinezhad et al. (2019) D. Hajinezhad, M. Hong, and A. Garcia.  Zone: Zeroth-order nonconvex multiagent optimization over networks.  _IEEE Transactions on Automatic Control_ , 64(10):3995–4010, 2019.  doi: 10.1109/TAC.2019.2896025.
  * Jamieson et al. (2012) K. G. Jamieson, R. D. Nowak, and B. Recht.  Query complexity of derivative-free optimization.  In _NIPS_ , 2012.
  * Jiang et al. (2022) X. Jiang, X. Zeng, J. Sun, and J. Chen.  Distributed stochastic gradient tracking algorithm with variance reduction for non-convex optimization.  _IEEE Transactions on Neural Networks and Learning Systems_ , pages 1–12, 2022.  doi: 10.1109/TNNLS.2022.3170944.
  * Kumar Sahu et al. (2018) A. Kumar Sahu, D. Jakovetic, D. Bajovic, and S. Kar.  Distributed zeroth order optimization over random networks: A kiefer-wolfowitz stochastic approximation approach.  In _2018 IEEE Conference on Decision and Control (CDC)_ , pages 4951–4958, 2018.  doi: 10.1109/CDC.2018.8619044.
  * LeCun and Cortes (2005) Y. LeCun and C. Cortes.  The mnist database of handwritten digits.  2005\.
  * Li et al. (2022) H. Li, L. Zheng, Z. Wang, Y. Yan, L. Feng, and J. Guo.  S-diging: A stochastic gradient tracking algorithm for distributed optimization.  _IEEE Transactions on Emerging Topics in Computational Intelligence_ , 6(1):53–65, 2022.  doi: 10.1109/TETCI.2020.3017242.
  * Li and Assaad (2021) W. Li and M. Assaad.  Distributed stochastic optimization in networks with low informational exchange.  _IEEE Transactions on Information Theory_ , 67(5):2989–3008, 2021.  doi: 10.1109/TIT.2021.3064925.
  * Liu et al. (2019) S. Liu, P.-Y. Chen, X. Chen, and M. Hong.  signsgd via zeroth-order oracle.  In _ICLR_ , 2019.
  * Liu et al. (2020) S. Liu, P.-Y. Chen, B. Kailkhura, G. Zhang, A. O. Hero III, and P. K. Varshney.  A primer on zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications.  _IEEE Signal Processing Magazine_ , 37(5):43–54, 2020.  doi: 10.1109/MSP.2020.3003837.
  * Lorenzo and Scutari (2016) P. D. Lorenzo and G. Scutari.  Next: In-network nonconvex optimization.  _IEEE Transactions on Signal and Information Processing over Networks_ , 2(2):120–136, 2016.  doi: 10.1109/TSIPN.2016.2524588.
  * Lu et al. (2019) S. Lu, X. Zhang, H. Sun, and M. Hong.  Gnsd: a gradient-tracking based nonconvex stochastic algorithm for decentralized optimization.  In _2019 IEEE Data Science Workshop (DSW)_ , pages 315–321, 2019\.  doi: 10.1109/DSW.2019.8755807.
  * McMahan et al. (2017) B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y. Arcas.  Communication-Efficient Learning of Deep Networks from Decentralized Data.  In A. Singh and J. Zhu, editors, _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , volume 54 of _Proceedings of Machine Learning Research_ , pages 1273–1282. PMLR, 20–22 Apr 2017.  URL https://proceedings.mlr.press/v54/mcmahan17a.html.
  * Nedić et al. (2017) A. Nedić, A. Olshevsky, and W. Shi.  Achieving geometric convergence for distributed optimization over time-varying graphs.  _SIAM Journal on Optimization_ , 27(4):2597–2633, 2017.  doi: 10.1137/16M1084316.  URL https://doi.org/10.1137/16M1084316.
  * Nesterov and Spokoiny (2017) Y. Nesterov and V. G. Spokoiny.  Random gradient-free minimization of convex functions.  _Foundations of Computational Mathematics_ , 17:527–566, 2017.
  * Pu (2020) S. Pu.  A robust gradient tracking method for distributed optimization over directed networks.  In _2020 59th IEEE Conference on Decision and Control (CDC)_ , pages 2335–2341, 2020.  doi: 10.1109/CDC42340.2020.9303917.
  * Pu and Nedić (2018) S. Pu and A. Nedić.  Distributed stochastic gradient tracking methods, 2018.  URL https://arxiv.org/abs/1805.11454.
  * Qu and Li (2018) G. Qu and N. Li.  Harnessing smoothness to accelerate distributed optimization.  _IEEE Transactions on Control of Network Systems_ , 5(3):1245–1260, 2018.  doi: 10.1109/TCNS.2017.2698261.
  * Shamir (2013) O. Shamir.  On the complexity of bandit and derivative-free stochastic convex optimization.  In S. Shalev-Shwartz and I. Steinwart, editors, _Proceedings of the 26th Annual Conference on Learning Theory_ , volume 30 of _Proceedings of Machine Learning Research_ , pages 3–24, Princeton, NJ, USA, 12–14 Jun 2013. PMLR.  URL https://proceedings.mlr.press/v30/Shamir13.html.
  * Shi et al. (2015) W. Shi, Q. Ling, G. Wu, and W. Yin.  Extra: An exact first-order algorithm for decentralized consensus optimization.  _SIAM Journal on Optimization_ , 25(2):944–966, 2015.  doi: 10.1137/14096668X.  URL https://doi.org/10.1137/14096668X.
  * Shi and Eryilmaz (2020) Z. Shi and A. Eryilmaz.  A zeroth-order admm algorithm for stochastic optimization over distributed processing networks.  In _IEEE INFOCOM 2020 - IEEE Conference on Computer Communications_ , pages 726–735, 2020.  doi: 10.1109/INFOCOM41043.2020.9155520.
  * Tang et al. (2021) Y. Tang, J. Zhang, and N. Li.  Distributed zero-order algorithms for nonconvex multiagent optimization.  _IEEE Transactions on Control of Network Systems_ , 8(1):269–281, 2021.  doi: 10.1109/TCNS.2020.3024321.
  * Vemula et al. (2019) A. Vemula, W. Sun, and J. A. Bagnell.  Contrasting exploration in parameter and action space: A zeroth-order optimization perspective.  _ArXiv_ , abs/1901.11503, 2019.
  * Xin et al. (2019) R. Xin, A. K. Sahu, U. A. Khan, and S. Kar.  Distributed stochastic optimization with gradient tracking over strongly-connected networks.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , pages 8353–8358, 2019.  doi: 10.1109/CDC40024.2019.9029217.