Variance-Reduced Gradient Estimator for Nonconvex Zeroth-Order Distributed Optimization

Huaiyi Mu¹, Yujie Tang¹, and Zhongkui Li¹ ¹The authors are with the College of Engineering, Peking University, China (e-mail: [email protected], [email protected], [email protected]).

Abstract

This paper investigates distributed zeroth-order optimization for smooth nonconvex problems. We propose a novel variance-reduced gradient estimator, which randomly renovates one orthogonal direction of the true gradient in each iteration while leveraging historical snapshots for variance correction. By integrating this estimator with gradient tracking mechanism, we address the trade-off between convergence rate and sampling cost per zeroth-order gradient estimation that exists in current zeroth-order distributed optimization algorithms, which rely on either the 2-point or $2d$ -point gradient estimators. We derive a convergence rate of $\mathcal{O}(d^{\frac{5}{2}}/m)$ for smooth nonconvex functions in terms of sampling number $m$ and problem dimension $d$ . Numerical simulations comparing our algorithm with existing methods confirm the effectiveness and efficiency of the proposed gradient estimator.

I INTRODUCTION

We consider a multi-agent system with $N$ agents, where the agents are connected by a communication network that allows them to exchange information for decentralized decision-making. The goal of this group of agents is to collaboratively solve the following consensus optimization problem in a decentralized manner:

\displaystyle\min_{x\in\mathbb{R}^{d}}\ f(x)\coloneqq\frac{1}{N}\sum_{i=1}^{N}f_{i}(x).

(1)

Here $x\in\mathbb{R}^{d}$ is the global decision variable. Each function $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ represents the local objective function for agent $i$ , known only to the agent itself; $f_{i}$ is assumed to be smooth but may be nonconvex. Each agent can only use zeroth-order information of $f_{i}$ during the optimization procedure.

Decentralized optimization has gained significant attention due to its wide range of applications including multi-agent system coordination [1], power systems [2], communication networks [3], etc. For smooth and convex objective functions, the decentralized gradient descent (DGD) algorithm with a decreasing step-size achieved a convergence rate of $\mathcal{O}(\frac{\log t}{\sqrt{t}})$ [4], while [5] and subsequent works proposed gradient tracking (GT) mechanisms with a fixed step-size, achieving a sublinear convergence rate of $\mathcal{O}(\frac{1}{t})$ , which is comparable to centralized gradient descent method. In many real-world applications, the objective functions can be nonconvex, and distributed nonconvex optimization has applications in areas such as machine learning [6] and robotics control [7]. For smooth nonconvex functions, DGD achieves convergence to a stationary point with a rate of $\mathcal{O}(\frac{1}{\sqrt{t}})$ [8, 9], while various GT methods can achieve convergence to a stationary point with a rate of $\mathcal{O}(\frac{1}{t})$ [10].

The aforementioned optimization algorithms rely on first-order information. However, in some scenarios, the gradient is unavailable or expensive to obtain, and agents can only access zeroth-order information of the objective functions. Such scenarios arise in optimization with black-box models [11], optimization with bandit feedback [12], fine-tuning language models [13], etc. To address this issue, various gradient-free methods, particularly algorithms based on zeroth-order gradient estimators, have attracted considerable attention [14]. The work [15] investigated the 2-point zeroth-order gradient estimator, which produces a biased stochastic gradient by using the function values of two randomly sampled points. In [16], the $2d$ -point gradient estimator was proposed, where $d$ is the dimension of the state variable for each agent. [17] combined the 2-point gradient estimator with DGD and the $2d$ -point gradient estimator with GT for nonconvex multi-agent optimization, which lead to convergence rates that are comparable with their first-order counterparts. However, [17] also argued that there appears to be a trade-off between the convergence rate and the sampling cost per zeroth-order gradient estimation when combining zeroth-order gradient estimation techniques with distributed optimization frameworks. This trade-off arises from the high sampling burden of the $2d$ -point estimator and the inherent variance of the 2-point estimator in distributed settings.

To address this trade-off, we aim to design a variance-reduced zeroth-order gradient estimator with a scalable sampling number of function values that is independent of the dimension $d$ . Variance reduction is widely applied in machine learning [18] and stochastic optimization [19]. In [20], variance reduction was used for centralized stochastic gradient descent with strongly convex objectives, achieving a linear convergence rate. [21] applied a 2-point gradient estimator and employed variance reduction for zeroth-order nonconvex centralized stochastic optimization, achieving sublinear convergence. Note that these works only focused on centralized problems. Decentralized finite-sum minimization problems were considered in [22], where variance reduction was used to accelerate convergence; we note that in this work, the variance reduction technique was employed to reduce the variance caused by the finite-sum structure, rather than to address the inherent variance of the 2-point zeroth-order gradient estimator.

In this paper, we propose a new distributed zeroth-order optimization method that incorporates variance reduction techniques as well as the gradient tracking framework, to address the trade-off between convergence rate and sampling cost per zeroth-order gradient estimation in existing zeroth-order distributed optimization algorithms. Specifically, We employ the variance reduction (VR) mechanism to design a novel variance-reduced gradient estimator for distributed nonconvex zeroth-order optimization problems, as formulated in (1). We then combine this new zeroth-order gradient estimation method with the gradient tracking framework, and the resulting algorithm is able to achieve both fast convergence and low sampling cost per zeroth-order gradient estimation. To the best of the authors’ knowledge, this is the first work that attempts to address the aforementioned trade-off for general zeroth-order distributed optimization problems. We also provide rigorous convergence analysis of our proposed algorithm under the smoothness assumption. Although the derived oracle complexities have a higher dependence on the dimension $d$ compared to the algorithms in [17], numerical experiments demonstrate that our proposed algorithm enjoys superior convergence speed and accuracy, reaching lower optimization errors with the same number of samples compared to existing zeroth-order distributed optimization algorithms [17, 23].

Notations: The set of positive integers up to $m$ is denoted as $[m]=\{1,2,\cdots,m\}$ . The $i$ -th component of a vector $x$ is denoted as $[x]_{i}$ . The spectral norm and spectral radius of a matrix $A$ are represented by $\sigma(A)$ and $\rho(A)$ , respectively. For a vector $x\in\mathbb{R}^{d}$ , $\|x\|$ refers to the Euclidean norm, and $\|x\|_{\infty}^{\pi}=\max_{i}\left|[x]_{i}\right|/[\pi]_{i}$ refers to the weighted infinity norm, where $\pi\in\mathbb{R}^{d}$ has all positive components. For a matrix $A$ , $\|A\|_{2}$ represents the matrix norm induced by $\|\cdot\|$ , and $\interleave A\interleave_{\infty}^{\pi}$ represents the matrix norm induced by $\|\cdot\|_{\infty}^{\pi}$ . For two matrices $M$ and $N$ , $M\otimes N$ denotes the Kronecker product. We denote $\mathbb{B}_{d}$ as the closed unit ball in $\mathbb{R}^{d}$ , and $\mathbb{S}_{d-1}=\{x\in\mathbb{R}^{d}:\|x\|=1\}$ as the unit sphere. $\mathcal{U}(\cdot)$ denotes the uniform distribution.

II Formulation And Preliminaries

II-A Problem Formulation

We consider a network consisting of $N$ agents connected via an undirected communication network. The topology of the network is represented by the graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$ , where $\mathcal{N}$ and $\mathcal{E}$ represent the set of agents and communication links, respectively. The distributed consensus optimization problem (1) can be equivalently reformulated as follows:

	$\displaystyle\min_{x_{1},\ldots,x_{N}\in\mathbb{R}^{d}}$	$\displaystyle\frac{1}{N}\sum_{i=1}^{N}f_{i}(x_{i})$		(2)
	s.t.	$\displaystyle x_{1}=x_{2}=\cdots=x_{N},$		(2)

where $x_{i}\in\mathbb{R}^{d}$ now represents the local decision variable of agent $i$ , and the constraint $x_{1}=\cdots=x_{N}$ requires the agents to achieve global consensus for the final decision. During the optimization procedure, each agent may obtain other agents’ information only via exchanging messages with their neighbors in the communication network. We further impose the restriction that only zeroth-order information of the local objective function is available to each agent. In other words, in each iteration, agent $i$ can query the function values of $f_{i}$ at finitely many points.

The following assumption will be employed later in this paper.

Assumption 1.

Each $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ is $L$ -smooth, i.e., we have

\|\nabla f_{i}(x)-\nabla f_{i}(y)\|\leq L\|x-y\|

(3)

for all $x,y\in\mathbb{R}^{d}$ and $i=1,\ldots,N$ . Furthermore, $f^{*}=\inf_{x\in\mathbb{R}^{d}}f(x)>-\infty$ .

II-B Preliminaries on Distributed Zeroth-Order Optimization

When gradient information of the objective function is unavailable, one may construct gradient estimators by sampling the function values at a finite number of points, which has been shown to be a very effective approach by existing literature. We first briefly introduce two types of gradient estimators [17] that are commonly used in noiseless distributed optimization.

Let $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a continuously differentiable function. One version of the 2-point zeroth-order gradient estimator for $\nabla h(x)$ has the following form:

\displaystyle G_{h}^{(2)}(x,u,z)=d\cdot\frac{h(x+uz)-h(x-uz)}{2u}z,

(4)

where $u$ is a positive scalar called the smoothing radius and $z$ is a random vector sampled from the distribution $\mathcal{U}(\mathbb{S}_{d-1})$ . One can show that the expectation of the 2-point gradient estimator is the gradient of a smoothed version of the original function [24, 25], i.e.,

\mathbb{E}_{z\thicksim\mathcal{U}(\mathbb{S}_{d-1})}[G_{h}^{(2)}(x,u,z)]=\nabla h^{u}(x),

where $h^{u}(x)=\mathbb{E}_{y\thicksim\mathcal{U}(\mathbb{B}_{d})}[h(x+uy)]$ . As the smoothing radius $u$ tends to zero, the expectation of the 2-point gradient estimator approaches to the true gradient $\nabla h(x)$ .

By combining the simple 2-point gradient estimator (4) with the decentralized gradient descent framework, one obtains the following algorithm for distributed zeroth-order consensus optimization (2):

x_{i}^{k+1}=\sum_{j=1}^{N}W_{ij}\!\left(x^{k}_{j}-\eta_{t}\,G_{f_{j}}^{(2)}(x^{k}_{j},u^{k},z^{k}_{j})\right),

(5)

which we shall call DGD-2p in this paper. Here $x^{k}_{i}$ denotes the local decision variable of agent $i$ at the $k$ -th iteration, $W\in\mathbb{R}^{N\times N}$ is a weight matrix that is taken to be doubly stochastic, and $\eta_{k}$ is the step-size at iteration $k$ . Since each construction of the 2-point gradient estimator (4) requires sampling only two function values, we can see that DGD-2p can achieves low sampling cost per zeroth-order gradient estimation. However, as shown by [17], DGD-2p achieves a relatively slow convergence rate $O(\sqrt{d/m}\log m)$ , where $m$ denotes the number of function value queries. [17] argued that this slow convergence rate is mainly due to the inherent variance of the 2-point gradient estimator, bounded by

\mathbb{E}_{z\thicksim\mathcal{U}(\mathbb{S}_{d-1})}\!\left[\|G_{h}^{(2)}(x,u,z)\|^{2}\right]\lesssim d\|\nabla h(x)\|^{2}+u^{2}L^{2}d^{2}

under the assumption that function $h$ is $L$ -smooth. In a distributed optimization algorithm, each agent’s local gradient $\nabla f_{i}(x^{k}_{i})$ does not vanish to zero even if the system has reached consensus and optimality. Consequently, the inherent variance of 2-point gradient estimator is inevitable and will considerably slow down the convergence rate.

To achieve higher accuracy for zeroth-order gradient estimation, existing literature has also proposed the $2d$ -point gradient estimator:

\displaystyle G_{h}^{(2d)}(x,u)=\sum_{l=1}^{d}\frac{h(x+ue_{l})-h(x-ue_{l})}{2u}e_{l}.

(6)

Here $e_{l}\in\mathbb{R}^{d}$ is the $l$ -th standard basis vector such that $[e_{l}]_{j}=1$ when $j=l$ and $[e_{l}]_{j}=0$ otherwise. It can be shown that $\|G_{h}^{(2d)}(x,u)-\nabla h(x)\|\leq\frac{1}{2}uL\sqrt{d}$ when $h$ is $L$ -smooth (see, e.g., [17]). Consequently, if we assume the function values of $h$ can be obtained accurately and machine precision issues in numerical computations are ignored, then the $2d$ -point gradient estimator can achieve arbitrarily high accuracy when approximating the true gradient. By combining (6) with the distributed gradient tracking method, one obtains the following algorithm:

	$\displaystyle x^{k+1}_{i}={}$	$\displaystyle\sum_{j=1}^{N}W_{ij}\!\left(x^{k}_{j}-\eta\,s^{k}_{j}\right),$		(7)
	$\displaystyle s^{k+1}_{i}={}$	$\displaystyle\sum_{j=1}^{N}W_{ij}\!\left(s^{k}_{j}\!+\!G_{f_{j}}^{(2d)}(x^{k+1}_{j},u^{k+1})\!-\!G_{f_{j}}^{(2d)}(x^{k}_{j},u^{k})\right),$		(7)

which we shall call GT- $2d$ . Here the auxiliary state variable $s^{k}_{j}$ in (7) tracks the global gradient across iterations. Distributed zeroth-order optimization algorithms that utilize the $2d$ -point gradient estimator, such as GT- $2d$ , can achieve faster convergence due to precise estimation of the true gradients that allows further incorporation of gradient tracking techniques. However, GT- $2d$ has higher sampling cost per gradient estimation compared to DGD-2p: As shown in (6), $2d$ points need to be sampled for each construction of the gradient estimator. This high sampling cost may lead to poor scalability when the dimension $d$ is large.

We remark that the $2d$ -point gradient estimator (6) can also be interpreted as the expectation of the following coordinate-wise gradient estimator:

\displaystyle G_{h}^{(c)}(x,u,l)=d\cdot\frac{h(x+ue_{l})-h(x-ue_{l})}{2u}e_{l},\ \ l\in[d],

(8)

and we have

\displaystyle G_{h}^{(2d)}(x,u)=\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[G_{f}^{(c)}(x,u,l)\right],

(9)

where $\mathcal{U}[d]$ denotes the discrete uniform distribution over the set $\{1,\ldots,d\}$ . The coordinate-wise gradient estimator in (8) shares the similar structure with the 2-point gradient estimator in (4). The key difference is that in (8), we restrict the perturbation direction $ue_{l}$ to lie in the $d$ orthogonal directions associated with the standard basis, instead of uniformly sampled from the unit sphere.

III Our Algorithm

To address the trade-off between convergence rate and sampling cost per gradient estimation in zeroth-order distributed optimization, we employ a variance reduction mechanism [20] to design an improved gradient estimator. The intuition is to combine the best of both worlds, i.e., the precise approximation feature of the $2d$ -point gradient estimator and the low-sampling feature of the 2-point gradient estimator.

Let $k$ denote the iteration number. For each agent, we keep a snapshot point $\tilde{x}_{i}^{k}$ together with its associate smoothing radius $\tilde{u}_{i}^{k}$ at which we have conducted relatively accurate gradient estimation. We then use a random variable $\zeta_{i}^{k}$ generated from the Bernoulli distribution $\mathrm{Ber}(p)$ as an activation indicator for the update of the snapshot point $\tilde{x}_{i}^{k}$ : When $\zeta_{i}^{k}=1$ , agent $i$ updates the snapshot point to be the current iterate, i.e., $\tilde{x}_{i}^{k}=x_{i}^{k}$ and $\tilde{u}_{i}^{k}=u_{i}^{k}$ , and takes a snapshot of the $2d$ -point gradient estimator $G_{f_{i}}^{(2d)}(\tilde{x}_{i}^{k},\tilde{u}_{i}^{k})$ which provides a more accurate gradient estimation at iteration $k$ . When $\zeta_{i}^{k}=0$ , the snapshot point $\tilde{x}_{i}^{k}$ together with the smoothing radius $\tilde{u}^{i}_{k}$ remain unchanged from the previous iteration, i.e., $\tilde{x}_{i}^{k}=\tilde{x}_{i}^{k-1},\tilde{u}_{i}^{k}=\tilde{u}_{i}^{k-1}$ .

In each iteration, agent $i$ constructs the coordinate-wise gradient estimators $G_{f_{i}}^{(c)}(x_{i}^{k},u_{i}^{k},l_{i}^{k})$ and $G_{f_{i}}^{(c)}(\tilde{x}_{i}^{k},\tilde{u}_{i}^{k},l_{i}^{k})$ at the two points $x^{k}_{i}$ and $\tilde{x}^{k}_{i}$ . We then propose the following variance-reduced gradient estimator (VR-GE):

	$\displaystyle g_{f_{i}}^{k}={}$	$\displaystyle G_{f_{i}}^{(c)}(x_{i}^{k},u_{i}^{k},l_{i}^{k})-G_{f_{i}}^{(c)}(\tilde{x}_{i}^{k},\tilde{u}_{i}^{k},l_{i}^{k})$		(10)
		$\displaystyle+G_{f_{i}}^{(2d)}(\tilde{x}_{i}^{k},\tilde{u}_{i}^{k}),$		(10)

where $l_{i}^{k}\in[d]$ is randomly selected by each agent $i$ at each iteration $k$ . By (9), we can see that

	$\displaystyle\mathbb{E}\left[g_{f_{i}}^{k}\mkern 2.0mu\middle\|\mkern 2.0mux^{k}_{i},\tilde{x}^{k}_{i}\right]$
$\displaystyle={}$	$\displaystyle\mathbb{E}\left[G_{f_{i}}^{(c)}(x_{i}^{k},u_{i}^{k},l_{i}^{k})\mkern 2.0mu\middle\|\mkern 2.0mux^{k}_{i},\tilde{x}^{k}_{i}\right]-\mathbb{E}\left[G_{f_{i}}^{(c)}(\tilde{x}_{i}^{k},\tilde{u}_{i}^{k},l_{i}^{k})\mkern 2.0mu\middle\|\mkern 2.0mux^{k}_{i},\tilde{x}^{k}_{i}\right]$
	$\displaystyle+G_{f_{i}}^{(2d)}(\tilde{x}_{i}^{k},\tilde{u}_{i}^{k})$
$\displaystyle={}$	$\displaystyle G_{f_{i}}^{(2d)}(x^{k}_{i},u^{k}_{i}),$	(11)

i.e., VR-GE provides a stochastic gradient of $f_{i}$ with a bias $\|G_{f_{i}}^{(2d)}(x^{k}_{i},u^{k}_{i})-\nabla f_{i}(x^{k}_{i})\|\leq\frac{1}{2}u^{k}_{i}L\sqrt{d}$ , assuming $f_{i}$ is $L$ -smooth.

The expected number of function value samples required per construction of VR-GE is $4+2dp$ . When $p=\frac{C}{2d}$ for some absolute constant $C$ , this becomes $4+C$ which is independent of the dimension $d$ . This gives VR-GE the potential to decrease the sampling cost in large-dimensional zeroth-order distributed optimization by appropriately adjusting the probability $p$ . In the following section, specifically in Lemma 1, we will rigorously analyze the variance of VR-GE and demonstrate its variance reduction property.

In designing our distributed zeroth-order optimization algorithm, we further leverage the gradient tracking mechanisms. Existing literature (including [5, 26, 27], etc.) has demonstrated that gradient tracking mechanisms help mitigate the gap in the convergence rates between distributed optimization and centralized optimization when the objective function is smooth. Drawing inspiration from this advantage, we incorporate the variance-reduced gradient estimator with gradient tracking mechanism to design our algorithm.

The details of the proposed algorithm are outlined in Algorithm 1. Here $\alpha>0$ is the step-size; Steps 1 and 6 implement the gradient tracking mechanism, while Steps 2–5 implement our proposed variance-reduced gradient estimator (10). The convergence guarantees of Algorithm 1 will be provided and discussed in the next section.

Algorithm 1 Distributed Zeroth-Order Optimization Algorithm with Variance Reduced Gradient Tracking Estimator

Initialization :

\tilde{x}_{i}^{0}=x_{i}^{0}\in\mathbb{R}^{d},s_{i}^{0}=g_{f_{i}}^{0}=0

for

k=0,1,2,\cdots

for each

i\in[N]

1. Update

x_{i}^{k+1}

x_{i}^{k+1}=\sum_{j=1}^{N}W_{ij}(x_{j}^{k}-\alpha s_{j}^{k}).

2. Select

l_{i}^{k+1}

uniformly at random from

[d]

3. Generate

\xi_{i}^{k+1}\sim\mathrm{Ber}(p).

4. If

\zeta_{i}^{k+1}=1

, update

\tilde{x}_{i}^{k+1}=x_{i}^{k+1},\,\tilde{u}_{i}^{k+1}=u_{i}^{k+1}

; 4. If

\zeta_{i}^{k+1}=0

, update

\tilde{x}_{i}^{k+1}=\tilde{x}_{i}^{k},\,\tilde{u}_{i}^{k+1}=\tilde{u}_{i}^{k}

5. Construct the VR-GE

g_{f_{i}}^{k+1}

	$\displaystyle g_{f_{i}}^{k+1}={}$	$\displaystyle G_{f_{i}}^{(c)}(x_{i}^{k+1},u_{i}^{k+1},l_{i}^{k+1})$
		$\displaystyle-G_{f_{i}}^{(c)}(\tilde{x}_{i}^{k+1},\tilde{u}_{i}^{k+1},l_{i}^{k+1})$
		$\displaystyle+G_{f_{i}}^{(2d)}(\tilde{x}_{i}^{k+1},\tilde{u}_{i}^{k+1}).$

6. Update

s^{k+1}_{i}

\displaystyle s_{i}^{k+1}={}

\displaystyle\sum_{j=1}^{N}W_{ij}(s_{j}^{k}+g_{f_{j}}^{k+1}-g_{f_{j}}^{k}).

end

IV Main Results

In this section, we present the convergence result of Algorithm 1 under Assumption 1.

For the subsequent analysis, we denote

x^{k}=\begin{bmatrix}x^{k}_{1}\\ \vdots\\ x^{k}_{N}\end{bmatrix},\quad\tilde{x}^{k}=\begin{bmatrix}\tilde{x}^{k}_{1}\\ \vdots\\ \tilde{x}^{k}_{N}\end{bmatrix},\quad s^{k}=\begin{bmatrix}s^{k}_{1}\\ \vdots\\ s^{k}_{N}\end{bmatrix},

and define the following quantities:

	$\displaystyle\delta^{k}$	$\displaystyle=\mathbb{E}[f(\bar{x}^{k})]-f^{*},$	$\displaystyle E_{x}^{k}$	$\displaystyle=\mathbb{E}[\\|x^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k}\\|^{2}],$
	$\displaystyle E_{\tilde{x}}^{k}$	$\displaystyle=\mathbb{E}[\\|\tilde{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k}\\|^{2}],$	$\displaystyle E_{s}^{k}$	$\displaystyle=\mathbb{E}[\\|s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k}\\|^{2}],$

where $\bar{x}^{k}=\frac{1}{N}\sum_{i=1}^{N}x_{i}^{k}$ and $\bar{g}^{k}=\frac{1}{N}\sum_{i=1}^{N}g_{f_{i}}^{k}$ . Here, $\delta^{k}$ quantifies the optimality gap in terms of the objective value, $E_{x}^{k}$ and $E_{\tilde{x}}^{k}$ characterize the consensus errors, and $E_{s}^{k}$ characterizes the tracking error. We also denote $\sigma=\|W-\frac{1}{N}\mathbf{1}_{N}\mathbf{1}_{N}^{T}\|_{2}$ .

Theorem 1.

Under Assumption 1, suppose the parameters of Algorithm 1 satisfy $p\in\big{(}\frac{1-\sigma^{2}}{d},1\big{]}$ , $\sum_{\tau=0}^{\infty}(du_{i}^{\tau})^{2}<\infty$ , $u_{i}^{k}$ is non-increasing, and

\displaystyle\alpha L\leq\min\!\left\{\mfrac{1}{6\sqrt{d}},\mfrac{1}{12\sqrt{d}}\bigg{(}\!\sqrt{\mfrac{2\!+\!\sigma^{2}}{1\!-\!p}}-\sqrt{3}\bigg{)},\mfrac{(1\!-\!\sigma^{2})^{3}}{216\sqrt{29}d^{\frac{5}{2}}}\right\}

Then we have

\displaystyle\frac{1}{k}\sum_{\tau=0}^{k-1}\mathbb{E}[\|\nabla f(\bar{x}^{\tau})\|^{2}]\leq\frac{1}{k}\!\left(\frac{6}{\alpha}\delta^{0}+\frac{6LR_{0}}{\sqrt{29d}N\alpha}+36L^{2}R_{u}\right)\!,

and

		$\displaystyle\frac{1}{k}\sum_{\tau=0}^{k-1}\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}[\\|x_{i}^{\tau}-\bar{x}^{\tau}\\|^{2}]$
	$\displaystyle\leq{}$	$\displaystyle\frac{1}{k}\left(\frac{4\sqrt{29}\delta^{0}}{(1-\sigma^{2})L\sqrt{d}}+\frac{8R_{0}}{(1-\sigma^{2})Nd}+\frac{7R_{u}}{54d^{2}}\right),$
		$\displaystyle\frac{1}{k}\sum_{\tau=0}^{k-1}\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}[\\|s_{i}^{\tau}-\nabla f(\bar{x}^{\tau})\\|^{2}]$
	$\displaystyle\leq{}$	$\displaystyle\frac{1}{k}\left(\frac{580\delta^{0}}{(1-\sigma^{2})^{2}\alpha}+\frac{40\sqrt{29}LR_{0}}{(1-\sigma^{2})N\alpha\sqrt{d}}+\frac{\sqrt{29}LR_{u}}{10\alpha d^{\frac{3}{2}}}\right),$

where $R_{0}=\frac{d}{1-\sigma^{2}}E_{x}^{0}$ , $R_{u}=\frac{1}{p}(du_{i}^{0})^{2}+\sum_{\tau=1}^{\infty}(du_{i}^{\tau})^{2}$ .

Proof.

See Section V. ∎

Remark 1.

The convergence rate of Algorithm 1 under Assumption 1 is $\mathcal{O}(\frac{1}{k})$ , which aligns with the rate achieved for distributed nonconvex optimization with gradient tracking using first-order information [10]. In addition, each iteration of VR-GE requires $4+2dp$ function value queries on average. As long as $p<1-\frac{2}{d}$ , the averaged sampling number for VR-GE is less then that for the $2d$ -point gradient estimator.

Remark 2.

Assuming $p\propto\frac{1}{d}$ and that all agents start from the same initial point, the convergence rate becomes $\mathcal{O}\big{(}d^{\frac{5}{2}}/m\big{)}$ with respect to the number of function value queries $m$ , which can be justified by simple algebraic calculation. Although this rate is not as favorable as the $\mathcal{O}(d/m)$ rate of GT- $2d$ in terms of the dependence on the problem dimension $d$ , we suspect that this discrepancy is primarily due to the complicated analysis procedure; the conditions in Theorem 1 serve as sufficient but not necessary conditions for convergence, meaning that the $\mathcal{O}\big{(}d^{\frac{5}{2}}/m\big{)}$ rate may not be a theoretically tight bound. As demonstrated in the simulation section, Algorithm 1 converges faster than GT- $2d$ and achieves higher accuracy with the same number of samples. Deriving a tighter convergence result remains a topic for future work.

Next, we present the convergence rate when the probability $p$ is fixed as $1/d$ in the following corollary.

Corollary 1.

Under Assumption 1, suppose the parameters of Algorithm 1 satisfy $p=1/d$ , $\sum_{\tau=0}^{\infty}(du_{i}^{\tau})^{2}<\infty$ , $u_{i}^{k}$ is non-increasing, and

\alpha L\leq\min\left\{\mfrac{1}{6\sqrt{d}},\mfrac{\sqrt{3}}{12\sqrt{d}}\bigg{(}\!{-1}+\sqrt{1\!+\!\mfrac{\sigma^{2}}{d\!-\!1}}\bigg{)},\mfrac{(1-\sigma^{2})^{3}}{264\sqrt{29}d^{\frac{3}{2}}}\right\}.

Then we have

\displaystyle\frac{1}{k}\sum_{\tau=0}^{k-1}\mathbb{E}[\|\nabla f(\bar{x}^{\tau})\|^{2}]\leq{}\frac{1}{k}\!\left(\frac{6}{\alpha}\delta^{0}+\frac{6\sqrt{d}L\tilde{R}_{0}}{\sqrt{29}N\alpha}+36L^{2}R_{u}\right)\!,

(12)

where $\tilde{R}_{0}=\frac{1}{1-\sigma^{2}}E_{x}^{0}$ .

Proof.

See Appendix F. ∎

Remark 3.

With some standard algebraic calculation, it can be shown that the convergence rate in (12) is $\mathcal{O}\big{(}d^{\frac{3}{2}}/m\big{)}$ in terms of the number of function value queries $m$ , which improves upon the $\mathcal{O}\big{(}d^{\frac{5}{2}}/m\big{)}$ rate in Theorem 1. Roughly speaking, this improvement is due to that, by fixing $p=1/d$ , the step-size’s dependency on the dimension can be decreased from $1/d^{\frac{5}{2}}$ to $1/d^{\frac{3}{2}}$ . In light of this improvement, setting $p=1/d$ offers practical guidance for selecting the probability $p$ .

V Outline of Convergence Analysis

V-A Bounding the Variance of VR-GE

The variance of VR-GE is essential for convergence proof of Algorithm 1 and we provide analysis details in this subsection. We first rewrite Algorithm 1 as follows:


	$\displaystyle x^{k+1}=(W\otimes I_{d})(x^{k}-\alpha s^{k}),$		(13a)
	$\displaystyle s^{k+1}=(W\otimes I_{d})(s^{k}+g^{k+1}-g^{k}),$		(13b)

where $g^{k}=[(g_{f_{1}}^{k})^{T},(g_{f_{2}}^{k})^{T},\cdots,(g_{f_{N}}^{k})^{T}]^{T}$ .

We now derive a bound on the expected difference between variance-reduced gradient estimator and the true gradient in the following lemma.

Lemma 1.

Let $h:\mathbb{R}^{d}\to\mathbb{R}$ be $L$ -smooth. Then for any $x,\tilde{x},y\in\mathbb{R}^{d}$ and $0<u\leq\tilde{u}$ , it holds that

		$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|g-\nabla h(x)\\|^{2}\right]$		(14)
	$\displaystyle\leq{}$	$\displaystyle 12dL^{2}\\|x-y\\|^{2}+12dL^{2}\\|\tilde{x}-y\\|^{2}+\frac{7}{2}\tilde{u}^{2}L^{2}d^{2},$		(14)

where $g=G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)+G_{h}^{(2d)}(\tilde{x},\tilde{u})$ .

Proof.

Using $\|a-b\|^{2}\leq 2\|a\|^{2}+2\|b\|^{2}$ , we derive that

		$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|g-\nabla h(x)\\|^{2}\right]$
	$\displaystyle={}$	$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\Big{[}\Big{\\|}\left(G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)\right)$
		$\displaystyle\!-\!\Big{(}G_{h}^{(2d)}\!(x,u)\!-\!G_{h}^{(2d)}\!(\tilde{x},\tilde{u})\Big{)}\!-\!\left(\nabla h(x)\!-\!G_{h}^{(2d)}\!(x,u)\!\right)\!\!\Big{\\|}^{2}\Big{]}$
	$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}_{l\sim\mathcal{U}[d]}\!\Big{[}\Big{\\|}\left(G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)\right)$
		$\displaystyle-\left(G_{h}^{(2d)}(x,u)-G_{h}^{(2d)}(\tilde{x},\tilde{u})\right)\Big{\\|}^{2}\Big{]}$
		$\displaystyle+2\mathbb{E}[\\|\nabla h(x)-G_{h}^{(2d)}(x,u)\\|^{2}].$

Then, by $\mathbb{E}[\|\zeta-\mathbb{E}[\zeta]\|^{2}]\leq\mathbb{E}[\|\zeta\|^{2}]$ for an arbitrary random vector $\zeta$ , we can derive from the above inequality that

		$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|g-\nabla h(x)\\|^{2}\right]$
	$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)\\|^{2}\right]+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle={}$	$\displaystyle 2d\\|G_{h}^{(2d)}(x,u)-G_{h}^{(2d)}(\tilde{x},\tilde{u})\\|^{2}+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle={}$	$\displaystyle 2d\\|(G_{h}^{(2d)}(x,u)-\nabla h(x))-(G_{h}^{(2d)}(\tilde{x},\tilde{u})-\nabla h(\tilde{x}))$
		$\displaystyle+(\nabla h(x)-\nabla h(\tilde{x}))\\|^{2}+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle\leq{}$	$\displaystyle 6dL^{2}\\|x-\tilde{x}\\|^{2}+\frac{3}{2}u^{2}L^{2}d^{2}+\frac{3}{2}\tilde{u}^{2}L^{2}d^{2}+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle\leq{}$	$\displaystyle 12dL^{2}\\|x-y\\|^{2}+12dL^{2}\\|\tilde{x}-y\\|^{2}+\frac{7}{2}\tilde{u}^{2}L^{2}d^{2},$

where we have used (3) and inequality $\|G_{h}^{(2d)}(x,u)-\nabla h(x)\|\leq\frac{1}{2}uL\sqrt{d}$ in the second inequality, and the fact that $u\leq\tilde{u}$ in the last inequality. ∎

Remark 4.

The estimation error of VR-GE (i.e., the left-hand side of (14)) comprises two parts: consensus error and the smoothing radius effect. Consensus error is inherent due to the distributed nature of the algorithm, while the smoothing radius can be tuned properly. As Algorithm 1 approaches consensus and the smoothing radius approaches zero, the estimation error between the VR-GE and the true gradient diminishes. Thus, VR-GE offers reduced variance compared to the 2-point gradient estimator while requiring fewer samples than the $2d$ -point gradient estimator.

V-B Proof Sketch of Theorem 1

The proof relies on four lemmas. The first lemma analyzes the evolution of function value $f(\bar{x}^{k})$ using $L$ -smooth property. We denote $\tilde{u}^{k}=\max_{i}\{\mathbb{E}[\tilde{u}_{i}^{k}]\},i\in[N]$ and assume $d\geq 3$ in the following analysis.

Lemma 2.

Suppose $\alpha L\leq\frac{1}{6\sqrt{d}}$ , we have

	$\displaystyle\delta^{k+1}\leq{}$	$\displaystyle\delta^{k}-\frac{\alpha}{3}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+\frac{18\alpha dL^{2}}{N}E_{x}^{k}$
		$\displaystyle+\frac{16\alpha dL^{2}}{N}E_{\tilde{x}}^{k}+5\alpha L^{2}d^{2}(\tilde{u}^{k})^{2}.$		(15)

Proof.

See Appendix B. ∎

Lemma 2 derives a bound for optimization error $\delta^{k}$ . We need to further bound consensus errors: $E_{x}^{k}$ and $E_{\tilde{x}}^{k}$ . This is tackled by the following lemma.

Lemma 3.

Suppose we choose $p\in\big{(}\frac{1-\sigma^{2}}{d},1\big{]}$ and $\alpha L\leq\min\left\{\frac{(1-\sigma^{2})^{2}}{12\sigma^{2}(1+2\sigma^{2})\sqrt{d(1+p)}},\\ \frac{1}{6\sqrt{d}},\frac{\sqrt{3}}{12d}\!\left(-\sqrt{d}\!+\!\sqrt{\frac{d-1+\sigma^{2}}{1-p}}\right)\!\right\}$ . Then we have the following inequality:

v^{k+1}\leq Av^{k}+B_{k},

(16)

where

	$\displaystyle v^{k}={}$	$\displaystyle\!\left[E_{x}^{k},E_{\tilde{x}}^{k},\mfrac{\alpha}{3\sqrt{29}L\sqrt{d}}E_{s}^{k}\right]^{T},$
	$\displaystyle A={}$	$\displaystyle\!\begin{bmatrix}\frac{1+2\sigma^{2}}{3}&0&\frac{9\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L\\ 17\sqrt{d}\alpha L+\frac{1+2\sigma^{2}}{3}&1-\frac{1-\sigma^{2}}{d}&\frac{9\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L\\ \frac{9\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L&\frac{3\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L&\frac{2+\sigma^{2}}{3}\end{bmatrix},$
	and	$\displaystyle\!\!\!B_{k}\!=\!\begin{bmatrix}0\\ \frac{N\alpha}{L\sqrt{d}}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+5Nd^{\frac{3}{2}}\alpha L(\tilde{u}^{k})^{2}\\ \frac{\sqrt{29}N\alpha}{3(1-\sigma^{2})L\sqrt{d}}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]\!+\!\frac{2\sqrt{29}Nd^{\frac{3}{2}}\alpha L}{1-\sigma^{2}}(\tilde{u}^{k})^{2}\\ \end{bmatrix}\!\!.$

Proof.

See Appendix C. ∎

For subsequent analysis, we need to bound the induced weighted matrix norm of the matrix $A$ in the inequality (16).

Lemma 4.

Consider the matrix $A$ defined in Lemma 3, and suppose that $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{12\sqrt{29}d^{\frac{5}{2}}}$ . Then

\rho(A)\leq\interleave A\interleave_{\infty}^{\pi}\leq 1-\frac{1-\sigma^{2}}{2d},

(17)

where $\pi=[\frac{1-\sigma^{2}}{d},3,1]^{T}$ .

Proof.

See Appendix D. ∎

Next, we derive a bound on the accumulated consensus errors over iterations using Lemma 3.

Lemma 5.

Suppose we choose $p\in(\frac{1-\sigma^{2}}{d},1]$ and $\alpha L\leq\min\!\left\{\frac{1}{6\sqrt{d}},\frac{(1-\sigma^{2})^{3}}{12\sqrt{29}d^{\frac{5}{2}}},\frac{\sqrt{3}}{12d}\left(-\sqrt{d}+\sqrt{\frac{d-1+\sigma^{2}}{1-p}}\right)\!\right\}$ . We have

	$\displaystyle\max\left\{\frac{d}{1-\sigma^{2}}\sum_{\tau=0}^{k}E_{x}^{\tau},\frac{1}{3}\sum_{\tau=0}^{k}E_{\tilde{x}}^{\tau},\frac{\alpha}{3\sqrt{29}L\sqrt{d}}\sum_{\tau=0}^{k}E_{s}^{\tau}\right\}$
$\displaystyle\leq{}$	$\displaystyle\frac{4d}{1-\sigma^{2}}R_{0}+\frac{2\sqrt{29}\sqrt{d}N\alpha}{3(1-\sigma^{2})^{2}L}\sum_{m=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{m})\\|^{2}]$
	$\displaystyle+\frac{4\sqrt{29}Nd^{\frac{5}{2}}\alpha L}{(1-\sigma^{2})^{2}}\sum_{m=0}^{k-1}(\tilde{u}^{m})^{2},$	(18)

where $R_{0}=\frac{d}{1-\sigma^{2}}E_{x}^{0}$ . We also derive that $\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}\leq\frac{1}{p}(u_{i}^{0})^{2}+\sum_{\tau=1}^{\infty}(u_{i}^{\tau})^{2}$ .

Proof.

See Appendix E. ∎

Based on Lemmas 2 and 5, we are now ready to prove Theorem 1. Plugging the inequality (18) into (15) leads to

$\displaystyle 0\leq{}$	$\displaystyle\delta^{0}-\frac{\alpha}{3}\sum_{\tau=0}^{k}\mathbb{E}[\\|\nabla f(\bar{x}^{\tau})\\|^{2}]+5\alpha L^{2}d^{2}\sum_{\tau=0}^{k}(\tilde{u}^{\tau})^{2}$
	$\displaystyle+\frac{54\alpha dL^{2}}{N}(\frac{4d}{1-\sigma^{2}}R_{0}+\frac{4\sqrt{29}Nd^{\frac{5}{2}}\alpha L}{(1-\sigma^{2})^{2}}\sum_{m=0}^{k-1}(\tilde{u}^{m})^{2}$
	$\displaystyle+\frac{2\sqrt{29}\sqrt{d}N\alpha}{3(1-\sigma^{2})^{2}L}\sum_{m=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{m})\\|^{2}])$
$\displaystyle\leq{}$	$\displaystyle\delta^{0}+\frac{L}{\sqrt{29}N\sqrt{d}}R_{0}-\frac{\alpha}{6}\sum_{\tau=0}^{k}\mathbb{E}[\\|\nabla f(\bar{x}^{\tau})\\|^{2}]$
	$\displaystyle+6\alpha L^{2}d^{2}\sum_{\tau=0}^{k}(\tilde{u}^{\tau})^{2},$	(19)

where we have used $\frac{18\alpha dL^{2}}{N}\cdot\frac{1-\sigma^{2}}{d}+\frac{16\alpha dL^{2}}{N}\cdot 3<\frac{54\alpha dL^{2}}{N}$ in the first inequality and $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{216\sqrt{29}d^{\frac{5}{2}}}\leq\frac{1-\sigma^{2}}{216\sqrt{29}d^{\frac{3}{2}}}$ in the last inequality.

Using the bound $\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}$ in Lemma 5 and $(u_{i}^{k})^{2}$ is summable, we see that $(\tilde{u}^{k})^{2}$ is also summable. Consequently, we have $\mathbb{E}[\|\nabla f(\bar{x}^{k})\|^{2}]\rightarrow 0$ as $k\rightarrow\infty$ by (19). We can further derive from the inequality (19) that

		$\displaystyle\frac{1}{k}\sum_{\tau=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{\tau})\\|^{2}]$		(20)
	$\displaystyle\leq{}$	$\displaystyle\frac{1}{k}\left[\frac{6}{\alpha}\delta^{0}+\frac{6L}{\sqrt{29}\sqrt{d}N\alpha}R_{0}+36L^{2}d^{2}\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}\right].$		(20)

Combining (20) with (18), we have

\frac{1}{N}\sum_{\tau=0}^{\infty}E_{x}^{\tau}\leq\frac{4\sqrt{29}\delta^{0}}{(1\!-\!\sigma^{2})L\sqrt{d}}+\frac{8R_{0}}{(1\!-\!\sigma^{2})Nd}+\frac{7}{54}\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2},

(21)

where we have used $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{216\sqrt{29}d^{\frac{5}{2}}}$ . Similarly,

		$\displaystyle\frac{1}{N}\sum_{\tau=0}^{\infty}\mathbb{E}[\\|s^{\tau}-\mathbf{1}_{N}\otimes\nabla f(\bar{x}^{\tau})\\|^{2}]$
	$\displaystyle\leq{}$	$\displaystyle\frac{3}{2N}\sum_{\tau=0}^{\infty}\mathbb{E}[\\|s^{\tau}\!-\!\mathbf{1}_{N}\otimes\bar{g}^{\tau}\\|^{2}]+3\sum_{\tau=0}^{\infty}\mathbb{E}[\\|\bar{g}^{\tau}\!-\!\nabla f(\bar{x}^{\tau})\\|^{2}]$
	$\displaystyle\leq{}$	$\displaystyle\frac{3}{2N}\!\sum_{\tau=0}^{\infty}\!E_{s}^{\tau}\!+\!3\sum_{\tau=0}^{\infty}\!\left[\frac{26dL^{2}}{N}E_{x}^{\tau}\!+\!\frac{24dL^{2}}{N}E_{\tilde{x}}^{\tau}\!+\!7(\tilde{u}^{\tau})^{2}L^{2}d^{2}\right]$
	$\displaystyle\leq{}$	$\displaystyle\frac{5\sqrt{29}\sqrt{d}L}{N\alpha}\bigg{(}\frac{4d}{1-\sigma^{2}}R_{0}+\frac{4\sqrt{29}Nd^{\frac{5}{2}}\alpha L}{(1-\sigma^{2})^{2}}\sum_{m=0}^{k-1}(\tilde{u}^{m})^{2}$
		$\displaystyle\!\!\!+\!\frac{2\sqrt{29}dN\alpha}{3(1\!-\!\sigma^{2})^{2}L\sqrt{d}}\!\sum_{m=0}^{k-1}\!\mathbb{E}[\\|\nabla f(\bar{x}^{m})\\|^{2}]\bigg{)}\!+\!21d^{2}L^{2}\!\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}$
	$\displaystyle\leq{}$	$\displaystyle\frac{580\delta^{0}}{(1-\sigma^{2})^{2}\alpha}+\frac{40\sqrt{29}LR_{0}}{(1\!-\!\sigma^{2})N\alpha\sqrt{d}}+\frac{\sqrt{29}\sqrt{d}L}{10\alpha}\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2},$

where we have used inequality (18) and $\frac{3}{2N}\cdot\frac{3\sqrt{29}L\sqrt{d}}{\alpha}+\frac{78dL^{2}}{N}\cdot\frac{1-\sigma^{2}}{d}+\frac{72dL^{2}}{N}\cdot 3\leq\frac{5\sqrt{29}\sqrt{d}L}{N\alpha}$ in the third inequality.

The proof of Theorem 1 can now be concluded.

VI Simulation

We consider a multi-agent nonconvex optimization problem adapted from [17] with $N=50$ agents in the network, and the objective function of each agent is given as follows:

\displaystyle f_{i}(x)=\frac{\alpha_{i}}{1+e^{-\zeta_{i}^{T}x-v_{i}}}+\beta_{i}\ln(1+\|x\|^{2}),

(22)

where $\alpha_{i},\beta_{i},v_{i}\in\mathbb{R}$ are randomly generated parameters satisfying $\frac{1}{N}\sum_{i}\beta_{i}=1$ , each $\zeta_{i}\in\mathbb{R}^{d}$ is also randomly generated, and the dimension $d$ is set to $64$ .

For the following numerical simulation of Algorithm 1, we set the step-size $\eta=0.02$ and the smoothing radius $u_{i}^{k}=3/k^{\frac{3}{4}}$ . All agents start from the same initial points to ensure consistency in the initial conditions across the network.

VI-A Comparison with Other Algorithms

Fig. 1 compares Algorithm 1 with DGD-2p, GT- $2d$ , and ZONE-M (with $J$ =100). In the figure, the probability used for Algorithm 1 is $p=0.1$ . The horizontal axis is normalized and represents the sampling number $m$ (i.e., the number of zeroth-order queries). The two sub-figures illustrate the stationarity gap $\|\nabla f(\bar{x}^{k})\|^{2}$ and the consensus error $\frac{1}{N}\sum_{i}\|x_{i}^{k}-\bar{x}^{k}\|^{2}$ , respectively.

By inspecting Fig. 1, we first see that DGD-2p converges faster than ZONE-M with $J=100$ . When comparing DGD-2p and GT- $2d$ , we can see a clear difference between their convergence behavior: DGD-2p achieves fast convergence initially but slows down afterwards due to the inherent variance of the 2-point gradient estimator, whereas GT- $2d$ achieves higher eventual accuracy but slower initial convergence before approximately $1.5\times 10^{4}$ zeroth-order queries due to the higher sampling burden of the $2d$ -point gradient estimator.

As demonstrated in Fig. 1, Algorithm 1 offers both high eventual accuracy and a fast convergence rate in terms of stationarity gap and consensus error. This improvement is attributed to the variance reduction mechanism employed in designing VR-GE, which effectively balances the sampling number and expected variance, thereby addressing the trade-off between convergence rate and sampling cost per zeroth-order gradient estimation that exists in current zeroth-order distributed optimization algorithms.

Refer to caption — Figure 1: Convergence of Algorithm 1, ZONE-M with J =100, DGD-2p, GT- $2d$ .

VI-B Comparison of Algorithm 1 under Different Probabilities

Fig. 2 compares the convergence of Algorithm 1 under different choices of the probability $p$ , which reflects the frequency with which each agent takes snapshots. The three sub-figures illustrate the stationarity gap $\|\nabla f(\bar{x}^{k})\|^{2}$ , the consensus error $\frac{1}{N}\sum_{i}\|x_{i}^{k}-\bar{x}^{k}\|^{2}$ , and the tracking error $\frac{1}{N}\sum_{i}\|s_{i}^{k}-\nabla f(\bar{x}^{k})\|^{2}$ , respectively.

The results demonstrate that Algorithm 1 with a lower probability achieves better accuracy with fewer sampling numbers. However, a lower probability also results in more fluctuation during convergence. This is expected because, with a lower probability, the snapshot variables are updated less frequently, leading to a greater deviation from the true gradient as iterations progress.

Two notable cases are $p=0$ and $p=1$ . When $p=1$ , VR-GE behaves similarly to GT- $2d$ , taking snapshots at every iteration, which results in poorer convergence performance compared to cases where $p\in(0,1)$ . Conversely, when $p=0$ , agents do not take any snapshots, and $\tilde{x}_{i}^{k}$ remains its initial value $x_{i}^{0}$ . In this scenario, VR-GE uses solely the initial iterate for gradient correction, which introduces variance and reduces convergence accuracy, as illustrated in Fig. 2.

VI-C Comparison of Algorithm 1 under Different Dimensions

Fig. 3 compares the convergence of Algorithm 1 across different agent dimensions, alongside varying probabilities for taking snapshots within the algorithm. The results show that Algorithm 1 can effectively handle different scenarios, such as when $d=300$ , achieving stationarity gaps that are below $10^{-6}$ .

As the dimension increases, VR-GE requires more samples to accurately estimate the gradient. To maintain similar convergence performance across higher dimensions, the probability $p$ for taking snapshots can be adjusted to lower values. As shown in Fig. 3, decreasing the probability as the dimension grows allows Algorithm 1 to achieve a convergence rate and optimization accuracy that are comparable to cases with lower dimensions. However, this adjustment also leads to increased fluctuation during the convergence process. This fluctuation is a result of the randomness introduced by the snapshot mechanism.

VII Conclusion

In this paper, we proposed an improved variance-reduced gradient estimator and integrated it with gradient tracking mechanism for nonconvex distributed zeroth-order optimization problems. Through rigorous analysis, we demonstrated that our algorithm achieves sublinear convergence for smooth nonconvex functions that is comparable with first-order gradient tracking algorithms, while maintaining relatively low sampling cost per gradient estimation. Comparative evaluations with existing distributed zeroth-order optimization algorithms verified the effectiveness of the proposed gradient estimator. Future work will focus on refining the step-size bounds of the algorithm and reducing the dependence of the convergence rate on the dimension $d$ .

Appendix A. Auxiliary Lemmas

This section summarizes some auxiliary lemmas for convergence analysis.

Lemma 6.

Suppose $f:\mathbb{R}^{d}\to\mathbb{R}$ is $L$ -smooth and $\inf_{x\in\mathbb{R}^{d}}f(x)=f^{*}>-\infty$ . Then

\|\nabla f(x)\|^{2}\leq 2L(f(x)-f(x^{*})).

(23)

Proof.

This lemma is standard in optimization theory, and can be shown via

f^{\ast}\leq f\!\left(x-\mfrac{1}{L}\nabla f(x)\right)\leq f(x)-\mfrac{1}{2L}\|\nabla f(x)\|^{2},

where the last step follows from the $L$ -smoothness of $f$ . ∎

Lemma 7 ([17]).

Let $f:\mathbb{R}^{d}\to\mathbb{R}$ be $L$ -smooth. Then for any $x\in\mathbb{R}^{d}$ ,

\|G_{f}^{(2d)}(x,u)-\nabla f(x)\|\leq\frac{1}{2}uL\sqrt{d}.

(24)

Lemma 8 ([26]).

Let $\sigma\triangleq\|W-\frac{1}{N}\mathbf{1}_{N}\mathbf{1}_{N}^{T}\|_{2}<1$ . For any $x_{1},\cdots,x_{N}\in\mathbb{R}^{d}$ , we have

\|(W\otimes I_{d})(x-\mathbf{1}_{N}\otimes\bar{x})\|\leq\sigma\|x-\mathbf{1}_{N}\otimes\bar{x}\|,

where we denote $x=\begin{bmatrix}x_{1}\\ \vdots\\ x_{N}\end{bmatrix},\bar{x}=\dfrac{1}{N}\sum_{i=1}^{N}x_{i}$ .

Lemma 9 ([28]).

Let $A\in\mathbb{R}^{n\times n}$ be non-negative and $\pi\in\mathbb{R}^{n}$ be positive. If $A\pi\leq\beta\pi$ for $\beta>0$ , then $\rho(A)\leq\interleave A\interleave_{\infty}^{\pi}\leq\beta$ .

Appendix B. Proof of Lemma 2

First, by left multiplying $\mathbf{1}_{N}^{T}\otimes T_{d}$ on both sides of equation (13b), and using the initialization $s^{0}=g^{0}$ , we derive that

\bar{s}^{k}=\bar{g}^{k},

where $\bar{s}^{k}=\frac{1}{N}\sum_{i=1}^{N}s_{i}^{k}$ and $\bar{g}^{k}=\frac{1}{N}\sum_{i=1}^{N}g_{f_{i}}^{k}$ .

Then, from (13a), we have

\bar{x}_{k+1}=\bar{x}^{k}-\alpha\bar{g}^{k}.

(25)

Leveraging the $L$ -smoothness of the function $f$ , we have

$\displaystyle\mathbb{E}[f(\bar{x}_{k+1})]\leq{}$	$\displaystyle\mathbb{E}[f(\bar{x}^{k})]\!-\!\alpha\mathbb{E}[\langle\nabla f(\bar{x}^{k}),\bar{g}^{k}\rangle]\!+\!\frac{\alpha^{2}\!L}{2}\mathbb{E}[\\|\bar{g}^{k}\\|^{2}]$	(26)
$\displaystyle={}$	$\displaystyle\mathbb{E}[f(\bar{x}^{k})]-\alpha\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+\frac{\alpha^{2}L}{2}\mathbb{E}[\\|\bar{g}^{k}\\|^{2}]$
	$\displaystyle-\alpha\mathbb{E}[\langle\nabla f(\bar{x}^{k}),\frac{1}{N}\sum\limits_{i=1}^{N}(g_{f_{i}}^{k}-\nabla f_{i}(\bar{x}^{k}))\rangle]$
$\displaystyle\leq{}$	$\displaystyle\mathbb{E}[f(\bar{x}^{k})]\!-\!\frac{\alpha}{2}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]\!+\!\frac{\alpha^{2}L}{2}\mathbb{E}[\\|\bar{g}^{k}\\|^{2}]$
	$\displaystyle+\frac{\alpha}{2}\mathbb{E}\bigg{[}\bigg{\\|}\frac{1}{N}\sum\limits_{i=1}^{N}(g_{f_{i}}^{k}-\nabla f_{i}(\bar{x}^{k}))\bigg{\\|}^{2}\bigg{]}.$

For the last term on the right-hand side, we can derive by using Lemma 1 that

	$\displaystyle\mathbb{E}\bigg{[}\bigg{\\|}\frac{1}{N}\sum\limits_{i=1}^{N}(g_{f_{i}}^{k}-\nabla f_{i}(\bar{x}^{k}))\bigg{\\|}^{2}\bigg{]}$	(27)
$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}\bigg{[}\bigg{\\|}\frac{1}{N}\sum_{i=1}^{N}\left(\nabla f_{i}(x_{i}^{k})-\nabla f_{i}(\bar{x}^{k})\right)\bigg{\\|}^{2}\bigg{]}$
	$\displaystyle+\frac{2}{N}\mathbb{E}\bigg{[}\sum_{i=1}^{N}\\|g_{f_{i}}^{k}-\nabla f_{i}(x_{i}^{k})\\|^{2}\bigg{]}$
$\displaystyle\leq{}$	$\displaystyle\frac{2L^{2}}{N}E_{x}^{k}+\frac{2}{N}(12dL^{2}E_{x}^{k}+12dL^{2}E_{\tilde{x}}^{k}+\frac{7}{2}N(\tilde{u}^{k})^{2}L^{2}d^{2})$
$\displaystyle={}$	$\displaystyle\frac{2L^{2}(1+12d)}{N}E_{x}^{k}+\frac{24dL^{2}}{N}E_{\tilde{x}}^{k}+7(\tilde{u}^{k})^{2}L^{2}d^{2}.$

Next, we bound the term $\mathbb{E}[\|\bar{g}^{k}\|^{2}]$ as follows:

	$\displaystyle\mathbb{E}[\\|\bar{g}^{k}\\|^{2}]$
$\displaystyle={}$	$\displaystyle\mathbb{E}[\\|\frac{1}{N}\sum_{i=1}^{N}g_{f_{i}}^{k}\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+2\mathbb{E}[\\|\frac{1}{N}\sum_{i=1}^{N}\left(g_{f_{i}}^{k}-\nabla f_{i}(\bar{x}^{k})\right)\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+\frac{4L^{2}(1+12d)}{N}E_{x}^{k}+\frac{48dL^{2}}{N}E_{\tilde{x}}^{k}$
	$\displaystyle+14(\tilde{u}^{k})^{2}L^{2}d^{2},$	(28)

where we have used (27) in the second inequality.

Hence, plugging inequalities (27) and (28) into (26), we have

		$\displaystyle\mathbb{E}[f(\bar{x}_{k+1})]$
	$\displaystyle\leq{}$	$\displaystyle\mathbb{E}[f(\bar{x}^{k})]-\frac{\alpha}{2}(1-2\alpha L)\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]$
		$\displaystyle+\frac{7\alpha(\tilde{u}^{k})^{2}L^{2}d^{2}}{2}(1\!+\!2\alpha L)+\frac{\alpha L^{2}(1\!+\!12d)}{N}(1\!+\!2\alpha L)E_{x}^{k}$
		$\displaystyle+\frac{12\alpha dL^{2}}{N}(1+2\alpha L)E_{\tilde{x}}^{k}$
	$\displaystyle\leq{}$	$\displaystyle\mathbb{E}[f(\bar{x}^{k})]-\frac{\alpha}{3}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+\frac{18\alpha dL^{2}}{N}E_{x}^{k}$
		$\displaystyle+\frac{16\alpha dL^{2}}{N}E_{\tilde{x}}^{k}+5\alpha(\tilde{u}^{k})^{2}L^{2}d^{2},$

where we have used $\alpha L\leq\frac{1}{6}$ in the second inequality.

Appendix C. Proof of Lemma 3

First, following from (13) and (25), we derive that

$\displaystyle E_{x}^{k+1}=$	$\displaystyle\mathbb{E}[\\|x^{k+1}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}]$	(29)
$\displaystyle={}$	$\displaystyle\mathbb{E}[\\|(W\otimes I_{d})[x^{k}\!-\!\mathbf{1}_{N}\otimes\bar{x}^{k}\!-\!\alpha(s^{k}\!-\!\mathbf{1}_{N}\otimes\bar{g}^{k})]\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle\sigma^{2}\left(1+\frac{1-\sigma^{2}}{3\sigma^{2}}\right)\mathbb{E}[\\|x^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k}\\|^{2}]$
	$\displaystyle+\sigma^{2}\left(1+\frac{3\sigma^{2}}{1-\sigma^{2}}\right)\alpha^{2}\mathbb{E}[\\|s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k}\\|^{2}]$
$\displaystyle={}$	$\displaystyle\frac{1+2\sigma^{2}}{3}E_{x}^{k}+\frac{\sigma^{2}(1+2\sigma^{2})}{1-\sigma^{2}}\alpha^{2}E_{s}^{k},$

where we have used Lemma 8 in the inequality.

Second, we bound the consensus error $E_{\tilde{x}}^{k}$ . We have

$\displaystyle E_{\tilde{x}}^{k+1}={}$	$\displaystyle\mathbb{E}[\\|\tilde{x}_{k+1}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}]$
$\displaystyle={}$	$\displaystyle p\mathbb{E}[\\|x^{k+1}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}]$
	$\displaystyle+(1-p)\mathbb{E}[\\|\tilde{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}]$
$\displaystyle={}$	$\displaystyle pE_{x}^{k+1}+(1-p)\mathbb{E}[\\|\tilde{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}].$	(30)

For the second term of (30), we have

	$\displaystyle\mathbb{E}[\\|\tilde{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}]$	(31)
$\displaystyle={}$	$\displaystyle\mathbb{E}[\\|\tilde{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k}+\mathbf{1}_{N}\otimes\bar{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}_{k+1}\\|^{2}]$
$\displaystyle={}$	$\displaystyle\mathbb{E}[\\|\tilde{x}^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k}+\alpha\mathbf{1}_{N}\otimes\bar{g}^{k}\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle(1+\beta)E_{\tilde{x}}^{k}+\left(1+\mfrac{1}{\beta}\right)N\alpha^{2}\mathbb{E}[\\|\bar{g}^{k}\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle\!\!\!\left[1\!+\!\beta\!+\!\left(\!1\!+\!\mfrac{1}{\beta}\!\right)\!48\alpha^{2}L^{2}d\right]\!E_{\tilde{x}}^{k}\!+\!\left(\!1\!+\!\mfrac{1}{\beta}\!\right)\!4\alpha^{2}L^{2}(1\!+\!12d)E_{x}^{k}$
	$\displaystyle\!\!\!+\!\left(\!1\!+\!\mfrac{1}{\beta}\!\right)\!2N\alpha^{2}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]\!+\!\left(\!1\!+\!\mfrac{1}{\beta}\!\right)\!14N\alpha^{2}(\tilde{u}^{k})^{2}L^{2}d^{2}\!,$

where $\beta>0$ is an auxiliary parameter that will be determined later.

Now, we can plug (29) and (31) into (30) to obtain

$\displaystyle E_{\tilde{x}}^{k+1}\leq{}$	$\displaystyle(1-p)\left[1+\beta+\left(1+\mfrac{1}{\beta}\right)48\alpha^{2}L^{2}d\right]E_{\tilde{x}}^{k}$	(32)
	$\displaystyle+\!\left[(1\!-\!p)\!\left(1\!+\!\mfrac{1}{\beta}\right)4\alpha^{2}L^{2}(1\!+\!12d)+p\frac{1\!+\!2\sigma^{2}}{3}\right]\!E_{x}^{k}$
	$\displaystyle+\frac{\sigma^{2}(1+2\sigma^{2})}{1-\sigma^{2}}\alpha^{2}pE_{s}^{k}$
	$\displaystyle+(1-p)\left(1+\mfrac{1}{\beta}\right)2N\alpha^{2}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]$
	$\displaystyle+(1-p)\left(1+\mfrac{1}{\beta}\right)14N\alpha^{2}(\tilde{u}^{k})^{2}L^{2}d^{2}.$

Letting $\beta=4\sqrt{3}\alpha L\sqrt{d}$ and using $\alpha L\leq\frac{1}{6\sqrt{d}}$ leads to

	$\displaystyle E_{\tilde{x}}^{k+1}\leq{}$	$\displaystyle(1-p)(1+48\alpha^{2}L^{2}d+8\sqrt{3}\alpha L\sqrt{d})E_{\tilde{x}}^{k}$
		$\displaystyle+\left(17\alpha L\sqrt{d}+\frac{1\!+\!2\sigma^{2}}{3}\right)E_{x}^{k}+\frac{\sigma^{2}(1\!+\!2\sigma^{2})}{1\!-\!\sigma^{2}}\alpha^{2}E_{s}^{k}$
		$\displaystyle+\frac{N\alpha}{L\sqrt{d}}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+5N\alpha Ld^{\frac{3}{2}}(\tilde{u}^{k})^{2}.$

Using $\alpha L\leq\frac{\sqrt{3}}{12d}\left(-\sqrt{d}+\sqrt{\frac{d-1+\sigma^{2}}{1-p}}\right)$ and under the condition that $p\in(\frac{1-\sigma^{2}}{d},1]$ , we get

(1-p)(1+48\alpha^{2}L^{2}d+8\sqrt{3}\alpha L\sqrt{d})\leq 1-\frac{1-\sigma^{2}}{d}.

Consequently, we derive that

$\displaystyle E_{\tilde{x}}^{k+1}\leq{}$	$\displaystyle\left(1-\frac{1\!-\!\sigma^{2}}{d}\right)E_{\tilde{x}}^{k}+\left(17\alpha L\sqrt{d}+\frac{1\!+\!2\sigma^{2}}{3}\right)E_{x}^{k}$	(33)
	$\displaystyle+\frac{\sigma^{2}(1+2\sigma^{2})}{1-\sigma^{2}}\alpha^{2}E_{s}^{k}+\frac{N\alpha}{L\sqrt{d}}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]$
	$\displaystyle+5N\alpha Ld^{\frac{3}{2}}(\tilde{u}^{k})^{2}.$

Third, we bound the consensus error $E_{s}^{k}$ . Note that

$\displaystyle E_{s}^{k+1}={}$	$\displaystyle\mathbb{E}[\\|s^{k+1}-\mathbf{1}_{N}\otimes\bar{g}^{k+1}\\|^{2}]$	(34)
$\displaystyle={}$	$\displaystyle\mathbb{E}[\\|(W\otimes I_{d})(s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k}+g^{k+1}-g^{k}$
	$\displaystyle-\mathbf{1}_{N}\otimes\bar{g}^{k+1}+\mathbf{1}_{N}\otimes\bar{g}^{k})\\|^{2}].$

For the term $g^{k+1}-g^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k+1}+\mathbf{1}_{N}\otimes\bar{g}^{k}$ in (34), we have

	$\displaystyle\\|g^{k+1}-g^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k+1}+\mathbf{1}_{N}\otimes\bar{g}^{k}\\|^{2}$
$\displaystyle={}$	$\displaystyle\\|g^{k+1}-g^{k}\\|^{2}+N\\|\bar{g}^{k+1}-\bar{g}^{k}\\|^{2}$
	$\displaystyle-2\sum_{i=1}^{N}\langle g_{f_{i}}^{k+1}-g_{f_{i}}^{k},\bar{g}^{k+1}-\bar{g}^{k}\rangle$
$\displaystyle={}$	$\displaystyle\\|g^{k+1}-g^{k}\\|^{2}-N\\|\bar{g}^{k+1}-\bar{g}^{k}\\|^{2}$
$\displaystyle\leq{}$	$\displaystyle\\|g^{k+1}-g^{k}\\|^{2}.$	(35)

Combining (35) with (34), we derive that

	$\displaystyle E_{s}^{k+1}=\mathbb{E}[\\|s^{k+1}-\mathbf{1}_{N}\otimes\bar{g}^{k+1}\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle\sigma^{2}\mathbb{E}[(\\|s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k}\\|+\\|g^{k+1}-g^{k}\\|)^{2}]$
$\displaystyle\leq{}$	$\displaystyle\frac{1+2\sigma^{2}}{3}E_{s}^{k}+\frac{\sigma^{2}(1+2\sigma^{2})}{1-\sigma^{2}}\sum_{i=1}^{N}\mathbb{E}[\\|g_{f_{i}}^{k+1}-g_{f_{i}}^{k}\\|^{2}],$	(36)

where we have used inequality $(a+b)^{2}\leq(1+\varpi)a^{2}+(1+\frac{1}{\varpi})b^{2}$ for all $a,b\in\mathbb{R}$ and $\varpi>0$ in the second inequality.

Now we consider the second item in (36):

	$\displaystyle\sum_{i=1}^{N}\mathbb{E}[\\|g_{f_{i}}^{k+1}-g_{f_{i}}^{k}\\|^{2}]$	(37)
$\displaystyle={}$	$\displaystyle\mathbb{E}\bigg{[}\sum_{i=1}^{N}\\|(g_{f_{i}}^{k+1}-\nabla f_{i}(x_{i}^{k+1}))-(g_{f_{i}}^{k}-\nabla f_{i}(x_{i}^{k}))$
	$\displaystyle+(\nabla f_{i}(x_{i}^{k+1})-\nabla f_{i}(x_{i}^{k}))\\|^{2}\bigg{]}$
$\displaystyle\leq{}$	$\displaystyle 3\sum_{i=1}^{N}\!\mathbb{E}[\\|g_{f_{i}}^{k+1}\!-\!\nabla f_{i}(x_{i}^{k+1})\\|^{2}]\!+\!3\sum_{i=1}^{N}\!\mathbb{E}[\\|g_{f_{i}}^{k}\!-\!\nabla f_{i}(x_{i}^{k})\\|^{2}]$
	$\displaystyle+3\sum_{i=1}^{N}\mathbb{E}[\\|\nabla f_{i}(x_{i}^{k+1})-\nabla f_{i}(x_{i}^{k})\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle 3L^{2}\mathbb{E}[\\|x^{k+1}-x^{k}\\|^{2}]+36dL^{2}E_{x}^{k}+36dL^{2}E_{\tilde{x}}^{k}$
	$\displaystyle+36dL^{2}E_{x}^{k+1}+36dL^{2}E_{\tilde{x}}^{k+1}+21N(\tilde{u}^{k})^{2}L^{2}d^{2},$

where we have used Lemma 1 in the second inequality.

According to (13), we have

	$\displaystyle x^{k+1}-x^{k}={}$	$\displaystyle(W\otimes I_{d}-I_{Nd})x^{k}-\alpha(W\otimes I_{d})s^{k}$
	$\displaystyle={}$	$\displaystyle(W\otimes I_{d}-I_{Nd})(x^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k})$
		$\displaystyle-\alpha(W\otimes I_{d})(s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k})-\alpha\mathbf{1}_{N}\otimes\bar{g}^{k}.$

Therefore, we can get

	$\displaystyle\mathbb{E}[\\|x^{k+1}-x^{k}\\|^{2}]$
$\displaystyle={}$	$\displaystyle\mathbb{E}[\\|(W\otimes I_{d}-I_{Nd})(x^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k})$
	$\displaystyle-\alpha(W\otimes I_{d})(s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k})-\alpha\mathbf{1}_{N}\otimes\bar{g}^{k}\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle 3\mathbb{E}[\\|\alpha(W\otimes I_{d})(s^{k}-\mathbf{1}_{N}\otimes\bar{g}^{k})\\|^{2}]+3\mathbb{E}[\\|\alpha\mathbf{1}_{N}\otimes\bar{g}^{k}\\|^{2}]$
	$\displaystyle+3\mathbb{E}[\\|(W\otimes I_{d}-I_{Nd})(x^{k}-\mathbf{1}_{N}\otimes\bar{x}^{k})\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle 3\alpha^{2}\sigma^{2}E_{s}^{k}+3\alpha^{2}N\mathbb{E}[\\|\bar{g}^{k}\\|^{2}]+12E_{x}^{k}$
$\displaystyle\leq{}$	$\displaystyle(12\!+\!12\alpha^{2}L^{2}(1\!+\!12d))E_{x}^{k}+144\alpha^{2}L^{2}dE_{\tilde{x}}^{k}+3\alpha^{2}\sigma^{2}E_{s}^{k}$
	$\displaystyle+6\alpha^{2}N\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+42N\alpha^{2}(\tilde{u}^{k})^{2}L^{2}d^{2}$
$\displaystyle\leq{}$	$\displaystyle 17E_{x}^{k}+4E_{\tilde{x}}^{k}+3\alpha^{2}\sigma^{2}E_{s}^{k}$
	$\displaystyle+6\alpha^{2}N\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+2N(\tilde{u}^{k})^{2}d,$	(38)

where we used $\|W\otimes I_{d}-I_{Nd}\|\leq 2$ in the second inequality, used (28) in the third inequality, and used $\alpha L\leq\frac{1}{6\sqrt{d}}$ in the last inequality.

Combine (29), (33), (37), and (38) with (36), and we get

$\displaystyle E_{s}^{k+1}\leq{}$	$\displaystyle\bigg{[}\frac{1+2\sigma^{2}}{3}+36dL^{2}\alpha^{2}(1+p)\frac{\sigma^{4}(1+2\sigma^{2})^{2}}{(1-\sigma^{2})^{2}}$	(39)
	$\displaystyle+9\alpha^{2}L^{2}\frac{\sigma^{4}(1+2\sigma^{2})}{1-\sigma^{2}}\bigg{]}E_{s}^{k}+\frac{9L^{2}}{1-\sigma^{2}}(70d\!+\!17)E_{x}^{k}$
	$\displaystyle+\frac{252dL^{2}}{1-\sigma^{2}}E_{\tilde{x}}^{k}+\frac{20N}{1-\sigma^{2}}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]$
	$\displaystyle+\frac{171NL^{2}d^{2}}{1-\sigma^{2}}(\tilde{u}^{k})^{2}.$

Using $\alpha L\leq\frac{(1-\sigma^{2})^{2}}{12\sigma^{2}(1+2\sigma^{2})\sqrt{d(1+p)}}$ , we have

		$\displaystyle 36dL^{2}\alpha^{2}(1+p)\frac{\sigma^{4}(1+2\sigma^{2})^{2}}{(1-\sigma^{2})^{2}}+9\alpha^{2}L^{2}\frac{\sigma^{4}(1+2\sigma^{2})}{1-\sigma^{2}}$
	$\displaystyle<{}$	$\displaystyle 45dL^{2}\alpha^{2}(1+p)\frac{\sigma^{4}(1+2\sigma^{2})^{2}}{(1-\sigma^{2})^{2}}\leq\frac{1-\sigma^{2}}{3}.$

We can further derive from inequality (39) that

	$\displaystyle E_{s}^{k+1}\leq{}$	$\displaystyle\frac{2+\sigma^{2}}{3}E_{s}^{k}+\frac{783dL^{2}}{1-\sigma^{2}}E_{x}^{k}+\frac{252dL^{2}}{1-\sigma^{2}}E_{\tilde{x}}^{k}$		(40)
		$\displaystyle+\frac{20N}{1-\sigma^{2}}\mathbb{E}[\\|\nabla f(\bar{x}^{k})\\|^{2}]+\frac{171NL^{2}d^{2}}{1-\sigma^{2}}(\tilde{u}^{k})^{2}.$		(40)

Combining (29), (33) and (40), we complete the proof.

Appendix D. Proof of Lemma 4

Using Lemma 9, we derive a range of $\alpha$ and a positive vector $\pi=[\pi_{1},\pi_{2},\pi_{3}]$ such that the following inequality holds:

A\pi\leq(1-\frac{1-\sigma^{2}}{2d})\pi,

which can be equivalently written as


	$\displaystyle c_{1}\alpha L\frac{\pi_{3}}{\pi_{1}}\leq\left(\frac{2}{3}-\frac{1}{2d}\right)(1-\sigma^{2}),$		(41a)
	$\displaystyle\!\!\left(\!c_{2}\alpha L\!+\!\frac{1\!+\!2\sigma^{2}}{3}\!\right)\!\frac{\pi_{1}}{\pi_{2}}+c_{1}\alpha L\frac{\pi_{3}}{\pi_{2}}\leq\left(\frac{1}{d}\!-\!\frac{1}{2d}\right)\!(1\!-\!\sigma^{2}),$		(41b)
	$\displaystyle c_{1}\alpha L\frac{\pi_{1}}{\pi_{3}}+\frac{1}{3}c_{1}\alpha L\frac{\pi_{2}}{\pi_{3}}\leq\left(\frac{1}{3}-\frac{1}{2d}\right)(1-\sigma^{2}),$		(41c)

where $c_{1}=\frac{9\sqrt{29}\sqrt{d}}{1-\sigma^{2}}$ and $c_{2}=17\sqrt{d}$ . Under the condition $d\geq 3$ , to show the inequalities (41), it is sufficient to check the following inequalities hold:

	$\displaystyle\frac{1}{3}c_{1}\alpha L\frac{\pi_{3}}{\pi_{1}}\leq\frac{1-\sigma^{2}}{2d},$
	$\displaystyle(c_{2}\alpha L+\frac{1+2\sigma^{2}}{3})\frac{\pi_{1}}{\pi_{2}}+c_{1}\alpha L\frac{\pi_{3}}{\pi_{2}}\leq\frac{1-\sigma^{2}}{2d},$
	$\displaystyle c_{1}\alpha L\frac{\pi_{1}}{\pi_{3}}+\frac{1}{3}c_{1}\alpha L\frac{\pi_{2}}{\pi_{3}}\leq\frac{1-\sigma^{2}}{2d},$

and by further letting $\frac{\pi_{1}}{\pi_{2}}=\frac{1-\sigma^{2}}{3d}$ and $\frac{\pi_{2}}{\pi_{3}}=3$ , it is sufficient to show the following inequalities:

		$\displaystyle\frac{1}{3}c_{1}\alpha L\frac{d}{1-\sigma^{2}}\leq\frac{1-\sigma^{2}}{2d},$
		$\displaystyle c_{2}\alpha L\frac{1-\sigma^{2}}{d}+c_{1}\alpha L\leq\frac{1-\sigma^{2}}{2d},$
		$\displaystyle c_{1}\alpha L\frac{1-\sigma^{2}}{d}+c_{1}\alpha L\leq\frac{1-\sigma^{2}}{2d}.$

Using inequality $c_{2}<c_{1}$ , we have the following inequality:

	$\displaystyle\max\Bigg{\{}$	$\displaystyle c_{2}\alpha L\frac{1\!-\!\sigma^{2}}{d}+c_{1}\alpha L,c_{1}\alpha L\frac{1\!-\!\sigma^{2}}{d}+c_{1}\alpha L,$
		$\displaystyle\frac{1}{3}c_{1}\alpha L\frac{d}{1\!-\!\sigma^{2}}\Bigg{\}}<\alpha L\left(c_{1}\frac{1\!-\!\sigma^{2}}{d}+c_{1}\frac{d}{3(1\!-\!\sigma^{2})}\right).$

Consequently, it suffices to ensure that

\alpha L\left(c_{1}\frac{1-\sigma^{2}}{d}+c_{1}\frac{d}{3(1-\sigma^{2})}\right)\leq\frac{1-\sigma^{2}}{2d},

which can now be justified by the condition $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{12\sqrt{29}d^{\frac{5}{2}}}$ . It is straightforward to verify that $\pi=[\frac{1-\sigma^{2}}{d},3,1]^{T}$ satisfies $\frac{\pi_{1}}{\pi_{2}}=\frac{1-\sigma^{2}}{3d}$ and $\frac{\pi_{2}}{\pi_{3}}=3$ , which completes the proof.

Appendix E. Proof of Lemma 5

Based on Lemmas 3 and 4 , we are ready to prove Lemma 5.

It is straightforward to verify that

\alpha L\leq\min\!\left\{\frac{1}{6\sqrt{d}},\frac{\sqrt{3}}{12d}\!\left(\!-\sqrt{d}\!+\!\sqrt{\mfrac{d\!-\!1\!+\!\sigma^{2}}{1-p}}\right)\!,\frac{(1-\sigma^{2})^{3}}{12\sqrt{29}d^{\frac{5}{2}}}\right\}

is a sufficient (and unnecessary) condition for

	$\displaystyle\alpha L\leq\min\Bigg{\{}$	$\displaystyle\frac{1}{6\sqrt{d}},\frac{(1-\sigma^{2})^{2}}{12\sigma^{2}(1+2\sigma^{2})\sqrt{d(1+p)}},$
		$\displaystyle\frac{\sqrt{3}}{12d}\left(-\sqrt{d}+\sqrt{\mfrac{d-1+\sigma^{2}}{1-p}}\right)\Bigg{\}}.$

Thus, Lemma 3 holds.

Note that the matrix norm $\interleave\cdot\interleave_{\infty}^{\pi}$ is induced by the vector norm $\|\cdot\|_{\infty}^{\pi}$ . From the inequality (16), we derive that

	$\displaystyle\\|v^{k+1}\\|_{\infty}^{\pi}$	$\displaystyle\leq\\|Av^{k}\\|_{\infty}^{\pi}+\\|B_{k}\\|_{\infty}^{\pi}$		(42)
		$\displaystyle\leq\interleave A\interleave_{\infty}^{\pi}\\|v^{k}\\|_{\infty}^{\pi}+\\|B_{k}\\|_{\infty}^{\pi},$		(42)

and by induction on (42), we get

\displaystyle\|v^{k}\|_{\infty}^{\pi}\leq(\interleave A\interleave_{\infty}^{\pi})^{k}\|v_{0}\|_{\infty}^{\pi}+\sum_{m=0}^{k-1}(\interleave A\interleave_{\infty}^{\pi})^{m}\|B_{k-m-1}\|_{\infty}^{\pi}.

Using the definition of the vector norm $\|x\|_{\infty}^{w}=\max_{i}|x_{i}|/w_{i}$ , we derive that

	$\displaystyle\max\left\{\frac{d}{1-\sigma^{2}}E_{x}^{k},\frac{1}{3}E_{\tilde{x}}^{k},\frac{\alpha}{3\sqrt{29}L\sqrt{d}}E_{s}^{k}\right\}$	(43)
$\displaystyle\leq{}$	$\displaystyle\sum_{m=0}^{k-1}(\interleave A\interleave_{\infty}^{\pi})^{m}\bigg{(}\frac{\sqrt{29}N\alpha}{3(1-\sigma^{2})L\sqrt{d}}\mathbb{E}[\\|\nabla f(\bar{x}^{k-m-1})\\|^{2}]$
	$\displaystyle+\frac{2\sqrt{29}}{1-\sigma^{2}}Nd^{\frac{3}{2}}\alpha L(\tilde{u}^{k-m-1})^{2}\bigg{)}+(\interleave A\interleave_{\infty}^{\pi})^{k}\\|v_{0}\\|_{\infty}^{\pi}.$

Note that for any nonnegative sequence $(a_{m})_{m\in\mathbb{N}}$ and $\rho\in(0,1)$ , we have

$\displaystyle\sum_{\tau=1}^{k}\sum_{m=0}^{\tau-1}\rho^{m}a_{\tau-m-1}$	$\displaystyle=\sum_{\tau=1}^{k}\sum_{m=0}^{\tau-1}\rho^{\tau-m-1}a_{m}$	(44)
	$\displaystyle=\sum_{m=0}^{k-1}\rho^{-m-1}a_{m}\sum_{\tau=m+1}^{k}\rho^{\tau}$
	$\displaystyle\leq\frac{1}{1-\rho}\sum_{m=0}^{k-1}a_{m}.$

Consequently, we obtain

		$\displaystyle\max\left\{\sum_{\tau=0}^{k}\frac{d}{1-\sigma^{2}}E_{x}^{\tau},\sum_{\tau=0}^{k}\frac{1}{3}E_{\tilde{x}}^{\tau},\sum_{\tau=0}^{k}\frac{\alpha}{3\sqrt{29}L\sqrt{d}}E_{s}^{\tau}\right\}$
	$\displaystyle\leq{}$	$\displaystyle\frac{4d}{1-\sigma^{2}}R_{0}+\frac{2\sqrt{29}\sqrt{d}N\alpha}{3(1-\sigma^{2})^{2}L}\sum_{m=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{m})\\|^{2}]$
		$\displaystyle+\frac{4\sqrt{29}Nd^{\frac{5}{2}}\alpha L}{(1-\sigma^{2})^{2}}\sum_{m=0}^{k-1}(\tilde{u}^{m})^{2},$

where $R_{0}=\|v_{0}\|_{\infty}^{\pi}=\frac{d}{1-\sigma^{2}}E_{x}^{0}$ and we have used $\interleave A\interleave_{\infty}^{\pi}\leq 1-\frac{1-\sigma^{2}}{2d}$ in the inequality.

For the bound on $\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}$ , we first derive that

	$\displaystyle(\tilde{u}^{k})^{2}$	$\displaystyle=\max_{i}\{\mathbb{E}[(\tilde{u}_{i}^{k})^{2}]\}$
		$\displaystyle=(1-p)\max_{i}\{\mathbb{E}[(\tilde{u}_{i}^{k-1})^{2}]\}+p\max_{i}\{(u_{i}^{k})^{2}\}$
		$\displaystyle=(1-p)^{k}(u_{i}^{0})^{2}+p\sum_{m=0}^{k-1}(1-p)^{m}(u_{i}^{k-m})^{2},$

where we have used the condition that all the $u_{i}^{k}$ follow from the same iterative equation in the third equality. Then, using inequality (44), we further get

	$\displaystyle\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}={}$	$\displaystyle\sum_{\tau=0}^{\infty}(1\!-\!p)^{\tau}(u_{i}^{0})^{2}+\sum_{\tau=1}^{\infty}\sum_{m=0}^{\tau-1}p(1\!-\!p)^{m}(u_{i}^{\tau-m})^{2}$
	$\displaystyle\leq{}$	$\displaystyle\frac{1}{p}(u_{i}^{0})^{2}+\sum_{m=1}^{\infty}(u_{i}^{m})^{2}.$

Since $(u_{i}^{k})^{2}$ is summable, we derive that $(\tilde{u}^{k})^{2}$ is also summable, which completes the proof.

Appendix F. Proof of Corollary 1

Let $p=\frac{1}{d}$ and $\beta=4\sqrt{3}\alpha L\sqrt{d}$ . Using $\alpha L\leq\frac{\sqrt{3}}{12\sqrt{d}}\left(-1+\sqrt{1+\frac{\sigma^{2}}{d-1}}\right)$ , we derive from (32) that

$\displaystyle E_{\tilde{x}}^{k+1}\leq{}$	$\displaystyle(1-\frac{1-\sigma^{2}}{d})E_{\tilde{x}}^{k}+(17\alpha L\sqrt{d}+\frac{1+2\sigma^{2}}{3d})E_{x}^{k}$	(45)
	$\displaystyle+\frac{\sigma^{2}(1+2\sigma^{2})}{(1-\sigma^{2})d}\alpha^{2}E_{s}^{k}+\frac{N\alpha}{L\sqrt{d}}\\|\nabla f(\bar{x}^{k})\\|^{2}$
	$\displaystyle+5N\alpha Ld^{\frac{3}{2}}(\tilde{u}^{k})^{2}.$

Combining (29), (40) and (45), we get

\displaystyle v^{k+1}\leq A^{\prime}v^{k}+B_{k},

where

A^{\prime}=\begin{bmatrix}\frac{1+2\sigma^{2}}{3}&0&\frac{9\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L\\ 17\sqrt{d}\alpha L+\frac{1+2\sigma^{2}}{3d}&1-\frac{1-\sigma^{2}}{d}&\frac{9\sqrt{29}}{(1-\sigma^{2})\sqrt{d}}\alpha L\\ \frac{9\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L&\frac{3\sqrt{29}\sqrt{d}}{1-\sigma^{2}}\alpha L&\frac{2+\sigma^{2}}{3}\end{bmatrix}.

We then demonstrate that the vector $\pi^{\prime}=[1-\sigma^{2},3,1]$ and the step-size condition $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{45\sqrt{29}d^{\frac{3}{2}}}$ together ensure the inequality

A^{\prime}\pi^{\prime}\leq(1-\frac{1-\sigma^{2}}{2d})\pi^{\prime}.

Similar to the proof of Lemma 4 , it suffices to verify

		$\displaystyle\frac{1}{3(1-\sigma^{2})}c_{1}\alpha L\leq\frac{1-\sigma^{2}}{2d},$		(46)
		$\displaystyle(1-\sigma^{2})c_{2}\alpha L+\frac{c_{1}\alpha L}{d}\leq\frac{1-\sigma^{2}}{2d},$
		$\displaystyle(1-\sigma^{2})c_{1}\alpha L+c_{1}\alpha L\leq\frac{1-\sigma^{2}}{2d}.$

It is straightforward to verify that $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{45\sqrt{29}d^{\frac{3}{2}}}$ is a sufficient condition of inequality (46). Using Lemma 9, we then obtain

\rho(A^{\prime})\leq\interleave A^{\prime}\interleave_{\infty}^{\pi^{\prime}}\leq 1-\frac{1-\sigma^{2}}{2d}.

(47)

Next, similar to Lemma 5, we can show that when

\alpha L\leq\min\!\left\{\frac{1}{6\sqrt{d}},\frac{\sqrt{3}}{12\sqrt{d}}\!\left(-1\!+\!\sqrt{1\!+\!\mfrac{\sigma^{2}}{d\!-\!1}}\right),\frac{(1-\sigma^{2})^{3}}{45\sqrt{29}d^{\frac{3}{2}}}\right\},

the following inequality holds:

	$\displaystyle\max\left\{\frac{1}{1-\sigma^{2}}\sum_{\tau=0}^{k}E_{x}^{\tau},\frac{1}{3}\sum_{\tau=0}^{k}E_{\tilde{x}}^{\tau},\frac{\alpha}{3\sqrt{29}L\sqrt{d}}\sum_{\tau=0}^{k}E_{s}^{\tau}\right\}$	(48)
$\displaystyle\leq{}$	$\displaystyle\frac{4d}{1-\sigma^{2}}\tilde{R}_{0}+\frac{2\sqrt{29}\sqrt{d}N\alpha}{3(1-\sigma^{2})^{2}L}\sum_{m=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{m})\\|^{2}]$
	$\displaystyle+\frac{4\sqrt{29}Nd^{\frac{5}{2}}\alpha L}{(1-\sigma^{2})^{2}}\sum_{m=0}^{k-1}(\tilde{u}^{m})^{2},$

where $\tilde{R}_{0}=\frac{1}{1-\sigma^{2}}E_{x}^{0}$ .

Now, we combine (48) with (15) to derive

$\displaystyle 0\leq{}$	$\displaystyle\delta^{0}-\frac{\alpha}{3}\sum_{\tau=0}^{k}\mathbb{E}[\\|\nabla f(\bar{x}^{\tau})\\|^{2}]+5\alpha L^{2}d^{2}\sum_{\tau=0}^{k}(\tilde{u}^{\tau})^{2}$	(49)
	$\displaystyle+\frac{66\alpha dL^{2}}{N}(\frac{4d}{1-\sigma^{2}}\tilde{R}_{0}+\frac{4\sqrt{29}Nd^{\frac{5}{2}}\alpha L}{(1-\sigma^{2})^{2}}\sum_{m=0}^{k-1}(\tilde{u}^{m})^{2}$
	$\displaystyle+\frac{2\sqrt{29}\sqrt{d}N\alpha}{3(1-\sigma^{2})^{2}L}\sum_{m=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{m})\\|^{2}])$
$\displaystyle\leq{}$	$\displaystyle\delta^{0}+\frac{\sqrt{d}L}{\sqrt{29}N}\tilde{R}_{0}-\frac{\alpha}{6}\sum_{\tau=0}^{k}\mathbb{E}[\\|\nabla f(\bar{x}^{\tau})\\|^{2}]$
	$\displaystyle+6\alpha L^{2}d^{2}\sum_{\tau=0}^{k}(\tilde{u}^{\tau})^{2},$

where we have used $\frac{18\alpha dL^{2}}{N}\cdot(1-\sigma^{2})+\frac{16\alpha dL^{2}}{N}\cdot 3<\frac{66\alpha dL^{2}}{N}$ in the first inequality and $\alpha L\leq\frac{(1-\sigma^{2})^{3}}{264\sqrt{29}d^{\frac{3}{2}}}\leq\frac{1-\sigma^{2}}{264\sqrt{29}d^{\frac{3}{2}}}$ in the last inequality. We then derive from (49) that

		$\displaystyle\frac{1}{k}\sum_{\tau=0}^{k-1}\mathbb{E}[\\|\nabla f(\bar{x}^{\tau})\\|^{2}]$
	$\displaystyle\leq{}$	$\displaystyle\frac{1}{k}\left[\frac{6}{\alpha}\delta^{0}+\frac{6\sqrt{d}L}{\sqrt{29}N\alpha}\tilde{R}_{0}+36L^{2}d^{2}\sum_{\tau=0}^{\infty}(\tilde{u}^{\tau})^{2}\right].$

The proof of Corollary 1 can now be concluded.

References

[1] Shu Liang, Le Yi Wang, and George Yin. Distributed smooth convex optimization with coupled constraints. IEEE Transactions on Automatic Control, 65(1):347–353, 2020.
[2] Xuan Zhang, Antonis Papachristodoulou, and Na Li. Distributed control for reaching optimal steady state in network systems: An optimization approach. IEEE Transactions on Automatic Control, 63(3):864–871, 2018.
[3] Jie Liu, Daniel W. C. Ho, and Lulu Li. A generic algorithm framework for distributed optimization over the time-varying network with communication delays. IEEE Transactions on Automatic Control, 69(1):371–378, 2024.
[4] Angelia Nedić and Asuman Ozdaglar. Distributed subgradient methods for multi-agent optimization. IEEE Transactions on Automatic Control, 54(1):48–61, 2009.
[5] Wei Shi, Qing Ling, Gang Wu, and Wotao Yin. EXTRA: An exact first-order algorithm for decentralized consensus optimization. SIAM Journal on Optimization, 25(2):944–966, 2015.
[6] Angelia Nedić. Distributed gradient methods for convex machine learning problems in networks: Distributed optimization. IEEE Signal Processing Magazine, 37(3):92–101, 2020.
[7] Guido Carnevale, Nicola Mimmo, and Giuseppe Notarstefano. Nonconvex distributed feedback optimization for aggregative cooperative robotics. Automatica, 167:111767, 2024.
[8] Xiangru Lian, Ce Zhang, Huan Zhang, Cho-Jui Hsieh, Wei Zhang, and Ji Liu. Can decentralized algorithms outperform centralized algorithms? A case study for decentralized parallel stochastic gradient descent. In Advances in Neural Information Processing Systems, volume 30, 2017.
[9] Jinshan Zeng and Wotao Yin. On nonconvex decentralized gradient descent. IEEE Transactions on signal processing, 66(11):2834–2848, 2018.
[10] Gesualdo Scutari and Ying Sun. Distributed nonconvex constrained optimization over time-varying digraphs. Mathematical Programming, 176:497–544, 2019.
[11] Zhongguo Li, Zhen Dong, Zhongchao Liang, and Zhengtao Ding. Surrogate-based distributed optimisation for expensive black-box functions. Automatica, 125:109407, 2021.
[12] Sébastien Bubeck and Nicolo Cesa-Bianchi. Regret analysis of stochastic and nonstochastic multi-armed bandit problems. Foundations and Trends® in Machine Learning, 5(1):1–122, 2012.
[13] Sadhika Malladi, Tianyu Gao, Eshaan Nichani, Alex Damian, Jason D. Lee, Danqi Chen, and Sanjeev Arora. Fine-tuning language models with just forward passes. In Advances in Neural Information Processing Systems, volume 36, pages 53038–53075, 2023.
[14] Anit Kumar Sahu, Dusan Jakovetic, Dragana Bajovic, and Soummya Kar. Communication-efficient distributed strongly convex stochastic optimization: Non-asymptotic rates. arXiv preprint arXiv:1809.02920, 2018.
[15] Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization of convex functions. Foundations of Computational Mathematics, 17(2):527–566, 2017.
[16] Jack Kiefer and Jacob Wolfowitz. Stochastic estimation of the maximum of a regression function. The Annals of Mathematical Statistics, pages 462–466, 1952.
[17] Yujie Tang, Junshan Zhang, and Na Li. Distributed zero-order algorithms for nonconvex multiagent optimization. IEEE Transactions on Control of Network Systems, 8(1):269–281, 2021.
[18] Yongyi Guo, Dominic Coey, Mikael Konutgan, Wenting Li, Chris Schoener, and Matt Goldman. Machine learning for variance reduction in online experiments. In Advances in Neural Information Processing Systems, volume 34, pages 8637–8648, 2021.
[19] Chong Wang, Xi Chen, Alexander J. Smola, and Eric P. Xing. Variance reduction for stochastic gradient optimization. In Advances in Neural Information Processing Systems, volume 26, 2013.
[20] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems, volume 26, 2013.
[21] Sijia Liu, Bhavya Kailkhura, Pin-Yu Chen, Paishun Ting, Shiyu Chang, and Lisa Amini. Zeroth-order stochastic variance reduction for nonconvex optimization. In Advances in Neural Information Processing Systems, volume 31, 2018.
[22] Ran Xin, Usman A Khan, and Soummya Kar. Variance-reduced decentralized stochastic optimization with accelerated convergence. IEEE Transactions on Signal Processing, 68:6255–6271, 2020.
[23] Davood Hajinezhad, Mingyi Hong, and Alfredo Garcia. ZONE: Zeroth-order nonconvex multiagent optimization over networks. IEEE Transactions on Automatic Control, 64(10):3995–4010, 2019.
[24] Arkadij Semenovič Nemirovskij and David Borisovich Yudin. Problem Complexity and Method Efficiency in Optimization. Wiley-Interscience, 1983.
[25] Abraham D Flaxman, Adam Tauman Kalai, and H Brendan McMahan. Online convex optimization in the bandit setting: gradient descent without a gradient. In Proceedings of the Sixteenth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 385–394, 2005.
[26] Guannan Qu and Na Li. Harnessing smoothness to accelerate distributed optimization. IEEE Transactions on Control of Network Systems, 5(3):1245–1260, 2017.
[27] Angelia Nedić, Alex Olshevsky, and Wei Shi. Achieving geometric convergence for distributed optimization over time-varying graphs. SIAM Journal on Optimization, 27(4):2597–2633, 2017.
[28] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 2nd edition, 2012.

		$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|g-\nabla h(x)\\|^{2}\right]$
	$\displaystyle={}$	$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\Big{[}\Big{\\|}\left(G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)\right)$
		$\displaystyle\!-\!\Big{(}G_{h}^{(2d)}\!(x,u)\!-\!G_{h}^{(2d)}\!(\tilde{x},\tilde{u})\Big{)}\!-\!\left(\nabla h(x)\!-\!G_{h}^{(2d)}\!(x,u)\!\right)\!\!\Big{\\|}^{2}\Big{]}$
	$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}_{l\sim\mathcal{U}[d]}\!\Big{[}\Big{\\|}\left(G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)\right)$
		$\displaystyle-\left(G_{h}^{(2d)}(x,u)-G_{h}^{(2d)}(\tilde{x},\tilde{u})\right)\Big{\\|}^{2}\Big{]}$
		$\displaystyle+2\mathbb{E}[\\|\nabla h(x)-G_{h}^{(2d)}(x,u)\\|^{2}].$

		$\displaystyle\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|g-\nabla h(x)\\|^{2}\right]$
	$\displaystyle\leq{}$	$\displaystyle 2\mathbb{E}_{l\sim\mathcal{U}[d]}\!\left[\\|G_{h}^{(c)}(x,u,l)-G_{h}^{(c)}(\tilde{x},\tilde{u},l)\\|^{2}\right]+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle={}$	$\displaystyle 2d\\|G_{h}^{(2d)}(x,u)-G_{h}^{(2d)}(\tilde{x},\tilde{u})\\|^{2}+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle={}$	$\displaystyle 2d\\|(G_{h}^{(2d)}(x,u)-\nabla h(x))-(G_{h}^{(2d)}(\tilde{x},\tilde{u})-\nabla h(\tilde{x}))$
		$\displaystyle+(\nabla h(x)-\nabla h(\tilde{x}))\\|^{2}+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle\leq{}$	$\displaystyle 6dL^{2}\\|x-\tilde{x}\\|^{2}+\frac{3}{2}u^{2}L^{2}d^{2}+\frac{3}{2}\tilde{u}^{2}L^{2}d^{2}+\frac{1}{2}u^{2}L^{2}d$
	$\displaystyle\leq{}$	$\displaystyle 12dL^{2}\\|x-y\\|^{2}+12dL^{2}\\|\tilde{x}-y\\|^{2}+\frac{7}{2}\tilde{u}^{2}L^{2}d^{2},$

	$\displaystyle\sum_{i=1}^{N}\mathbb{E}[\\|g_{f_{i}}^{k+1}-g_{f_{i}}^{k}\\|^{2}]$	(37)
$\displaystyle={}$	$\displaystyle\mathbb{E}\bigg{[}\sum_{i=1}^{N}\\|(g_{f_{i}}^{k+1}-\nabla f_{i}(x_{i}^{k+1}))-(g_{f_{i}}^{k}-\nabla f_{i}(x_{i}^{k}))$
	$\displaystyle+(\nabla f_{i}(x_{i}^{k+1})-\nabla f_{i}(x_{i}^{k}))\\|^{2}\bigg{]}$
$\displaystyle\leq{}$	$\displaystyle 3\sum_{i=1}^{N}\!\mathbb{E}[\\|g_{f_{i}}^{k+1}\!-\!\nabla f_{i}(x_{i}^{k+1})\\|^{2}]\!+\!3\sum_{i=1}^{N}\!\mathbb{E}[\\|g_{f_{i}}^{k}\!-\!\nabla f_{i}(x_{i}^{k})\\|^{2}]$
	$\displaystyle+3\sum_{i=1}^{N}\mathbb{E}[\\|\nabla f_{i}(x_{i}^{k+1})-\nabla f_{i}(x_{i}^{k})\\|^{2}]$
$\displaystyle\leq{}$	$\displaystyle 3L^{2}\mathbb{E}[\\|x^{k+1}-x^{k}\\|^{2}]+36dL^{2}E_{x}^{k}+36dL^{2}E_{\tilde{x}}^{k}$
	$\displaystyle+36dL^{2}E_{x}^{k+1}+36dL^{2}E_{\tilde{x}}^{k+1}+21N(\tilde{u}^{k})^{2}L^{2}d^{2},$

	$\displaystyle\\|v^{k+1}\\|_{\infty}^{\pi}$	$\displaystyle\leq\\|Av^{k}\\|_{\infty}^{\pi}+\\|B_{k}\\|_{\infty}^{\pi}$		(42)
		$\displaystyle\leq\interleave A\interleave_{\infty}^{\pi}\\|v^{k}\\|_{\infty}^{\pi}+\\|B_{k}\\|_{\infty}^{\pi},$		(42)