An Optimal Structured Zeroth-order Algorithm for Non-smooth Optimization
Authors: Marco Rando
Cesare Molinari
Lorenzo Rosasco
Silvia Villa
Venue: NeurIPS 2023 poster (2023)
Keywords: nonsmooth optimization;zeroth order optimization;nonsmooth zeroth-order
TL;DR: Zeroth order descent with structured directions for non-smooth optimization

Abstract: Finite-difference methods are a class of algorithms designed to solve black-box optimization problems by approximating a gradient of the target function on a set of directions. In black-box optimization, the non-smooth setting is particularly relevant since, in practice, differentiability and smoothness assumptions cannot be verified. To cope with nonsmoothness, several authors use a smooth approximation of the target function and show that finite difference methods approximate its gradient. Recently, it has been proved that imposing a structure in the directions allows improving performance. However, only the smooth setting was considered. To close this gap, we introduce and analyze O-ZD, the first structured finite-difference algorithm for non-smooth black-box optimization. Our method exploits a smooth approximation of the target function and we prove that it approximates its gradient on a subset of random {\em orthogonal} directions. We analyze the convergence of O-ZD under different assumptions.  
For non-smooth convex functions, we obtain the optimal complexity. In the non-smooth non-convex setting, we characterize the number of iterations needed to bound the expected norm of the smoothed gradient. For smooth functions, our analysis recovers existing results for structured zeroth-order methods for the convex case and extends them to the non-convex setting. We conclude with numerical simulations where assumptions are satisfied, observing that our algorithm has very good practical performances.

An Optimal Structured Zeroth-order Algorithm for Non-smooth Optimization
Marco Rando*
Cesare Molinarit
Silvia Villat
Lorenzo Rosasco *
Abstract
diffe-edifeerenc methods are class of algorithmss designed to solve black-box optimization problems by approximating gradient of the target function set of directions. black-box optimization, the non-smooth setting is particularly relevant since, in practice, differentiability and smoothness assumptions cannot be verified. To cope with nonsmoothness, several authors use a smooth approximation of the target function and show that finite difference methods approximate its gradient. Recently, has been proved that imposing a structure in the directions allows improving performance. However, only the smooth setting considered. To close this gap, we introduce and analyze O-ZD, the first structured finite- difference algorithm for non-smooth black-box optimization. Our method exploits a smooth approximation of the target function and we prove that approximates its gradient subset of random orthogonal directions. We analyze the convergence O-ZD under different assumptions. For non-smooth convex functions, we obtain the optimal complexity. In the non-smooth non-convex setting, characterize the number of iterations needed to bound the expected norm of the smoothed gradient. For smooth functions, zeroth-order methods for the convex case and extends them to the non-convex setting. We conclude with numerical simulations where assumptions are satisfied, observing that our algorithm has very good practical performances.
analysis recovers existing results for structured
Introduction
Black-box optimization problems are a class of problems for which only values of the target function are available and no first-order information is provided. Typically, these problems arise when the evaluation of the objective function is based on a simulation and no analytical form of the gradient
accessible or its explicit calculation is too expensive [14
36 35].
To face these problems, different methods that do require first-order information have been proposed see for instance [39 47
30 33 and references therein. These techniques are called derivative-free methods and a wide class of these is the class of finite-difference algorithms
19 26 20
31 39 17]. These iterative procedures mimic first-orde optimization strategies by replacing the gradient of the objective function with an approximation built through finite differences in random directions.
Two types of fifte-diifeernnc methods can be identified: unstructured and structured ones, depending on the way in which the random directions are generated. In the former, directions are sampled i.i.d. while in the latter, directions have satisfy some structural Several authors
from some distribution 39 constraints, e.g. orthogonality
17
10
46
32
42
(32
42 5 theoretically and empirically
Ga-.GGaDIIRRS.
University
of Genova,
IT
(marco. rrwddoeeddd unige itt
lorenzo. rossasoo
1 MaL.Ga DIMA, University Genova, Italy (molinari indiiiiiiddmaa...m unige. it, silvia. villa@unige..s it). Istituto Italiano di Tecnologia, Genova, Italy and CBMM MIT, Cambridge, MA, USA
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

observed that imposing orthogonality among the directions provides better performance than using unstructured directions. Intuitively, imposing orthogonality allows us to avoid cases in which the gradient approximation is built using similar or redundant directions.
Notably, methods based on structured finite differences have been analyzed only for smooth and convex (or specific non-convex) target functions. This represents a strong limitation, since smoothness and convexity are in practice hardly verifiable, due to the nature of black-box optimization problems. The aim of this paper is the analysis of structured finite difference algorithms dropping smoothness and convexity assumptions.
For unstructured ducing smooth approximation of the target function, also known "smoothing" [6 and using as a surrogate of the target. Although the gradient of the smoothing is not computable, different authors showed that for certain types of smoothing the unstructured finite-difference approximation provides an unbiased estimation of such a gradient see, for example, (39, 46 18] 22 21 24]. This
methods, common way to face nonsmoothness consists in intro-
key observation allows prove that unstructured finite-difference methods approximate solution in
different nonsmooth settings [17 39, 46
18
22 34].
For structured finteeddifference methods, the analysis in the nonsmooth setting is not available. A key step which is missing is the proof of the fact that the surrogate of the gradient built with structured directions is estimation of the gradient of suitable smoothing.
In this paper, we close the gap, and introduce O-ZD, a structured finite-difference algorithm in the non-smooth setting. The algorithm builds an approximation of the gradient of smoothed target using a set of & S d orthogonal directions. We analyze the procedure proving that our finite-difference surrogate is an unbiased estimation of the gradient of smoothing. We provide convergence rates for non-smooth convex functions with optimal dependence on the dimension
[17] and rates for the non-smooth non-convex setting (for which lower bounds are not known 34). Moreover, for non-smooth convex functions we provide the first proof of convergence of the iterates for structured finite differences in this setting. For smooth convex functions, we recover the standard results for structured zeroth-order methods [32 42). We conclude with numerical illustrations. To the best of our knowledge, this is the first work on nonsmooth lnite-difference method with structured directions.
The paper is organized as follows. Section 2 we introduce the problem and the algorithm. In Section 3 we state and discuss the main results. In Section 4 we provide some experiments and in some final remarks.
Problem Setting & Algorithm
Given a function f Rdd + I and assuming that f has at least a minimizer in IRd, we consider the problem to find
x² E arg min f(z).
DER
In particular, we consider the non-smooth setting where f might be non-differentiable. To solve problem ① we propose zeroth-order algorithm, namely iterative procedure that uses only function values. At every iteration E N, a first-order information of f is approximated with
finite-differences using set of l < d random orthogonal directions (pk ())ge ) Such orthogonal directions are represented as rotations (and reflections) of the first vectors the canonical basis
(t) (e:) We set Pk Ggee, where Gş belongs to the orthogonal group defined as
(dd) [G e [x det G + G GT}.
Methods for generating orthogonal matrices are discussed in Appendix D Given G E O()), 0 << E and f we consider the following central finite-difference surrogate of the gradient information
(x) d f(z + +hGG,) - f(x hGes) Ge
(2)
2h
Then, introduce the following algorithm.

---

Algorithm O-ZD: Orthogonal Zeroth-order Descent
Input: Zo E R (ou)KeN C R4, (hE)KEN C R4, &EN.t. 15 3e
for k 1, do
sample G1 i.i.d. from O(d) 2k++1 2i h2 end for
Starting from an initial guess To E Rd at every iteration k E N. the algorithm samples an orthog- onal matrix G i.i.d. from the orthogonal group O(d) and computes a surrogate I(G,, of the gradient at the current iterate Then, it computes r&+1 by moving in the opposite direction of G(G. (xx). This approach belongs to the class of structured finte-ddifference methods, where bunch of orthogonal directions is used to approximate a gradient of the target function [32 42 Such algorithms have been proposed and analyzed only for smooth functions. To cope with this limitation, and extend the analysis to the nonsmooth setting, we exploit smoothing technique. For a fixed probability measure p on Id and a positive parameter h > 0 called smoothing parameter, we define the following smooth surrogate of? f
fn.p(x) := f(x + ha) dp(u).
As shown in (6. Sh.p is differentiable even when f' is not. In the literature, different authors have used this strategy to face non-smooth zeroth-order optimization with random finite-difference methods
18 39 49 (2) fixing specific smoothing distribution, but one applied and analyze it for structured methods.
In this work, is the uniform distribution over the lg unit ball Bd, defining the smooth surrogate
1 = [
fz(2))
f(x + hu) du,
vol(Ba)
where VBI( denotes the volume of idd One of our main contributions is the following Lemma which shows that the surrogate proposed in (2) is unbiased estimator of the gradient of the smoothing in (4).
Lemma 1 (Smoothing Lemma). Given a probability space (Ω, F, P), ler G Old) be a random variable where ((d) is the orthogonal group endowed with the Borel a-algebra. Assume that the probability distribution of G is the (normalized) Haar measure. Let h > 0 and let g be the gradient
surrogate defined in eq. (2). Then,
(Vz E Rd)
EclgGN))<>)] = (n()))
where fn is the smooth approximation of the target function f defined in eq.
The proof of Lemma T is provided in Appendix A.
Remark 1. Note that Lemma holds also using a or gaN defined as
f(x + hGe,) f(x) Ges (x) = d 7
f(x + hGee)
and
(. (x)
Ge
h
1=1
Bally (x)). Despite these two estimators being com- putationally cheaper than the proposed one, we use central finite differences since they allow to derive a better bound for BGllo(c.))) as observed [46] for the case 1
since Ealaa =
Thanks to Lemma [ we can interpret each step of Algorithm| D as a Stochastic Gradient Descent (SGD) the smoothed function fhg. But the analysis of the proposed algorithm does not follow from the SGD one for two reasons. First, the smoothing parameter h2 changes along the iterations; second (and more importantly), the set of minimizers of fr. and f are different in general. However, we will take advantage of the fact that fa can be seen as an approximation of f. The relationship between ff and its approximation fn depends the properties of f see Proposition T in Appendix
Our main contributions are the theoretical and numerical analysis of Algorithm I under different choices of the free parameters on, he, namely the stepsize and the smoothing parameter. To the best of our knowledge, Algorithm [ ii the first structured zeroth-order method for non- smothh optimization.
3

---

2.1 Related Work
In practice, the advantage of the use of structured directions in finite-difference methods has been observed in several applications [12] and motivated their study. In [.) the authors theoretically and empirically observed that structured directions provide better approximation of the gradient with respect to unstructured (Gaussian and spherical) ones. Next, we review the most related works, showing the differences with our algorithm.
Unstructured Finite-differences. Most of existing works focused on the theoretical analysis of methods using single direction to approximate the gradient see e.g. [39, 46
22
45
18
24
34
The main results existing so far analyze the convergence of the function values in expectation. They provide convergence rates in terms of the number of function evaluations and explicitly characterize the dependence on the dimension of the ambient space.
10].
Smooth setting: In (39 [7]. a finite difference algorithm with single direction is analyzed. Rates on the expected function values are shown. In [39] the dependence on the dimension in the rates is not optimal. In [[7]. both single and multiple direction cases are analyzed and lower bounds are derived. However, only the convex setting is considered. In [24], both convex and non-convex settings are analyzed. They obtain optimal dependence dimension in the complexity for convex functions. However, only the single-direction case is considered and only the smooth setting is analyzed. Comparing the result with rates, Algorithm [ achieves the same dependence on the dimension in the complexity taking as a fraction of d. Note that, by parallelizing the computation of the function evaluations, we can get a better result. Non-smooth setting: In [39 [7 the non-smooth setting is also analyzed. However, only the single direction case has been considered and both algorithms do not achieve the lower bound. More precisely. for convex functions, a complexity of O(de is achieved by [39] and (oglo by [17 while our algorithm gets the optimal dependence. Moreover, note that in [17 the strategy adopted (also called "double smoothing") requires tuning one more sequence of parameters, which challenging problem practice, and only the convex setting considered. non-convex setting is analyzed by bounding the expected norm of the smoothed gradient. However, they obtain a complexity of (dp while our algorithm gets a better dependence on dimension. In [46), the optimal complexity is obtained for convex functions. However, the author does not analyze the non-convex setting. Moreover, note that, despite the complexity in terms of function evaluations being the same, our algorithm gets a better complexity in terms of the number of iterations since uses multiple directions (and this is an advantage we can parallelize the function evaluations). Furthermore, note that the method proposed in [46] can be seen the special case of Algorithm [ with & = 1. In (34] the single direction case is analyzed only in the non-convex setting. The dependence the dimension of the complexity in the number of function evaluations achieved matches our result in this setting (again, in the number of iterations our method obtains a better dependence).
39 also the
authors analyze structured finite differences in both deterministic and stochastic settings. However, only the smooth convex setting is considered. In [12 orthogonal matrices used to build directions but analysis is provided. In (25). finite-iffferens with coordinate directions is analyzed. At every iteration, d + function evaluations are performed to compute the estimator and only the smooth setting is considered.
Structured finit-difffrrnne. In 332 42
Main Results
In this section, we analyze Algorithm considering both non-smooth and smooth problems. We present the rates obtained by our method for convex and non-convex settings and compare them with those obtained by ofate-of-the-art methods. Proofs are provided in Appendix B In the following, we call complexity in the number of iterations / function evaluations, respectively, the number of iterations/ function evaluations required to achieve an accuracy E e (0, 1).
3.1 Non-smooth Convex Setting
In this section, we provide the main results for non-smooth convex functions. In particular, we will assume that the target function is convex and satisfy the following hypothesis.

---

Assumption (Lo-Lipschitz continuous). The function f is Lo-Lipschitz continuous; i.e., for some Lo
(Vx, E 1R
(f(x) f(vil & Lol - vll-
Note that this assumption implies that also Jn is Lo-Lipschitz continuous see Proposition Moreover, to analyze the algorithm, we will consider the following parameter setting.
Assumption (Zeroth-oree non non-smooth convergence conditions). The step-size sequence ay and the sequence of smoothing parameters h& satisfy the following conditions:
ag & C'
a E and anhy E Ct
The assumption above is required to guarantee convergence to a solution. In particular, the first two conditions common for subgradient method and stochastic gradient descent, while the third condition was already used in structured zeroth-order methods [32] 42 and links the decay of the smoothing parameter with the stepsize's s An example of ork, hg that satisfy Assumption 2 is Ogk k and h4 -P with 0 E (1/2, 1) and s.t. + 1.
We state now the main theorem for non- non-smooth convex functions.
Theorem (Non-smooth convex). Under Assumption assume that f is convex and let (Ik)NEN
be sequence generated by Algorithm 1 Ex agui/AN. Ak. Then
For every l N. denote Ak EE= ai and set Is
|1:00
x=1
dLg
[Tf( min fl s Sk/ Ax
with
S :=
Σ
a
Lo
oghi,
2
i=0
t=0
where > is a constant independent of the dimension and x² is any solution in arg min f. Moreover, under Assumption we have
lim f(xx) = min f a.s.
and that there exists a random variable taking values in arg min f s.t. 2% > a.s.
In the next corollary, we derive explicit rates for specific choices of the parameters. Corollary 1. Under the assumptions of Thehormm let x E arg min f. Then, the following hold:
(i) Let 0 (1/2, 1) and p E such that p 1. For every N, let ag a(k + 1) and h4 = h(k + 1)-P with a > O and h > 0. Then
C
1
[[(( min f akk + 0
k11-0)
for some constant C. provided the proof.
(ii) For every le N, let a4 and h2 with a, h > Then
x1
cdl2
Bff())) min fl 2
+
e Loh,
where is constant independent of the dimension.
(iii) Fix an accuracy E E (0,1) and let K 8(cT2) ===)a=m0==())(()===="")) Set = / ||00='1| and hy s for every K. Then
V2nk
Ef min fl S
and the complexity in terms of number of iterations is ((d/f)==
Discussion. The bound in Theorem 1 depends the initialization and additional quantity that can be interpreted as an approximation error. The latter is composed of two parts. The first one is generated by the approximation of the gradient of the smoothed function; while the second, involving lg. generated by the smoothing. Since the rate depends on 1/ / Ois, we would like to choose the stepsize large as possible. However, to get convergence, need make the approximation errors

---

vanish sufficiently fast. To guarantee this, as can observe from Theorem we need to impose some conditions the stepsize ag and on the smoothing parameter h2 (i.e. Assumption 2), that will slow down the decay of the first term. In Corollary we provide two choices of parameters: the first one satisfies Assumption and ensures convergence; the second one corresponds to constant stepsize and smoothing parameter. For the first choice, we recover the rate of the subgradient method in terms of k. In particular, for 0 approaching 1/2, the convergence rate is arbitrarily close to O(k-1/2) and similar the one derived in [77 Theorem 2]. The dependence on the dimension depends on the choice of the constant a. The optimal dependence is obtained with the choice a vt/dd. Indeed, in that case, the complexity in the number of iterations of the order ((d/()= 2), which is better than the one achieved by [17] Theorem 2] and [39]. Note that also 46 Corollary 1] and [22] Theorem 2.4] obtain worse complexity terms of the number of iterations (since they use a single direction), but the same complexity in the number of function evaluations. Clearly, since multiple directions are used, single iteration O-ZD will be more expensive than one iteration of [46] 22 However, our algorithm is more efficient if the i' function evaluations required at each iteration can be parallelized. On the more theoretical side, we observe that the advantage of multiple orthogonal directions is in the tighter bounds for the variance of the estimator, namely for Bl]l ( (zx)
see [46] Lemma
5] and Lemma
3.2 Non-smooth Non-convex Setting
To analyze the non-convex setting, following 139. we provide bound the averaged expected square norm of the gradient of the smoothed target. In particular, we use the following notation:
(h) Mi
a E (())))1
/AR.
where A
ag.
Next, we state the main theorem for the non-convex non-smooth setting. Theorem (Non-smooth non-convex). Under Assumption [ let (Ik)ken be sequence generated by Algorithm ! with, for every & E N, h2 = for some h > O. Then
St/AA with S : fa( min f + Lidvd Uk
(b)
h
In the next corollary, we derive the rates for specific choices of the parameters. Corollary 2. Under the assumptions of Theorem the following statements hold. (i) aa a(k + 1)-0 with a > and E (1/2, 1). then
mh (b) < C Filzo) - min f
aki-6
k:1-6
where C is constant independent of the dimension.
(ii) Ifax with > for every k E N. then
The = [n(ro)- - min
S min
cLādva
th a,
ak
where is constant independent of the dimension.
(iii) Let E (0,1), let K 4(fn(xu) min f) de-2/(Eh)) and choose a
V (Ullaa) (aa((aa)nnnnaaaa fjch Then have that (h) S Thus, the number of function evalua- Keljavd tions required to get precision vi < s is of the order Odd
Relying on the results in (34 we show that this is related to a precise notion of approximate stationarity. To do so, we need to introduce a definition of subdifferential which is suitable to this setting. shown in 34 the Clarke subdifferential not the right notion, and the approximate Goldstein subdifferential should be used instead.

---

Definition 1 (Goldstein subdifferential and stationary point). Under Assumption 1 let a E R and O. The h-Goldstein subdifferential of f at fzf(z) == conv(U ff()) where of is the Clarke subdisferential [73] and B8, (x) is the ball centered in x with radius h. For E (0, 1), point E TR& is (h, Gooossetii stationary point for the Sunction J min و E 0nf(2)) S Corollary 3. Under the same assumptions of Theorem fix K N and let I be random variable taking values in {o................... K --11 } such that, for all i, P a] = os/Ak-1. Let also Sy be defined as in Theorem Then
E, [min 17 E (f(((r))) - Sk/Ak.
In the setting of Corollary (iii, we have also that E, [mn (Inall E 0fu(zi))] S..
Discussion. In Theorem we fix the smoothing of the target, i.e. we consider h&g constant, and we analyze the non-smooth non-convex setting providing a rate on the expected norm of the smoothed gradient. The resulting bound is composed of two parts. The first part is very natural, and due to the functional value at the initialization. The second part is the approximation error. Recall that Assumption holds, and therefore fn < f + Loh due to Proposition I This suggests taking h as small as possible in order to reduce the gap between fn and f. However, taking h too small would make the approximation error very big. In our analysis, we consider the case with h constant. Moreover, as for the convex case, the speed of the rate depends on Ag and so we would like to take the stepsize large as possible. But to control the approximation error, need assume of < L!. In Corollary 2 we consider two choices of stepsize. The first choice satisfies the property of og E él while the second one analyzes the case of constant step-size. Comparing our rate to the one in [39] we see that we obtain a better dependence iterations and function evaluations. Our results match the one of (34 Theorem 3.2] in terms rate and in terms of function evaluations. We get better dependence on the dimension in the number of iterations. Note again that, despite the complexity in terms of the number of function evaluations being the same, the possibility of parallelization for the function evaluations yields better result for our method. As for the convex setting, we have tighter upper-bound on the variance of the estimator of the smoothed gradient with respect to the dimension has been used to assess the approximate stationarity for first-order methods as well, see [15]. The latter work shows that cutting plane algorithm achieves rate of Ode-)) for Lipschitz functions.
the dimension in the complexity, both in terms of
34 Lemma D. Goldstein stationarity
3.3 Smooth Convex setting
We consider now the smooth setting, i.e. we assume that the target function satisfies the following hypothesis.
Assumption 3 (Lo-Smooth), The function f is L1 smooth; i.e. the function is differentiable and, for some L > 0,
(Vz,yeR")
[Vf(2) - f s 1110 - y|l.
This is the standard assumption for analyzing first-order methods and has been used in many other works in the literature for zeroth-order algorithms- see e.g. [7]. As shown in previous works, if f satisfies Assumption 3 then also ʃu satisfies it see Proposition [ We will consider also the following assumptions on the stepsize and the smoothing in order to guarantee convergence.
Assumption (Smooth zeroth-order convergence conditions). The stepsize sequence (ak)ReN and the smoothing sequence (hy) satisfy the following conditions:
Ogk & gl and aghy E El'.
Moreover, Og S ā < E/dL, for every N.
Note that this a weaker version of Assumption 2 Next, we state the main theorem for convex smooth functions.
Theorem (Smooth convex). Under Assumptions 3 and let (#K)keN be a sequence generated by
Algorithm and x¹ E arg min f. For every le E set A Ete ai and Ts = Σ
a,5.// Ag. Then,
for every EE N.
(+x( min fl <
Dz A
with Dr
LA + dā
Sk +
Σ ps ( + ++:)) i=0
2PA
7

---

where S = ||20 x112+
1.3d 2
a??h, Pk Lidayhu and A (t+ {a)
Corollary 4. Under the same Assumptions of Theorem 3 the following hold.
(i) for every & E IN we set a₁ = =a =>>0 and hyg h(k + 1) for h > 0 and > then
C
(I((E)) min S|
okk'
where C a constant provided in the proof. Moreover, if < f/(2dL1). lim [( =
k+200
min f a.s. and there exists random variable & taking values in arg min f s.I. IR + a.s.
>>0 =>> (ii) for every & E N we set as = =a
and < h& S h, then
C1 Czah + Cya Caa2
B[f(+)) - min fl S k
where C1 C2 C3 and Cs are non-negative constants.
Discussion. the previous cases, the bound in Theorem composed by two terms: the error due to the initialization and the one due to the approximation. An important difference with the results in the non-smooth setting is that every term in the approximation error decreasing with respect to the smoothing parameter f4. This allows obtaining convergence also with the constant step-size scheme, taking h& E C². In Corollary (i), we recover the result [2] Theorem 5,4] with specific choice of parameters a, f (up to constants). The complexity depends on the choice of Note that by Assumption 4 a < {/[]d) thus the dependence on the dimension in the rate will be at least d/í. In particular, taking a f/(2d[1). terms of function evaluations. This result has a better dependence on the dimension than (3). In Corollary 4 (ii), the dependence
obtain the optimal complexity of (de () in
the dimension in the complexity depends the choice of and h. Moreover, the rate obtained is equal (up to constants) to the rate obtained [32 in the same setting, i.e. ( /k) (in which we hide the dependence d and f). As for (32], for the first setting we can prove the almost sure convergence of the iterates.
3.4 Smooth Non-Convex setting
To analyze the smooth non-convex setting, we introduce the following notation:
(Ak E A
oz E(IV)())))))) /As-
Oki
TJk: =
1=0
Note that, in comparison with the quantity defined in Section 3.2 here 7i is related the exact objective function f and not to its smoothed version fn- Next, we state the main result for smooth non-convex functions.
Theorem 4 (Smooth non-convex). Suppose that Assumption 3 holds and assume that, for every k a á (/(2dL.). Let (Ix) be a sequence generated by Algorithm 7 Then
Th AA
((x0) min
Li 8
Σ
ogh
Lid
a
Lid
48
Corollary 5. Under the assumptions of Theorem the following hold.
(i) af S â and h4 hk-0 with h and 1, then
f(x)) min f h2
C2ah2d²
m <
+ Al
Aa
A
k'
where C1 and C.2 constants provided in the proof.
(ii) Ifaz S ã and hy = h > O, then
f(xo) min ff Cidh22
Czah2
mh
Aak
A
Al
where C₁ and C2 constants provided in the proof.
8

---

Discussion. As for the convex case, every term in the approximation enror depends on the smoothing
parameter hg- In Corollary 5 (i), we take constant step-size and hg E e.. With this choice of parameters, we get rate of On1 /k) which matches with the result obtained by B9 The dependence on the dimension depends on the choice of a and h. Note that a < </(2aaa) thus taking h = O /d).
in the rate we get a dependence on the dimension of d/e. Taking for instance = f/(3d.,) and
(1/d), we get a complexity of ()d=-) in terms of function evaluations.
4 Numerical Results
In this section, we provide some numerical experiments to assess the performances of our algorithm. We consider two target functions: a convex smooth one and a convex non-smooth one. Details on target functions and parameters of the algorithms reported in Appendix C To report our findings, we run the experiments 10 times and provide the mean and standard deviation of the results.
How to choose the number of directions? In these experiments, we set a fixed budget of 4000 function evaluations and we consider d 50. We investigate how the performance of Algorithm T changes as the value of increases. In Figure each function evaluation. If f > 1, then the target function values are repeated 2 times, since we need to perform 2 function evaluations to do one iteration. For a sufficiently large budget, increasing the number of directions l leads to better results compared to using single direction both smooth and non- mon-smooth settings.
1 we observe the mean sequence / [xxx) f(x") after
Smooth function
NonSmooth Convex
P1u) 0x 3
10
function wosluations
function evsiuations
Figure 1: From left to right, function values per function evaluation in optimizing smooth and non-smooth target functions with different numbers of directions.
Comparison with finite-differrnce methods. Now, compare Algorithm with other finite- difference methods. More precisely, we consider finite differences with single (and multiple) Gaussian (and spherical) directions. The budget of function evaluations is 1000 and the ambient dimension is = 10. For multiple direction methods, we fix the number of directions = d. Further experiments provided in Appendix F
Smooth Convex Target
NonSmooth Convex Target
Sogg Saher
Maz Sphetial
100
100
evaluations
evolustions
Figure 2: From left to right, function values per function evaluation in optimizing smooth and non-smooth convex functions with different difite-difference algorithms.

---

In Figure 2 we plot the sequence f(xx) f(=") with respect to the number of function evaluations. While terms of rates and complexity the different algorithms are the same, Algorithm shows better performances than random directions approaches, and we believe this is due to the use of structured (i.e. orthogonal) directions. Indeed, orthogonal directions yield better approximation of first-order information with respect to other methods. The practical advantages of structured directions were already observed [32 42 [ [2] and these experiments confirm that the good practical behavior holds even in the nonsmooth setting.
Conclusion
We introduced and analyzed O-ZD zeroth-order algorithm for non-smooth zeroth-order optimization. We analyzed the algorithm and derived rates for non-smooth and smooth functions. This work opens different research directions. An interesting would be the introduction of learning procedure for the orthogonal directions. Such an approach could have significant practical applications.
Acknowledgments and Disclosure of Funding
This project has been supported by the TraDE-OPT project, which received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Cun grant agreement No 137. L. R. and M. R. acknowledge the financial support of the European Research Council (grant SLING 819789), AFOSR projects FA9550-18 1-7000 (European Office of Aerospace Research and Development), the EU H2020-MSCA-RISE. project NoMADS DLV- 777826, and the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216. S. V. and L. R. acknowledge the support of the AFOSR project FA8655-22-1- 1-7-44 The research by S. V. and C. M. has been supported by the MIUR Excellence Department Project awarded to Dipartimento di Matematica, Università di Genova, CUP D33C23001 and C. M. are members of the Gruppo Nazionale per I' Analisii Matematica, 100000 10001. S. V.
Probabilità e le loro Applicazioni (GNAMPA) of the Istituto Nazionale di Alta Matematica (INdAM). This work represents only the view of the authors. The European Commission and the other organizations are responsible for any use that may be made of the information it contains.
References
[1] Hans Wilhelm Alt. Linear Functional AunctionaAnnalysiss An Applicaatonn-OOiinteed Introduction. Springer London, London, 2016.
[2] T. W. Anderson, I. Olkin, and L. G. Underhill Generation of random orthogonal matrices. SIAM Journal on Scientific and Statistical Computing. 84)::625-629, 1987.
[3] J.B. Baillon and G. Haddad. Quelques propriétés des opérateurs boréés etn-cycliquement monotones. Israel Journal of Mathematics, 26(2):137-150, Jun 1977.
[4] A. Barvinok. Approximating orthogonal matrices by permutation matrices. Pure and applied mathematics quarterly, 2, 11 2005.
[5] A. Berahas, Cao, K. Choromanski, and Scheinberg. theoretical and empirical compar- ison of gradient approximations in derivative-free optimization. Foundations of Computational Mathematics, 22(2):507-560, Apr 2022.
[6] D. P. Bertsekas. Stochastic optimization problems with nondifferentiable cost functionals. Journal of Optimization Theory and Applications, 12(2):218-22 2-11, Aug 1973.
[7] Á. Björck. Numerics of gram-schmidt orthogonalization. Linear Algebra and its Applications, 197-198:297-316. 1994.
[8] C. Boutsidis and A. Gittens. Improved matrix algorithms via the subsampled randomized hadamard transform. SIAM Journal on Matrix Analysis and Applications, 3403)1301111340, 2013.
10

---

[9] Lars Buitinck, Gilles Louppe, Mathieu Blondel, Fabian Pedregosa, Andreas Mueller, Olivier Grisel, Vlad Niculae, Peter Prettenhofer, Alexandre Gramfort, Jaques Grobler, Robert Layton. Jake VanderPlas, Arnaud Joly, Brian Holt, and Gaël Varoquaux. API design for machine learning software: experiences from the scikit-leam project. ECML. PKDD Workshop: Languages for Data Mining and Machine Learning, pages 108-122, 2013.
[10] R. Chen and S. Wild. Randomized derivative-free optimization of noisy convex functions. arXiv preprint arXiv: 2015.
[11] K. Choromanski, M. Rowland, W. Chen, and A. Weller. Unifying orthogonal monte carlo methods. In International Conference Machine Learning, pages
PMLR, 2019.
[12] K. Choromanski, M. Rowland, V. Sindhwani, R. Turner, and A. Weller. Structured evolution with compact architectures for scalable policy optimization. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 970-978, PMLR, 10-15 Jul 2018.
[13] F. H. Clarke. Optimization and Nonsmooth Analysis. Society for Industrial and Applied Mathematics, 1990.
[14] A. R. Conn, K. Scheinberg, and L. N. Vicente. Introduction Derivative-Free Optimization. Society for Industrial and Applied Mathematics, 2009.
[15] D. Davis, D. Drusvyatskiy, Y. T. Lee, S. Padmanabhan, and G. Ye. A gradient sampling method with complexity guarantees for lipschitz functions high and low dimensions. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 6692-6703. Curran Associates, Inc., 2022.
[16] J. C. Duchi, P. L. Bartlett, L. Peter, and M. J. Wainwright. Randomized smoothing for stochastic optimization. SIAM Journal on Optimization, 22(2):674-701, 2012.
[17] J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono. Optimal rates for zero-order convex optimization: The power of two function evaluations. IEEE Transactions on Information Theory, 61 (5):78888-28 2015.
[18] A. Flaxman, A. Tauman Kalai, and B. McMahan. Online convex optimization in the bandit setting: Gradient descent without gradient. In SODA '05 Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithms,, pages 385-394, January 2005.
[19] M. Fornasier, T. Klock, and K. Riedl. Consensus-hased optimization methods converge globally mean-field law. arXiv 2103.15133 2021.
[20]
Frazier. Bayesian Optimization, chapter 11, pages 255-278. INFORMS, 2018
[21] X. Gao, B. Jiang, and S. Zhang. On the information-adaative variants of the admm: An iteration complexity perspective. Journal of Scientific Computing, 76(1):327-363, Jul 2018.
[22] A. Gasnikov, A. Novitskii, V. Novitskii, F. Abdukhakimov, D. Kamzolov, A. Beznosikov, M. Takac, P. Dvurechensky, and B. Gu. The power of first-order smooth optimization for black-box non-smooth problems. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Confer- ence on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 7241-7265. PMLR, 17-23 Jul 2022.
[23] A. Genz. Methods for generating random orthogonal matrices. In Monte-Cario and Quasi- Monte Carlo Methods 1998: Proceedings of a Conference held at the Claremont Graduate University, Claremont, California, USA, June 22-26. 1998, pages 199-213, Springer, 2000.
[24] S. Ghadimi and G. Lan. Stochastic first- and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, 23(4):2341 2013.
[25] G. N. Grapiglia. Worst-case evaluation complexity derivative-free quadratic regularization method. Optimization Letters, Feb 2023.
11

---

[26] N. Hansen. The CMA Evolution Strategy: A Comparing Review, pages 75-102. Springer Berlin Heidelberg, Berlin, Heidelberg, 2006.
[27] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk,
M. Brett, A. Haldane, J. Fernández del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, Reddy, W. Weekesser. H. Abbasi, C. Gohlke, and T. Oliphant. Array programming with NumPy. Nature, 585(7825):357-362 September 2020.
[28] A. Hedayat and W. D. Wallis. Hadamard matrices and their applications. The Annals of Statistics, 6(6): 14-1238, 1978.
[29] J. D. Hunter. Matplotlib: A 2d graphics environment. Computing in Science & Engineering, 933:900-55, 2007.
[30] N. K. Jain, U. Nangia, and J. Jain. A review of particle swarm optimization. Journal of The Institution of Engineers (India): Series B, 99(4):407-411 Aug 2018.
[31] J. Kiefer and J. Wolfowitz. Stochastic estimation of the maximum of a regression function. Annals of Mathematical Statistics, 23:462-466, 1952.
[32] D. Kozak, C. Molinari, L. Rosasco, L. Tenorio, and S. Villa. Zeroth order optimization with orthogonal random directions, 2021.
[33] R. M. Lewis, V. Torezon, and M. W. Trosset. Direct search methods: then and now. Journal of Computational and Applied Mathematics, 124(1):191-22)) 207, 2000. Numerical Analysis 2000. Vol. IV: Optimization and Nonlinear Equations.
[34] T. Lin, Z. Zheng, and M. Jordan. Gradient-free methods for deterministic and stochastic nonsmooth nonconvex optimization. S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 26160-26175. Curran Associates, Inc., 2022.
[35] S. Liu, P. Y. Chen, B. Kailkhura, G. Zhang, A. O. Hero III, and P. K. Varshney. A primer zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications. IEEE Signal Processing Magazine, 37(5):43 54, 2020.
[36] H. Mania, A. Guy, and B. Recht. Simple random search of static linear policies is competitive for reinforcement learning. S. Bengio, H. Wallach, H. Larochelle, Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018
[37] P. Mattila. Geometry of Sets and Measures in Euclidean Spaces: Fractals and Rectifiability Cambridge Studies in Advanced Mathematics. Cambridge University Press, 1995.
[38] F. Mezzadri. How to generate random matrices from the classical compact groups. arXiv preprint math-ph/0609050, 2006.
[39] Y. Nesterov and V. Spokoiny. Random gradient-free minimization of convex functions. Founda- tions of Computational Mathematics, 17:527-566, 2017.
[40] Z. Opial. Weak convergence of the sequence of successive approximations for nonexpansive mappings. Bulletin of the American Mathematical Society, 73(4):591 - -997, 1967.
[41] B. T. Polyak. Introduction optimization. Optimization Software Inc., Publications Division, New York, 1:32, 1987.
[42] M. Rando, C. Molinari, S. Villa, and L. Rosasco. Stochastic zeroth order descent with structured directions, 2022.
[43] H. Robbins and D. Siegmund. A convergence theorem for non negative almost supermartingales and some applications. J. Rustagi, editor, Optimizing Methods Statistics, pages 233-257. Academic Press, 1971.
12

---

[44] C. Rusu and L. Rosasco. Fast approximation of orthogonal matrices and application to pea. Signal Processing, 194:10845,, 2022.
[45] T. Salimans, J. Ho, X. Chen, Sidor, and Sutskever. Evolution strategies as a scalable alternative to reinforcement learning, 2017.
[46] O. Shamir. An optimal algorithm for bandit and zero-order convex optimization with two-point feedback. The Journal of Machine Learning Research, 18(1):1703-1713, 2017.
[47] C. Totzeck. Trendy in Consensus-Basee Optimization, pages 201-226. Springer International Publishing, Cham, 2022.
[48] Trogdon. On spectral and numerical properties of random butterfly matrices. Applied Mathematics Letters, 95:48-58, 2019.
[49] F. Yousefian, A. Nedić, and U. V. Shanbhag, On stochastic gradient and subgradient methods with adaptive steplength sequences, Automatica, 48(1):56-67, 2012.
13

---

A Auxiliary Results
In this appendix, we state and collect lemmas and propositions required to prove the main results.
Notation. In the following sections, we denote with F& the filtration o(G1. Gk-1). Moreover, to simplify the notation, we define Ik as the gradient surrogate eq.2) time-step i.c. 9k D(G Auk) (Tx) and g(-) = ((G,,)()) for an arbitrary G O(d) and h > 0. We denote the normalized
Haar measure
37 by μ.. We define the unit ball I and the unit sphere Sd-1 as follow
1
R
||l|
1 and
sd-1 = {u R
||l|| = 1}.
We denote by a and ON the spherical measure and the normalized spherical measure on Sd-1 respectively. Moreover, we denote with Ia E RRdx the (truncated) identity matrix.
Lemma 2. Let B(Sa-1) be the surface area gdd and let E TRdxd be the identity matrix. Then,
3(8"-1) J8L To do(u) = d 1.
Proof. This result is proved in [21] Lemma 7.3, point (b)].
Lemma 3. Let Φ Rd + R be funcciiitt If u is uniformly distributed on St-1, then
(Blo(u) E[o())])) (Edd()) L
d
for some numerical constant 0.
Proof. The proof follows the same line as 46 Lemma 9].
A.1 Smoothing Lemma & Properties
In this appendix, provide the proof of the Smoothing Lemma (i.e. Lemma T
Proof of Smoothing Lemma.
By eq. 22.
Ealg(a,))())
Σ
8
/
f(x + hGe,) - f(x - hGee)
Ge (G(G).
1=1 0(1)
2h
By 37 Theorem 3.7].
Σ f ==1 ga-1
Ea(a(,,) 2th
(f(x ho(l)) f(x -he do xo
Since (00) is uniformly distributed on the sphere, which is symmetric with respect to the origin, we have
d Eclura,))) (x)) th Σ gg
f( )o da (u()).
As consequence of Stokes' Theorem (details in [8 Lemma 1] and Theorem A8.8)), get
B(gi).)( 1 Σ fi(()) with fa(x)) = vol( I f(x hu) du.
Rearranging terms, we get the claim.
Proposition 1 (Smoothing properties). Let fm. be the smooth approximation of f defined in eq. Then the following hold:
f is convex then Sú. is convex and, for every x E RRd.
f(x) S fa(x).
14

---

I is Lo-Lisschitz continuous i.e. Vz, y E R If(x)) - f(u)] - Loll= - vll. then fn is Lo- Lo-Lippchitz
continuous, differentiable and for every x. y E RRd
[V) - Vfflw)ll & Lovd h
1 - vll and fa(z) S f(x) + Loh.
If f is L1 L1-smooth -ee. f is differentiable and Vx, R [V) < = L1ll yll then fn is Li-soooth and for every x E R4.
[V vf(x))l hdLi and fa(a) - f(x) + L h²a
Proof. These are standard results proposed and proved in different works see for example [16 Lemma 8) 21 Proposition 7.5], 34 Proposition 2.2]. 49
Lemma 4 (Approximation Error). Let g(-) be the surrogate defined eq. 2 for arbitrary and G E O(d). Then the following hold:
(i) I is Ln-Lipschitz (see Assumption then, for every x E TRd
Balla(>))1 s 2c dLB
where is numerical constant.
(it) If f is L1 -mooohh (see Assumption then, for every Rdd
Ealllg<))1 < 2d
L.d
[+(()]
h?.
28
Proof. Note that, since directions are orthogonal, we have
de
[G]lg(())1
Σ
Eollf(s+ hGe,) - f(x NGGe))
462
By 37 Theorem 3.7].
d²
Ballg(1)1
Evg (v hw/lt) f(x 110(0)))111
(5)
4
where each (8)) is uniformly distributed on gd-1
(i): Set Ey [f + hu(()] for every (this expectation does not depend on i). Then
d²
BG(lo(())
Byt ((f( + (0) f(x - -ho(()) + (10
2
d²
Σ
Buco [([[(x +o(0)) - 77 - (J(x hol()) - (1
d²
B,oo [[[[( (h00) - 77 + (f(x - -ho()) - -))3) 10
2021
d²
Σ
E.to ho( - (10)
20212
+ Eyco [f(x hv) 7) || [1
Since ((0) is uniformly distributed on gd-1, it satisfies (1)) 12 and by symmetry we have
d²
Ea [a]l]9(>))11]
Exco [[f(x + ho))) - -1.
22
15

---

The definition of y yields
de
2
Ec[lg(z)]l
E e [((f( ( - -7)
d²
Σ E [[[ + hv)()) - Exyco [f )))
The claim follows by Lemma and the fact that f(x - hwl()) is hLooLipschiz continuous w.r.t to
v/())
(it): Equation (5) yields
d2
[G]lg(())133
E [(f([(((xx
< d²
Σ
Epw ho)
2 12
im1
+ E, [[ - - -f(=)))(10000)11
-
d
2 Boc (ff hul0) f()))
22
where the last equation follows by symmetry. Adding and subtracting ho())) derive
d²
ffx)) -(()) (v))))(((()), (V ))
Ballg(=))1
Σ
Euto
22
z 2d
t
(Fiot holl) ff(x) (vf(t), )² 2]
Σ
Ent0
h2
i=1
((VJLC)))ool))) (ool)) )²))
+ Byco
Denote by B($ the surface area of gd-1
The Descent Lemma
41
implies
2d
L2
((TT(o)))ao)))
EG Ea]lg(<))13
Σ
221
2d
8
((V.(e)))oo0)) 2
122
+
E
20
12
1=1
d.
2d
h22 +
(0)
28
(33(Sd-1)
88-1
By Lemma 2 we get the claim. Indeed,
L7
2d
3(0d
Ea[lg(())111 <
112
{(()))
(23()
2d
13
122
28
A.2 Auxiliary results and proof' for the nonsmooth setting, convex, and nonconvex.
In this subsection, for every k. we will denote by Fu the o-algebra oGoo...., GR-1).
16

---

Lemma's 5. Let f I + I. be a lower semi-continuous function and denote with S = arg min f
and f min f. Then,
(A)
Vr* S, lim | x²Il lim inf f(xx) f
3o E S s. Ik Ioo0
(B)
k
Proof. Since (B) holds, we have that exists (Ir; JoEN subsequence of (TR)Ren such that f((kk)) - S Since S + 0 and (A) we have that
E and
B lim || - x² &
Thus, the sequence (NN) is bounded and, therefore, also (2TR,) jEN is bounded. Taking convergent subsequence (Ik,,, NAE of (I )jEN, we have that exists Too s.t.
Tkin Too.
Since is assumed be lower semi-continuous, we have that (xxx) S lim inf (Ix,,,)
= -
lim ((rky)
Thus, we have that Too E S which implies, by (A), that
Ilim || zooll and lim | Tooll =0.
Hence, since Thy, subsequence of Tk+
lim | Tooll O,
and, therefore, TA - Eee E S.
Lemma (Convergence: convex non-smooth). Assume that is convex and Lo Lipschitz continuous. Let (IK)eeN be the sequence generated by Algorithm 0 and let x² E arg min f. Then, for every N, the following inequality holds:
L&d 1 - x²1 (Fx)( S 2c og 2
E||aa+11 x² [1(]
where non-negative constant independent from the dimension. Moreover, the stepsizes satisfy Assumption]
2
we have
lim f(x))= f(x) a.s.
k+++0
and there exists a random variable & taking values in in arg min f such that Tkk ->
Proof. Let ls E N. By Algorithm
T
|| xx" =1
|| x²11 19l12111 20c (gk: Tk x²)
(6)
Since fh, is convex by Proposition [ and Flgg|Fu fhn (see Lemmal I we have
(V fn (Ix), 24 x²) S (x*) fhr (*x). Thus, taking the conditional expectation with respect to Fs. by Lemma we get,
1| - | 111 < 2c L3d a 2ax (fre (xx) - frne (x*)).
E - x
Then, by Proposition I
E|]|=*+1 x [|
|| x²122 < Cs 2ax f(x")) )) 2Loaxh
Next suppose that Assumption holds. Rearranging the terms,
E|] - x 1 || - 1= + 2o (s(x)) - f(x < C. 24gohhk
with Ck e and oncy et Therefore, Robbins-Siegmund Theorem 43 implies that ([<a
x² I[eeen is a.s. convergent, au (fx((())) - f(x*)) E e1 a.s. and thus, since dr g el
Inf [ f(x") a.s.
k1+00
We derive from 32 Lemma 9.9] and Lemma 5 that there exists a random variable taking values in arg min such that Ta & Finally, continuity of / yields that lim f(xx) f(x.) a.s.
17

---

In the next Lemma, to derive bounds on function values, we study the sequence (fhr (TI-1) fI, (ox) RENN is the difference between the smoothed function iteration & evaluated 2g and at 1-+1- It corresponds to the function value decrease between the iterations k + 1 and l if hg is constant.
Lemma 7 (Function Value decrease: nonconvex non-smooth setting). Under Assumption 7 let (oN)keN be the sequence generated by Algorithm 7
Then,
Lgdv d o l
Bi (o -nn (1x) S - oo.
hg
where с is a numerical constant.
Proof. By Lemma
have that Sh, is Lo dovd/h smoothh. Thus, by the Descent Lemma 41
Lovd 2h
Jhu (I++1) Ira (xx) S on (V (Ix)), g6) +
Taking the conditional expectation with respect to Fk.
Lovd ag
E (I)++)) - [ha = -ax V.n -axll (TW))1
Bllggn gxllgl11 11]]...
(8)
2/h
The claim follows from Lemma 4
A.3 Auxiliary results for smooth setting.
Lemma 8 (Function value decrease: convex smooth setting). Under Assumption 3 let (oR) be the sequence generated by Algorithm [
Then the following holds:
1[ - f(xx) - -ax -ax( ( Lid
Lidanhh
Lid
as
+
a
e
8
46
Proof. By the Descent Lemma 41 and Algorithm
!
4
f(*x+1) f(*)) S aa (Vf(z)). gx) +
2 ax 1
Taking the conditional expectation and by Lemma
L1 2d
L7d
B[s(T++))\FF) - f(xx) S -as (V/(x)), V (xx))
[((x))11
28 1
2 O
Adding and subtracting (xxx),
B[f(xx+11)7H) f(xx) S -0k (Vf(x)), Vhn (xx) (Nx)-)((x))) au aaal-aa/////i
I
2d
1..22
+ 2
a
[f(xx)11
hi
28
By Cauchy- Sauchy-Schwazz inequaliyyy and Proposition
Lid
E [(((()+)) f(xr)) S ax
2
Li
2d
+
2
k
28
1/
By Young's inequality,
B( f(xx) <
on
+ an
[[(xx)]1 a (v/>))
8
Ia
2d
T.d2
0
[/(xx))1
h
28
1
Lid
13daash
Ldd
ak
a#
[V5(2=))11
0 oi
2
8
4
This concludes the proof.
18

---

Lemma (Convergence in smooth setting), Let (3% KEN be the sequence generated by Algorithm 7 and let arg min fxr) Then, under Assumption the following inequality holds =ER
2d a (f(x)))] + L2d a hi
0 x [1 100 x=1 =
2
+ Lydayhu 1|= xx'1 || 20x (VI(z)), II - x²)
Moreover, if f is convex, Assumption holds and as S a < e/(2dL1). Then
(=xIlV((==)11) MEN ll a.s.
(| INEN is a.s. convergent.
(on(f(zz)) - f(z") MREN E ll a.s.
there exists random variable à taking values in arg min f such that 2s + it a.s. and lim f(xx) min
k0000
Proof. We have
| 2=1 11. - 2²1 = 11 20 (gk, IK x²).
Taking the conditional expectation,
[|| - x= [R - |20 - x= =f | - 20x (V She (2x), Ix - x²).
For every k, set ux = 1, - x=1 ||- By Lemma 4
2d
L2
Eluh+1 ui - a (
h 2ax (Tx) -x-x''
2
Note that
-20x Sh ( -x*) x²) = 20 (i. (xx) Vf(Ix), - -xa) - 200 (V [ - =x))
Thus,
u < 2d
f. a V L7dL
E U uk4+1
Fu]
28 ag
20x (V (Ixx) (f(x)), -3x) -20R (VI(I)), x²)
By the Cauchy-Schwerr inequality,
2d u < 7
2 ( + L.d
E 1TR
of
2
+ 2a V fn. (Tx) - vf(xx)]] 20x (Vf(I)), 2x -x)).
The first claim follows from Proposition 1 By Proposition land Young's 's inequality with parameter
Th 0g hgg, we get
2d
Ldd fxx))
2
FR - u 5
2
+
21 a
Lid
2 + LidTA ui - 20y (VI(z)), - =x)).
27k
2
2d a V (xx))1
L3dd
(9)
28
Lid
Lid
04 h4
ag u
20 (V)())), xz - x"
Since f is convex, by Baillon-Hadder Theorem D. we derive that
1 d ur 5 -2 ( TM az
LddZ
E
2
28
+ hh + Lid
Lid
2 ach
19

---

By Assumption
4
Ind
E
[ic u I -2
11 -
V/
oxhi uf
2
=iA
pk
L,d
T.dd
oxhx +
a
28
=:Cx
Note that A 0. Thus, rearranging the terms
E Elu + [Fa] (1 + p=)u 2Aox
/ [+x)) Cs.
Since Pk. Ck E e by Assumption Robbins-Siegmund Theorem [43 ensures that (ut) is convergent and [ ((x) kEN E e a.s. Since f convex, it follows from (9) that
2d
Fi (1 pa)u a ( (fx(f(zx) f(x") Ck.
Goobins-Siegmund Theorem implies that Og E (t therefore
43 implies that (fx()( [(x))))EEN E ll a.s. Assumption
lim inf f(xx) f(x") 0a.s.
(10)
By Lemma 8 and Assumption we have that the sequence [ -((((())))) -f(xx))(((x)) upper-bounded by a sequence in é1, Thus, by Gobbins--iegminn Theorem [43] lime (f(I)) f(x ()))) exists a.s. Then, it follows from (TO) that
(I
lim f(zx) = f(x") a.s. k+000
Moreover, saw before, ([] x=" k)N is convergent for every x² arg min f. Then, by Opial's Lemma [40]. there exists a random variable à taking values in arg min f such that 2y +
Lemma 10 (Gradient bound: convex smooth setting). Suppose that Assumptions and hold, and assume J to be convex. Let (Ix (EE))ENN be the sequence generated by Algorithm Then, for every l E M and every x² E arg min S,
a ([)
Ss +
Σ
pi - 1²1
24
and
/
/E[l| - 10 V + Σ
Pii
1=0
where
A
a
S = ||:0 x111 |||}
Σ
C.
L1
1=0
L7d
C
and
Pk = Lidaphx.
28
Proof. By Lemma 9 we derive
2d
E G[| x² [1]]]
||
- 2212 <
ag 2 [ C
pk|lax x= 2ax (VI(x)), vs x²).
By Theorem and Assumption
E|]| 2*[[[[Fa] || 201 -2 [[(x))] + + + pk|1 =='|-
20

---

Let up = z
x² 12. Taking the full expectation, by Jensen inequality we have
"ht o S 22 V Pauk + Ck.
Summing the previous inequality from 0, k, we get
24 #k+
04 E V o5 + E Σ
P:ut.
(11)
=
Since Uy is non-negative, the first claim of the lemma follows. Since A >> 0, Ph 2 0. Sq. is non decreasing, and S 2 ug in (T), then
4 S S +
P.u
i=0
Thus, the (discrete) Bihari's Lemma 32 Lemma 9.8] yields
Uk+1 S
Pi + Sk +
1 (
Σ Pi
}
2/2
S
k
Pi
2 1=0
1=0
t=0
concluding the proof.
B Proofs of Main Results
B.1 Proof of Theorem 1
By Lemma 6
l&d
E El x [1 - || x²11 20x(f(x))))))))))))))))))))))))) f(x))) s 2c
i a 2Loaghkk
Rearranging the terms, taking the full expectation, and summing the first iterations
||a0 11
dL
Σ a [[((x.)) - -((((x)) f(x))))
Σ a + Lo
Σ ashs
Let T2 = Σ (uri:/( / Σ oz). Dividing by Σ az and observing that by convexity we have
1=0 ==
==0
Σ a
E((((x;) - f(x)))
[fx( min f]
a
t=0
we get the first claim. Under Assumption 2 the second claim holds by Lemma 6
B.2 Proof of Corollary 1
By Theorem
1
dL3
[ff( min fl
+
a
Lu
Σ
aghs
2
Σ
1=0
1=0
1=0
Replacing ag and h4 with the sequences in the statement,
Ci
C2 ht dCss ko
Elf(Ex) - f(x))]
akl-o
+
ko a.
21

---

with
C1 (1 12 x²1
Lo 8)
cLi 8)
C2
and C3 (1 20)
(1 1 p)
The second point of the corollary can be proved replacing at a and hy h. Now, to prove the third point, fix & E (0, 1). Since we want [[f(Tx) - f(x" S E, we impose
20 x=12
adT
+ Loh S.
20
Choosing hig
h S e
21.0 get the previous inequality it sufficient to impose
I00 12
cdLo
2ock
2
We fix a priori number of iterations and we minimize the left handside with respect to a, obtaining
100 - x"1l V2cKLn Lo
Thus, for hy S zi.. as above and
20
8
2 Lied
K
we have ([f((x) S E. Note that, since the computation of the surrogate requires 2 function evaluations, ensure error of we need to perform number of function evaluations of the order
O(de 3).
This concludes the proof.
B.3
Proof of Theorem
2
By Lemma 7
Lğdvd oj.
fl (*+41)) fa((x)) < -ax -allVV//<>>>
l
Taking the full expectation and rearranging the terms,
B/ < (()) - fn((x+1)] + Lÿdv dd
h
Next sum from = 0 to
By definition of fho we have fn( 2 min f for every x E Rd. thus,
a V 1 (zo) min fl + Lddvd
(12)
h
0m0
The claim follows.
B.4 Proof of Corollaries and 3
By Theorem
(M) ((in(50)- ((x))) Lodvd
0s
Due to the choice of or = aKk: + -0 with 0 (1/2,11 1) and a
0, we get
(h) < a(k + 1)1-0
C₁
+ Cadvda th (k: 13
22

---

where
C1 = 100 x x111 0) and C2 : cL 0)
20)
we choose ag = a, we derive
fi((z)) min f eLIdvdo
(h) Tk 1
(13)
th If fix a priori a number of iteration K and minimize the right handside with respect to a, get
ak
(a((o))-((x1))
KeLidVd
Let E (0, 1). Choosing a á, we get 77 S = for
K (fn(zo) - f(xx' )LL
(14)
Lh
This concludes the proof of Corollary To prove Corollary we a maximum number of iterations K E N and consider the random variable T of the statement. Let 0. be the h-Goldstein subdifferential defined in Definition I almost surely, therefore
follows from [34] Theorem 3.1] that (x( E fx1(
E, min[|nnl?
n E 0nf((r)) s E, ([IV\(=>>))))
In addition, Theorem 2
yields
E, Vcl\V/(())1 = Σ aj B VcllV(((>>)111
Σ
aj
S (zfa min fI + Lādva
&
l
h
==
Thus,
Ey[l||| n E (zr)) < E, E fo (z1))
m/ (K)
Hence, for a = ā and K chosen s.t. inequality [4] holds, we have
Br[|n| 1 (xr))
This concludes the proof.
B.5
Proof of Theorem 3
By Lemma 9
2d - 101e11=2111115 21119FF || x²12 S
1 a
v (zx)]1 + +aa (Vf(xx), x²
xx)
17
+
Lodog 110 xxll +
2
28
=
By convexity,
2d
B B [R) - || - x x=112 <
2 [ - (fa f(x)))
-
{ ax
+ >M||x>* - xxll + Ck-
Rearranging the terms and taking the full expectation,
2d
2 2Elau - f(x")) < E[|lz* - x 2=112 8 || - x"11 a Elll V((x))))
+ ps E|||=*
xx + Ch.
23

---

Since Bi|l2" - xxil] BIv x- 24 ). Jensen's inequality implies that
2d
2 2E(ax(((x)) (x)))] = B|| xx 112 1 x01
1 o
E
(fxx))
+ Pk |||| - x& 123 + Ck.
Denoting with un = =Bll2x - xx [|a] and taking the sum from to
2d
2 Σ ai Bls(s)) [(x n
C.
a
V(
Bll
Σ pine
i=U
S:Sk
2d
k
S
S +
[[ +
p:us:
1=0
1=0
where the last inequality holds by Assumption 4 have
Let A = (1/L1 (d/E)a). By Lemma 10
k
dā
o. [[f())
f(x'''))))))))))))))
-
Sh +
Sy +
Σ
Piu
+
Σ
Psul:)
2
At
i=0
0=0
AA dā
Sk + Σ
pilli)
21A
1=0
A dā
Sk + Σ
p.VVS, S. + Σ E.))
21A
t=0
/=0
k Let 8x Σ a ( Σ o:). Dividing both sides by
Oi convexity yields
Σ a
- -((())))(())))
[[f(fx) min f
a
B.6 Proof of Corollary 4
In this proof, we use the same notation as the one in the proof of Theorem parameters, we have
3 By the choices of the
k
1
Σ
ps S
Cidah
with C1 =
0 0-1 1
1=0
+ C2 d
with C3 = L
S < [00 x 2²11
2222
20 1
Thus, using these inequalities in Theorem 3
we get
(A + dā 1100-x*11200 da212
D&
+ VIz0 x= PPdaah
2
C4 dah
Vi Coda2t2).
with
130
L 2
Cs = I,70
C = 1-1
C =
- (1)
Dividing by Σ
as we get
1=0
Eff(Ex) min fl < C
ak'
24

---

Note that by Assumption 4 a < f/(dL1), thus 1/a > (dL1)/E. The algorithm performs 2 function evaluations at each iteration. Thus, to guarantee [[f( min fl S E for (0, 1), the algorithm has to perform number of function evaluations in the order of
O(de
Assuming, instead, ag S â (2dL))) by Lemma 9 we get the last claim; i., there exists random variable & taking values in arg min f s.t. 24 + 2 a.s.
B.7 Proof of Theorem
Set C1 = (dL1)/2. It follows from Lemma 8 that
[(f(*++1) f(x)) < - _ Lid a ) ( + Gauh
Lide
48 agl
2
Taking the full expectation and rearranging the terms, and recalling the definition of A
B[f(zk) f(xx+1)]+ Caahh Summing for i = 0, ..... k and observing that min f S ((x) for every x,
L7d af
Aag E|
4l
à
Gaah
Lid k
a El Vf < (zo) min f+ Σ
Σ n?
2
4
0=
1
i=0
Divinding by
we get the claim.
B.8 Proof of Corollary S
(i): From the choice of og and hgs we have
k
20ch²
k
20
oah? z
Σ
2 <
20
20 -
1=
It follows from Theorem 4 that
f(x0) - min f &ah + 1222
S Aakk
l
with C1 4 and C2
1
1.0
2(20 1)
(Wi): It follows directly from Theorem taking into account that
Σ a kah2.
Σ af ka²
and setting C1 1./8 and C2 /./4..
C Experimental Details
In this appendix, we report details the experiments performed. We implemented every script in Python3 (version 3.9.11) 1) and used numpy (version 1.22.2) 27 and matplotlib (version 3.5.1) 29 libraries.
Machine used perform the experiments. In the following table, we describe the features the machine used to perform the experiments in Section 4
Table 1: Machine used perform the experiments
Feature
OS
Debian GNU/Linux
CPU(s)
4 Intel(R) Core(TM) i7- 165G7 11th Gen @ 2.80GHz
CPU Core(s)
4
RAM
8 GB
25

---

Target Functions. We considered two synthetic target functions: a convex smooth function fi and convex non-smooth function f2 defined as follows
(Convex Smooth) [1(x) - | with A E Rdxd
(Convex Non-smooth) [2(x) = || T||
where A is a random Gaussian matrix (i.e. At, ~ (O, 1)) and ;= [0, 1. 1JT.
Choice of the number of directions. We report here the details of the first experiment of Section 4 For these experiments, we consider 50 and we use, for the smooth convex case, the following parameters
10
Og 0.99
and
dL,
++1
The constant L1 is computed as the maximum eigenvalue of the matrix AT Note that this parameter choice satisfies Assumption 4 For the non-smooth target, we used
0+ V 1/2-10 and hy =
10-7
k +
Note that this parameter configuration satisfies Assumption 2 evaluations considered is 4000. The direction matrices G are generated with the QR method see Appendix D
The maximum number of function
Comparison with Finite-difference methods. In Section we compare finite-difference method with different choice of directions. In order to make fair comparison consider only central finite- differences. However, note that Algorithm| [ can be modified (in practice) considering computationally cheaper gradient estimators see Remark [ For these experiments, we consider for methods with multiple directions. The maximum number of function evaluations is 1000 for both smooth and non-smooth targets and the direction matrices G2 for Algorithm are generated with the QR method see Appendix D To solve the smooth problem we consider the following parameter choice for every method
10 and l
E
10-7
c dL,
and
h4
d2 + +1)
where L, is computed taking the maximum eigenvalue of ATA. For Algorithm and finite-difference with single and multiple spherical directions c = 0.99 while it is equal to = 0.11 for Liite-difference with single and multiple Gaussian directions. We made this choice since for finite- differeeffe methods with Gaussian directions we observed divergence for larger choices of c see Figure
Smooth Convex Target
Smooth Convex Target
397
Single Spherical
Sayle Gaupiin
399
Mus Sphericl
Single Countin
Seyle Sphenoo
1.2
Spherica
Our
10-
10'
zic
600
60
1306
60o
1006
function availstionn
function evalustions
Figure 3: From left to right, comparison of finite-difference methods for smooth convex target with é = 0.99 and c = 0.2 for methods with Gaussian directions.
For the non-smooth convex target, we considered the following parameter choice
a2 = 7 k-12 1/2-100-
and hg = [2(k+1)
26

---

For every method, we selected = 0.65 except for the method with multiple Gaussian directions in which we selected 0.08 since it provided better performances see Figure
Figure 4: From left to right and up to down, comparison of finite difference method with different directions and different values of for multiple Gaussian directions. The values of considered are the following [0.085, 0.089, 0.09, 0.1, 0.2, 0.3, 0.5, 0.65)
D Techniques to Generate Orthogonal Direction Matrices
In the literature, different algorithms were proposed generate orthogonal matrices see for instance 23 38
2 4. 44 8 and references therein. Such methods can be used to generate the direction matrices C required for the iteration proposed in Algorithm (T) In this appendix, briefly discuss three of them.
1 28
QR factorization. As observed in 332 42). way generate orthogonal consists in generating a random Gaussian matrix A E idddd with A.j ~ (000, 1) and perform the QR factorization i.c. QR. Then, the direction matrix the truncation of the matrix i.c. Qla,e-
Householder Reflection. To obtain a direction matrix, we can use a Householder reflector. This can be done by sampling vector from the unit sphere gd-1, The direction matrix defined as Householder reflector, given by
G == 200
with Rdxd identity matrix. To obtain the desired matrix, we compute the product of G with Ia.e- ie.. we take the first é columns. The (truncated) identity matrix can be generated and stored offline (note that since it is very sparse, can be stored using a sparse format (e.g. the COO format proposed in scikit-learn library(9).. In this way, we can save resources in high -igh-dimensionll settings. In order to quantify the time-cost of this procedure, we compared the time of generating this kind of matrix with random matrices with different dimensions. For this experiment, we consider the { = d case i.e. the most expensive. Matrices are computed in CPU and the details of the machine used are described in Appendix C We report the mean and standard deviation of the time using 500 repetitions. In Figure 5) we compare the time-cost of generating orthogonal matrices with this procedure against generating random matrices while in Table 2 we report the mean and standard deviation of the results.
27

---

Computational Time to generate matrices
random Geussian
10"²
random spherical
Householder
10aa-
10"3
187
10
256
Figure 5: Time comparison in CPU of different methods to generate direction matrices.
observe that using this strategy limit the cost of generating random orthogonal matrices. particular, for dimensions larger than 32, our method is faster than random gaussian and spherical directions.
In Figure 5]
Table Comparison of the time-cost (seconds) of generating random and orthogonal matrices with different dimensions
d
Random Gaussian
Random Spherical
Householder
9.27
10-7 + 7.96
10-7
5.49 x 10-6 2.05 x 10-6
9.32 10 3.34
10-6
4
1.30 10 6 7.21
10-7
6.56 x 10 2.63 x 10
1.12 x 5
10
8
2.18 x 20-6 x 10-7
8.01
10-645..33 x 10
1.11 x 5
x 10
8
2.18 10-6 10- 6.06 x 10-7
8.01 10-5532 x 10-6
1,11 x 10 10
10
16
5.69 x 10 #1.11.61
1.61
10
1.15 155510044.100 10-6
1.18 10-5
10
32
1.78 x 1 6.42 10-6
2.49 149x10-5+1.33 10
1.16 x 5
10
64
6.58 10-5 7.03 10-6
7.74x106+11.5xx-
1.62 10 4
10-6
128
2.73 4 10-5
2.98 108x100-+2455x000-
3.32 x 10-5 +.002 10-6
256
1.26 x 2
10-5 10-5
x 1.36 10-3 10-5
1.20 x 10-4+1.04 x 10
512
5.50 x 10-3 10
5.91 x 10-3 10
1,22 10-3 10
1024
2.16 x 10-2
10
2.41 x 10 10
4.83 x 10
10
2048
8.92 x 10- 8.19
x 10-2
1.04 x 10-1 1.03 x 10-1
2.40
3.87 10-2
Moreover, if more computational resources available, can build Householder reflectors
GI.... Gm using m random vectors U1, "., Um sampled i.i from Sd-1 and define the direction matrix as
G G2 GmI
is important to note that when m d, this procedure is equivalent to using the QR factorization.
Haar Butterfly matrices. We can build orthogonal matrices using Butterfly matrices [88 Let
( [1], build orthogonal matrix of dimension d = 2" with the following recursion
((a) cos con (0n)) sin (in cos(O,,) G(n
(0, (n-11
where 0,, sampled uniformly in [0. 2T)) Then we compute GTa. (we take first columns). The construction of Haar butterfly matrices is faster than previous methods because it only requires simple operations. However, this procedure allows build only matrices with 2n for n 2 In literature, different methods were proposed to cope with this limitation e.g. [23
28

---

E Limitations
In this appendix, we discuss the main practical limitations of Algorithm Like all finite-difference methods with multiple directions, O-ZD requires multiple function evaluations execute a single step. In many practical upplications, function evaluations can be conesmonsu,nng, leading to the use of small number of directions É. This may result in poor performance as observed in numerical experiments. As for the subgradient method, in O-ZD the step size significantly affects performance, and tuning it can be challenging. To address this limitation, an adaptive stepsize selection method could be proposed. Furthermore, decreasing the sequence hgg too quickly can lead to numerical instability, as noted
42
F Other Experiments
We performed other experiments in minimizing convex functions. We considered the targets defined and, for each experiment, we reported the mean and standard deviation using 20 repetitions.
in Table and,
Elestic Net
Total Variation
Huber Loss
Sigo wkl
5 Spheicaal
Saherol
Sio 2n1 BIl 1100 1ite 2000
256
1900
200 USM
function Evalual ona
Infinity Norm
Sparse Group Lasso
L1 Norm
Spherca
8.127
Siga CAALE
3/17
Figure 6: Function values per function evaluation in optimizing functions with different algorithms.
In Figure 6 methods.
observe that structured fiite-ddifference performs better than unstructured
Table 3: Functions used and relative dimension and number of directions considered.
Name
Definition
d
Sparse Group Lasso
f(x) Σ
(000
50
25
0.5 1
|| S 8
Huber Loss
f(x) 8||x||12 0.55
for
8 >
0
50
25
otherwise
Elastic Net
f(x) = ad |
50
25
LI
f(x) [x]]
50
25
Infinity Norm
f(z) = r
50
20
Total Variation
f(z) TV
50
25
In Table 3
we define the function used for the experiments. In particular:
29

---

Sparse Group Lasso: is set to and given an æ E R ((1) is vector obtained by taking
3 entries of
Huber Loss: 0 is set to 0.5.
Elastic Net: 3 are set to 0.5.
30


Reviews:
[{"review_id": "FmG4GJlXQ9", "reviewer": "Reviewer_eK7Y", "summary": "This paper introduces and analyzes  a structured finite difference algorithm for non-smooth optimization problems. The algorithm is built on a smooth approximation of the nonsmooth loss function and a structured finite difference approximation. The convergence of the proposed algorithm in non-smooth convex, non-smooth non-convex, smooth convex and smooth non-convex cases are studied.", "questions": "1. In Theorem 1 and Corollary 1, the error is linear in \\ell. The convergence rate seems only depends on \\theta. Does \\ell affect the convergence rate?\n\n2. From the theories in this paper, setting \\theta towards 1/2 provides better convergence rate. But what is the hurt by doing that?", "limitations": "The limitation is not discussed.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. A simple algorithm using structured finite difference approximation of the gradient is proposed.\n\n2. The convergence behavior of the proposed algorithm in four cases are studied.", "weaknesses": "I think the relation between the convergence of the proposed method and the number of directions chosen is not clearly analyzed. The error in Theorem 1 is only linear in \\ell. While numerical experiments suggests \\ell has an impact on the convergence rate.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her useful comments. We answer here to his/her questions.\n\n## Questions\n**Dependence on $\\ell$ in the rate:** Yes, the number of directions $\\ell$ affects the convergence rate but only in the constants. Indeed,  Theorem $1$ states that\n\n\\begin{equation}\n\\\\mathbb{E}[f(\\bar{x}_k) - \\min f] \\leq S_k / A_k = \\frac{1}{\\sum\\_{i=0}^k \\alpha_i } \\Big( \\frac{\\|\\| x_0 - x^*\\| \\|^2}{2} + c\\frac{d L\\_0^2}{\\ell} \\sum\\limits\\_{i=0}^k \\alpha_i^2 + L\\_0 \\sum\\limits\\_{i = 0}^k \\alpha_i h_i \\Big).\n\\end{equation}\n\nNote that the impact of the dimension in the rates depends on the choice of the stepsize $\\alpha_k$ and the discretization parameter $h_k$. Since the second term on the right part of the inequality depends on $\\frac{d}{\\ell}$, in order to reduce the impact of the dimension in the rate and obtain the optimal complexity, we need to include a $\\sqrt{\\frac{\\ell}{d}}$ in the stepsize $\\alpha_k$. The same observation holds for Corollary $1$ (specifically in points $(i)$ and $(ii)$). We will include these observations in the discussion paragraph below Corollary $1$.\n\n**Choice of $\\theta$ in the stepsize:** As for classic methods like subgradient descent and stochastic gradient descent, the choice of $\\theta$ towards $1/2$ is the best choice we can provide. Indeed, note that the second term in the right part of the inequality of Theorem $1$ (see the inequality above) is an error term which does not depend on the discretization sequence $h_k$ and thus it is bounded if and only if $\\alpha_k^2$ goes to $0$. Thus, ideally, we would choose a sequence that decreases as fast as possible. However, if $\\alpha_k$ goes to $0$ too fast, we might have no convergence, for this reason we need to take $\\alpha_k \\not\\in \\ell^{1}$ and $\\alpha_k^2 \\in \\ell^{1}$. The result is not surprising and is in line with the classical results obtained by studying algorithms with errors or based on subgradients.\n\n## Limitations\nLimitations of the algorithm are discussed in Appendix E. We will add a reference in the main paper. \n"}, {"review_id": "MzeN2ibXby", "reviewer": "Reviewer_VJ1k", "summary": "Zeroth-order optimization is the sub-field of optimization concerned with solving \n$$ x^{\\star} \\in \\operatorname*{argmin}_{x \\in \\mathbb{R}^d}$$ \n_without_ using gradient information. This paper studies the most general case, where $f$ is assumed neither smooth nor convex. They propose a structured approach to the standard finite-difference gradient approximation, which means that the sampling directions are taken to be orthogonal.\n\nThe main contributions of this paper are theoretical; it provides a fine-grained analysis of the convergence rate of the proposed algorithm under a variety of assumptions. The theoretical claims are supported by two numerical experiments.", "questions": "1. In the related work section, you state that your method achieves the optimal dependence on dimension and then cite [1]. But this paper deals with stochastic zeroth-order optimization with a particular \"two-point\" query model, whereas in your paper you consider the non-stochastic/noise-free setting. \n\nI did look at [1] but could not find a clear statement of the optimal dependence on dimension in the noise-free setting (I believe it is $\\mathcal{O}(d)$). Could you make this a bit clearer in your exposition? Also, what is the authoritative reference for the optimal $d$-dependence for noise-free DFO?\n\n[1] J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono. Optimal rates for zero-order convex optimization: The power of two function evaluations.", "limitations": "Yes.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- The first 4 pages are easy to read, and position the problem under consideration (orthogonal sample directions for non-smooth optimization) nicely with respect to prior work\n- This paper is definitely a technical advancement on prior work.\n- The mathematical technique displayed in this paper is impressive. As an example, I enjoyed the use of the Goldstein subdifferential. The correct notion of \"approximate stationarity\" in this setting is quite subtle but you handle it very well.\n - In Corollaries 1, 2, and 4 I appreciated that you provided convergence guarantees in both the constant step-size and decreasing step-size settings. Both are useful\n", "weaknesses": "- Overall I found this paper to be a little too heavy on theory and a little too light on experiment. I would suggest relegating some of the parameter configurations of Corollaries 2, 4, and 5 to an appendix and adding an additional experiment (see next point).\n - I'd strongly suggest adding more experiments; in my opinion the two simple functions you tested do not provide enough data to draw any conclusions. It would be particularly interesting to see more non-smooth experiments. \n\nMinor points:\n - In line 235 \"the ball centered in 0\" should be \"the ball centered at x\"\n - In line 239 you write $\\partial f_h(x_I)$ but I think you mean $\\partial_h f(x_I)$. \n - Why is there a staircase pattern in the purple graph in Figure 1?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her observations and comments. We answer here to his/her questions\n\n## Weakness\nThe main goal of this article is to provide the first analysis of a structured zeroth-order algorithm in non-smooth setting providing the mathematical tools (i.e. the Smoothing Lemma) required to analyze its convergence, and hopefully future extensions (e.g. stochastic setting). The purpose of the numerical experiments section is to empirically show the properties indicated by the theory. \n\nAn exhaustive empirical comparison is out of the scope of this work but we are considering it as a research direction.\n\nWe performed other experiments that we will include in the paper (some of them are included in the global answer - see Figure 1). However, we want to underline that in order to understand the practical behavior of the algorithm, an exhaustive empirical analysis should be performed keeping in account also the other parameters (e.g. stepsize choice or the discretization choice) and this is out of the scope of this work.\n\n\n**Minor:** in figure 1 we analyze the impact of the number of directions $\\ell$. In particular, we plotted in x-axis the number of function evaluations. Now, since to compute the gradient estimator we need to perform $2\\ell$ function evaluations, we repeated $2\\ell$ times the target function values (this is explained in line 323-324).\n\n## Questions\nThe reference is the same, note that by replacing the stochastic oracle with a noise-free oracle the proofs (of the propositions) follow the same line. \n"}, {"review_id": "DmiwQEqBLW", "reviewer": "Reviewer_KbSU", "summary": "This paper proposed a structured zeroth-order estimator for non-smooth optimization. The proposed algorithm using this estimator could achieve optimal convergence rate for non-smooth convex optimization, also achieve a convergence rate in terms of Goldstein stationarity for non-smooth non-convex optimization. Numerical experiments are provided to show the efficiency of the proposed algorithm.", "questions": "Major:\n1. Could the authors talk about possible applications of the proposed method? For example one could argue that zeroth-order smooth optimization can be applied to black-box attack of neural nets. But I’m not aware of the application of non-smooth situation.\n\n2. The design of the zeroth-order estimator may cost more time, due to the fact that it requires sampling orthogonal matrices (either using QR or other matrix operations). More specifically, I’m wondering how the proposed method behaves comparing to simply sample $\\ell$ vectors (with/without replacement) from the canonical orthogonal basis.\n\n3. Regarding the numerical experiments. The choice of $\\ell$ seem to be arbitrary. One could imagine that with a larger $\\ell$, the per-iteration convergence would be faster, but it would cost more time to construct the estimator. Can the author give some insight on how to choose $\\ell$ in an efficient manner? Another question related is that would the convergence and numerical behavior behave better if we vary $\\ell$ in each iteration?\n\n4. For the numerical experiment. I understand that due to page limit the authors moved a lot of details into the appendix, but it should be necessary to at least include the problem on which the experiments are conducted. Also since all the plots in the numerical experiments use “function evaluation” as ax-axis. I’m wondering what would the plot look like if we could the CPU time, since the proposed method could consume more time to construct the orthogonal matrix.\n\n5. This is a personal comment and the authors may consider this for future works. The proposed method and analysis are all for deterministic optimization, whereas modern machine learning is about stochastic optimization. It would be interesting to see the convergence behavior of the proposed method for stochastic optimization problems, where the smooth stochastic case has been studied[1]. For zeroth-order non-smooth stochastic optimization, Goldstein stationarity seems to be necessary.\n\n\nMinor:\n1. Please add some reference in the abstract, e.g. line 7 “Recently,… improve performance”.\n\nReferences:\n[1] Balasubramanian, Krishnakumar, and Saeed Ghadimi. \"Zeroth-order nonconvex stochastic optimization: Handling constraints, high dimensionality, and saddle points.\" Foundations of Computational Mathematics (2022): 1-42.", "limitations": "The limitation is well stated in weakness and question sections.\n\nI’m not aware of any potential negative social impact of this work.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The theory is well-rounded with both theoretical and numerical evidences. The proposed structured zeroth-order estimator is the first one for non-smooth optimization. The use of the Goldstein stationary is pretty novel and an interesting direction to further explore.", "weaknesses": "(Please respond to the Questions section directly) The design of the zeroth-order estimator may cost more time, due to the fact that it requires sampling orthogonal matrices; The problem studied doesn’t cover stochastic situations; The numerical experiments are not adequate.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her comments and his/her suggestions. We answer here to his/her questions\n\n## Questions \n1. An example of an application of our algorithm we are working on is gain tuning in robotics.\n2. Also coordinate directions (i.e. random canonical bases) provide good results. However, if the target function is \"sparse\" in the sense that the intrisic dimensionaly is smaller than the dimension (e.g. only some dimensions are relevant), sampling random canonical bases can provide no improvement while with random rotation we should be able to cover such settings. \n3. In practice, the choice of $\\ell$ depends on the function evaluation budget i.e. the maximum number of function evaluations that we can spend to face the optimization problem. Such budget depends on the time-cost of the single function evaluation and on the amount of time we want to spend to solve such a problem. According to the theory, it is better to perform less iterations with large $\\ell$ instead of many iteration with small $\\ell$. Thus a possible criterion can be to choose $\\ell$ as large as possible (according to time constraints and computational resources). We thank the reviewer for the interesting idea on varying the number of directions $\\ell$ per iteration. To the best of our knowledge there is no finite-difference method that analyze such a setting and it can be an interesting research direction.\n4. We agree with the reviewer that it could be interesting to plot the CPU time vs function values. We repeated the experiments and we will include such plots either in Appendix or in the Numerical Experiments section. We included such plots in the global answer (i.e. the pdf file) - see Figure 2. In that case the orthogonal matrix is generated as a single Householder reflector (see Appendix D), in that way, we can precompute the truncated identity $I \\in \\mathbb{R}^{d \\times \\ell}$ and at every iteration we just have to generate a vector $v \\in \\mathbb{S}^{d - 1}$ and compute the outer product. The parameters are chosen as indicated in the paper (see Appendix C). Note that in order to make the comparison fair, the number of iterations performed by methods with multiple directions is smaller than the number of iterations performed by single direction methods (more precisely, given a budget of $T = 1000$ function evaluations, the number of iterations performed is $T/(2\\ell)$). According to these experiments, we can still observe an advantage in using orthogonal directions. Moreover, we want to underline that in the literature there are different (and faster) methods to generate orthogonal matrices (see references in Appendix D). We want also to underline that in order to understand this kind of properties, we should perform an exhaustive empirical analysis which is out-of-the-scope of this work. However, note that previous works confirm that structured methods provide better practical performance (see e.g. [1]).\n5. We thank the reviewer for the suggestion and we confirm that a research direction that we are considering consists in extending these results to the stochastic setting (considering different noise models). In this work, we wanted to introduce the first analysis for non-smooth structured zeroth and in particular the Smoothing Lemma for structured directions providing the basis to analysis structured zeroth order methods in non-smooth setting. \n\n**Minor:** we will add some references.\n\n**Reference**\n1. A. S. Berahas, L. Cao, K. Choromanski, and K. Scheinberg. A theoretical and empirical comparison of gradient approximations in derivative-free optimization. Foundations of Computational Mathematics, 22(2):507–560, Apr 2022\n"}, {"review_id": "9gQzWLVWOs", "reviewer": "Reviewer_77xV", "summary": "This paper analyzes the convergence rate of structured zeroth order optimization, Whose iterations’ descent direction Is chosen by random orthogonal group. It applies to the most general non smooth setting, and this paper gives the convergence rate of all specific settings of interest.", "questions": "In Eq(2), is there a motivation why you want to fix $l$ and vary $G$, is $h$ fixed or random? \n\nIn Algo. 1, I think it means sample $G_k$ i.i.d. from “a uniform distribution” on $O(d)$, is that correct?\n\nWhat structure do the Algo’s chosen direction satisfy? For example, the paper mentions “orthogonality [30,40]” as a type of structure, but are the directions in Algo. 1 orthogonal? What is the difference between 1) “sampling $G_k$ randomly” and 2) zeroth order GD when you just randomly choose a point nearby and estimate $f(x+dx) – f(x-dx)$? With this part clarified, I would be willing to increase the score.\n\nThere is another type of zeroth order method – although not always applicable in practice but worth mention. When the objective function is analytic, one can use complex number to make estimation variance smaller by estimating the gradient by \n$$ Im( f(x+yi) – f(x-yi) ) / 2y $$\nwhere $y$ is small. For example, if $f(x) = x^3$, then \n$$ Im( f(x+yi) – f(x-yi) ) / 2y = 3x^2 + O(y^2)$$\n$$ ( f(x+y) – f(x-y) ) / 2y = 3x^2 + O(y)$$\nWith complex number, the noise or variance on the real part can be ignored, with only high order noise kept.\n\n======\n\nRaised to 7.\n", "limitations": "N/A", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper analyzes an important optimization method, and provides concrete mathematical analysis to prove the theories. The result is novel and comprehensive. This paper has a very detailed discussion of many settings, each of them has the convergence rate in its case. ", "weaknesses": "No major weaknesses, questions below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her comments. We answer here to his/her questions.\n\n## Question \n\n**Eq. 2 and gradient estimator.** In eq.2, we introduce our gradient estimator. The parameter $h$ controls the smoothness of the smoothed target (see Proposition 1) and, in the smooth setting, the quality of the estimator (see Lemma 4). It is not random, but it is fixed (and the results depend on the choice of this parameter). Specifically, in Eq. 2 $h$ is fixed. In Algorithm 1 we consider a sequence of $h_k$ where $k$ is the iteration counter, and for every $k=0, 1, 2, \\cdots$,  we compute\n\\begin{equation}\n\tg_k(x_k) := \\frac{d}{\\ell} \\sum\\limits\\_{i = 1}^\\ell \\frac{f(x_k + h_k G_k e_i) - f(x_k - h_k G_k e_i)}{2h_k} G_k e_i.\n\\end{equation}\nAccording to the theoretical result, the best choice of $\\ell$  is $\\ell =d$ because it let reduce or remove (in non-smooth setting) the dependence on the dimension in the variance upper bound (see Lemma 4 in Appendix). However, a sequence of $\\ell$ can be considered and it can be useful in pratical scenarios where a budget of function evaluations is provided (and it can be \"derived\" by considering the cost in time of the single function evaluation and how much time we want to spend to solve the optimization problem).\n\n**$G_k$ sampling.** In Algorithm 1, $G_k$ is i.i.d. uniformly sampled from $O(d)$.\n\n**Structured directions and differences with random directions.** Directions in Algorithm 1 are orthogonal in the sense that, since $G_k \\in \\mathcal{O}(d)$, it implies that $G_k^\\intercal G_k = G_k G_k^\\intercal = I$. Of course, since $\\ell \\leq d$ we have to truncate it. The main difference between using orthogonal directions and random directions is in the variance of the estimator obtained as we indicated in line 210-212. Moreover, in previous works (see e.g. [1, 2, 3, 4]), authors showed that gradient estimators built  with structured directions (in particular orthogonal directions) provide better performance than the ones built using random directions. In particular, in [2] the authors empirically show that to obtain a gradient accuracy comparable to methods that use orthogonal directions, other methods (e.g. random Gaussian or spherical) can require significantly more samples. Informally, such an improvement can be justified by noting that structured directions provide a better local exploration of the space than non-structured ones, reducing the probability to generate bad directions.\n\n**Other zeroth-order method.** We thank the reviewer to indicate these zeroth-order methods. We will include some references in the \"Related Work\" section.\n\n**References**\n1. K. Choromanski, M. Rowland, V. Sindhwani, R. Turner, and A. Weller. Structured evolution with compact architectures for scalable policy optimization.\n2. A. S. Berahas, L. Cao, K. Choromanski, and K. Scheinberg. A theoretical and empirical comparison of gradient approximations in derivative-free optimization.\n3. M. Rando, C. Molinari, S. Villa, and L. Rosasco. Stochastic zeroth order descent with structured directions.\n4. D. Kozak, C. Molinari, L. Rosasco, L. Tenorio, and S. Villa. Zeroth order optimization with orthogonal random directions\n"}, {"review_id": "thbSC5W9I8", "reviewer": "Reviewer_RtjB", "summary": "This paper presents a study on a structured finite-difference algorithm for non-smooth black box optimization. The authors successfully demonstrate that their finite-difference surrogate serves as an unbiased estimator of the gradient for the smoothing approximation of the target function. The proposed O-ZD method's convergence analysis is established under different assumptions.", "questions": "See above", "limitations": "not applicable", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Strengths And Weaknesses:\nOverall, the paper is well-structured and provides valuable insights into the structured finite-difference method for non-smooth black box optimization. The authors establish the optimal complexity in the non-smooth convex case and demonstrate convergence rates in the non-smooth nonconvex and smoothing settings. However, it is unclear to me what is the exact theoretical complexity improvement of this method compared to state-of-the-art methods in different scenarios. Additionally, the objective function used in the numerical experiments appears to be relatively simple.\n\nFurthermore, I kindly request the authors to address the following questions for better understanding:\n\n1. Could the authors please provide the exact formula of C in the Corollary 1? It seems that C should be somehow related to \\theta.\n\n2. In line 205, the authors claim that the complexity in terms of the number of iterations is better than [15] and [37]. However, I believe it would be more appropriate to compare the number of function evaluations since the iterations of O-ZD are more computationally expensive.\n\n3. I suggest the authors include a table explicitly listing the number of function evaluations' complexity of this work and the related works, rather than just providing a discussion for each case separately.", "weaknesses": "See above", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her observations and questions. We answer here to his/her questions.\n\n## Weakness \nIn the literature, different works (empirically) showed that imposing a structure on the direction matrices let zeroth order method obtain better performance. However, no non-smooth analysis was provided. The main goal of this work was to provide analysis of structured zeroth order method in non-smooth setting. Indeed, one of our main results is the Smoothing Lemma for direction matrices in $\\mathcal{O}(d)$. In terms of advantages w.r.t. the state-of-the-art is in the variance of the gradient estimator (as indicated in lines 210 to 212). Moreover, we agree with the reviewer that the objective functions considered in the numerical experiments are relatively simple. Indeed, the main contributions of this work are theoretical and the experimental part is included just to confirm the theoretical results. We will extend the experimental part by repeating the experiments on different functions (some of these experiments are included in the global response).\n\n## Questions \n* Q1: We provided the explicit formulation of the constant $C$ of Corollary $1$ (point $(i)$) in the proof (in Appendix B.2). We will write it in the main paper.\n* Q2 and Q3: In line 205 we discuss the result obtained for non-smooth convex functions. In [30] and [15, Theorem 2], the authors consider estimators built using a single directions and thus the complexity in terms of number of iterations coincide with the complexity in terms of function evaluations (i.e. we can obtain such a complexity by multiplying the complexity in terms of number of iterations by $2$ since a single iteration costs $2$ function evaluations). In [30] (more precisely in Theorem 6), to derive the complexity, authors consider a constant step-size scheme (see [30] eq. 46) and they obtain a complexity of $\\mathcal{O}(d^2 \\varepsilon^{-2})$. For our algorithm, the constant step-size scheme is considered in Corollary $1$ point $(iii)$. The complexity in terms of function evaluations is $\\mathcal{O}(d \\varepsilon^{-2})$ (since every iteration costs $2\\ell$ function evaluations, we can compute the complexity in terms of number of function evaluations by multiplying the complexity in terms of number of iterations by $2\\ell$). In [15, Theorem 2], the authors consider a stepsize sequence $\\alpha_k$ that satisfies the Monroe conditions. More precisely, they consider a stepsize sequence $\\alpha_k$ as $\\alpha (1/\\sqrt{k})$ with $\\alpha$ constant (note that it does not satisfy the Monroe conditions, however the choice $1/(k^{1/2 + \\delta})$ satisfies it with $\\delta$ arbitrary near to $0$). Using their double smoothing scheme they obtain a complexity in terms of function evaluations in the order of $\\mathcal{O}(d \\log d \\varepsilon^{-2})$ (you can obtain it upper-bounding eq. $18$ by $\\varepsilon$ and solving the inequality for $k$ i.e. searching $k$ s.t. the right part of eq. $18$ is smaller than $\\varepsilon$). To make a fair comparison, we consider the choice of parameters proposed in Corollary 1 point $(i)$ i.e. $\\alpha_k = \\alpha k^{-\\theta}$. As indicated in the discussion, in order to obtain the optimal dependence on the dimension we need to include $\\sqrt{\\ell/d}$ in the stepsize e.g. by taking $\\alpha = \\sqrt{\\ell/d}$ (line 204-205). Again the complexity in terms of function evaluations is $\\mathcal{O}(d \\varepsilon^{-2})$. \nWe will include an Appendix \"Expanded Discussion\" in which we extend the discussions below the corollaries including also the table with complexities in terms of the number of function evaluations of this work and the related works.\n\n\n"}]

Author rebuttal:
Several reviewers have requested more experiments, including nonsmooth functions, and comparison in terms of CPU time. We include some of these in the attached pdf. We run these experiments $20$ times and provide the mean and standard deviation of the results.


Meta-review: The paper achieves the optimal complexity in the non-smooth convex case and provided convergence rates for other settings. The zeroth-order estimator is novel in the context of non-smooth optimization. Some reviewers have requested more experiments, which the authors have provided during the rebuttal. I encourage the authors include these new experiments in the appendix of the updated version.
Decision: Accept (poster)