Year: 2021

Abstract:
Zeroth-order optimization has become increasingly important in complex optimization and machine learning when cost functions are impossible to be described in closed analytical forms. The key idea of zeroth-order optimization lies in the ability for a learner to build gradient estimates by queries sent to the cost function, and then traditional gradient descent algorithms can be executed replacing gradients by the estimates. For optimization over large-scale multi-agent systems with decentralized data and costs, zeroth-order optimization can continue to be utilized to develop scalable and distributed algorithms. In this paper, we aim at understanding the trend in performance transitioning from centralized to distributed zeroth-order algorithms in terms of convergence rates, and focus on multi-agent systems with time-varying communication networks. We establish a series of convergence rates for distributed zeroth-order subgradient algorithms under both one-point and two-point zeroth-order oracles. Apart from the additional node-to-node communication cost due to the distributed nature of algorithms, the established rates in convergence are shown to match their centralized counterpart. We also propose a multi-stage distributed zeroth-order algorithm that better utilizes the learning rates, reduces the computational complexity, and attains even faster convergence rates for compact decision set.

text: Zeroth-order optimization has become increasingly important in complex optimization and machine learning when cost functions are impossible to be described in closed analytical forms. The key idea of zeroth-order optimization lies in the ability for a learner to build gradient estimates by queries sent to the cost function, and then traditional gradient descent algorithms can be executed replacing gradients by the estimates. For optimization over large-scale multi-agent systems with decentralized data and costs, zeroth-order optimization can continue to be utilized to develop scalable and distributed algorithms. In this paper, we aim at understanding the trend in performance transitioning from centralized to distributed zeroth-order algorithms in terms of convergence rates, and focus on multi-agent systems with time-varying communication networks. We establish a series of convergence rates for distributed zeroth-order subgradient algorithms under both one-point and two-point zeroth-order oracles. Apart from the additional node-to-node communication cost due to the distributed nature of algorithms, the established rates in convergence are shown to match their centralized counterpart. We also propose a multi-stage distributed zeroth-order algorithm that better utilizes the learning rates, reduces the computational complexity, and attains even faster convergence rates for compact decision set.

heading: 1 INTRODUCTION
text: Various machine learning tasks ultimately boil down to solving optimization problems of different forms, where the cost functions are formed jointly by the data accumulated in experiences and the model used in representing the learning framework. Gradient descent algorithms have been playing a foundational role in practically solving such optimization problems. However, for learning tasks with high-dimensional data and involved learning representations, access to the gradient of the cost function may turn out not possible: the cost function supporting the learning may not have a closed analytical form; or it is simply too computationally costly to be properly differentiated. Zeroth-order optimization provides a systemic way of facilitating gradient descent without direct access to gradient information, where oracles query the cost function values and generate gradient estimates. Zeroth-order methods have shown a number of successful applications, e.g., searching for adversarial attacks in deep learning Chen et al. (2019); Liu et al. (2019) and policy search in reinforcement learning Vemula et al. (2019).
The literature has also explored the potential in extending the standard (centralized) zeroth-order optimization to distributed settings over multi-agent systems, where the data and cost functions are scattered across a network of decentralized agents. With the help of a communication network, the agents may collaboratively solve the network-level optimization task by iteratively exchanging decisions obtained from local zeroth-order descent. The rates of convergence of centralized zerothorder optimization algorithms are now well understood for several sub-classes of convex functions. We are interested in systematically investigating these convergence rates scale for the corresponding distributed algorithms, and focus on the case of time-varying communication networks.

heading: 1.1 PROBLEM DEFINITION
text: Consider a network of agents (nodes) V = {1, . . . , N}. The agents aim to collectively solve the following distributed optimization problem
minimize f(x) := N∑ i=1 fi(x)
subject to x ∈ X. (1)
Here x ∈ Rd is the decision variable, X ⊆ Rd is a convex decision space, and fi : Rd → R is a private convex objective function associated with agent i.
The communication network connecting the nodes is described by a time-varying graph G(t) = (V,E(t)), where E(t) is the set of activated links at time t. Let A(t) be a weight matrix at time t for the graph G(t): for each link (i, j) ∈ E(t), a weight [A(t)]ij > 0 is assigned, and [A(t)]ij = 0 for (i, j) /∈ E(t). We impose the following assumption on the communication network E(t) and the weight matrix A(t).
Assumption 1 (i) There exists a positive integer B such that the union graph (V,E(kB+ 1)∪ · · ·∪ E((k+1)B)) is strongly connected for all k ≥ 0; (ii) A(t) is doubly stochastic, i.e., ∑N i=1[A(t)]ij =
1 and ∑N j=1[A(t)]ij = 1; (iii) [A(t)]ii ≥ ξ for all i, and [A(t)]ij ≥ ξ if (j, i) ∈ E(t), where ξ > 0.

heading: 1.2 FUNCTION CLASSES
text: Let Fcvx denote the set of all convex functions on Rd. We define the following three classes of convex functions in Fcvx.
• The Lipschitz continuous class Flip(Lf ,X) contains the functions in Fcvx that admit a finite Lipschitz constant Lf over X, i.e.,
Flip(Lf ,X) := {g ∈ Fcvx : ∀x,x′ ∈ X, |g(x)− g(x′)| ≤ Lf‖x− x′‖}.
• The smooth class Fsmo(sf ,X) contains the functions that admit a sf -Lipschitz continuous gradient over X, i.e.,
Fsmo(sf ,X) = {g ∈ Fcvx : ∀x,x′ ∈ X, ‖∇g(x)−∇g(x′)‖ ≤ sf‖x− x′‖}.
• The strongly convex class Fsc(µf ,X) contains the functions that are µf -strongly convex, i.e.,
Fsc(µf ,X) = { g ∈ Fcvx : ∀x,x′ ∈ X, g(x) ≥ g(x′) + 〈∇g(x′),x− x′〉+
µf 2 ‖x− x′‖2
} .

heading: 1.3 CONTRIBUTIONS AND RELATED WORK
text: Contributions. We first present MAZOPA, a multi-agent zeroth-order projection averaging algorithm. In MAZOPA, the agents iteratively carry out local zeroth-order descents for their private costs to generate intermediate decisions, send these intermediate decisions to their neighbors over the graph G(t), and then update their decisions by projecting the average neighboring intermediate decisions onto X. For distributed zeroth-order oracles based on one-point or two-point estimates, a series of convergence rate results are established for the three basic function classes. Remarkably, the convergence rates for distributed algorithms are found to be matching their centralized counterpart, and sometimes even tighter rates are obtained, as summarized in Table 1. These results show that by paying the price of node-to-node communication, distributed zeroth-order optimization provides equal performance guarantees as those of centralized approaches. Next, we generalize the MAZOPA to a multi-stage setting, where the local zeroth-order descents take place for multiple steps before the projected averaging in a sequence of epochs. Such multi-stage MAZOPA is shown to be able to reduce the computational complexity, while providing improved convergences rates compared to MAZOPA when the decision set is compact.
Related Work. Recently, many types of centralized zeroth-order optimization algorithms have been studied, and their convergence rates (and the way they depend on the dimension) have been established in different settings. For unconstrained convex optimization, Nesterov & Spokoiny (2017) develops several types of two-point gradient estimators and achieves convergence rates that scale with dimension as O(d2). For constrained stochastic optimization, Duchi et al. (2015) establishes that the convergence rates are sharp up to factors at most logarithmic in the dimension. Zeroth-order optimization has a natural connection to bandit online optimization, where the latter focuses on dynamic environment where the objective functions are varying over time (see, e.g., Flaxman et al. (2005); Agarwal et al. (2010); Shamir (2013; 2017); Bubeck et al. (2017); Lattimore (2020); Hazan & Levy (2014)). In particular, the seminal work Flaxman et al. (2005) constructs a one-point gradient estimator (or one-point bandit feedback model) and achieves an O(d/T 1/4) average regret. For two-point gradient estimator, Shamir (2017) establishes the tightness of the dimension-dependent factor O( √ d) in the framework of zeroth-order stochastic mirror descent.
It is worth zooming into the literature on distributed zeroth-order/bandit online optimization. Due to the absence of a central coordinator, the algorithms developed should always rely on local computations and communications (e.g., Yuan & Ho (2015); Yi et al. (2020); Jakovetic et al. (2018); Hajinezhad et al. (2019); Wang et al. (2019); Pang & Hu (2019); Hajinezhad & Zavlanos (2018); Wan et al. (2020)). This makes the convergence analysis of the distributed zeroth-order/bandit online optimization algorithms more challenging. In Yuan & Ho (2015), the authors develop a class of distributed zeroth-order optimization algorithms that require two functional evaluations at each iteration, and establishes asymptotic convergence of the algorithm. Non-asymptotic convergence is established in Jakovetic et al. (2018); Hajinezhad et al. (2019); Wang et al. (2019); Pang & Hu (2019); Hajinezhad & Zavlanos (2018), but the dimension-dependence factors are either O(d2) or far from optimal. The work Yi et al. (2020) considers distributed online optimization with long-term constraints and establishes bounds on regret as well as constraint violations. To avoid Euclidean projection onto the constraint set, Wan et al. (2020) develops a distributed bandit online optimization algorithm based on conditional gradient descent and one-point bandit feedback, and achieves a regret scaling of O(T 3/4 √ lnT ).
2 THE MAZOPA ALGORITHM AND ITS CONVERGENCE RATES
In this section, we present the MAZOPA algorithm and establish the convergence rates for the three function classes.

heading: 2.1 DISTRIBUTED ZEROTH-ORDER ORACLES
text: Let n be a random vector in Rd drawn from some probability distribution. Then
f̂i(x; δ) := En [fi(x + δn)] (2)
is a smoothed function for fi. Here δ > 0 is a parameter setting the level of the smoothing. We introduce the following definition on distributed zeroth-order oracles (DistZOO).
Definition 1 (DistZOO) A vector g̃i(x; δ) ∈ Rd is called a distributed zeroth-order oracle at node i if the following conditions hold:
(i) E [g̃i(x; δ)] = ∇f̂i(x; δ) for all x ∈ Rd; (ii) If fi ∈ Flip(Lf ), then f̂i ∈ Flip(Lf ) as well, and there holds ∣∣f̂i(x; δ) − fi(x)∣∣ ≤ pdLfδ,
with pd being some positive constant;
(iii) If fi ∈ Fsmo(sf ), then ∣∣f̂i(x; δ)− fi(x)∣∣ ≤ 12 p̃dsfδ2 with p̃d being some positive constant.
A number of DistZOO satisfying Definition 1 can be obtained using existing gradient estimators, see, e.g., Liu et al. (2020). In the paper, we provide two representative gradient estimators that are commonly adopted in the literature. Let ui be a random vector independently generated from a unit sphere B1 in Rd. Then (e.g., Flaxman et al. (2005))
g̃OPi (x; δ) := fi(x + δtui)uid/δ (3)
is a one-point DistZOO satisfying Definition 1. Moreover,
g̃TPi (x; δ) := d
2δ
( fi(x + δu)− fi(x− δu) ) u (4)
is a two-point DistZOO satisfying Definition 1 (e.g., Shamir (2017)).
2.2 THE MAZOPA ALGORITHM
We present the following Multi-Agent Zeroth-Order Projection Averaging (MAZOPA) algorithm, which consists of two steps, a local zeroth-order optimization step and a distributed averaging step. MAZOPA, whose pseudo-code is presented in Algorithm 1, is a variation of the multi-agent subgradient averaging algorithm proposed in Nedic et al. (2008); Nedic & Ozdaglar (2009); Nedic et al. (2010), where the local optimization step is executed by sub-gradient descent.
Algorithm 1 MAZOPA: x̂i(T ) = MAZOPA (xi(1), ηt, δt,X) Require: step size ηt, DistZOO g̃i(x; δt) with exploration parameter δt for all i ∈ V Ensure: xi(1) ∈ X, ∀i ∈ V
1: for t = 1 to T do 2: Node i queries the DistZOO at point xi(t) and receives g̃i(xi(t); δt) 3: Node i computes
vi(t) = xi(t)− ηt · g̃i(xi(t); δt)
4: Node i updates its state by using the information received from its instant neighbors
xi(t+ 1) = projX ( N∑ j=1 [A(t)]ijvj(t) )
5: end for Output: x̂i(T ) = 1T ∑T t=1 xi(t)

heading: 2.3 MAIN RESULTS
text: Let x̂i(T ) be the output of Algorithm 1 at agent i. We denote the optimal solution of problem (1) by x? = arg minx∈X f(x). Defining X◦ := {x + u : x ∈ X,u ∈ B1}, we present the following results on the convergence rate of the MAZOPA algorithm.
Theorem 1 Let Assumption 1 hold. Let DistZOO take the form of g̃OPi (·). Further assume that |fi(xi(t) + δtui(t))| ≤ C for all i ∈ V. We have the following convergence results for every i ∈ V and all T ≥ 1.
(i) Consider fi ∈ Flip(Lf ,X◦) for all i ∈ V. Setting ηt = 1dT 3/4 and δt = 1 t1/4 , t = 1, . . . , T , it holds that E [ f(x̂i(T )) ] − f(x?) = O ( d
T 1/4
) .
(ii) Consider fi ∈ Flip(Lf ,X◦) ∩ Fsmo(sf ,X◦). Setting ηt = 1dT 2/3 and δt = 1 t1/6 , t = 1, . . . , T , it holds that E [ f(x̂i(T )) ] − f(x?) = O ( d
T 1/3
) .
Theorem 2 Let Assumption 1 hold. Let DistZOO take the form of g̃TPi (·). Set ηt = 1√dT and δt =
1√ t , t = 1, . . . , T . Consider fi ∈ Flip(Lf ,X◦), i ∈ V. Then, for every i ∈ V and all T ≥ 1, we have E [ f(x̂i(T )) ] − f(x?) = O (√ d T ) .
With strong convexity, the convergence rates established above can be further strengthened.
Theorem 3 Let Assumption 1 hold. Let DistZOO take the form of g̃OPi (·). Further assume that |fi(xi(t) + δtui(t))| ≤ C for all i ∈ V. We have the following convergence results for every i ∈ V and all T ≥ 1.
(i) Consider fi ∈ Flip(Lf ,X◦) ∩ Fsc(µf ,X◦) for all i ∈ V. Setting ηt = 1µf t and δt = 1 t1/3 , t = 1, . . . , T , it holds that E [ f(x̂i(T )) ] − f(x?) = O ( d2
T 1/3
) .
(ii) Consider fi ∈ Flip(Lf ,X◦)∩Fsmo(sf ,X◦)∩Fsc(µf ,X◦). Setting ηt = 1µf t and δt = 1 t1/4 , t = 1, . . . , T , it holds that E [ f(x̂i(T )) ] − f(x?) = O ( d2√ T ) .
Theorem 4 Let Assumption 1 hold. Let DistZOO take the form of g̃TPi (·). Set ηt = 1µf t and δt = 1 t , t = 1, . . . , T . Consider fi ∈ Flip(Lf ,X◦)∩Fsc(µf ,X◦), i ∈ V. Then, for every i ∈ V and all T ≥ 1, we have E [ f(x̂i(T )) ] − f(x?) = O (d ln(T ) T ) .
3 MULTISTAGE MAZOPA: ADAPTIVE LOCAL DESCENT
We now propose a multi-stage variant of Algorithm 1. We impose the following compactness assumption on the constraint set X.
Assumption 2 There exists 0 < RX <∞ such that ‖x‖ ≤ RX for all x ∈ X.

heading: 3.1 THE ALGORITHM
text: The basic idea is to divide the optimization process into a sequence of epochs, each of which has an exponentially decreasing step size and an exponentially increasing iteration number. The updates in the inner loop of each stage are just made according to Algorithm 1 with fixed step size. In each stage only the average point is maintained and used as the starting point of the next stage. This idea of setting up multi-stage optimization algorithms was originally explored in Hazan & Kale (2011).
Take positive integersm ≥ 1 and a ≥ 2. Let k\ = ⌊ loga ( T m + 1 )⌋ , where bxc represents the largest inter with value no greater than x ∈ R. We divide the T time steps into k\ epochs by
Epoch 1 : 1, . . . , T (1);
Epoch 2 : T (1) + 1, . . . , T (2);
...
Epoch k\ : T (k \−1) + 1, . . . , T (k \).
Here T (1) = m, T (2) = am, . . . , T (k \) = ak \
m. For the jth epoch, all agents will run the MAZOPA algorithm, and denote output of the j-th epoch at agent i by x̂(j)i (T
(j)). The pseudo-code of the resulting multi-stage MAZOPA is presented in Algorithm 2. Compared to the MAZOPA algorithm, the multi-stage MAZOPA has the following advantages:
(i) Multistage MAZOPA only requires each node projects its estimates onto the ball BRX , rather than the constraint set X in each epoch. In particular, multistage MAZOPA algorithm significantly reduces the number of Euclidean projections onto the constraint set X from T to k\. This makes the algorithm more computationally efficient.
(ii) Multistage MAZOPA better utilizes the step size rules, in the sense that at earlier epochs of the algorithm, larger step sizes are adopted to facilitate convergence, while smaller step sizes are adopted to achieve better accuracy at later epochs.

heading: 3.2 OPTIMAL CONVERGENCE RATES
text: We now modify the definitions of Flip, Fsmo and Fsc by replacing X with BRX , with slight abuse of notation. As it turns out, the multistage MAZOPA enjoys refined convergence rates.
Algorithm 2 Multistage MAZOPA Require: exploration parameter δ(1), step size η(1), T (1) = m, total number of iterations T , integer
a ≥ 2, and scalar b > 1 Ensure: x(1)i (1) ∈ X for all i ∈ V, and set k = 1
1: while j = 1, . . . , k\ do 2: Call Algorithm 1 to obtain
x̂ (j) i (T (j)) = MAZOPA ( x (j) i (1), η (j), δ(j),BRX )
3: Compute x(j+1)i (1) = projX ( x̂ (j) i (T (j)) ) 4: Update η(k+1) = 1aη (k) and δ(j+1) = 1b δ (j) 5: Update T (j+1) = aT (j) 6: Update j = j + 1 7: end while
Output: x̄i(T ) = projX ( x̂ (k\) i (T (k\)) )
Theorem 5 Let Assumptions 1 and 2 hold. Let DistZOO take the form of g̃TPi (·). Set a = b, T (1) = m = 1, η(1) = 4a3µf and δ (1) = 1. Consider fi ∈ Flip(Lf ,BRX) ∩ Fsc(µf ,BRX), i ∈ V. We
have E [ maxi∈V { ‖x̄i(T )− x?‖2 }] = O ( d
T+1
) .
The idea of the analysis leading to Theorem 5 can also be extended to one-point oracles. If DistZOO takes the form of g̃OPi (·), one needs to impose the following assumption on the objective functions, that is, |fi(xi(t) + δtui(t))| ≤ C for all i ∈ V. For fi ∈ Flip(Lf ,BRX) ∩ Fsc(µf ,BRX), setting a = b3, the final estimates enjoy a convergence rate of E [ maxi∈V { ‖x̄i(T ) − x?‖2 }] =
O ( d2
(T+1)1/3
) . For fi ∈ Flip(Lf ,BRX)∩Fsmo(sf ,BRX)∩Fsc(µf ,BRX), setting a = b4, there holds
E [ maxi∈V { ‖x̄i(T )− x?‖2 }] = O ( d2√ T+1 ) .

heading: 4 NUMERICAL EXAMPLES
text: In this section, we evaluate the performance of the proposed algorithms on a distributed ridge regression problem.
System setup. The optimization problem has the following form:
minimize f(x) = N∑ i=1 ( 1 2 (a T i x− bi)2 + ρ‖x‖2 ) subject to ‖x‖1 ≤ k
(5)
where x ∈ Rd is the optimization variable, the data pair (ai, bi) ∈ Rd × R is only known to node i with ai and bi being generated uniformly from the unit normal distribution.
Network setup. We implement the proposed algorithms over a randomly generated network that consists of N = 50 nodes, which is shown in Fig. 1. In the simulations, we set d = 10, k = 3/4, ρ = 1/2, and RW = 3/4. We evaluate the performance of the algorithms via the average of 10 implementations. The weight matrix associated with the graph is generated according to the maximum-degree weights:
[A(t)]ij =  1 1+dmax , (j, i) ∈ Et
1− di1+dmax , i = j 0, (j, i) /∈ Et
where dmax = maxi∈V{di} is the maximum degree of Gt (di denotes the degree of node i). Results. The performance of algorithms MAZOPA and multistage MAZOPA is illustrated via plotting the maximum function errors, maxi∈V f(x̂i(T )) and maxi∈V f(x̄i(T )), as a function of the number of iterations T in Fig. 2. As a benchmark, the convergence performance of the gradient
algorithm is displayed in Fig. 2 as well. From the numerical results it is clear that the maximum function errors are vanishing for all zerroth-order algorithms. In fact, the convergence performance of two-point MAZOPA is even comparable to the gradient method. Moreover, the multistage variants in general exhibit better convergence performance, and this is more obvious for the case of two-point MAZOPA. These numerical results are in compliance with the theoretical findings in the paper.
Reproduction of the results. The code used for producing this numerical example is provided in the suplementary material.

heading: 5 CONCLUSIONS
text: We have established a series of convergence rates for distributed zeroth-order subgradient algorithms that match their centralized counterpart for Lipschitz, smooth, and strongly convex function classes. These results provided the theoretical benchmarks for zeroth-order approaches over complex dynamic networks. We also proposed a multi-stage variant of the algorithm that better utilizes the learning rates and attains even improved convergence rates. In future work, it is worth exploring the connection between the convergence rates and the underlying communication complexity for distributed zeroth-order algorithms.

heading: A KEY LEMMAS
text: We first establish the basic convergence result for Algorithm 1 that is based on DistZOO g̃i(x; δt), which plays a crucial role in subsequent analyses. We will sometimes use i• to denote a node in V just to highlight the focus on a given node (but i• indeed may take any value in V and therefore it is a generic node).
Lemma 1 Let Assumption 1 hold. Let g̃i(xi(t); δt) be a DistZOO that satisfies Definition 1 (i) and (ii). Then, for any i• ∈ V and all T ≥ 1, there holds
T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ T∑ t=1 N∑ i=1 E [∣∣fi(x?)− f̂i(x?; δt)∣∣]
+ T∑ t=1 N∑ i=1 E [∣∣fi(xi•(t))− f̂i(xi•(t); δt)∣∣]
+ T∑ t=1 E [Λ?(t)]− E [Λ?(t+ 1)] 2ηt + p1Lf
+ 1
2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] + p2Lf
T−1∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖ ] where p1 = 2N maxi∈V{‖xavg(1)−xi(1)‖}+ 2Nαβ1−β (∑N i=1 ‖xi(1)‖ ) , p2 = 2N ( 3α 1−β + 4 ) ,
and Λ?(t) = ∑N i=1 ‖xi(t) − x?‖2 with xavg(1) = 1 N ∑N i=1 xi(1), α = ( 1− ξ4N2 )−2 , and
β = ( 1− ξ4N2 )1/B .
Before presenting the proof of Lemma 1, we provide the following two supporting lemmas. The first lemma characterizes the convergence property of the transition matrix induced by weight matrix A(t) (see Nedic et al. (2008)).
Lemma 2 Define the transition matrix as A(t : `) = A(t)A(t − 1) · · ·A(` + 1)A(`) for all t ≥ ` ≥ 1, and write A(t : t) = A(t). A(t : `) satisfies∣∣∣∣[A(t : `)]ij − 1N ∣∣∣∣ ≤ αβt−`+1 where α = ( 1− ξ4N2 )−2 and β = ( 1− ξ4N2 )1/B .
The second lemma establishes the accumulated disagreement for every node in the network.
Lemma 3 (Disagreement) Let Assumption 1 hold. For every node i ∈ V, we have T∑ t=1 ‖xavg(t)− xi(t)‖ ≤ ‖xavg(1)− xi(1)‖+ αβ 1− β ( N∑ i=1 ‖xi(1)‖ )
+
( 3α
1− β + 4 ) T−1∑ t=1 ηt N∑ i=1 ‖g̃i(xi(t); δt)‖
where xavg(t) = 1N ∑N i=1 xi(t).
Proof. To simplify the presentation, we denote
ṽi(t) = N∑ j=1 [A(t)]ijvj(t)
si(t) = projX (ṽi(t))− ṽi(t).
Step 3 in Algorithm 1 can be rewritten as
xi(t+ 1) = ṽi(t) + si(t).
Our analysis relies on the estimate of ‖si(t)‖, which can be bounded as follows:
‖si(t)‖ ≤ ∥∥∥projX (ṽi(t))− N∑
j=1
[A(t)]ijxj(t) ∥∥∥+ N∑
j=1
[A(t)]ij ‖ηtg̃j(xj(t); δt)‖
where the inequality is based on Step 3 in Algorithm 1 and A(t) is double stochastic (cf. Assumption 1). Using the non-expansiveness of the Euclidean projection projX(·) and the fact that∑N j=1[A(t)]ijxj(t) ∈ X, we have
‖si(t)‖ ≤ 2 N∑ j=1 [A(t)]ij ‖ηtg̃j(xj(t); δt)‖ . (6)
We now derive the general expressions for xavg(t+ 1) and xi(t+ 1), respectively. For xavg(t+ 1), we have
xavg(t+ 1) = xavg(t)− 1
N N∑ i=1 ηtg̃i(xi(t); δt) + 1 N N∑ i=1 si(t).
Applying the preceding inequality recursively, we get
xavg(t+ 1) = xavg(1)− t∑ `=1 1 N N∑ i=1 η`g̃i(xi(`); δ`) + t∑ `=1 1 N N∑ i=1 si(`). (7)
Similarly, for xi(t+ 1), we have
xi(t+ 1) = N∑ j=1 [A(t : 1)]ijxj(1)− t∑ `=1 N∑ j=1 [A(t : `)]ijη`g̃j(xj(`); δ`) + t−1∑ `=1 N∑ j=1 [A(t : `+ 1)]ijsj(`) + si(t). (8) Combining (7) and (8), gives ‖xavg(t+ 1)− xi(t+ 1)‖ ≤ N∑ j=1 ∣∣∣∣[A(t : 1)]ij − 1N ∣∣∣∣ ‖xj(1)‖+ t∑ `=1 N∑ j=1 ∣∣∣∣[A(t : `)]ij − 1N ∣∣∣∣ η`‖g̃j(xj(`); δ`)‖
+ t−1∑ `=1 N∑ j=1 ∣∣∣∣[A(t : `+ 1)]ij − 1N ∣∣∣∣ ‖sj(`)‖+ ‖si(t)‖+ 1N N∑ i=1 ‖si(t)‖.
(9) Combining the results in (6), (9) and Lemma 2, leads to
‖xavg(t+ 1)− xi(t+ 1)‖ ≤ αβt (
N∑ i=1 ‖xi(1)‖
) + 3α
t∑ `=1 βt−`η` N∑ i=1 ‖g̃i(xi(`); δ`)‖+ 4ηt N∑ i=1 ‖g̃i(xi(t); δt)‖
where we used the following relation, based on (6):
‖si(t)‖+ 1
N N∑ i=1 ‖si(t)‖ ≤ 2 N∑ i=1 ‖si(t)‖ ≤ 4 N∑ i=1 N∑ j=1 [A(t)]ij ‖ηtg̃j(xj(t); δt)‖ ≤ 4ηt N∑ i=1 ‖g̃i(xi(t); δt)‖.
This implies that
T∑ t=1 ‖xavg(t)− xi(t)‖ ≤ ‖xavg(1)− xi(1)‖+ α ( N∑ i=1 ‖xi(1)‖ ) T−1∑ t=1 βt
+ 3α T−1∑ t=1 t∑ `=1 βt−`η` N∑ i=1 ‖g̃i(xi(`); δ`)‖+ 4 T−1∑ t=1 ηt N∑ i=1 ‖g̃i(xi(t); δt)‖.
(10)
This, in combination with (10), leads to the final bound.
[Proof of Lemma 1]. Denote
Λ(t) = N∑ i=1 ‖xi(t)− x‖2, ∀x ∈ X, t ≥ 1. (11)
We follow the standard analysis by deriving the general evolution of ∆(t),
Λ(t+ 1) = N∑ i=1 ∥∥∥∥∥projX  N∑ j=1 [A(t)]ijvj(t) − x∥∥∥∥∥ 2 ≤ N∑ i=1 ‖vi(t)− x‖2 (12)
where the inequality follows from the non-expansiveness of the Euclidean projection and the convexity of norm square function. Expanding the term further gives
Λ(t+ 1) = Λ(t) + N∑ i=1 ‖ηtg̃i(xi(t); δt)‖2 − 2ηt N∑ i=1 〈g̃i(xi(t); δt),xi(t)− x〉 (13)
Taking the expectation on both sides and using the following property of DistZOO (cf. Definition 1(i)):
E [g̃i(xi(t); δt)] = ∇f̂i(xi(t); δt) we further obtain
E [Λ(t+ 1)] = E [Λ(t)] + η2t N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] − 2ηt N∑ i=1 ( E [ f̂i(xi(t); δt) ] − f̂i(x; δt) ) (14)
which implies T∑ t=1 N∑ i=1 E [ f̂i(xi(t); δt) ] − f̂(x; δt) ≤ T∑ t=1 E [Λ(t)]− E [Λ(t+ 1)] 2ηt + 1 2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] . (15) We turn our attention to the left-hand side of (15). By adding and subtracting the term f̂i(xi•(t); δt) and using the Lipschitz continuity of f̂i (cf. Definition 1(ii)), it follows that
T∑ t=1 N∑ i=1 f̂i(xi(t); δt) ≥ T∑ t=1 f̂(xi•(t); δt)− Lf T∑ t=1 N∑ i=1 ‖xi(t)− xi•(t)‖
which, together with (15), yields
T∑ t=1 E [ f̂(xi•(t); δt) ] − f̂(x; δt) ≤ T∑ t=1 E [Λ(t)]− E [Λ(t+ 1)] 2ηt
+ 1
2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] + Lf T∑ t=1 N∑ i=1 ‖xi(t)− xi•(t)‖.
(16)
The desired result follows by relating the left-hand side to the original function f , using the disagreement estimate in Lemma 3, and setting x = x?.

heading: B PROOFS OF THEOREMS 1 AND 2
text: We first provide the the following lemma that characterizes the properties of the DistZOOs in (3) and (4). Its proof can be derived by resorting to Flaxman et al. (2005); Shamir (2017), which we omit here to save space.
Lemma 4 Suppose that fi ∈ Flip(Lf ) for all i ∈ V. We have the following.
(i) For g̃OPi (·), there hold pd = 1 and
E [ ‖g̃OPi (xi(t); δt)‖2 ] ≤ ( Cd
δt )2 where C = maxi∈V |fi(xi(t) + δtui(t))| with xi(t) ∈ X and ui(t) uniformly drawn from B1. We have p̃d = 1 when fi ∈ Fsmo(sf ).
(ii) For g̃TPi (·), there hold pd = 1 and E [ ‖g̃TPi (xi(t); δt)‖2 ] ≤ cL2fd
where c is some universal constant. In addition, p̃d = 1 when fi ∈ Fsmo(sf ).
Now we are ready to prove Theorems 1 and 2.
(i) First, using the property of the DistZOO g̃OPi (·) (cf. Definition 1(iii)), it follows that N∑ i=1 ∣∣fi(xi•(t))− f̂i(xi•(t); δt)∣∣ ≤ NLfδt N∑ i=1 ∣∣fi(x?)− f̂i(x?; δt)∣∣ ≤ NLfδt (17)
We now focus on the case of fi ∈ Flip(Lf ,X◦). Combining with inequality (17), the results in Lemmas 4(i) and 1, gives
T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ p1Lf + 2NLf T∑ t=1 δt + 1 2ηt E [Λ?(1)]
+ 1
2 NC2d2 T∑ t=1 ηt δ2t + p2NLfCd T∑ t=1 ηt δt
(18)
where we used the fact that ηt is a function of T and Jensen’s inequality, i.e., E [ ∥∥g̃OPi (xi(t); δt)∥∥ ] ≤ (E[ ∥∥g̃OPi (xi(t); δt)∥∥2 ])1/2 ≤ Cdδt .
Substituting the explicit expressions of η = 1 dT 3/4 and δt = 1t1/4 into (18) and dividing both sides by T , we find that
1
T T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) = O ( d T 1/4 ) (19)
where we used the following inequality thet ∑T t=1 t
a = O(T 1+a),∀a 6= −1. The desired result follows by using the convexity of function f , i.e., 1T ∑T t=1 f(xi•(t)) ≥ f(x̂i•(T )). When fi ∈ Flip(Lf ,X◦) ∩ Fsmo(sf ,X◦), it follows from the property of the DistZOO g̃OPi (·) (cf. Definition 1(iii)) that
N∑ i=1 ∣∣fi(xi•(t))− f̂i(xi•(t); δt)∣∣ ≤ 1 2 Nsfδ 2 t
N∑ i=1 ∣∣fi(x?)− f̂i(x?; δt)∣∣ ≤ 1 2 Nsfδ 2 t .
(20)
We then combine the preceding inequality and the results in Lemmas 4(i) and 1 to get T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ p1Lf +Nsf T∑ t=1 δ2t + 1 2ηt E [Λ?(1)]
+ 1
2 NC2d2 T∑ t=1 ηt δ2t + p2NLfCd T∑ t=1 ηt δt .
(21)
In contrast with the bound in (18), the second term on the right-hand side of (21) now becomes Nsf ∑T t=1 δ 2 t , which gives us much more space when choosing δt; we can show that the choices of η = 1 dT 2/3 and δt = 1t1/6 yield the optimal convergence rate O ( d T 1/3 ) .
(ii) When fi ∈ Flip(Lf ,X◦), we have the following result for Algorithm 1 running with DistZOO g̃TPi (·):
T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ p1Lf + 2NLf T∑ t=1 δt + 1 2ηt E [Λ?(1)]
+
( 1
2 cd+ p2
√ c √ d ) NL2fηtT
(22)
where we have used the bounds in (17), Lemma 4(ii) and Lemma 1. Then we can deduce from the terms 1ηt and ηtT that the optimal choice of ηt is 1√ T
. Hence, substituting the explicit expressions for ηt = 1√dT and δt = 1√ t
into (22), dividing both sides with T , and using the convexity of function F , the desired bound can be concluded.
When fi ∈ Flip(Lf ,X◦) ∩ Fsmo(sf ,X◦), it can be shown that the term (2NLf ∑T t=1 δt) on the
right-hand side of (22) is replaced by (Nsf ∑T t=1 δ 2 t ), because of (20). As we discussed in the case of fi ∈ Flip(Lf ,X◦), the convergence rate is determined by the terms involving 1ηt and ηtT . Hence, the convergence rate is the same as that of the case when fi in only Lipschitz continuous. The proof is complete.

heading: C PROOFS OF THEOREM 3 AND 4
text: First, we claim that for the DistZOOs g̃OPi (·) and g̃TPi (·), the strongly convexity of fi implies the strongly convexity of its smoothed variant f̂i, and its proof is straightforward. We now establish the basic convergence results for Algorithm 1 running with g̃OPi (·) and g̃TPi (·). It follows from (13) that
E [Λ(t+ 1)] = E [Λ(t)] + η2t N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] − 2ηt N∑ i=1 E [〈 ∇f̂i(xi(t); δt),xi(t)− x 〉] ≤ E [Λ(t)] + η2t
N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] − 2ηt N∑ i=1 ( E [ f̂i(xi(t); δt) ] − f̂i(x) ) − µfηtE [Λ(t)]
(23) where in the equality we used the relation E [g̃i(xi(t); δt)] = ∇f̂i(xi(t); δt) (cf. Definition 1(i)) and in the inequality we used the strongly convexity of function f̂i. Summing the inequalities in (23) over t = 1 to t = T and regrouping the terms, we obtain T∑ t=1 N∑ i=1 ( E [ f̂i(xi(t); δt) ] − f̂i(x) ) ≤ 1 2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2
] + 1
2 T∑ t=1 ( 1 ηt (E [Λ(t)]− E [Λ(t+ 1)])− µfE [Λ(t)] )
= 1
2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] + 1 2 ( 1 η1 −µf ) E [Λ(1)]
+ 1
2 T∑ t=2 ( 1 ηt − 1 ηt−1 −µf ) E [Λ(t)]− 1 2ηT E [Λ(T + 1)] .
(24) By substituting the expression for ηt = 1µf t into (24) and dropping the negative term, it follows
T∑ t=1 N∑ i=1 ( E [ f̂i(xi(t); δt) ] − f̂i(x) ) ≤ 1 2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] . (25)
Then, following the same lines as that of the proof of Lemma 1, we have that, for DistZOOs g̃OPi (·) and g̃TPi (·), T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ T∑ t=1 E [∣∣f(x?)− f̂(x?; δt)∣∣]+ T∑ t=1 E [∣∣f(xi•(t))− f̂(xi•(t); δt)∣∣]
+ p1Lf + 1
2 T∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖2 ] + p2Lf T−1∑ t=1 ηt N∑ i=1 E [ ‖g̃i(xi(t); δt)‖ ] .
(26)
(i) We now derive the convergence rate results for Algorithm 1 running with DistZOOs g̃OPi (·). When fi ∈ Flip(Lf ,X◦) ∩ Fsc(µf ,X◦), we combine the results in (17), (26) and Lemma 4(i) to get T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ p1Lf + 2NLf T∑ t=1 δt + 1 2 NC2d2 T∑ t=1 ηt δ2t + p2NLfCd T∑ t=1 ηt δt .
By substituting ηt = 1µf t and δt = 1 t1/3 into the preceding inequality and using the convexity of F , it yields the following optimal bound
1
T T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) = O ( d2 T 1/3 ) .
When fi ∈ Flip(Lf ,X◦) ∩ Fsmo(sf ,X◦) ∩ Fsc(µf ,X◦), it follows from an argument similar to that of (21) that
T∑ t=1 ( E [ f(xi•(t)) ] − f(x?) ) ≤ p1Lf +Nsf T∑ t=1 δ2t + 1 2 NC2d2 T∑ t=1 ηt δ2t + p2NLfCd T∑ t=1 ηt δt .
It can be proven that the choice of δt = 1t1/4 yields the optimal convergence rate, that is, O ( d2√ T ) .
(ii) The proof for Algorithm 1 running with DistZOO g̃TPi (·) can be obtained in a similar way by exploiting the properties of the DistZOO g̃TPi (·).

heading: D PROOF OF THEOREM 5
text: We provide the basic convergence result for each stage k, and we start by deriving a similar bound as that of Lemma 1 as follows: T (k)∑ t=1 ( E [ f(x (k) i• (t)) ] − f(x?) ) ≤ T (k)∑ t=1 E [∣∣f(x?)− f̂(x?; δ(k))∣∣]+ T (k)∑ t=1 E [∣∣f(x(k)i• (t))− f̂(x(k)i• (t); δ(k))∣∣]
+ T (k)∑ t=1
E [ Λ(k),?(t) ] − E [ Λ(k),?(t+ 1) ] 2η(k) + p (k) 1 Lf + 1 2 T (k)∑ t=1 η(k) N∑ i=1 E [∥∥g̃i(x(k)i (t); δ(k))∥∥2]
+ p2Lf T (k)∑ t=1 η(k) N∑ i=1 E [∥∥g̃i(x(k)i (t); δ(k))∥∥]
(27) where Λ(k),?(t) = ∑N i=1 ‖x (k) i (t) − x?‖2 and p (k) 1 satisfies the following bound, according to compactness of the set X,
p (k) 1 = 2N max i∈V {‖x(k)avg(1)− x (k) i (1)‖}+
2Nαβ
1− β ( N∑ i=1 ‖x(k)i (1)‖ ) ≤ ( 4N + 2αβ 1− β N2 ) RX.
(28)
The left-hand side of (27) can be further bounded by using the strongly convexity of F , that is,
1
T (k) T (k)∑ t=1 ( f(x (k) i• (t))− f(x ?) ) ≥ 〈 ∇f(x?), x̂(k)i• (T (k))− x? 〉 + Nµf 2 ‖x̂(k)i• (T (k))− x?‖2
Applying the first-order optimality condition to the preceding inequality, i.e., 〈 ∇F (x?),x−x? 〉 ≥ 0 for any x ∈ X, yields
1
T (k) T (k)∑ t=1 ( f(x (k) i• (t))− f(x ?) ) ≥ Nµf 2 ‖x(k+1)i• (1)− x ?‖2 (29)
where the second inequality follows from the non-expansiveness of the Euclidean projection projX(·), and the last equality from Step 3 in Algorithm 2. Combining the inequalities (27), (28) and (29), we have for any i• ∈ V,
Nµf 2
E [ ‖x(k+1)i• (1)− x ?‖2 ] ≤ ∑N i=1 E [ ‖x(k)i (1)− x?‖2 ] 2η(k)T (k) + ( 4N + 2αβ 1− β N2 ) LfRX 1 T (k)
+ 1
T (k) T (k)∑ t=1 E [∣∣f(x?)− f̂(x?; δ(k))∣∣]+ 1 T (k) T∑ t=1 E [∣∣f(x(k)i• (t))− f̂(x(k)i• (t); δ(k))∣∣]
+ 1
2T (k) T (k)∑ t=1 η(k) N∑ i=1 E [∥∥g̃i(x(k)i (t); δ(k))∥∥2]+ p2Lf 1T (k) T (k)∑ t=1 η(k) N∑ i=1 E [∥∥g̃i(x(k)i (t); δ(k))∥∥].
(30) From (30) we find that the convergence depends on the properties of the DistZOOs, and we first derive the dimension-dependence error bounds for DistZOO g̃OPi (·) and the bounds for g̃TPi (·) naturally follows from the derivations.
For DistZOO g̃OPi (·), when fi ∈ Flip(Lf ,BRX) ∩ Fsc(µf ,BRX), it follows from (17), (30) and Lemma 4(i) that
E [
max i∈V
{ ‖x(k+1)i (1)− x ?‖2 }] ≤ E [ max i∈V { ‖x(k)i (1)− x ?‖2 }] 1 µfη(k)T (k)
+ 4 Lf µf δ(k) + 4
( 2 + αβ
1− β N ) Lf µf RX 1 T (k) + 2p2 Lf µf Cd η(k) δ(k) + 1 µf C2d2 η(k) (δ(k))2 . (31)
On the other hand, we have
T (k) = ak−1T (1), η(k) = 1
ak−1 η(1), δ(k) =
1
bk−1 δ(1) (32)
This, together with inequality (31), gives E [
max i∈V
{ ‖x(k+1)i (1)− x ?‖2 }] ≤ E [ max i∈V { ‖x(k)i (1)− x ?‖2 }] 1 µfη(1)T (1)
+ 1( min { a, b, ab , a b2
})k−1 × (
4 Lf µf δ(1) + 4
( 2 + αβ
1− β N ) Lf µf RX 1 T (1)
+ 2p2 Lf µf Cd
η(1) δ(1) + 1 µf C2d2 η(1) (δ(1))2
) .
(33)
By substituting T (1) = m = 1, η(1) = 4min{b, a b2 }
3µf and δ(1) = 1 into the preceding relation, we
arrive at E [
max i∈V
{ ‖x(k+1)i (1)− x ?‖2 }] ≤ E [ max i∈V { ‖x(k)i (1)− x ?‖2 }] 3 4 min { b, ab2 } + R1( min { b, ab2
})k−1 (34)
where we have used the fact that min { a, b, ab , a b2 } = min { b, ab2 } , due to a > b2 (because of a b2 > 1) and b > 1, and
R1 = 4 Lf µf δ(1) + 4
( 2 + αβ
1− β N ) Lf µf RX 1 T (1) + 2p2 Lf µf Cd η(1) δ(1) + 1 µf C2d2 η(1) (δ(1))2 .
We next show by induction that E [
max i∈V
{ ‖x(k)i (1)− x ?‖2 }] ≤ 4 max{R1, R 2 X}
hk−2 (35) where h = min { b, ab2 } . For k = 1, we can use the following bound to deduce that inequality (35)
holds, maxi∈V { ‖x(1)i (1) − x?‖2 } ≤ 4R2X ≤ 4hmax{R1, R2X}. We then assume that inequality (35) holds for k and show it holds for k + 1 as well,
E [
max i∈V
{ ‖x(k+1)i (1)− x ?‖2 }] ≤ 3 4h E [ max i∈V { ‖x(k)i (1)− x ?‖2 }] + R1 hk−1
≤ 3 max{R1, R 2 X}
hk−1 + max{R1, R2X} hk−1 ≤ 4 max{R1, R 2 X} hk−1
which leads to the conclusion in (35). It is easy to verify that the total number of stages in Algorithm 2 is k\ = ⌊ loga ( T m + 1 )⌋ , and the final estimates returned by Algorithm 2 are x̄i(T ), i ∈ V. Hence, applying k\ + 1 to (35), we have
E [
max i∈V
{ ‖x̄i(T )− x?‖2 }] ≤ 4 max{R1, R 2 X}
hk\+1−2 ≤ 4h 2 max{R1, R2X} hloga( T m+1) = 4h2 max{R1, R2X}(
T m + 1
) 1 logh(a)
(36) where we used the inequality that k\ ≥ loga ( T m + 1 ) − 1. We are left to find the minimum of logh(a) = logmin{b, a b2 }(a), which is achieved when b = a b2 . This yields the following final conver-
gence rate, that is, E [ maxi∈V { ‖x̄i(T ) − x?‖2 }] ≤ 4b
2 max{R1,R2X} (T+1)1/3
= O ( d2
(T+1)1/3
) , where the
equality follows from R1 = O ( d2 ) .
When fi ∈ Flip(Lf ,BRX) ∩ Fsmo(sf ,BRX) ∩ Fsc(µf ,BRX), the dimension-dependence bound can be obtained in a similar fashion.
For DistZOO g̃TPi (·), when fi ∈ Flip(Lf ,BRX) ∩ Fsc(µf ,BRX), it follows from (17), (30) and Lemma 4(ii) that
E [
max i∈V
{ ‖x(k+1)i (1)− x ?‖2 }] ≤ E [ max i∈V { ‖x(k)i (1)− x ?‖2 }] 1 µfη(k)T (k)
+ 4 Lf µf δ(k) + 4
( 2 + αβ
1− β N ) Lf µf RX 1 T (k) + L2f µf ( cd+ 2p2 √ c √ d ) η(k). (37)
By setting b = a and then substituting T (1) = 1, η(1) = 4a3µf and δ (1) = 1 into the preceding inequality it follows that
E [
max i∈V
{ ‖x(k+1)i (1)− x ?‖2 }] ≤ E [ max i∈V { ‖x(k)i (1)− x ?‖2 }] 3 4a + R2 ak−1
(38)
where R2 = 4 Lf µf δ(1) + 4 ( 2 + αβ1−βN ) Lf µf RX 1 T (1) + L2f µf ( cd+ 2p2 √ c √ d ) η(1). Then, following an argument similar to that of part (i), we obtain
E [
max i∈V
{ ‖x̄i(T )− x?‖2 }] ≤ 4a
2 max{R2, R2X} T + 1 = O
( d
T + 1
) .
Similarly, when fi ∈ Flip(Lf ,BRX) ∩ Fsmo(sf ,BRX) ∩ Fsc(µf ,BRX), we have E [
max i∈V
{ ‖x̄i(T )− x?‖2 }] ≤ 4a
2 max{R′2, R2X} T + 1 = O
( d
T + 1 ) where R′2 = 2 sf µf (δ(1))2 + 4 ( 2 + αβ1−βN ) Lf µf RX 1 T (1) + L2f µf ( cd+ 2p2 √ c √ d ) η(1). The proof is complete.

Review summaries:
This paper proposes a zeroth-order optimization algorithm for distributed, multi-agent systems with time-varying communication networks. The authors show that their presented multi-agent zeroth-order projection averaging algorithm (and its improved multi-stage version) has a convergence rate that matches the centralized counterpart algorithms under different assumptions. A small numerical experiment is also conducted to illustrate their theoretical findings.