Distributed Zero-Order Optimization under Adversarial Noise
Year: 2021

Abstract:
We study the problem of distributed zero-order optimization for a class of strongly convex functions. They are formed by the average of local objectives, associated to different nodes in a prescribed network. We propose a distributed zero-order projected gradient descent algorithm to solve the problem. Exchange of information within the network is permitted only between neighbouring nodes. An important feature of our procedure is that it can query only function values, subject to a general noise model, that does not require zero mean or independent errors. We derive upper bounds for the average cumulative regret and optimization error of the algorithm which highlight the role played by a network connectivity parameter, the number of variables, the noise level, the strong convexity parameter, and smoothness properties of the local objectives. The bounds indicate some key improvements of our method over the state-of-the-art, both in the distributed and standard zero-order optimization settings. We also comment on lower bounds and observe that the dependency over certain function parameters in the bound is nearly optimal.

heading: 1 Introduction
text: We study the problem of distributed optimization where each node (or agent) has an objective function fi : Rd → R and exchange of information is limited between neighbouring agents within a prescribed network of connections. The goal is to minimize the average of these objectives on a closed bounded convex set Θ ⊂ Rd,
min x∈Θ
f(x) where f(x) = 1
n n∑ i=1 fi(x). (1)
Distributed optimization has been widely studied in the literature, we refer to Tsitsiklis et al. [1986], Nedic and Ozdaglar [2009], Nedic et al. [2010], Boyd et al. [2011], Duchi et al. [2012], Jakovetić et al. [2014], Lobel et al. [2011], Kia et al. [2015], Shi et al. [2014], Jakovetić [2019], Scaman et al. [2019], Pu et al. [2021] and references therein. This problem has broad applications such as multi-agent target seeking Liu et al. [2017], distributed learning Kraska et al. [2013], and wireless networks Park et al. [2020], among others.
We address problem (1) from the perspective of zero-order distributed optimization. That is we assume that only function values can be queried by the algorithm, subject to measurement noise.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
During the optimization procedure, each agent maintains a local copy of the variables which are sequentially updated based on local and neighboring functions’ queries. We wish to devise such optimization procedures which are efficient in bounding the average optimization error and cumulative regret in terms of the functions’ properties and network topology.
Contributions Our principal contribution is a distributed zero-order optimization algorithm, introduced in Section 2, which we show to achieve tight rates of convergence under certain assumptions on the objective functions, outlined in Section 3. Specifically, we consider that the local objectives fi are β-Hölder and the average objective f is α-strongly convex. The algorithm relies on a novel zero-order gradient estimator, presented in Section 4. Although conceptually very simple, this estimator, when employed within our algorithm, allows us to obtain an O(d2) computational gain as well as improved error rates than previous state-of-the-art zero-order optimization procedures Akhavan et al. [2020], Bach and Perchet [2016], in the special case of standard (undistributed) setting. Another key advantage of our approach is due to the general noise model presented in Section 5, under which function values are queried. The noise variables do not need to be zero mean or independently sampled, and thus they include “adversarial” noise. In Section 6, we derive the rates of convergence for the cumulative regret and the optimization error of the proposed algorithm, and in Section 7 we consider the special case of 2-smooth functions. The rates highlight the dependency with respect to the number of variables d, the number of function queries T , the spectral gap of the network matrix 1− ρ, and the parameters n, α and β. The bounds enjoy a better dependency on 1− ρ than previous bounds on zero-order distributed optimization Qu and Li [2018], Yu et al. [2019], Tang et al. [2019]. We also compare our bounds to related lower bounds in Akhavan et al. [2020] for undistributed setting, observing that our rates are optimal either with respect to T and α, or with respect to T and d.
Previous Work We briefly comment on previous related work and defer to Section 8 for a more in depth discussion and comparison. For both deterministic and stochastic scenarios of problem (1), a large body of literature is devoted to first-order gradient based methods with a consensus scheme (see the papers cited above and references therein). On the other hand, the study of zero-order methods was started only recently Qu and Li [2018], Sahu et al. [2018b,a], Hajinezhad et al. [2019], Yu et al. [2019], Tang et al. [2019]. The works Qu and Li [2018], Yu et al. [2019], Tang et al. [2019] are dealing with zero-order distributed methods in noise-free settings while the noisy setting is developed in Hajinezhad et al. [2019], Sahu et al. [2018b,a]. Namely, Hajinezhad et al. [2019] considers 2-point zero-order methods with stochastic queries for non-convex optimization but assume that the noise is the same for both queries, which makes the problem analogous to noise-free scenario in terms of optimization rates. Papers Sahu et al. [2018b,a] study zero-order distributed optimization for strongly convex and β-smooth functions fi with β ∈ {2, 3}. They derive bounds on the optimization error, though without providing closed form expressions.
Notation Throughout we denote by 〈·, ·〉 and ‖ · ‖ be the standard inner product and Euclidean norm on Rd, respectively, and by ‖ · ‖∗ the spectral norm of a matrix. The notation I is used for the n-dimensional identity matrix and 1 for the vector in Rn with all entries equal to 1. We denote by ej the j-th canonical basis vector in Rd. For any set A, the number of elements in A is denoted by |A|. For x ∈ R, the value bxc is the maximal integer less than x. For every closed convex set Θ ⊂ Rd and x ∈ Rd we denote by ProjΘ(x) = argmin{‖z − x‖ : z ∈ Θ} the Euclidean projection of x onto Θ. We denote by diam(Θ) the Euclidean diameter of Θ. Finally we let U [−1, 1] be the uniform distribution on [−1, 1].

heading: 2 The Problem
text: Let n be the number of agents and let G = (V,E) be an undirected graph, where V = {1, . . . , n} is the set of nodes and E ⊆ V × V is the set of edges. The adjacency matrix of G is the symmetric matrix (Aij)ni,j=1 defined as Aij = 1, if (i, j) ∈ E and zero otherwise. We consider the following sequential learning framework, where each agent i gets values of function fi corrupted by noise and shares information with other agents. At step t, agent i acts as follows:
• makes queries and gets noisy values of fi, • provides a local output ui(t) based on these queries and on the past information, • broadcasts ui(t) to neighboring agents,
Algorithm 1 Distributed Zero-Order Gradient
Input Communication matrix (Wij)ni,j=1, step sizes (ηt > 0) T0−1 t=1 Initialization Choose initial vectors x1(1) = · · · = xn(1) ∈ Rd For t = 1, . . . , T0 − 1
For i = 1, . . . , n 1. Build an estimate gi(t) of the gradient∇fi(xi(t)) using noisy evaluations of fi 2. Update xi(t+1)= ∑n k=1WikProjΘ(x
k(t)− ηtgk(t)) End
End Output Approximate minimizer x̄(T0) = 1n ∑n i=1 x i(T0) of the average objective f = 1n ∑n i=1 fi
• updates its local variable using information from other agents as follows:
xi(t+1)= n∑ j=1 Wiju j(t),
where W = (Wij)ni,j=1 is a given matrix called the consensus matrix.
Below we use the following condition on the consensus matrix. Assumption A. Matrix W is symmetric, doubly stochastic, and ρ := ∥∥W − n−111>∥∥∗ < 1. Matrix W accounts for the connectivity properties of the network. If Wij = 0 the agents i and j are not connected (do not exchange information). Often W is defined as a doubly stochastic matrix function of the adjacency matrix A of the graph. One popular example is as follows:
Wij =  Aij γmax{d(i),d(j)} if i 6= j,
1− ∑ k:k 6=i Aki γmax{d(i),d(k)} if i = j,
where d(i) = ∑n j=1Aij is the degree of node i and γ > 0 is a constant. Then, clearly, W = (Wij) is a symmetric and doubly stochastic matrix, andWij = 0 if agents i and j are not connected. Moreover, we have ρ < 1 − c/n2 for a constant c > 0 (see Qu and Li [2018], Olshevsky [2014]). Values of spectral gaps ρ for some other W reflecting different network topologies can be found in Duchi et al. [2012]. Typically, ρ < 1− an, where an = Ω(n−1) or an = Ω(n−2). Parameter ρ can be viewed as a measure of difference between the distributed problem and a standard optimization problem. If the graph of communication is a complete graph a natural choice is W = n−11n1>n and then ρ = 0. For more examples of consensus matrices W , see Olshevsky and Tsitsiklis [2009], Duchi et al. [2012] and references therein.
The local outputs ui can be defined in different ways. Our approach is outlined in Algorithm 1. At Step 1, an estimate of the gradient of the local objective fi at xi(t) is constructed. This involves a randomized procedure that we describe and justify in Section 4. The local output ui is defined as an update of the projected gradient algorithm with such an estimated gradient. At Step 2 of the algorithm, each agent computes the next point by a local consensus gradient descent step, which uses local and neighbor information. Step 2 of the algorithm is known as gossip method, see e.g., Boyd et al. [2006]), which was initially introduced as an approach for the networks with the imposed connection between the nodes changing by time. We also refer to Sayin et al. [2017] for similar algorithms in the context of distributed stochastic first-order gradient methods.

heading: 3 Assumptions on Local Objectives
text: In this section, we give some definitions and introduce our assumptions on the local objective functions f1, . . . , fn. Definition 1. Denote by Fβ(L) the set of all functions f : Rd → R that are ` = bβc times differentiable and satisfy, for all x, z ∈ Rd the Hölder-type condition∣∣∣∣f(z)− ∑
0≤|m|≤`
1
m! Dmf(x)(z − x)m
∣∣∣∣ ≤ L‖z − x‖β , (2)
Algorithm 2 Gradient Estimator with 2d Queries
Input Function F : Rd → R and point x ∈ Rd Requires Kernel K : [−1, 1]→ R, parameter h > 0 Initialization Generate random r from uniform distribution on [−1, 1] For j = 1, . . . , d
1. Obtain noisy values yj = F (x+ hrej) + ξj and y′j = F (x− hrej) + ξ′j 2. Compute gj = 12h (yj − y ′ j)K(r)
End Output g = (gj)dj=1 ∈ Rd estimator of∇F (x)
where L > 0, the sum is over the multi-index m = (m1, ...,md) ∈ Nd, we used the notation m! = m1! · · ·md!, |m| = m1 + · · ·+md, and we defined, for every ν = (ν1, . . . , νd) ∈ Rd,
Dmf(x)νm = ∂|m|f(x)
∂m1x1 · · · ∂mdxd νm11 · · · ν md d .
Elements of the class Fβ(L) are referred to as β-Hölder functions. Definition 2. Function f : Rd → R is called 2-smooth if it is differentiable on Rd and there exists L̄ > 0 such that, for every (x, x′) ∈ Rd × Rd, it holds that
‖∇f(x)−∇f(x′)‖ ≤ L̄‖x− x′‖.
Definition 3. Let α > 0. Function f : Rd → R is called α-strongly convex if f is differentiable on Rd and
f(x)− f(x′) ≥ 〈∇f(x′), x− x′〉+ α 2 ‖x− x′‖2 , ∀x, x′ ∈ Rd.
Assumption B. Functions f1, . . . , fn: (i) belong to the class Fβ(L), for some β ≥ 2, and (ii) are 2-smooth.
In Section 6 we will analyse the convergence properties of Algorithm 1 when the objective function f in 1 is α-strongly convex. We stress that we do not need the functions f1, . . . , fn, to be as well α-strongly convex. It is enough to make such an assumption on the compound function f , while the local functions fi only need to satisfy the smoothness conditions stated in Assumption B above.

heading: 4 Gradient Estimator
text: In this section, we detail our choice of gradient estimators gi(t) used at Step 1 of Algorithm 1. We consider Algorithm 2. For any function F : Rd → R and any point x, the vector g returned by Algorithm 2 is an estimate of∇F (x) based on noisy observations of F at randomized points. The estimator is computed for every node i at each step t, thus giving the vectors g = gi(t) in Algorithm 1. The gradient estimator crucially requires a kernel function K : [−1, 1] → R that allows us to take advantage of possible higher order smoothness properties of f . Specifically, in what follows we assume that∫
uK(u)du = 1, ∫ ujK(u)du = 0, j = 0, 2, 3, . . . , `, and κβ ≡ ∫ |u|β |K(u)|du <∞, (3)
for given β ≥ 2 and ` = bβc. In Polyak and Tsybakov [1990] such kernels can be constructed as weighted sums of Legendre polynomials, in which case κβ ≤ 2 √ 2β with β ≥ 1; see also Appendix A.3 in Bach and Perchet [2016] for a derivation.
The gradient estimator in Algorithm 2 differs from the standard 2d-point Kiefer-Wolfowitz type estimator in that it uses multiplication by a random variable K(r) with a well-chosen kernel K. On the other hand, it is also different from the previous kernel-based estimators in zero-order optimization literature Polyak and Tsybakov [1990], Bach and Perchet [2016], Akhavan et al. [2020] in that it needs 2d function queries per step, whereas those estimators require only one or two queries; see, in particular, Algorithm 1 in Akhavan et al. [2020] for a comparison. At first sight, this seems a big drawback of the estimator proposed here, however we will show below that thanks to this estimator
we achieve both a more efficient optimization procedure and better rate of convergences for the optimization error.
When the estimator in Algorithm 2 is used at the t-th outer step of Algorithm 1, it should be intended as a random variable that depends on the randomization used during the current estimation at the given node, as well as on the randomness of the past iterations, inducing the σ-algebra Ft (see Section 5 for the definition). Bounds for the bias of this estimator conditional on the past and for its second moment play an important role below, in our analysis of the convergence rates. These bounds are presented in the next two lemmas, whose proofs are presented in Appendix B. We state them in the simpler setting of Algorithm 2, with no reference to the filtration (Ft)t∈N. Lemma 1. Let f : Rd → R be a function in Fβ(L), β ≥ 2, and let the random variables ξ1, . . ., ξd and ξ′1, . . ., ξ ′ d be independent of r and satisfy E[|ξj |] < ∞, E[|ξ′j |] < ∞, for j = 1, . . ., d. Let the kernel satisfy conditions (3). If the gradient estimator g of f given by Algorithm 2 then, for all x ∈ Rd,
‖E[g]−∇f(x)‖ ≤ Lκβ √ dhβ−1.
It is straightforward to see that the bound of Lemma 1 holds when the estimators are build recursively during the execution of Algorithm 1 and the expectation is taken conditionally on Ft. This will be used in the proofs.
Lemma 2. Let f : Rd → R be 2-smooth and let maxx∈Θ ‖∇f(x)‖ ≤ G, κ ≡ ∫ K2(u)du <∞. Let the random variables ξ1, . . . , ξd and ξ′1, . . . , ξ ′ d be independent of r and E[ξ 2 j ] ≤ σ2, E[(ξ′j)2] ≤ σ2 for j = 1, . . . , d. If g is defined by Algorithm 2, where x is a random variable with values in Θ independent of r and depending on ξ1, . . . , ξd and ξ′1, . . . , ξ ′ d in an arbitrary way, then
E‖g‖2 ≤ 3dκ 2
( σ2
h2 +
3L̄2 4 h2 ) + 9G2κ.

heading: 5 Noise Model
text: Algorithm 2 is called to compute estimators of gradients of the local functions fi, i = 1, . . . n, at each iteration t of Algorithm 1. Thus, we assume that agent i at iteration t generates a uniform random variable ri(t) ∼ U [−1, 1] and gets 2d noisy observations, defined, for j = 1, . . . , d
yi,j(t) = f(x i(t) + htri(t)ej) + ξi,j(t) y′i,j(t) = f(x i(t) + htri(t)ej) + ξ ′ i,j(t)
where the parameters ht > 0 will be specified later.
In what follows, we denote by Ft the σ-algebra generated by the random variables xi(t), for i = 1, . . . , n. In order to meet the conditions of Lemmas 1 and 2 for each (i, t), we impose the following assumption on the collection of random variables (ri(t), ξi,j(t), ξ′i,j(t)).
Assumption C. For all integers t and i ∈ {1, . . . , n} the following properties hold.
(i) The random variables ri(t) ∼ U [−1, 1] are independent of ξi,1(t), . . . ξi,d(t), ξ′i,1(t), . . . , ξ ′ i,d(t) and from the σ-algebra Ft,
(ii) E[(ξi,j(t))2] ≤ σ2, E[(ξ′i,j(t))2] ≤ σ2 for j = 1, . . . , d, and some σ ≥ 0.
Assumption C is very mild. Indeed, its part (i) occurs as a matter of course since it is unnatural to assume dependence between the random environment noise and artificial random variables ri(t) generated by the agents. We state (i) only for the purpose of formal rigor. Remarkably, we do not assume the noises ξi,j(t) and ξ′i,j(t) to have zero mean. What is more, these variables can be deterministic and no independence between them for different i, j, t is required, so we consider an adversarial environment. Having such a relaxed assumption on the noise is possible because of the multiplication by the zero-mean variable K(r) in Algorithm 2. This and the fact that all components of the vectors are treated separately allows the proofs go through without the zero-mean assumption and under arbitrary dependence between the noises.

heading: 6 Main Results
text: In this section, we provide upper bounds on the performance of the proposed algorithms. Recall that T0 is the number of outer iterations in Algorithm 1. Let T be the total number of times that we observed noisy values of each fi. At each iteration of Algorithm 2 we make 2d queries. Thus, to keep the total budget equal to T we need to make T0 = T/(2d) steps of Algorithm 1 (assuming that T/(2d) is an integer). We compare our results to lower bounds for any algorithm with the total budget of T queries.
For given β ≥ 2, we choose the tuning parameters ηt and h = ht in Algorithms 1 and 2 as
ηt = 2 αt , and ht = t− 1 2β . (4)
Inspection of the proofs in Appendix C shows that these values of ηt and ht lead to the best rates minimizing the bounds. As one can expect, there are two contributions to the bounds, one representing the usual stochastic optimization error, while the second one accounts for the distributed character of the problem. This second contribution to the bounds is driven by the following quantity that we call the mean discrepancy: ∆(t) ≡ n−1 ∑n i=1 E[
∥∥xi(t)− x̄(t)∥∥2]. It plays an important role in our argument and may be of interest by itself, cf. Tang et al. [2019]. The next lemma gives a control of the mean discrepancy. Lemma 3. Let Assumptions A, B, and C hold. Let Θ be a convex compact subset of Rd. Assume that diam(Θ) ≤ K and maxx∈Θ ‖∇f(x)‖ ≤ G. If the updates xi(t), x̄(t) are defined by Algorithm 1, in which the gradient estimators for i-th agent are defined by Algorithm 2 with F = fi, i = 1, . . . , n, and parameters (4) then
∆(t) ≤ A ( ρ
1− ρ
)2 d
α2 t− 2β−1 β , (5)
where A is a constant independent of t, d, α, n, ρ. The explicit value of A can be found in the proof. Proof Sketch. Let V (t) = ∑n i=1 ∥∥xi(t)−x̄(t)∥∥2, and zi(t) = ProjΘ(xi(t)− ηtgi(t))− (xi(t) − ηtg
i(t)). The first step is to show that, due to the definition of the algorithm and Assumptions A on matrix W , we have
V (t+ 1) ≤ ρ2 n∑ i=1 ∥∥xi(t)− x̄(t)− ηt(gi(t)− ḡ(t)) + zi(t)− z̄(t)∥∥2 , (6) where ḡ(t) and z̄(t) denote the averages of gi(t)’s and zi(t)’s over the agents i. From (6), by using the fact that ‖zi(t)‖ ≤ ηt‖gi(t)‖, applying Lemma 1 conditionally on Ft, taking expectations and then applying Lemma 2 we deduce the recursion
∆(t+ 1) ≤ ρ∆(t) +A1 ρ2 1− ρ · d α2 t− 2β−1 β ,
where A1 > 0 is a constant. The initialization of Algorithm 1 is chosen so that ∆(1) = 0. It follows that ∆(t) is bounded by a discrete convolution that can be carefully evaluated leading to (5).
Using Lemma 3 we obtain the following theorem. Theorem 4. Let f be an α-strongly convex function and let the assumptions of Lemma 3 be satisfied. Then for any x ∈ Θ the cumulative regret satisfies
T0∑ t=1 E [ f(x̄(t))− f(x) ] ≤ d α(1− ρ) T 1 β 0 ( B1 + B2ρ2 ) + B3 α(1− ρ) (log(T0) + 1),
where the positive constants Bi are independent of T0, d, α, n, ρ. The explicit values of these constants can be found in the proof. Furthermore, if x∗ is the minimizer of f over Θ the optimization error of the averaged estimator x̂(T0) = 1T0 ∑T0 t=1 x̄(t) satisfies
E[f(x̂(T0))− f(x∗)] ≤ d α(1− ρ) T − β−1β 0
( B1 + B2ρ2 ) +
B3 α(1− ρ) ( log(T0) + 1 T0 ) . (7)
Proof sketch. Note first that, due to the definition of Algorithm 1 and to the properties of matrix W we have x̄(t+1) = x̄(t)−ηtḡ(t)+ z̄(t). This resembles the usual recursion of the gradient algorithm with an additional term z̄(t) = n−1 ∑n i=1 z
i(t), where ‖zi(t)‖ ≤ ηt‖gi(t)‖. Using this bound and α-strong convexity of f , analyzing the recursion in the standard way and taking conditional expectations we obtain that, for any x ∈ Θ,
f(x̄(t))− f(x) ≤ 1 2ηt
E [ at − at+1|Ft ] − αat
2 + 2ηt n n∑ i=1 E [ ∥∥gi(t)∥∥2 |Ft]
+ ∥∥E[ḡ(t)|Ft]−∇f(x̄(t))∥∥ ‖x̄(t)− x‖︸ ︷︷ ︸
Bias1
+ 1 ηt E [ 〈z̄(t), x̄(t)− x〉|Ft ] ︸ ︷︷ ︸
Bias2
, (8)
where at = ‖x̄(t)− x‖2. Here, the term Bias2 is entirely induced by the distributed nature of the problem. Using the properties of Euclidean projection and some algebra, it can be bounded as
Bias2 ≤ 3ηt 2(1− ρ)n n∑ i=1 E [ ∥∥gi(t)∥∥2 |Ft]+ 1− ρ 2nηt n∑ i=1 ∥∥xi(t)− x̄(t)∥∥2 . (9) On the other hand, Bias1 accumulates two contributions, the first due to the gradient approximation (cf. Lemma 1) and the second due to the distributed nature of the problem:
Bias1 ≤ κβL √ dhβ−1t ‖x̄(t)− x‖+ L̄
n n∑ i=1 ∥∥xi(t)− x̄(t)∥∥ ‖x̄(t)− x‖ ≤ ( (κβL)2
α dh
2(β−1) t + αat 4
) +
( L̄tα(1− ρ)
n
n∑ i=1 ∥∥xi(t)− x̄∥∥2 + L̄K2 4tα(1− ρ)
) . (10)
Next, we combine inequalities (8)–(10), take expectations of both sides of the resulting inequality, and use Lemmas 2 and 3 to bound the second moments E [ ∥∥gi(t)∥∥2 ] and the mean discrepancy. The final result is obtained by summing up from t = 1 to t = T0 and recalling that ηt = 2αt , ht = t − 12β .
Due to α-strong convexity of f , Theorem 4 immediately implies a bound on the estimation error E[‖x̂(T0)− x∗‖2]. The bound is of the order of the right-hand side of (7) divided by α. Furthermore, we get the following result about local estimators, which follows from a slight modification of Lemma 3 and Theorem 4. Corollary 5. Let Assumptions A, B, and C hold. Let Θ be a convex compact subset of Rd. Assume that diam(Θ) ≤ K and maxx∈Θ ‖∇f(x)‖ ≤ G. If the updates xi(t) are defined by Algorithm 1, in which the gradient estimators for i-th agent are defined by Algorithm 2 with F = fi, i = 1, . . . , n, and parameters ηt = 4α(t+1) , ht = t
− 12β then the local average estimator x̂i(T0) = 2
T0(T0+1)
∑T0 t=1 tx i(t) satisfies
E[‖x̂i(T0)− x∗‖2] ≤ Cmin { 1, d
α2(1− ρ) T − β−1β 0
( 1 +
nρ2
(1− ρ)T0
)} , i = 1, . . . , n,
where C > 0 is a positive constant independent of T0, d, α, n, ρ.
We now state a corollary of Theorem 4 for an algorithm with total budget of T queries. Assume that T0 = T/(2d) is an integer. As our algorithm makes 2d queries per step the estimator x̂(T/(2d)) uses the total budget of T queries. Combining Theorem 4 with the trivial bound E[f(x̂(T/(2d))−f(x∗)] ≤ GK we get the following result. Corollary 6. Let T ≥ 2d and let the assumptions of Theorem 4 be satisfied. Then we have
E[f(x̂(T/(2d))− f(x∗)] ≤ Cmin { 1, d2−1/β
α(1− ρ) T−
β−1 β } ,
where C > 0 is a positive constant independent of T, d, α, n, ρ.
We now state several important implications of our results.
Remark 1. Previous bounds on zero-order distributed optimization Qu and Li [2018], Yu et al. [2019], Tang et al. [2019] contain a dependency of (1 − ρ)−2 in the "connectivity" parameter ρ. While Theorem 4 covers a more difficult noisy setting, our bound displays a better dependency of (1 − ρ)−1. Since most common values of 1 − ρ are of the order n−2 (or n−1), this represents a substantial gain. Remark 2. The case n = 1, ρ = 0 corresponds to usual (undistributed) zero-order stochastic optimization. Then Corollary 6 gives a bound of order min ( 1, d 2−1/β
α T − β−1β
) . This improves upon
the bound1 min ( 1, d 2
α T − β−1β
) obtained under the same assumptions in Akhavan et al. [2020]. Still
our bound does not match the minimax lower bound established in Akhavan et al. [2020] and equal to min (
max(α, T−1/2+1/β), d√ T , d α T− β−1 β
) . (11)
For α 1 the lower bound (11) scales as min ( 1, dαT − β−1β ) . It has the same behavior in the interesting regime of α not too small (α ≥ T−1/2+1/β) and T ≥ d. Note, however, that the lower bound (11) is obtained for the setting with i.i.d. noise, while our upper bound is valid under adversarial noise. Therefore, it may seem rather surprising that the ratio is only d1−1/β . Remark 3. With the same budget of queries T , the 2d-point method in Algorithm 2 is computationally simpler than the methods with one or two queries per step Polyak and Tsybakov [1990], Bach and Perchet [2016], Akhavan et al. [2020] previously suggested for the same setting. For example, the method in Bach and Perchet [2016], Akhavan et al. [2020] prescribes, at each step t = 1, . . . , T , to generate a random variable uniformly distributed on the unit sphere in Rd. This requires of order d calls of one-dimensional random variable generator. Overall, in T steps, the number of calls is of order dT . For our method with the same budget T , we make of order T0 = T/(2d) steps and at each step we need to call the generator only once in order to get r ∼ U [−1, 1]. Thus, with the same budget of queries, Algorithm 2 needs ∼ 1/d2 less calls of random variable generator than the gradient estimator in Bach and Perchet [2016], Akhavan et al. [2020].
In Appendix E, we present a numerical comparison between the algorithm proposed in this paper and that in Akhavan et al. [2020]. The results confirm our theoretical findings. The algorithm of this paper converges faster and the advantage is more pronounced as d increases.

heading: 7 Improved Bounds for β = 2
text: In this section we provide improved upper bounds for the case β = 2 in Corollary 6, where we relax the dependency over d, from d3/2 to d.
Following the literature on undistributed zero-order optimization, we use a standard 2-point method with elements of the analysis developed in Flaxman et al. [2005], Agarwal et al. [2010], Duchi et al. [2015], Shamir [2013, 2017], Akhavan et al. [2020] among others. Specifically, we define
gi(t) = d
2ht (yi(t)−y′i(t))ζi(t) (12)
where yi(t) = fi(xi(t)+htζi(t))+ξi(t), y′i(t) = fi(x i(t)−htζi(t))+ξ′i(t),
with the random variables ζi(t), 1 ≤ i ≤ n, 1 ≤ t ≤ T , that are i.i.d. uniformly distributed on the unit Euclidean sphere in Rd. We make the following assumption on the noise analogous to Assumption C. Assumption D. For all integers t and all i ∈ {1, . . . , n} the following properties hold.
(i) The random variables ζi(t) are independent of ξi(t), ξ′i(t) and from the σ-algebra Ft,
(ii) E[(ξi(t))2] ≤ σ2, E[(ξ′i(t))2] ≤ σ2 for some σ ≥ 0.
Theorem 7. Let f be an α-strongly convex function. Let Assumptions A, B, and D hold with β = 2. Let Θ be a convex compact subset of Rd, and assume that diam(Θ) ≤ K. Assume that
1The recent work Novitskii and Gasnikov [2021] obtains the same improvement using the gradient estimator of Akhavan et al. [2020]. However, as we explain in Remark 3 that estimator is less appealing from the computational point of view.
maxx∈Θ ‖∇fi(x)‖ ≤ G, for 1 ≤ i ≤ n. Let the updates xi(t), x̄(t) be defined by Algorithm 1, in which the gradient estimator for i-th agent is defined by (12), and ηt = 1αt , ht = ( 3d2σ2 2Lαt+9L2d2 )1/4 .
Then for the estimator x̃(T ) = 1T−bT/2c ∑T t=bT/2c+1 x̄(t) we have
E[f(x̃(T ))− f(x∗)] ≤ B 1− ρ ( d√ αT + d2 αT ) ,
where B > 0 is a constant independent of T, d, α, n, ρ.
The main idea of the proof is to use surrogate functions f̂ it (x), for 1 ≤ i ≤ n, defined, for every x ∈ Rd, as f̂ it (x) = Efi(x + htζ̃), where the expectation with respect to the random vector ζ̃ uniformly distributed on the unit ball Bd = {u ∈ Rd : ‖u‖ ≤ 1}. A result, which can be traced back to Nemirovsky and Yudin [1983] implies the fact that gi(t) is an unbiased estimator of the gradient of the surrogate function f̂ it at x
i(t). Thus, we can consider Algorithm 1 as a gradient descent for the surrogate function. Then replacing fi and f by the surrogate functions with the cost of the order h2t , we can recover the initial problem. This method does not work for β > 2 since the error of approximation by surrogate function becomes of bigger order than the optimal rate T− β−1 β . The results that we implement as tools for this section are given in Appendix D.
Combining Theorem 7 with the obvious bound E[f(x̃(T ))− f(x∗)] ≤ GK we obtain
E[f(x̃(T ))− f(x∗)] ≤ B ′
1− ρ min
( 1,
d√ αT
) , (13)
where B′ > 0 is a constant independent of T, d, α, n, ρ. By comparing this upper bound with the minimax lower bound (11) for β = 2, one can note that (13) is optimal with respect to the parameters T and d when α 1.

heading: 8 Discussion
text: We expand our discussion on previous related work, comparing our results to the state-of-the-art distributed and undistributed zero-order optimization settings, and highlight few key open problems.
Comparison to Zero-Order Distributed Settings Distributed opimization with noisy functions’ queries was considered in detail in Sahu et al. [2018b,a], where the setting differs from ours in some key aspects: the updates are obtained not as in Step 2 of Algorithm 1 but rather via decentralized techniques, matrix W is random, the noise is zero-mean random rather than adversarial, and 2-point gradient estimator is used. Papers Sahu et al. [2018b,a] provide, for β = 2 and β = 3, bounds on E[‖xi(T )− x∗‖2] of the order at least n 3/2
(1−ρ)2T −1/2 and n
3/2
(1−ρ)2T −2/3, respectively, as functions of
n, ρ and T . Their bounds contain uncontrolled terms of the form E[‖xi(k0)− x∗‖2] for some large enough k0 = k0(n, α, d) leaving unclear the resulting rate. Paper Hajinezhad et al. [2019] considers 2-point methods with stochastic queries but assume that the noise is the same for both queries and deal with non-convex optimization. Noisy-free zero-order distributed optimization is studied by Qu and Li [2018], Yu et al. [2019], Tang et al. [2019]. From these, Tang et al. [2019] is the closest to our work as it builds on the updates as at Step 2 of Algorithm 1 (though without projections). The bounds obtained therein are of the order (1− ρ)−2 considered as functions of ρ, although they hold for the larger class of gradient dominant functions. As noted in Remark 1 the bound of Theorem 4 scales only as (1− ρ)−1 and this bound holds true, in particular, for noisy-free setting, which is its special case corresponding to σ = 0. Since most common values of 1− ρ are of the order n−2 (or n−1), this represents a substantial gain. Moreover, Theorem 4 covers a difficult noise setting as we deal with adversarial noise. It is also worthwhile to note that the first-order distributed optimization exhibits much better dependency on ρ since bounds that scale as (1 − ρ)−1/2 can be achieved, see [Duchi et al., 2015, Scaman et al., 2019]. Some of the references mentioned above considered unconstrained optimization while our results deal with constrained optimization. Note that the only difference in the proofs of the upper bounds for constrained and unconstrained cases is in the presence of an additional term proportional to the second moment of the gradient at the update (see Lemma 2.4 in Akhavan et al. [2020] for a similar argument). Since this additional term is bounded independently of ρ the overall dependency on ρ remains the same.
Computational and Statistical Advantage of the Proposed Gradient Estimator As we highlighted in Section 6 the gradient estimator in Algorithm 2 requires 2d function queries. At first sight this seems problematic when the dimension d is high, as they need at least T = 2d queries. However, the lower bounds in Shamir [2013], Akhavan et al. [2020] reported in (11) above indicate that no estimator can achieve nontrivial convergence rate for zero-order optimization when T . d β β−1 . Thus, having the total budget of T d queries is a necessary condition for success of any zero-order stochastic optimization method. Algorithms with one or two queries per step can, of course, be realized for T . d but in this case they do not enjoy any nontrivial error behavior. Moreover, by Remark 3, with the same total budget of queries T , the gradient estimator from Algorithm 2 is computationally more efficient2 than the estimators in Polyak and Tsybakov [1990], Bach and Perchet [2016], Akhavan et al. [2020]. Indeed, with the same budget of queries, it needs 1/d2 less calls of random variable generator than it would be for the gradient estimator in Bach and Perchet [2016], Akhavan et al. [2020]. At the same time, as detailed in Remark 2 the proposed gradient estimator yields better rates on the optimization error. We conclude that the proposed zero-order optimization procedure provides both a computational and statistical improvement over the state-of-the-art methods in Akhavan et al. [2020].
Limitations and Future Work A main problem, which remains open, is to study whether the dependency of (1 − ρ)−1 in the upper bounds in Corollary 6 and Theorem 7 is minimax optimal. Moreover, in the standard (undistributed) setting it remains an open problem to design a zeroorder optimization procedure that meets the minimax lower bound (11) with respect to all problem parameters (T, d, β and α). Further directions of research include the analysis of disturbed zero-order algorithms for larger classes of functions, such as α-gradient dominant ones, as well as extension of our results to stochastic updates or asynchronous activation schemes.
Funding Disclosure
A. Akhavan and M. Pontil were partially supported by SAP SE. The research of A.B. Tsybakov is supported by a grant of the French National Research Agency (ANR), “Investissements d’Avenir” (LabEx Ecodec/ANR-11-LABX-0047).

Review summaries:
This paper studies zero-order optimization for strongly convex functions. The authors propose a novel 2d point zero-order gradient estimator, combing with distributed projected gradient descent, they provide upper bounds that are robust to adversary noises and have better dependence on (1 - \rho). In addition, they provide an improved analysis for 2-Holder functions using 2 point gradient estimator in that the dimension dependence is improved from d^1.5 to d. 

This paper considers distributed optimization with zero-order information and under rather general setting of noise in the evaluation of function values. The authors propose a zero-order projected gradient descent algorithm, in which a novel zero-order gradient estimation scheme is used. An important contribution of this paper is to provide explicit upper bounds for the average cumulative regret and optimization error of the proposed algorithm. The obtained bound is clearly presented and show the advantage of the proposed algorithm over existing ones. 

This paper considers the problem of distributed zeroth-order optimization for a class of strongly convex functions. The authors propose a distributed zeroth-order optimization algorithm that only requires only functional evaluations of the objective functions, subject to a general noise model. The error bounds of the algorithm have also been established and improve upon the previous works. Lower bounds have also been discussed and they are nearly optimal. 

This paper studies a distributed optimization problem where the nodes conform to a given adjacency graph. The prevailing formulation of this problem in the literature is based on first order (gradient based) algorithms. This paper however, following [11, 27, 28, 29, 35, 37] considers a zeroth order optimization problem where only direct observations from the function value are accessible (instead of gradient). The average objective function is assumed to be strongly convex and the local objectives are assumed to be Hölder continuous. A distributed optimization algorithm based on a novel gradient estimator is presented. It is shown that the proposed algorithm is favorable over the existing ones in both convergence rate and computational complexity.