Accelerating Model-Free Optimization via Averaging of Cost Samples
Authors: Guido Carnevale
Giuseppe Notarstefano
Venue: NeurIPS 2025 poster ()
Keywords: Optimization
Model-Free Optimization
Zeroth-Order Optimization
Derivative-Free Optimization
Primary area: optimization
TL;DR: Model-free optimization methods typically use only current cost samples (e.g., one per iteration) by discarding all the past cost samples. We introduce a simple yet memory mechanism to maintain and use them until new data become available.

Abstract: Model-free optimization methods typically rely on cost samples gathered by perturbing the current solution estimate along a finite and fixed set of directions. However, at each iteration, only the current cost samples are used, while potentially informative, previously collected samples are discarded. In this work, we challenge this conventional approach by introducing a simple yet effective memory mechanism that maintains an auxiliary vector of iteratively updated cost samples. By leveraging this stored information, our method estimates descent directions through an averaging of all perturbing directions weighted by the auxiliary vector components. This results in faster convergence without increasing the number of function queries. By interpreting the resulting algorithm as a time-varying dynamical system, we are able to establish its convergence properties in the strongly convex case. In particular, by using tools from system theory based on timescale separation, we are able to guarantee a linear convergence rate toward an arbitrarily small neighborhood of the optimal solution. Numerical simulations on regression problems demonstrate that the proposed approach significantly outperforms existing model-free optimization methods.

Accelerating Model-Free Optimization via Averaging of Cost Samples
Guido Carnevale
Department of Electrical, Electronic and Information Engineering
Alma Mater Studiorum Universita' di Bologna
Bologna, Italy
guido. carnevale@unibo
Giuseppe Notarstefano
Department of Electrical, Electronic and Information Engineering
Alma Mater Studiorum Universita" Bologna
Bologna, Italy
giuseppe.notarstee www.rtteaanonuuuutnoo...o io.
Abstract
Model-free optimization methods typically rely on cost samples gathered by per- turbing the current solution estimate along finite and fixed set of directions. However, at each iteration, only the current cost samples are used, while poten- tially informative, previously collected samples are discarded. challenge this conventional approach by introducing simple yet effective memory mechanism that maintains an auxiliary vector of iteratively updated cost samples. By leveraging this stored information, our method estimates descent directions through an averaging of all perturbing directions weighted by the auxiliary vector components. This results in faster convergence without increasing the number of function queries. By interpreting the resulting algorithm as a time-varying dy- namical system, are able to establish its convergence properties in the strongly convex case. In particular, by using tools from system theory based on timescale separation, are able to guarantee a linear convergence rate toward an arbitrarily small neighborhood of the optimal solution. Numerical simulations on regres- sion problems demonstrate that the proposed approach significantly outperforms existing model-free optimization methods.
this work, we
Introduction
At the core of many machine learning and deep learning tasks lie complex optimization problems that often cannot be solved analytically. In some cases, the objective function is not even known in closed form, allowing function evaluations but not gradient computations Liu et al. [2020). The class of optimization methods that operate without analytical knowledge of the objective function (and thus without using its gradient) referred to as model-free (or black-box) optimization. This branch of schemes Conn et al. (2009), [nook et al. (2122 Larson et al. [2019] is receiving increasing attention learning-oriin -oriine- and related domains, such as signal processing Turner and Rasmussen [2022). reinforeement learning Malik al. [2020). Qian and Yu (202 1oT management Chen and Giannakis [2018), and system theory Galarza- Jimenez et al. [0022] He et al. [2023]. The most popular model-free optimization methods are the so-called zeroth-order methods Liu et al. [2020), in which the function evaluations are typically used to approximate its gradient. A distinction within this field involves the way in which the cost function is sampled, which can be done by using points that
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

are selected in a random or deterministic fashion Kiefer and Wolfowitz [952 Nesterov and Spokoiny 2017
Another distinction is between 1-point Flaxman al. [2005] Saha and Tewan 2011 Dekel
et al. 2015]
Gasnikov et al. (2017] Roy et al. [2022] and mulli-poing zeroth-order methods Agarwal
et al. 2010 Duchi et al. 2015] Lian et al. 2016]. Shamir (0177] Baalaubraaaaiian and Ghadimi [2022], based on the number of points selected at each iteration for cost evaluation. Indeed, in this field, cost function evaluations are typically expensive and, therefore, there is interest in reducing their number without compromising the convergence properties of the algorithm. Further, multi-point methods, besides being expensive, are often not practical, as they require simultaneous evaluations of the cost function at different points that are not possible in many applications, such when the environment is non-stationary Hazan and Levy 20144), Bubeck et al. 2015), Yang and Mohri (0016), Zhao et al. [021]. On the other hand, while 1-point zeroth- methods are easily implementable and significantly reduce the number of eva evaluations, they typically suffer from slow convergence Liu et al. 2020]. To improve their performance, Chen et al. [2022] show the benefits of high-pass and low-pass filters. Xiao et al. 22023 introduce a suitable checking mechanism to possibly skip the cost
function evaluation at each iteration. The works wrans/Zhagg al. [2022] 2024) use the residual between two feedback points at consecutive iterations to improve the performance of a single-point method.
Our main contribution is novel algorithmic paradigm for model-free optimization. We assume that the unknown loss function can be evaluated by perturbing the current solution estimate along a finite, fixed set of directions. In this setting, associate auxiliary variable each direction, which stores the corresponding cost sample as soon as becomes available. Rather than updating the solution estimate the descent estimate solely based the cost samples gathered in the current iteration, direction by averaging all perturbation directions weighted by the auxiliary variables. This enables faster convergence without increasing the number of cost function queries. In the strongly convex case, by assuming given gradient estimation technique and regular use of each perturbation direction, we formally show that the proposed method linearly converges to a neighborhood of the optimal solution whose radius can be made arbitrarily small by tuning the hyperparameters. Our line of proof interprets the resulting method tools based on timescale separation between the memory mechanism and the solution update.
time-varying dynamical system and leverages jystem-theoretic
The paper unfolds follows. In Section 2 we introduce the problem setup. In Section 3 present the proposed algorithm and state its convergence guarantecs. In Section 4 we analyze the proposed algorithm. Finally, in Section 5
we provide some numerical simulations regression scenarios.
Notation
The identity matrix of order n is In. The all all-zero vector in R" is denoted as O,, The vertical concatenation of column vectors u1...., Uy is CoL(EI..... uy). The diagonal matrix with diagonal entries an E R is denoted as diag(ar...... aN} E R....
Problem Setup
We consider unconstrained optimization problems of the form
min {(),
OER
where e: _ - TR. denotes the cost function. We focus on the following class of functions.
Assumption 1. The cost function E is u-strongly convex for some μt 0. Moreover, E is differentiable and its gradient is L-Lipschitz continuous for some L > 0.
By Assumption I problem ① admits unique solution, denoted by 0, E R Our goal is to devise iterative methods to address problem ([) in a model-free scenario. In detail, we assume that the function is unknown but can be evaluated at each iteration (c.g., around the current solution estimate) to approximate its gradient according to a given, generic estimation technique introduced below.
2.1 Gradient Estimation Technique
Existing derivative-frie free methods in the literature properly sample the cost function é to approximate the unavailable gradient VL. large part of these methods perform this sampling phase according to fiiite fixed set of additive perturbing (or dither) directions ed! ...., ed E R², where > 0 is a parameter tuning the perturbation amplitude, while each d E R characterizes the perturbation direction. More in detail, given generic point 0 E R", the corresponding gradient Vi(0) is estimated

---

by suitably manipulating D E IN cost samples E(O + +dd). .add)....,((6 + ed") gathered by perturbing 0 along these directions. We formalize follows the availability of a generie gradient estimation technique based on this sampling approach and we characterize its "ideal" estimation capabilities.
Assumption 2. There exist D N directions d' ,..., dd E [K" and a function of IR R" I parametrized in O such that, for all compact sets S [R,, there exists I. > 0 such that
EP=19 (e(0 + ed"), d VE(O)]I S seL,
for all 0 E S andd e
Assumption 2 characterizes the estimation capabilities of Σド=1 I (E(0 ed'), d)) in the ideal setting where D queries of can be performed simultaneously. In detail, in such an ideal case, Assumption 2 ensures that the estimation technique yields an estimate of V&(0) (see (2) whose accuracy can be made arbitrarily high by tuning the amplitude parameter c, which remains a free design choice in the algorithm. As one may expect, the perturbation directions d², their number D. the constant L, and the explicit structure of g depend on the specific gradient estimation method that is employed.
2.1.1 Examples of Gradient Estimation Techniques
We provide some explicit examples about estimation techniques widely used in the literature that satisfy Assumption 2 In the schemes based on 2n-point gradient estimators, the perturbation directions d' ae are given by the canonical basis (C (1,0,... ".....,,).,,........ 0, in I and the corresponding estimation technique g explicitly reads as
(2, = =d4))))=d/22)),,
(3)
for some R (ideally, E(O + ed), Assumption 2 If the function has Lipschitz
continuous gradient, it can be shown that this technique satisfies Assumption 2 (see, e.g.. Kiefer and Wolfowitz [1952] Tang al. (0200). Extremun-seeking methods (see, e.g., the recent review Scheinker [2024] and the popular works Wittenmark and Urquhart [1995), Teel and Popovic [2001), Choi et al. 20022, Ariyur and Krstic [2003] J, Krstić and Wang [2000]. Tan et al. 2006 use sinusoidal perturbation directions, For instance, in the scalar case n = 1, we have d = sin
Φ)
for all E (L...., D} where ó E R is phase shift. Thus, in these methods, the number of directions D E N corresponds to the period of the sinusoidal signal, while the estimation function g reads as
g (2, d) 2zd/(eD))
(4)
for some 2 e R. If é is C and with a proper choice of the sinusoidal functions' frequencies, can be shown that this technique satisfies Assumption 2 (see, e.g., [Mimmo et al. 2024 Lemma IJ).
3 Averaged Model-Free Meta-Algorithm
In this section, we develop novel method to address problem ① in model-free fashion. We propose as a meta-algorithm in the sense that we show the updates of the solution estimate based on new memory mechhniim paradigm and generic gradient estimation oracle that can be made explicit in a specific algorithm (c.g., by using one of the schemes described in Section 2.1.I
3.1 Meta-Algorithm Design
We remark that, in our model-free setting, it is not possible to simultaneously use all the perturbation directions d² ..... dD each iteration. The idea behind the most popular 1-point zeroth-order methods (see (Ceen et al. (2022) ) is instead to use only a single sample per iteration, with the rationale that, on average, the resulting descent direction approximately corresponds to the one characterized by (2) in Assumption 2 Mathematically, given the iteration index f E N. this would result in updating the current estimate 2, IR" about the solution problem (I) according the time-varying law
Zi41 2t n (f(2, + edz), dz).
(5)
where d, E R is the direction used at iteration and y > is the step size. Instead, we pursue different paradigm based on a memory mechanism that stores the queries ( +dd ).... f(zz + ed?) whenever available. Namely, we use auxiliary variables _!...., ... 2
and update them as
{(x; + ed') z
if E(xre + edl) is gathered otherwise,
2 =
(6)
3

---

for j E fl...., D). By introducing selector signal st E [) for all j !) to model
whether (z, + edd) ) gathered not at iteration t, we can equivalently rewrite (6) as
* = # + $y (e(xx, + -
Then, these variables 2 are used to approximate Σ y (l(T, + ed'), d²) (which, turn, approximates the exact gradient V((x;), see Assumption 2 and update the solution estimate via
Σ g (#.11.d)
2++1 =
The overall model-free method arising from (6 and is reported in Algorithm T
We highlight that
Algorithm Averaged Model-Free Meta-Algorithm
Initialization: Ty IR", "....., =0 fort = 0, 0,1,... do for j = 1.2.......
do (x, + ed?) is gathered, i.., of = then * = f(Et + edl)
else
2++1 = z
end if
end for
Ft+11 = Ft 7 E}= 9f (+t41)
end for
Algorithm does not necessarily increase the number of cost queries compared 1-point approaches such as 5 Algorithm uses the proxies g (= d²) to mimic multiple queries at the price of an increased memory burden due to storing the D-dimensional vector 2 == co.,2)....,,2 E RDP
2.)
3.2 Meta-Algorithm Convergence Properties
In this section, we state the convergence properties of Algorithm| [ To this end, we first introduce the following condition on the sampling sequences (ee,,,,,.... ( HEN
Assumption There exists T Tmax E such that, for all N and j {L....,D} D). there exists 7 [t,t+ +Tmxx - 1] such that sf = 1.
Assumption 3 Dimposes an essentially cyclic selection of the directions (d' HIEN. Specifically, starting from all E N. it ensures the existence of an upper bound Tmux (independent of t) on the number of iterations required to select all directions at least once. We emphasize that Assumption 3 is very general and includes, as special cases, scenarios in which fixed number of samples is collected at each iteration, ranging from single sample all D samples without loss of generality. Moreover, our setting accommodates more challenging situations involving irregular number of cost samples across iterations, including extreme cases where some iterations may involve no cost samples all.
Now, we are ready to establish a linear convergence rate for Algorithm T towards a neighborhood of the optimal solution 9, with an arbitrarily small radius p >>
Theorem 1. Let Assumptions 1 2 and 3 Then, for all (zo, 20) R x RR,, RR.p>>
p>>>>> 0, K E (0, 1), and v E (0, 20L [ L)). there exist 7y and é such that, for all & E (0, 7) and e (0, 2), holds
||.; - 041l < 3111 min min(nu,
/
Timex
|[
[CoL(s)
E(xo cd' )..., zp E(xo + ca e
T0 - 0
Toax
for all L N.
The proof of Theorem is provided in Section 4.4
More in detail, the proof of Theorem is based on interpreting Algorithm [ varying dynamical system and, then, on using system theory tools to characterize its evolution through some preparatory steps carried out in Section 4

---

Theorem ensures that for any initial conditions (Zo, zo) and desired final accuracy p. it is possible to tune A to get the desired performance. From system theory perspective, this means that (0., cO + ed ..... COL, + e))))) is a semi-globally practically exponentially stable equilibrium point (see Definition 2 in Appendix A.T of the dynamical system describing Algorithm
4 Analysis
In this section, we interpret Algorithm 1 as a dynamical system and, then, we analyze it through system theory tools grounded on timescale separation arguments. Indeed, in Section 4.1 we can interpret the algorithm as a time-varying woo-tmme-scalee system (see Appendix (A.I). namely, the feedback interconnection between fast and slow subsystem. Hence, we separately analyze the identified subsystems by focusing on the so-called boundary-layer (cf. Section 4.2) and reduced (cf. Section (4.3) systems, two auxiliary schemes associated to the fast and slow subsystems, respectively. Finally, in Section 4.4 we provide the proof of Theorem by suitably combining the results obtained in the previous sections with timescale separation and syapunov-based arguments.
Assumptions 2 and hold true throughout the whole section.
Remark 1. The key tool employed in the proof of Theorem D is (deterministic) timescale separation. In the case of stochastic direction selection (rather than the deterministic imposed by Assumption 3
our proof can be adapted by relying on stochastic timescale separation (see, e.g., Carnevale and Notarstefano (2024)). In this case, the convergence results would be stated in an almost sure sense, provided that the expected value E[s[] of each sampling selector index s) is strictly positive.
Remark 2. Our proof technique bused on timescale separation enables an efficient extension of the meta-algorithm analysis more general scenarios (e.g., the nonconvex case). In particular, promotes modularity, as allows for modifying only the analysis of the reduced system (i.e., the gradient method, see Lemma 2) to adapt the specific optimization setting. For the same reasons, this modularity paves the way for extending our meta-al vetta-algoritmm to more advanced variants using, e.g., accelerated methods Lin al. 2020 baseline in place of the gradient descent method. D
4.1 Two-Time-Scale System Interpretation
sP) } DxD, Gi: : R" x R RR,, and L R" + R. be defined as Let S, == diag(s/...... .
()(0,)) := E + (/(0 + )d)) - ) () : [e(o + edt))
6(O + edPP]
With this notation at hand, we compactly rewrite Algorithm I as the time-varying dynamical system
x1+1 = It - (x,122) 2+11 24 + S((((()) 22)
(8a)
(8b)
System (8) exhibits the following four peculiar features. First, subsystem (bb) admits equilibrium
manifold of the form = L (I2). namely, that is parametrized by the other state variable xy. Second, for arbitrarily fixed J; x, it can be shown that the equilibrium manifold
L (ar) ) is globally A. Third, the variations of xr over I can be made arbitrarily small by reducing the parameter 7. Fourth, in light of Assumption [ we note that subsystem (8a), in the ideal case in which 2, = L (xx)) would approximate a gradient descent method applied to (I) with an error characterized in 2 Neglecting this error, we observe that systems with these four features are typically referred to as two-time-scale systems in the literature and are widely studied in system theory, see the survey Abdelgalil et al. [2023) and the dedicated Appendix A. The
exponentially stable for (8b) (see Definition 2 in Appendix[A.]
key idea is that subsystem (8a) can be made arbitrarily slow via y, allowing the other subsystem (8b) to stay close to the current equilibrium L (x,) and, in turn, allowing subsystem (8a) to evolve approximately as gradient descent method. Accordingly, we refer to (8a) as the slow subsystem, while (8b) is termed the fast subsystem. To conveniently exploit this two-time-seale interpretation, we temporarily disregard the approximation error by considering an auxiliary system that serves as "nominal" version of the original system (). Essentially, this nominal system, we explicitly include an additive perturbation in the slow dynamics (a) that allows for exactly recovering the gradient
5

---

descent method applied to problem when 2 = (zz). Hence, this nominal system is given by
I++1 z₁ YG] (2,, (x,))
(9a)
2+1 = 24 + S,(C 21),
(9b)
where (1x)) := (ar, (zx)) V((x)). Finally, since we want to base our analysis on Theorem 2
(cf. Appendix A.I) about generic wwo-imee-scabee systems, we rewrite (9 to shift the equilibrium
point (cf. Definition in Appendix AIII of subsystem (9a) and match condition (36) in Theorem
To this end, we introduce the novel coordinate == 0. R" and rewrite (9) as
E1+1 = u. -G, 2n) (
(10a)
1++1 24 + ( (5x) - =t),
(10b)
where G3 R" RDD + Rd 2 R DR and éi R² + R defined
(ā) = C 0), ()) (T 0.), [ = e 0.).
(11)
With this nominal system at hand, the following sections adopt the customary approach in the analysis of two-time-scale systems (see Theorem in Appendix A.I). In detail, we separately study the stability of the fast dynamics (1Ob) and the slow one (IOa) in two idcalized scenarios giving rise to Then,
the so-called boundary-layer and reduced systems, respectively (cf. Sections 4.2 and in Section 4.4 we leverage these steps to establish the stability and convergence properties of the interconnected system (TO), for sufficiently small values of that is, when there is a sufficiently large timescale separation between the identified subsystems. Finally, we conclude the proof of Theorem by showing that the original system (8) evolves closely to the nominal system (10). with discrepancy governed by et (Ix)) and, thus, tunable by the parameter e (cf. Assumption 2
4.2 Boundary-Layer System Analysis
Now, we analyze the so-called boundary-layer system associated to (10.. which obtained by considering the fast dynamics (10b) with an arbitrarily fixed slow state ū; I E R" for all t E N (or, equivalently, by setting 7 0 in (Taa)). Hence, by using the error variable ā := 2n (T) E R
with respect to the parametrized equilibrium L (i), the boundary-layer system reads
8141 = (Ip S,)a.
(12)
The next lemma ensures that the origin a globally expomentially stable equilibrium of (1).
Lemma 1, Consider 12). Then, there exists a continuous function U RR. x - IR such that
||ā|| < U(E, 1) < T Tmax ||=||²
(13a)
U(LLD - 1) (2, t) < - ||
(13b)
((²,t) U( t) < Tmax || - - + +I
(13c)
for all ž, 5 E RR. and N.
The proof of Lemma] is provided in Appendix A.2
4.3 Reduced System Analysis
In a mirrored way, we now analyze the so-called reduced system of the interconnection [0).. which is obtained by considering the slow dynamics (IOa) with the fast state in the equilibrium manifold at each iteration t, namely, with 2, (ā,) for all E N. Therefore, the reduced system reads as
Ht 3. 7G E (ā+)) (Gx))..
(14)
In light of the definitions of C; L. and 2 (cf. 11 system (14 can be equivalently expressed as क+ς1 2,
Vi(a. A.).
(15)
Namely, the reduced system corresponds the gradient descent method applied problem I In the next lemma, we ensure that the origin is a globally exponentially stable equilibrium of (15)
Lemma 2. Consider
15 Then, for all y (0, p it holds
- - + 0.)1 || < -y 2μL
je T. ||
(16)
for all ā E R".
The proof of Lemma 2 provided in Appendix
A.3

---

4.4 Proof of Theorem 1
Now, we combine the results obtained in the previous sections to prove Theorem ! Let us first introduce == 2 (i) e.. E
COL 2) IR (x++D), and compactly rewrite system (T) as Et+1 = F(G,, t),
(17)
where, by using hybrid notation & = CØ 2). F R (n D) x N + R (n+ D), reads as
ā - (G;
2 (
() [ID - S,) - ( ( (z) +
(18)
Analogously, by using this notation, we also rewrite the original system (8) in the compact form
E++1 F(E,, t) + (Et),
(19)
in which we model the gradient estimation error term through E : ((++D) + E (a+ D)) defined as ( E (G
F. [-)((=))=((()) ( - rū; + (T)
2 (Z)))
Then, given the Lyapunov function U characterized in Lemmal
we define V (n+ x N + R as
V(g, t) U' + ||ā]|
(20)
where we use again & = (ā, 2). In light of 3a for all & E R(++D) and t E N, we have
|| < ((E S Tmax |1E112
(21)
Hence, in light of (2I). the level set Ω2 == (E E R(++)) (gg.t) S c, Vt E N] of is compact
for all e > O. We want to invoke Theorem 2 (cf. Appendix A about generic time-varying
two-time-scale systems to characterize the increment of V along the trajectories of the nominal,
interconnected system (17) rather than in the auxiliary systems [12 and (15) (cf. Lemma [2 respectively). Theorem requires (i) Lyapunov function ensuring exponential stability of the origin
for the boundary-layer system (12) (see (40).. (ii) a Lyapunov function ensuring exponential stability
of the origin for the reduced system 15 (see (41)) and (iii) that the system dynamics and equilibrium
function admit the bounds (37). The first point is achieved by using U characterized in (13) (cf.
Lemma 1.. The second one follows by setting ( || using (16) (cf. Lemma 2). and noting that, being W a quadratic function, it trivially satisfies (a and (41c). The third point is achieved since & is continuous (cf. Assumption I and ef is bounded in compact sets (cf. Assumption 2
Hence, by Theorem 2 for all g = (z. 2) E R2" TR.,, u E (0, 2/I,/(1r + L)),, and k E (0, 1), there exists >> 0 such that, for all 7 E (0,5) 7). the increment of V along the trajectories of (17) satisfies
V(F(S, (), 1) ((, S - ||ā||
<|=|?,
(22)
for all / E N. Moreover, the definition of V (cf. (20)) and the bound (3e) lead to
VE,, - (E', < (Tmux + 1) ||5 - - (IIEll + IIE'I)),
(23)
for all L E N and E. E R(n+D), Now, we consider & E Ω. for some > 0 to be defined later and add (V 1),++) 1) to the increment AV () == ((F(E, + (E) - V(E,) of
V along the trajectories of the original system 199. thus obtaining
AV(S,t) ((F(E, f), 1) ((E + ((F(s.t) + FF(S), - -VF(S))))+)
(0) < ー |ل x || (Tmax + 1) (E) (1E + E (I) (24)
where in (a)
the bound (22) along the nominal trajectory and (23). Being Ω. compact, then 2 Hence, by denoting the Lipschitz constant of L 7 (G, (i, 2 + (a)) 25 (+))) in Pe with Le (which finite since é is differentiable, see Assumption we have [E)( < y2L2 with L2 == LV1++ L2. and, thus, we can bound (24) as
[+) < L for all & COL 2) E Ωi by Assumption
AV < min / T Tmmx V(S, 1) + eeyTmmm + 1) O) + 72
(25)
By recalling the definition of Ōū L. and ēú (cf. 1II II,, and the gradient reconstruction property note that (0,, (0.))) - (0.)) vi(0.) O,, since 0, is the optimal solution to problem
7

---

By plugging this result in the definition of F (cf. (88) ) we get F( t) = On+D for all & e N. This not surprising since Ç COL(Z- 0. 2 L (x)) by definition and, thus, 5 0,+D means that the solution estimate corresponds to the optimal solution 0.. while the auxiliary variables correspond to their equilibrium values. Further, since L (see (IT) is continuous (cf. Assumption I and E; is
bounded compact sets (cf. Assumption 2.
Lipschitz continuous in & within l2c, that is
[[(({, - F('' t)ll & Lr llE - E'll,
(26)
for all E. E Qc. E N, and some finite LF > 0. Thus, for all & See we can further bound (25) as AV(G, f) v min (7v, mmax ((C, 1) + 67(7max + +1) ||Ell + veli))
(27)
We thus set y E (0, 5) and, given the desired final radius S (without loss of generality), we define : min min K/0²/ /(TT/TTmm (7 + +)L2 Then, for all é E (0, we can further bound 27 as
(28)
AV( t) - - min (72, / max tV)
(29)
The inequality 29 ensures that the level set Pe is nerward-invarian for system 19 i.c., Eto E
Ωe = & She for all f. 2 to along the trajectories of (9). To summarize, recall that this
property has been obtained for generic c, p > and by satisfying the corresponding bounds defined by 7 and Therefore, by repeating all the above steps with > 0 such that fo COL(x0 0,, 20
C (zo)) E Pec by noting that | e. || = ||ق||| S (V( in light of the definition of V (cf. (20). and by iterating inequality (29) from = O, we ensure that the origin is practically exponentially stable equilibrium (see Definition 2 in Appendix A.I) of system (99 and the proof concludes.
5 Numerical Simulations
This section numerically validates our theoretical findings with Monte Carlo simulations of = 20 trials in logistic regression scenario (ef. Section 5.1 and a ridge regression one (cf. Section 5.2).
5.1 Logistic Regression
In this first case, we consider mn N points Pl,... Pm and we use them to train linear classifier by addressing the problem
Pm. E R" with binary labels / ........................... fm (-1,1},
min OER"
E k=1
log (0T + - 11012
(30)
where C >> 0 is regularization parameter. In each trial, randomly generate synthetic dataset composed by m = 1000 points. See Appendix [A.4] for details on the dataset generation. We perform two sets of comparisons. In the first one, we consider n 10, 25, 50. The estimation technique we consider is inspired by extremum-seeking methods and is defined in 4 In detail, we consider D directions d! ,..., dD R and we generate them according to d' CoL(Sin( J é'),..., sin
)) for all j E {1...., D}. where, for each & {l...., n). we set &² x/2. 7 = n if l is even and o² = 0, 7 = D 2((-R+1)/2) if is odd. is possible to show that Assumption 2 satisfied with D 11 when n 10, D 29 when
25, and D 53 when n = 50 (see, e.g., Mimmo et al. (224)). We run Algorithm (by cyclically selecting a single direction d per iteration. Namely, for all l E N. we set #f = for j = (1. mod D) + and s! 0 for all other directions. Thus, Assumption is satisfied with T max D. We consider the dimensions = 10, 25, 50 and compare our scheme with the stae-oof-the-as 1 1-point schemes described in Appendix A.5 namely, with methods (57) Choi al. 20022, (58) Chen et [2022). and (59 Zhang et al. 2022 2024). We empirically select = 0.1 in all the schemes. We manually select the step sizes y to achieve the fastest possible convergence for each method. Table reports the chosen values.
Table 1: Step size y used for each algorithm in 1-point (left) and 2-point (right) comparisons.
n = 10
n 25
n = 50
Algorithm
0.001
0.0003
0.0001
50
200
300
1-point 57
0.0003
0.0001
0.00004
Algorithm 2 2-point (60)
0.02
0.004
0.002
1- 1-point
(8
0.003
0.0025
0.002
0.02
0.004
0.002
1 1-point
(59
0.003
0.0025
0.002
8

---

Fig. reports the average and 1-standard deviation band across the trials of the evolution over t of the distance - l. where 0, IR" the optimal solution to the trial problem. Fig. shows that our algorithm significantly outperforms methods
57
58 and (59) in terms of convergence speed.
10P
102
10
107
10
10-2
10-1
10-1
10
10-2
1.2
Itaention
Iteration
Iteration
-10F
Figure 1: Problem (30) with 10 (left), 25 (middle), and 50 (right): Monte Carlo
comparison of Algorithm I with single cost sample per iteration (blue) and the 1-point methods (77 (red), 588 (yellow), and (59) (violet).
In the second comparison, we consider 50, 200, 300, the 2-point method (60) Agarwal er al. 2010). |uchi 2uuhh. et al. (0155) Shamir [017] (cf. Appendix A.5). and equip Algorithm with the estimation technique 33. Namely, we consider D) 2n directions d': the first n correspond to the canonical basis R" while the remaining n are their negatives. In this case, at each iteration t N, we set s/. s + = for j = (t mod n) + 1 and oy = 0 for all the other directions, so that and the comparison with (60) is fair both methods take 2
Assumption is satisfied with Tax cost samples per iteration. We empirically select 0.01 in both algorithms,, while, as before, we select to achieve the fastest possible convergence of both methods, and report the chosen values in Table 1 In Fig. [ we report the average and 1-standard deviation band across the trials of the evolution the iterations of the optimality distance || 0-ll in this second set of comparisons. respectively, where 0, E TR" is the optimal solution to the problem of the trial. As before, Fig. P shows that our algorithm significantly outperforms its counterpart (60) in terms of convergence speed.
Ire Ie-66ll 6.1
10
10-3
10
10-2
10
0.5
0.0
0.0
12
Tteration
Derationn
Ereeatonn
Figure 2: Problem (30): Monte Carlo comparison of Algorithm with cost samples per iteration (blue) and the 2-point method (60) (red) with 50 (left), = 200 (middle), and n = 300 (right).
5.2 Ridge Regression
In this second case, we consider ridge regression problems described as
min DER"
Z | y|l + + 11011
(31)
where X E Rmxn is a data matrix, E Rm response vector, and C > is regularization
parameter. The first term enforces data fidelity, while the second term controls overfitting by promoting smoothness in the solution. See Appendix fo4 details on the dataset generation. We consider the two sets of comparisons performed in Section [ Hence, we first consider our method equipped with the gradient estimation technique (4) and compare it with the three 1-point methods described Section A The parameters of the gradient estimation technique (4) are the same used in Section 5
We empirically select 0.1 in all the algorithms. We manually select the step sizes to achieve the fastest possible convergence for each method, and report the chosen values Table
Table 2: Step size y used for each algorithm in 1-point (left) and 2-point (right) comparisons.
n = 10
n = 25
n = 50
Algorithm
0.0005
0.0001
0.00007
50
200
n 300
1-point 57
0.0001
0.00005
0.00002
Algorithm 2-point 66
0.01
0.003 0.003
0.002
1 1-point
58
0.001
0.00125
0.001
0.01
0.002
1 1-point
59
0.001
0.00125
0.001

---

In Fig. B we report the average and 1-standard deviation band across the trials of the evolution over the iterations the optimality distance ||,, 8+ll in this first comparisons. where 0, R" is the optimal solution to the problem of the trial. Fig. 3 shows that our algorithm significantly outperforms the 1-point methods (57) 58) and terms of convergence speed.
10
[x 0,1
[0 0.1
10
10"4
10-1
10-3
10
10
0.5
1.5
30
104
-10"
Figure 3: Problem with n = 10 (left), n 25 (middle), and n = 50 (right); Monte Carlo
comparison of Algorithm I with single cost sample per iteration (blue) and the 1-point methods (57) (red), 58 (yellow), and (59) (violet).
In the second comparison, we consider n = 50, 200, 300, the 2-point method (60) Agarwal et al. [2010). Duchi et al.] [2015), Shamir [017] (cf. Appendix A.5 and equip Algorithm with the estimation technique (3). Namely, we consider D = 2n directions d': the first n correspond to the canonical basis of F". while the remaining are their negatives. In this case, each iteration for j = (t mod n) + and s for all the other directions, so that
t N, we set sy - Assumption is satisfied with Tmax and the comparison with (60) is fair both methods take 2 cost samples per iteration. We empirically select e = 0.1 in both algorithms, while, as before, we select to achieve the fastest possible convergence of both methods, and report the chosen values Table 2 In Fig. 4 we report the average and 1-standard deviation band across the trials of evolution over the iterations of the optimality distance ||a,, 0.|] in this second set of comparisons, respectively, where 0, E I" is the optimal solution to the problem of the trial. As before, Fig. 4 shows that our algorithm significantly outperforms its counterpart 60 in terms of convergence speed.
10
102
[2y=0.]
10-2
10"2
10-4
10-7
10-12
750
1,200
2,200
3,000
0.5
10-12
Dicnation
1004
Iteration
-10"
Figure 4: Problem (31) Monte Carlo comparison of Algorithm with cost samples per iteration
(blue) and the 2-point method (60) (red) with n 50 (left), = 200 (middle), and n = 300 (right).
Conclusions
In this paper, we introduced novel paradigm for model-free optimization based memory mecha- nisms. We started from a generic technique that estimates the gradient of the cost function at a given point using cost samples obtained by perturbing that point along finite, fixed set of directions. Our main contribution is the introduction of memory mechanism that enables the reuse of the most recent cost samples for each direction. This mechanism removes the need evaluate the cost in all directions at every iteration, while still allowing the estimation of descent directions through weighted average of all directions. As result, by using this approximate descent in the solution estimate update, the proposed approach achieves faster convergence without increasing the oracle complexity. We analyzed the algorithm using system-theoreed tools based on timescale separation and, for strongly convex problems, established convergence an arbitrarily small neighborhood of the optimal solution. Finally, numerical experiments on regression problems demonstrated that our method significantly outperforms date-of-the-att model-free optimization algorithms.
10

---

Acknowledgments and Disclosure of Funding
Work supported by Fondi PNRR Bando PE Progetto PEI1 3A-ITALY, "Made in Italy Circolare Sostenibile" Codice PE0000004, CUP: J33C220022900001
References
Mahmoud Abdelgalil, Daniel E Ochoa, and Jorge I Poveda. Multi-time scale control and optimization via averaging and singular perturbation theory: From ODEs to hybrid dynamical systems. Annual Reviews in Control, 56: 100926, 2023.
Alekh Agarwal, Ofer Dekel, and Lin Xiao. Optimal algorithms for online convex optimization with multi-point bandit feedback. In Colt, pages 28 40. Citeseer, 2010.
Kalyanmoy Ariyur and Miroslav Krstic. Real-time optimization by extremim seeking control. John Wiley & Sons, 2003.
Krishnakumar Balasubramanian and Saeed Ghadimi. Zeroth-order nonconvex stochastic optimization: Handling constraints, high dimensionality, and saddle points. Foundations of Computational Mathematics, 22(1):35-76 2022.
Sébastien Bubeck, Ofer Dekel, Tomer Koren, and Yuval Peres. Bandit convex optimization: T regret in dimension. In Conference Learning Theory, pages 266-278. PMLR, 2015.
Guido Carnevale and Giuseppe Notarstefano. Timescale separation for discrete-time nonlinear stochastic systems. IEEE Control Systems Letters, 2024.
Guido Carnevale, Nicola Bastianello, Giuseppe Notarstefano, and Ruggero Carli. ADMM-tracking gradient for distributed optimization over asynchronous and unreliable networks. IEEE Transac- tions Automatic Control, 2025.
Tianyi Chen and Georgios B Giannakis. Bandit convex optimization for scalable and dynamic iot management. IEEE Internet Things Journal, 6(1): 6(1):1776-1266, 2018.
Xin Chen, Yujie Tang, and Na Li. Improve single-poin zeroth-order optimization using high-pass and low-pass filters. International Conference on Machine Learning, pages 3603-3620, PMLR. 2022.
loon-Young Choi, Miroslav Kistic, Kartik B Ariyur, and Jin Soo Lee. Extremum seeking control for discrete-time systems. IEEE Transactions on automatic control, 47(2):318-33) -23, 2002.
Andrew Conn, Katya Scheinberg. and Luis N Vicente. Introduction to derivative-free optimization. SIAM, 2009.
Ofer Dekel, Ronen Eldan, and Tomer Koren. Bandit smooth convex optimization: Improving the bias-variance tradeoff. Advances in Neural Information Processing Systems, 28, 2015.
John C Duchi, Michael Jordan, Martin J Wainwright, and Andre Wibisono. Optimal rates for zero-order convex optimization: The power of two function evaluations. IEEE Transactions on Information Theory, 61(5):2788-280)) 2015.
Abraham D Flaxman, Adam Tauman Kalai, and H Brendan McMahan. Online convex optimization in the bandit setting: gradient descent without gradient. In Proceedings of the sixteenth annual ACM-SIAM symposium on Discrete algorithmss, pages 385-394, 2005.
Felipe Galarza-Jimenez, Jorge Poveda, and Emiliano Dall' Anese. Sliding-seeking control: Model- free optimization with safety constraints. In Learning for Dynamics and Control Conference, pages 1100-1111. PMLR, 2022.
Alexander Gasnikov, Ekaterina Krymova, Anastasia A Lagunovskaya, Iluua Usmanova, and Fedor A Fedorenko. Stochastic online optimization. single-point and multi-point non-linear multi- armed bandits. convex and strongly-convex case. Automation and remote control, 78:224-234, 2017.
11

---

Wassim M Haddad and VijaySekhar Chellaboina. Nonlinear dynamical systems and control: Lyapunov- based approach. Princeton university press, 2008.
Elad Hazan and Kfrr Levy. Bandit convex optimization: Towards tight bounds. Advances Neural Information Processing Systems, 27, 2014.
Zhiyu He, Saverio Bolognani, Jianping He, Florian Dörfler, and Xinping Guan. Model-free nonlinear feedback optimization. IEEE Transactions on Automatic Control, 69(7):455-4469, 2023.
Jack Kiefer and Jacob Wolfowitz. Stochastic estimation of the maximum regression function.
The Annals of Mathematical Statistics, pages 462-466, 1952.
Miroslav Krstié and Hsin- Wang- Stability of extremum seeking feedback for general nonlinear dynamic systems. Automatica, 36(4):595-601 601, 2000.
Jeffrey Larson, Matt Menickelly, and Stefan M Wild. Derivative-free optimization methods. Acta Numerica, 28:287-4000 2019.
Xiangru Lian, Huan Zhang, Cho-Jui Hsieh, Yijun Huang, and Ji Liu. A comprehensive linear speedup analysis for asynchronous stochastic parallel optimization from zeroth-order to first-order. Advances in neural information processing systems, 29, 2016.
Zhouchen Lin, Huan Li, and Cong Fang. Accelerated optimization for machine learning. Nature Singapore: Springer. 2020.
Sijia Liu, Pin-Yu Chen, Bhavya Kailkhura, Gaoyuan Zhang, Alfred O Hero III, and Pramod K Varshney. primer on zeroth- optimization signal processing and machine learning: Principals, recent advances, and applications. IEEE Signal Processing Magazine, 37(5):33554, 2020.
Dhruv Malik, Ashwin Pananjady, Kush Bhatia, Koulik Khamaru, Peter Bartlett, and Martin J Wainwright. Derivative-free methods for policy optimization: Guarantees for linear quadratic systems. Journal of Muchine Leurning Reseurch, 21(21):1-51, 2020.
Nicola Mimmo, Guido Carnevale, Andrea Testa, and Giuseppe Notarstefano. Extremum seeking tracking for derivative-free distributed optimization. IEEE Transactions Control of Network Systems, 2024.
Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization convex functions.
Foundations of Computational Mathematics, 17(2):227-566, 2017.
Hong Qian and Yang Yu. Derivative-free reinforcement learning: A review. Frontiers of Computer Science, 15(6):156336, 2021.
Abhishek Roy, Krishnakumar Balasubramanian, Saeed Ghadimi, and Prasant Mohapatra. Stochastic zeroth-order optimization under nonstationarity and nonconvexity. Journal of Machine Learning Research, 23(64):1 2022.
Ankan Saha and Ambuj Tewari. Improved regret guarantees for online smooth convex optimization with bandit feedback. In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 636-642. JMLR Workshop and Conference Proceedings. 2011.
Alexander Scheinker. 100 years of extremum seeking: A survey. Automatica, 161:111481, 2024.
Ohad Shamir. An optimal algorithm for bandit and zero-order convex optimization with two-point feedback. Journal of Machine Learning Research, 18(52) 1-11, 2017.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems, 25, 2012.
Ying Tan, Dragan Nešić, and Iven Mareels. On non-local stability properties of extremum seeking control. Automatica, 42(())8899990)) 903, 2006.
Yujie Tang, Junshan Zhang, and Na Li. Distributed zoro- algorithms for nonconvex multiagent optimization. IEEE Transactions Control of Network Systems, 8(1):269-228 2020.
12

---

Andrew R Teel and Dobrivoje Popovic. Solving smooth and nonsmooth multivariable extremum seeking problems by the methods of nonlinear programming. In Proceedings of the 2001 American Control Conference. (Cat. No. O1CH37148), volume 3, pages 2394-2399. IEEE, 2001.
Ryan Turner and Carl Edward Rasmussen. Model based learning of sigma points in unscented kalman filtering. Neurocomputinng, 00,47-53,20 2012.
Björn Wittemmark and Alexander Urquhart. Adaptive extremal control. In Proceedings of 1995 34th IEEE Conference Deciston and Control, volume 2, pages 1639-1644. IEEE, 1995.
Quan Xiao, Qing Ling, and Tianyi Chen. Lazy queries can reduce variance in zeroth-order optimiza- tion. IEEE Transactions on Signal Processing, 71:3695-3709, 2023.
Scott Yang and Mehryar Mohri. Optimistic bandit convex optimization. Advances in Neural Information Processing Systems, 29, 2016.
Yan Zhang, Yi Zhou, Kaiyi Ji, and Michael Zavlanos, new one-point residual-feedbacc oracle for black-box learning and control. Automatica, 136: 35:1100066 2022.
Yan Zhang, Yi Zhou, Kaiyi Ji, Yi Shen, and Michael M Zavlanos. Boosting one-point derivative-free
online optimization via residual feedback. IEEE Transactions on Automatic Control, 2024.
Peng Zhao, Guanghui Wang, Lijun Zhang, and Zhi-Hua Zhou. Bandit convex optimization in non-stationary environments. Journal of Machine Learning Research, 22(125): 2021.
13

---

A Technical Appendices and Supplementary Material
A.1 Preliminaries Two-Time-Scale Discrete-Time Systems
In this section, we provide a generic stability result for discrete-time two-time-scale systems. Before doing so, introduce some preliminaries about discrete- dynamical systems. Consider generic time-varying discrete-time dynamical system parametrized by a parameter 7 > 0 and described by
X++1 = (xt,)))), t),
(32)
where X E R is the state vector and G. R x N - R describes its time-varying dynamics. The next definition introduces the concept of equilibrium point.
Definition 1. A point X E _R" is an equilibrium point of system 32) if
X = (x.t).
for any t E N and y > O.
The following definition introduces the concepts of global and semi remi-global exponential (or geometric) stability of an equilibrium point. The key distinction between the two lies in the dependence of the bound on the parameter (see system (32) with respect to the initial conditions. In the global case, there exists a bound on y that is uniform over all initial conditions for which stability is achieved. In contrast, in the semi-global case, for each initial condition, there exists (possibly different) bound on that ensures stability. The definition also introduces the concept of practical stability, which allows for arbitrarily small bound on the asymptotic distance from a given point.
Definition 2. The equilibrium point X is said to be globally exponentially (or geometrically) stable for 32 if there exists 7 >> such that, for all 7 (0,5). the trajectories of (32) satisfy
Iixe xl S a1 Ilxo - xxl-xxll,s XIl
(33)
for all Xo E IR" and some a1 > 02 E (0, 1). The equilibrium point X is said to be semi aoni-globally exponentially (or geometrically) stable if. for all Xo E RR" there exist 5, bi, by > 0 such that, for
all n E (0.7). 7) the exponential decay (33) holds true for some a1 > O, a2 E (0, 1). Finally, a poins
Xp E IR" is said to be semi- mlmi-globllly exponentially (or geometrically) practically stable if for all
Xo E R" and > 0, there exist 7 > 0 such that, for all E (0, 7). it holds
1xx= - Xpll S an lixo Xpll +
(34)
for some a1 0. 02 (0.
A popular tool for establishing the stability properties of an equilibrium point is the so-called Lyapunov approach. Essentially, it consists in showing the decrease of a so-called Lyapunov function V: R. x N + IR along the trajectories the system, namely. (X,xX,,)) (), + V(x, () S 0 for all Rn and / E N (see, e.g.. (Haddad and Chellaboina] 2008 Ch. 13) for a detailed discussion). The following theorem applics the Lyapunov approach to establish semi-global exponential stability properties for generic discrete-time two-time two-time-scale system. Such stability result extends [Carnevale] ct al. 2025 Th3] as it weakens the assumption of global Lipschitz continuity on the subsystems vector fields and the equilibrium function, requiring only condition (37 to hold in compact sets. As consequence, the next theorem establishes semi-global exponential stability for the interconnected system (35), in contrast to the global result presented in [Carnevale al. 2025] Th. II.3). Namely, the bound on 7 under which the stability of the interconnected system is guaranteed is no longer uniform over all initial conditions (see Definition 2 To the best of the authors' knowledge, this result is novel and, thus, represents a side contribution of this paper.
Theorem Maniillobal exponential stability for Lime-varying wo-timee-cale systems). Consider the tame-varying interconnected system
T++1 Iz + (f 2t) 2e41 = g(zt, Ne. t),
(35a)
(35b)
with in E DCC R",
21 E R.", f: TM x
+ D. D.g. : Rm x Rd X Rm and y > 0. Assume
that there exists Zeq Rn + Rmm such that for all x D it holds
O, = [O0m Zeg (0n) Zoq (x) = g (Zeg (x)), x t),
(36a)
(36b)
14

---

all E N. Moreover, assume that any compact set S C D x R the constants L,, Log Leg > 0 defined as
V V1f(I, 2 Zeg (z).
L (2,E)ES,1EN
sup
(2f 2 (z)
(37a)
V1g(ā+ + 2z)
Lg =
sup
V2g Zep (x),
(37b)
(2,
Ley sup (2,2)E8.tEN
V (I + (z 2eg
(37c)
exist and are finite. Let be the reduced system and be the boundary layer system with &, E Rm
J(+1 = Ft + (3 Zeg (3r))))
(38)
44 = g( + Zeg (z) xr 1) - Zegl (x)
(39)
Assume that there exists a contimuous funtinuonsuuuctiion U : R x N- + R such that
bi || U(, () < b2 || Zay (x), U(,tt) < b ||
(40a)
U(g(2+
2ng2) (2)
(40b)
- U('', - bu 1| E'11 (1 + [['1I),
(40c)
for all , E R", x E RR", EE N, and some bj. by, b3. by > O. Further, assume there exist a continuous function W D R and 71- > such that, for all y E (0, 51), it holds
ci ||||| < (x, t) $ 02 ||=11
(41a)
W(x + 7f(a, (1) t) + ((x,t) < -yc3 ||x||
(41b)
((x ((z',t) s ca caIı== - -''1=111|ax) + ca || x'll |x'1,,
(41c)
for all x, x' DD, E N. and some C2, Cz >
Then, for all (I, D x R", 03 (0, c3), and by E (0, b3). there exists & E (0, 91) such that. for all (0, it holds
Ulg(ā - ( x, 1) - Eeg(x + +f((x, , 2 (z) (), W(x -U(2)))-W(C S --23 |x:12 ba ||
(x),
for all t E N.
Proof. Let us define the error coordinate & == 2t 20g( R and accordingly rewrite the interconnected system (35] as
1141 zt + (f Eo Zeq (zt),
(42a)
411 = g( + Zoq (2t), t) 20y( + Aic ( 21),
(42b)
where we introduce the drift function A (E = - Zey (11+1) + Zeg (xg())t. Given c >> 0, let Pe C D 1nth be the level set of the overall Lyapunov function ((,, f) +)++)
+W W (z, 1) (see 40 and 41
namely
22 = (x, E D x R U(ž,t)+ W(,, I c, V4 E).
Now, we take a generic pair (Ir 2) E Ωg. In light of (40a) and (41a). we note that Ω. is compact for all c > 0. The compactness of 92,. turn, allows to claim that, for all E N. the functions fxx ن + zey (z). t) and g( + Zeg (z), x. t) are Lipschitz continuous in their first two arguments on Qe, and that Zey + 7 (z. , 2 Zeq (z). t)) Lipschitz continuous on Ue with finite constants L,, Lg- and Leg defined as in (37) with S Ω2. With these constants at hand, we can follow the remaining steps in the proof of [Crnevale et al. 2025 Th.11.3] to show that the overall Lyapunov function U(E, () (W))W((a (ar, 1) is decreasing along the trajectories of (42) rather than along the trajectories of the boundary-layer and reduced systems, namely, the auxiliary dynamics (38) and (39). respectively. For the sake of completeness, we report all the steps of the proof as follows.
15

---

We now evaluate the increment (W ž, := W(x 7f (zg()).), 1) ( t) of
(az t) along the trajectories of subsystem 42a) and obtain
AW (z, 2,t) - ((2 +
(x,t)
+W(x + x// Aan (z), - W + -====
~ -703 ||:|| + W( (1 + 7 (1, 2 + 2 (z) t), - ( +1
ycs ||=|| + 2 y || ||||||
caac || [ 2 Zag 11 + + || (cg]((,,,
(43)
where in (a) we add and subtract the term (x + 7f Zoq (2), t), t + 1). in (b) we use 4bb to bound the difference of the first two terms, in (c) we use (clc). the Lipschitz continuity of f. and the triangle inequality. By recalling that f02,: Feg (0.). t) O, (cf. (36), ) we can write
[a 2 + Eey (x)())))) = [x, 2 + ()) s(On: zeq (O,),
Lg |||||| by || + zeq(xx) Zeq (0n)l
133
(, + Loq) ||ا||| + Ly ||
(44)
where in (a) we use the Lipschitz continuity of f and zeg- while in (b) we combine the Lipschitz continuity of Zey and the triangle inequality. With similar arguments, we can show the bound
[1, Zag (x), S L + Lan) ||=|||
(45)
By using the inequalities (44 and
45 we then bound the right-hand side of (43) according to 2g T,y || ||~||| + 7 ca 3 || + 27 ( + Trea) ||||| ||]|||
AW 2, () S 703
||:||
(a)
ca ||=||3+~* ka || + (7/1 2 k2) || ||=|||
(46)
where in (a) we introduce the constants 11 k2, k3 0 defined as
l = 20,by,
k2 = 2c, 1 (1 + Leqg); k3= cn..
We now evaluate the increment (U,E, = (g( ( + zeg (2), x, 1) - zeq (x) + Azzg (f -
Zeg ()), 1) - t) of the function U (cf. (40) along the trajectories of (4b)). obtaining
AU = U((( + (x). x. t) zeq(ar) + A A (x + 7f & (x) U(2,t)
s U(g(a + zeq (1), x, =<q(x)) + 1)
U(((= + (0), x. 1) - =rg(t), 1
Ul((= Zay (,,) =x)(s) +
+ 7 (xn (x)(s)))),,))) 1)
(
S - by || - U (g(3 - -
+ Ug(( + (zg 2,,) - (x) + Aiog (x +
1)
(9)
S by 2|
+ by A Aeeg (I + ff 2
20g (+)
+ A2n + + + z (x), t). x)
by A Tey (x + (0) [gg + (4), 2) (x)]
b || + 44
+ {
x. + 2b || (x + 7f + (z), [ (x), - Zzg (2)1
(47)
where in (a) we add and subtract the term (g (w), x. 1) ))==)(=)), + 1), in (b) we exploit (40b) to bound the first two terms, in (c) use 40c to bound the difference of the last two terms, and
16

---

in (d) we use the triangle inequality. The definition of A + (x, + zog(a)) and the Lipschitz continuity of Bey lead
| (r + (z, 2 + (x) 1), x)
x Teg | + 7f(,, 2 Rey (x), - xll
Log [ Zeq (z), t)ll
(0)
- y Log Il + Zeg 2og( t) f(0m= 2og (0n), +)
(c)
2 7 Le I. T |||||| + 7 7 T.g ||
(48)
where in (a) we use the update (42a). in (b) we add the term f(On. Feg (0,,), t) since this is zero in light of (36), and in (c) we use the triangle inequality and the Lipschitz continuity of the functions
and Zeq. Moreover, since ((2e( (z), ar t) 2eq Zeq(2) (I) for all x E 'D (cf. (36)). we can write
19(( Zeq (2), x. () Zeq (x)) 1 + Zeg 2x)
(x), x. !II
Lg |
(49)
where the inequality is due to the Lipschitz continuity of the function g. By using the inequalities and (49). we then bound (47) as
(48
AU(2,,,1)
by || + 2-yb Teyy Tog Fg Tng + Tenn) ||0||| || 2Tyb Teey T Tos ||
+ b 23 + Leg) 1 ||0|| + 27 2b L 13 - Leg) ||0||| ||
+ 7 ba 1 ||=1?
< (-b3 + 7ks + 722kr] || + 2ks ||:|| + (ryks + 72kks) ||=||| |I
(50)
where we introduce the constants
Res Leg + Leq),
kg = 2b 2by 1.3 (1 Irey);
ks = 2b 2b
Ky == by L2 Li3
lg 13 + Log)²
We now introduce the overall Lyapunov function V D R N F defined
Viz, 2,t) := (z, t) (Z, t).
By evaluating the increment (v 8, = (2,+2
x,,) =x))) (x)
A (I + (z 5 zzy (t),
(z, = AW E, + AU,,,, t) of along the
trajectories of the interconnected system (22 we can use the results (46) and (50) to write
||
A (z <
(51)
|
where we introduce the matrix 2(7) = 2(7) E IR.2 defined as
7003 - 22ks
(21(7) 0y -ko 7 (ka k7)
921 (7)
with (7) == (1 + k4) - (Ak2 + kt)).. Now, we consider the constants ča, by introduced in the theorem statement and we recall that 0अ E (0, ca) and by (0, by). By Sylvester criterion, we
(Ag
[703 0
know that Q >>>
if and only if
0
ba
7(33 - 23) > 22 7(s (bs ba) p(7),
(52)
where the polynomial p(>) is defined as
p(7): : =0)>===2=~~>>0 ="
1 23) (k3 + ko + (b3 - ba)ks - - -aackis-"" ~ka k7).
First, we note that the first condition in (2 is satisfied for all E (0, (ca kk)) Then, we note that p is continuous function 7 and 1 p(T)/7 = O. Hence, there exists some 7 E (0,0 7)) such that (2) is satisfied for all y (0,7) 7) and the proof concludes.
17

---

A.2 Proof of Lemmal 1
Consider the linear, time-vane time-vayying system [22 Given tu. E we introduce the so-called state transition matrix (tt, E IR px of system defined as
(o,, II (ID Sr).
(53)
We note that this transition matrix allows us to write ā, = 2o along the trajectories of system (12). Moreover, in light of the essentially-cyclic condition ensured by Assumption 3 that
we note
(to + Tmax + to) = ODxD,
(54)
for all l0,7 T E N, where ODXD denotes the zero matrix of dimension D x D. Then, we use & to explicitly define the (time-varying) Lyapunov function U F. N R as
4+
U(2,t) =
Σ
|| (TTT))=||
(55)
Being the chosen function quadratic, it trivially satisfies the conditions condition (13b), by definition (55). we note that
13a and 13c
As for
U(LD - St),,t+1)
Σ
[ f 1) - St)
H+
( + T max: 1 + 1)(TT S,) + Σ
|
-s
++Tinn
(b)
| + Tmax 1) +
Σ T=1+1
(T 1)11
(c) U(2,1)- ||
where in (a) we isolate the last term of the sum, (b) we use fact that Φ(T, / + I)
S.) = ((T,() () for all 2 + by definiton (53). in (c) we use (54) to cancel out the first term and the definiton of U (cf. (55) manipulate the second term 2rat+1 [T+(
Taa [T, t)]l? (t, E| (žītt) t) - || and the proof concludes.
Lrst
A.3 Proof of Lemma
The proof is straightforward since the reduced system 15 is the gradient method applied to prob- lem (T). We report here for completeness. We consider the increment || >V + - and, by expanding the square norm, we obtain
|| - -V 0.)11 ||*1|
= 2-V + 0.) + Vi( + 01)
22 (Vl(ā + 0.) - v((O.)) T (a 0. > +0.) V((0.)]1
(6) 2pL || ( p - VV +0,) ?
~ -
2
(56)
where in (a) we use the fact that V((O.) = On, while in (b) we use the fact that is ju-strongly convex and its gradient is I.-Lipschitz continuous (ef. Assumption Tp. The proof concludes by using the condition E (0, #+L to neglect the second term in (56).
A.4 Dataset Generation for the Tests in Section 5
In this section, we describe how we generate the datasets used in the numerical simulations presented in 5cction[5]
18

---

A.4.1 Logistic Regression Dataset Generation
In the logistic regression scenario (cf. Section 5.1 for each trial and fl....,,, m). generate an imput vector pk E R" by sampling from a standard normal distribution R" and we assign the label 4 sign ( Pk + me), where the ground truth vector 0 R" is sampled from standard normal distribution in R", 77 is a zero-mean Gaussian noise with standard deviation 0.1, while sign : R {-1, denotes the sign function, i.e., sign(1)
if 1j 2 and sign(1) otherwise. Finally, we set C = 1 and, in each trial, we randomly initialize the algorithms' solution estimates by sampling 20 from zero-mean Gaussian distribution with covariance matrix 10 I,-
A.4.2 Ridge Regression Dataset Generation
In the ridge regression scenario (cf. Section 5.2 in each Monte Carlo trial, the entries of X are sampled from standard Gaussian distribution, while the ground truth vector & E R" is drawn from standard Gaussian distribution R". The response vector y is then constructed as
Xe +
where 1) E IR is randomly generated by sampling it from a zero-mean Gaussian distribution with standard deviation 0.1. Finally, we C = 1 and, each trial, we randomly initialize the algorithms' solution estimates by sampling 2o from zero- cero-mean Gaussian distribution with covariance matrix 10 In-
A.5
Description of the Algorithms Used in the Comparisons in Section 5
In this section, we provide the description of the algorithms used in the comparisons with our method performed in Section 5 We recall that the choice of the parameters yy and e for these methods is provided in Section and Appendix 5.2
A.5.1 1-point methods
We consider three 1-point methods. The first one a discrete-time extremum-seekknn method (see, e.g., Choi et al. [2002) and consists in maintaining solution estimate E 12" and a low-pass filter
2, E IR and updating them according to
2(f(2x: + ed.) - 21)
Tt-1 xz
(57a)
24 + (e(xr + edg) - 21),
(57b)
where de E R.n sinusoidal signal, namely, d, = coi((sn((5)).
sin)...,,inn(:))))))) + x²)),, where 7' and ø² E [0, 2n) for all E (1......7 n). In the simulations, all the parameters 7¹ and øi are chosen in method (see Section 5p. Second, we consider the method proposed in Chen et al. [2022) that consists in maintaining solution estimate 2e E R and a filter variable 2 E and updating them according
2141 = (1 (1--B) B)2t + E(xe + edl) - f(1-11 + edi-1)
(58a)
F1+1 Ftt 2
(58b)
where {deheN are i.i.d. random directions uniformly sampled from the sphere in R" with unitary
radius, while 3 0 is an additional tuning parameters. In all the simulations, we run (58) by
empirically setting 3 = 0.9. Finally, consider the 1-point method proposed in Zhang et al. [022 2024 and consists maintaining solution estimate &e E IR" and updating it according to
(f(x+ + ede) - ((x-11 + ed4-1))
Et+1 It
(59a)
where, also in this case, (de HeN are i.i.d. random directions uniformly sampled from the sphere in 1" with unitary radius.
19

---

A.5.2 2-point method
We consider the 2-point method proposed in Agarwal et al. (2010) Duchi et al. [2015). Shamir (0017] that consists in maintaining solution estimate 34 R and updating according to
Ft+1 It - 7 (e(xr + edg) f(ae edg))d
(60)
where (de are i.i.d random directions uniformly sampled from the sphere R" with unitary radius.
20

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims made in the abstract and the introduction accurately reflect the contributions and scope of the paper. In particular, they emphasize our paradigm shift in the design of model-free optimization algorithms. In the paper, the effectiveness of our approach is theoretically guaranteed, and its advantages over existing methods are demonstrated through experimental results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The paper discusses the fact that the memory mechanism of our approach increases the memory burden of the algorithm.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and proofs
21

---

Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer:[Yes]
Justification: Theorem is the main result of the paper and its proof is provided in Sec- ti.4[4.4] while its assumptions are reported in the statement of the theorem. More in detail, the proof reported
(oction[4.4]seee Lemma [and Lemma 2 (whose statements are reported in Appendix while proof are given in Appendix A.2 and A.3 respectively) as preparatory steps, as well as Theorem 2 i.e., a generic result about two-time-scale systems reported and proved in Appendix A.I
Guidelines:
The answer NA means that the paper does not include theoretical results, All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
The numerical experiments reported in Section 5 are carefully detailed to make them fully reproducible. In particular, the algorithms' and problem description are carefully detailed.
Justification:
are
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
22

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: We prefer keep the code and data private for now, but open sharing them upon request.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ wttppssccpppbllc/////ccc
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All the training and test detail (e.g., the number of trials of the Monte Carlo simulation) carefully reported. In particular, the paper provides all the details about the dataset generation used and the parameters of the algorithms.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core paper to a level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: The numerical results provided Section are provided in terms of the average achieved over the trials and using the standard deviation as confidence interval.
23

---

Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [No]
Justification: We believe that this detail is not relevant for the scope of the paper and, thus, we prefer omit in the main body of the paper. However, clarify here that the experiments were run by using Matlab on a single machine with 8GB of RAM and CPU cores.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics htttps://neurre http:///eeeeeee.c/cc//cc/cttttte
Answer: [NA]
the main scope of our paper is to provide a theoretical analysis of the proposed algorithm, which not related to any ethical issues. However, we aware of the NeurIPS Code of Ethics and we will make sure to comply with it in our future work.
Justification:
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).

---

10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: the main scope of our paper is to provide a theoretical analysis of the proposed optimization algorithm, which is not strictly related to any societal impact, the aim of our algorithm the same as the one already existing optimization algorithms.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: the paper poses such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: the paper does not use existing assets.
25

---

Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: the paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: the paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
26

---

Justification:
the paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: we only used LLMS to check grammar and spelling mistakes in the paper. The core method development this research does not involve LLMs as any important, original, non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM for what should or should not be described.
27


Reviews:
[{"review_id": "3SSqFjJzm5", "reviewer": "Reviewer_xBVs", "summary": "The authors propose a zeroth-order gradient method for unconstrained convex optimization when the gradient of the objective (or cost) is unavailable, but can be estimated at a point $x$ using cost evaluations at perturbations of $x$ with respect to a finite and fixed set of directions. Precisely, they assume that the gradient of the cost, $\\nabla l(x) = \\sum_{j=1}^{D}g_{\\epsilon}(l(x+\\epsilon d^j),d^j) - e_{\\epsilon}(x)$ for known directions $d^1,d^2,\\cdots,d^D$. In the setting where it is not possible to simultaneously use all the perturbation directions in each optimization step, so that cost evaluations $l(x_t+\\epsilon d^j)$ at a current solution $x_t$ is not available for $j=1,\\cdots,D$, they propose a descent method that reuses cost evaluations from previous optimization steps to estimate the gradient $\\nabla l(x_t)$ in round $t$. Then, under three main assumptions: strong-convexity, smoothness and an essentially-cyclic condition, they prove linear convergence of the last-iterate to a unique optimal solution. Their proof follows from an involved analysis of the descent method as a two-timescale dynamical system.", "questions": "### Some questions\n1. Can the authors elaborate on what they mean by “model-free” in the context of their paper?\n2. Just for clarity, can the authors also explain the essentially-cyclic condition in Assumption 3 in lay terms, and how this helps in the analysis?\n3. The algorithm is easy to follow, but I had a conflicting opinion about the definition of the auxiliary variables in Equation 7. Consider the intended descent step $x_{t+1} = x_{t}-\\gamma\\nabla l(x_{t})$. Under Assumption 2, the previous expression would be equivalent to $x_{t+1} = x_{t}-\\gamma\\bigg(\\sum_{j=1}^{D}g_{\\epsilon}(l(x_t+\\epsilon d^j),d^j) - e_{\\epsilon}(x_t)\\bigg)$. However, since the error is unknown and cost evaluations may not be available, the authors propose the update $x_{t+1} = x_{t}-\\gamma\\sum_{j=1}^{D}g_{\\epsilon}(z_t^j,d^j)$ in Equation 8 where $\\\\{ z\\_t^{j} \\\\} \\_{j=1}^{D}$ are introduced as auxiliary variables with $z_t^j$ storing the cost evaluation along direction $d^j$ if available at step $t$. My confusion lies with the fact that $z_t^j$ is defined in Equation 7 to collect cost samples at $x_{t-1}$ i.e $z_t^j=l(x_{t-1}+\\epsilon d^j)$. Meanwhile, I expected that $z_t^j=l(x_t+\\epsilon d^j)$ so that in the best case: when all the cost samples are available and Assumption 2 holds exactly (with no errors), $g_{\\epsilon}(z_t^j,d^j)=g_{\\epsilon}(l(x_t+\\epsilon d^j),d^j)$ and Equation 8 simply executes a true gradient descent update.\n\n    It is interesting that the analysis still seems to check out nonetheless. Precisely, $z_t^j=l(x_t+\\epsilon d^j)$ is the equilibrium point of the dynamic system formed by definition of the auxiliary variables in Equation 10b.\n\n4. L260: Is this expression $F(0,t)=0$ obvious. Also what does $\\xi=0$ mean? It would help if there is some distinction between scalar and vector variables in the paper, for example the real value $0$ and vector of all zeros.\n5. L535: I assume the authors also use the expression in equation 14?\n\n### Some typos\n1. L109: is this missing an index? The LHS is a vector while the RHS is a scalar value.\n2. L132: Do you mean $l(x_t,\\epsilon d^j)$ instead of $\\nabla l(x_t,\\epsilon d^j)$.\n3. Equation 8: Should be $z_t^j$ instead of $z_t$.\n4. Equation 12b: Should be $\\tilde{x}_t$ instead of $\\tilde{x}$\n5. Equation 24: $V$ is a multi-variable function by definition but is defined here with single input variables.\n6. L255: $\\xi^t$ should be $\\xi’$.\n7. Equations 35 and 36: These are exactly the same expression.\n8. L540: you mean $L$-smooth?", "limitations": "NA", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strength\n\nThe authors address an interesting problem in unconstrained convex optimization where the cost function and its gradient are unknown, however cost evaluations along a finite and fixed set of directions, when available, can be used to estimate the gradient. Moreso, their two timescale interpretation of the resulting descent method, and its analysis as a fast vs slow dynamical system, are quite interesting. Finally, from a quick glance at the technical content in the main text seem accurate.\n\n### Weakness\nThe paper is motivated by claims that multi-point zeroth other methods are too expensive, while existing single-point methods suffer from slow convergence. However, there is no mention of what kind of rates have been achieved in prior works. To properly validate the authors claim of “faster convergence” with their method, I believe it will be helpful to reference notable convergence rate results in related literature especially from Chen et al. 2022, Xiao et al 2023 and Zhang et al 2024 – which reportedly already achieve improved performance in terms of convergence of single-point methods.", "was_revised": true, "final_justification": "I maintain my positive score of the paper due to the listed strengths and satisfactory discussion with the authors. In particular, during the discussion period the authors 1) proposed to address the listed weakness with experiments to empirically compare their method with relevant related works, and 2) discovered an alternative version of the algorithm with more intuitive descent steps at little-to-no cost to the analysis.", "rebuttal": "We thank you for this exhaustive and positive feedback on our work, we are happy that you found the problem we address interesting and that you appreciated the two-timescale interpretation of our method. We really thank you for all the useful suggestions. In the following, we provide a detailed account of the changes we would make if the paper is accepted, as well as those we have already included in the revised version.\n\nWe thank you for the suggestion about comparisons with state-of-the-art algorithms. In the revised version of the paper, we added numerical experiments that show that our method sensibly outperform those proposed by Chen et al. (2022) and Zhang et al. (2024) in terms of convergence speed. We omitted the comparison with Xiao et al. (2023), as their method allows for possibly discarding a new cost sample at each iteration, which may offer an advantage in reducing the number of queries rather than accelerating convergence and, thus, our method would have unfairly outperformed it in terms of convergence speed. Although the aforementioned works only establish $\\mathcal{O}(1/\\sqrt{t})$ convergence rates, we chose not to explicitly claim the superiority of our exponential rate, as we focus on strongly convex problems, in contrast to their merely convex settings.\n\nWe thank you for the question about \"model-free\" meaning. We refer to model-free (or black-box) optimization methods as those that do not require analytical knowledge of the objective function. If the paper is accepted, we will clarify this point in the revised version.\n\nWe thank you for the question about the essentially cyclic rule. Assumption 3 ensures the existence of an upper bound $T_{\\text{max}} \\in \\mathbb{N}$ (independent of $t$) on the number of iterations required to select all directions $d^j$ at least once in the interval $[t,t+T_{\\text{max}} - 1]$.  In the analysis, this essentially-cyclic rule allows for establishing the exponential stability properties of the memory mechanism. If the paper is accepted, we will better clarify these aspects in the revised version of the paper.\n\nWe sincerely thank you for your insightful observation regarding the best-case scenario where all samples are available. Your suggestion led us to consider an alternative version of our meta-algorithm, where the solution estimates are updated according to $x_{t+1} = x_t - \\gamma \\sum_{j=1}^{D}(s^j_t g_\\epsilon(\\ell(x_t + \\epsilon d^j),d^j) + (1-s_t^j) g_\\epsilon(z_t^j,d^j))$, which employs the same cost samples as the original version, but benefits from the improvement you pointed out. Importantly, the analysis remains unchanged, as it still fits within the timescale separation framework developed for the original version. Indeed, the memory mechanism is unaffected, and the reduced system continues to correspond to the standard gradient method.\nWe have tested this alternative formulation and observed a slight performance improvement compared to the original version. A third version of the algorithm could instead use samples of the form $\\ell(x_{t+1} + \\epsilon d^j)$ instead of $\\ell(x_t + \\epsilon d^j)$ in the update of the auxiliary variables $z_t^j$.  From a technical standpoint, also this change does not affect the convergence analysis. However, the other two formulations offer the advantage of enabling a parallel update of the solution estimate $x_t$ and the auxiliary variables $z_t^j$ rather than a sequential one, which is particularly beneficial in large-scale applications where the problem dimension $n$ and/or the number of perturbation directions $D$ is large. If the paper is accepted, we will modify the algorithm design with the update described above. Once again, thank you for this valuable suggestion!\n\nThank you for the questions about $F(0_{n + D},t) = 0_{n + D}$ and $\\xi = 0_{n +D}$. First of all, in the revised version of the paper, we already followed your suggestion and used the notation $0_n$ to denote the vector of all zeros of dimension $n$. Since $\\xi$ is the stack of the error between (i) the solution estimate and the problem solution $\\theta_\\star$ and (ii) the auxiliary variables and their equilibrium values, $\\xi = 0_{n + D}$ means that the solution estimate corresponds to the problem solution $\\theta_\\star$ and that all the auxiliary variables are at their equilibrium values. Hence, the origin must be an equilibrium of the \"nominal\" system described by Eq. (19), i.e., $F(0_{n + D},t) = 0_{n +D}$.  If the paper is accepted, we will better highlight those steps to help grasping this result.\n\nWe thank you for the question about L535. We use the steps $\\sum_{\\tau=t+1}^{t+T_{\\text{max}}-1}\\|\\Phi(\\tau,t)\\tilde{z}\\|^2 = \\sum_{\\tau=t}^{t+T_{\\text{max}}-1}\\|\\Phi(\\tau,t)\\tilde{z}\\|^2 - \\|\\Phi(t,t)\\tilde{z}\\|^2 = U(\\tilde{z},t) - \\|\\tilde{z}\\|^2$, where the last equality uses the definition of $U$.\n\nWe finally thank you for your careful reading and for spotting so many typos. We accordingly fixed all the spotted typos except for the first one. Indeed, in L109, we are considering the scalar case $n=1$ and, thus, also the left-hand-side is a scalar value.\n\nThank you again for your patience and effort in reviewing our paper!"}, {"review_id": "iARSNvEaF1", "reviewer": "Reviewer_p5uW", "summary": "This paper introduces a meta-algorithm for zeroth-order optimization designed to accelerate convergence without increasing the number of function queries per iteration. The key idea is to challenge the standard practice of discarding cost function samples from previous iterations. The proposed method maintains an auxiliary vector of stored cost samples, one for each of a fixed set of $D$ perturbation directions. When a cost is evaluated for a particular direction, its corresponding entry in the memory vector is updated; otherwise, the old value is retained. The descent direction at each step is then calculated by averaging all $D$ perturbation directions, weighted by their corresponding cost samples from the memory vector.\n\nThe authors provide a rigorous theoretical analysis for the strongly convex case by interpreting the algorithm as a time-varying two-time-scale dynamical system. They prove that the method achieves a linear convergence rate to an arbitrarily small neighborhood of the optimal solution.", "questions": "Thank you for this interesting and well-written paper. I have a few questions to better understand the nuances of your work:\n\n1. On the $T_{\\max}$ Trade-off: The convergence rate in Theorem 1 includes the term $\\min{\\gamma \\mu, \\kappa}/T_{\\max}$. This seems to imply that a larger $T_{\\max}$ (resulting from using more directions $D$ with single-sample-per-iteration cyclic updates) would slow down the linear convergence factor. This appears to be in tension with the need for a sufficiently rich set of directions to accurately approximate the gradient. Could you comment on this trade-off? Is there an optimal number of directions $D$ that balances the quality of the gradient estimate against the staleness of the cost samples?\n\n2. Stochastic Sampling Analysis: You mention that Assumption 3 could be extended to a stochastic setting. Could you elaborate on the main technical challenges in doing so? In your two-time-scale analysis, the fast-layer dynamics for $z$ converge due to the guaranteed periodic updates. If the selector $s_\\tau^j$ becomes a random variable, how would you ensure the fast subsystem $z_t$ remains sufficiently close to the equilibrium $L_\\epsilon(x_t)$ for the slow dynamics to be well-approximated by the reduced system?\n\n3. Comparison to Momentum: How do you see your memory mechanism relating to momentum-based methods in ZO optimization? Momentum methods also leverage past information by maintaining a velocity vector, which is an exponential moving average of past gradient estimates. Your method, in contrast, stores raw cost samples. What are the conceptual advantages of storing costs over storing aggregated gradient estimates?", "limitations": "I can think of a few limitations:\n\n1. Reliance on Strong Convexity and Deterministic Sampling: The theoretical analysis (Theorem 1) relies on two strong assumptions. First, the objective function $l$ is assumed to be $\\mu$-strongly convex (Assumption 1). This is crucial for establishing the stability of the reduced system (Lemma 2) and ensuring that the gradient descent dynamics converge to a unique optimum. Extending the analysis to the general convex or non-convex case is a non-trivial but important next step. For non-convex functions, it is unclear if the method would be guaranteed to converge to a stationary point or how the memory mechanism would interact with local minima and saddle points. Of course, lack of strong convexity will probably not guarantee linear convergence, which is expected.\n\nSecond, the proof of linear convergence hinges on the deterministic, essentially-cyclic sampling condition (Assumption 3), which guarantees that every direction is sampled within a finite window $T_{\\max}$. This is essential for the stability of the boundary-layer system (Lemma 1), as it ensures that no component of the memory vector $z_t$ can remain stale indefinitely. Many popular ZO methods employ random sampling of perturbation directions. Adapting the analysis to such a stochastic setting would require different analytical tools, likely establishing convergence in expectation or with high probability.\n\n2. Memory and Dimensionality Trade-off: The core of the method is a memory vector $z$ of size $D$, where $D$ is the number of perturbation directions. This introduces a memory cost of $O(D)$ that is not present in standard single-point methods. For many gradient estimation techniques, $D$ scales with the problem dimension $n$ (e.g., $D=2n$ for a two-sided gradient estimator). While this is acceptable for moderately-sized problems, the memory overhead could become a practical limitation in very high-dimensional settings where $n$ is in the millions. This presents a clear trade-off: the algorithm trades $O(D)$ memory for faster convergence.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strengths:\n\n1. The idea of maintaining and reusing the most recent cost sample for each perturbation direction is simple, intuitive, and, to the best of my knowledge, novel in this context.\n\n2. The convergence analysis is a major strength. Modeling the algorithm as a two-time-scale system is an appropriate choice. The resulting guarantee of linear convergence to a tunable neighborhood (Theorem 1) is a strong theoretical result for this class of algorithms. The inclusion of the detailed derivation in the appendix (especially Theorem 2) is a valuable contribution in its own right.\n\n3. The algorithm's main promise is faster convergence without increasing function queries. The method appears easy to implement and adds only a small memory overhead ($D$ scalars).\n\nWeaknesses:\n\n1. The core idea is not too different from coordinate descent, where descent is performed in one direction at a time. In fact, a very recent work by Karandikar et al. (Revisiting Stochastic Approximation and Stochastic Gradient Descent) has strong parallels with the proposed work in a more general setting, where they also derive sufficient conditions for the convergence of zero-order SGD, wherein the stochastic gradient is computed using $2D$ function evaluations, but no gradient computations. As such the core idea does not appear to be novel.\n\n2. Assumption 3 introduces $T_{\\max}$, the maximum time between sampling any given direction. The convergence rate in Theorem 1 depends on $1/T_{\\max}$. This suggests that using more directions $D$ (which would likely increase $T_{\\max}$ under cyclic sampling) could lead to a slower convergence rate. This creates an interesting trade-off: more directions might give a better gradient estimate, but the staleness of information could slow convergence. A more explicit discussion of this trade-off would improve the paper. Also, convergence of algorithms under assumptions of a constant such as $T_{\\max}$ is widely dealt in the distributed optimization literature wherein the network is assumed to be connected over any $T_{\\max}$ time (though it may be disconnected at a given instant).\n\n3. The introduction positions the work against methods that discard samples. However, it would be valuable to contrast this \"memory of costs\" approach with other memory-based ZO methods, such as those using momentum or filter-based techniques (like the cited Chen et al. [2022]). The proposed memory is on the raw function values $l(x_t + \\epsilon d^j)$, not on the aggregated gradient estimate or the iterates themselves. Highlighting this distinction more clearly would further sharpen the paper's contribution.\n\nA minor typo: (8) should have $z_t^j$.", "was_revised": true, "final_justification": "Based on the discussion with authors, I have suitably updated my score. They seem to have addressed most of my concerns. The generality of the results is the only limiting factor in my opinion.", "rebuttal": "We are glad that you found our idea both intuitive and novel, and we truly appreciate your recognition of the rigor of our theoretical analysis. We sincerely thank you for your valuable suggestions and observations. In the following, we provide a detailed account of the changes we would make if the paper is accepted, as well as those we have already included in the revised version.\n\nWe thank you for pointing out the interesting reference Karandikar et al. (2023) that we were not aware of and that we will include in our bibliography. As mentioned in your comment too, the method in the pointed reference needs $2D$ simultaneous function evaluations per iteration and, thus, it does not employ our memory paradigm to mitigate the number of these queries. Therefore, we really thank you for pointing out this reference as it shows that the idea of using a memory mechanism can be beneficial also in more general settings, as the stochastic approximation one considered in Karandikar et al. (2023). If the paper is accepted, we will explicitly discuss this aspect in the revised version of the paper.\n\nWe thank you for raising this point about the potential tradeoff involving $T_{\\text{max}}$. We clarify that $T_{\\text{max}}$ denotes an upper bound on the number of iterations required to sample all directions at least once in the window $[t, t + T_{\\text{max}} - 1]$. While a natural scenario may involve a cyclic and single sampling of the directions, $T_{\\text{max}}$ is not necessarily equal to the number of directions $D$. Moreover, we emphasize that the formalization of our general estimation scheme (see Assumption 2) does not imply that increasing the number of directions $D$ improves gradient estimation accuracy at a given point. For example, when using the canonical basis, once all basis directions have been sampled, adding more directions does not yield a better estimate. To conclude, there are not tradeoffs since, in the absence of specific constraints, the most effective strategy is to cyclically sample the available directions, thereby minimizing $T_{\\text{max}}$ subject to the scenario’s sampling budget (e.g., two points per iteration). If the paper is accepted, we will provide some insights about this aspect in the revised version of the paper.\n\nWe thank you for the suggestions about comparisons with filter-based techniques. In the numerical simulations, we accordingly performed novel tests to show that our algorithm outperforms the method by Chen et al. (2022) in terms of convergence speed. We already included these tests in the revised version of the paper.\n\nWe thank you carefully reading our paper and spotting the typo in Equation (8).\n\nWe thank you for the question about a stochastic sampling extension. Extending Assumption 3 to a stochastic setting technically translates in considering stochastic timescale separation rather than deterministic timescale separation. More in detail, the technical assumption is that we would need a strictly positive expected value $\\mathbb{E}[s_t^j]$ of each selector $s_t^j$. With this assumption at hand, we should be able to get a stochastic counterpart of the current main result with almost sure guarantees. In other words, we should rely on the concept of almost sure convergence rather than on the concept of asymptotic stability. If the paper is accepted, we will better clarify this point in the revised version of the paper.\n\nWe sincerely thank you for your question regarding momentum-based methods and extensions to other scenarios without strong convexity. It helped us to realize that the modularity enabled by our timescale separation perspective naturally paves the way for extending our meta-algorithm to (i) more advanced variants, for example, by using accelerated methods as a baseline instead of plain gradient descent, and (ii) problems with more class of functions. Indeed, timescale separation introduces modularity in the sense that only the reduced system analysis (in the current version, the gradient descent analysis in strongly convex setups) needs to be adapted to the new baseline (for instance, an accelerated method) and setup (for instance, nonconvex functions), while the rest of the framework remains unchanged. If the paper is accept, we will explicitly mention this observation in the revised version. We agree with your observation: in the nonconvex case, we should lose linear convergence (unless one considers Polyak–Łojasiewicz functions). More specifically, convergence should be analyzed by using LaSalle’s Invariance Principle to establish convergence to the set of stationary points, rather than relying on Lyapunov-based theorems, which are suited for the stability of a single equilibrium point. If the paper is accepted, we will explicitly discuss these aspects in the revised version of the paper.\n\nWe agree with your observation on the memory tradeoff. In one sentence, our method trades memory for convergence speed.\n\nThank you again for your patience and effort in reviewing our paper!"}, {"review_id": "2KSuJ2Pyke", "reviewer": "Reviewer_VWjt", "summary": "This paper introduces a memory-based paradigm for model-free (derivative-free) optimization. The key innovation is an auxiliary memory mechanism that stores function evaluations for each perturbation direction and reuses them across iterations. The algorithm maintains memory variables that are updated when new function evaluations become available and uses weighted averages of all stored samples to estimate descent directions. The authors also provide theoretical proof of this method, then use numerical experiments on logistic regression demonstrate its performance.", "questions": "1. The experiments are limited to 53-dimensional logistic regression. Could the authors provide results on different dimensional case?\n\n2. In chapter 5, the author mentioned γ is manually selected, and also different for two method. It would be better to show how different γ affect the result and compare the results with same γ. Could the authors provide sensitivity analysis of performance to parameter choices\n\n3. Since the method needs to store a vector, can you discuss how memory efficiency and computational cost in your method?", "limitations": "No. I think that the sensitivity to parameter choices is limited in discussion.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The memory mechanism is intuitive and makes sense. This paper also use the two-timescale system framework to provide rigorous theoretical analysis.\n\nThe idea of using historical information is well-established in optimization (momentum, Adam, etc.), though the application to zeroth-order methods is novel.\n\nThis method requires Lipschitz continuous gradients and strong convexity, limiting applicability to broader function classes.\nAnd the numerical experiment only tested on logistic regression, lacking evaluation on other problem types also didn't compare with SOTA derivative-free optimization algorithms.", "was_revised": true, "final_justification": "I will keep my score after discussion with authors and other reviewers.", "rebuttal": "We are happy that you appreciated our memory mechanism as well as our theoretical analysis and we thank you for recognizing the novelty of our idea in this context.  We thank you for your useful suggestions. In the following, we provide a detailed account of the changes we would make if the paper is accepted, as well as those we have already included in the revised version.\n\n\nWe thank you for your comment about the applicability of our analysis. The main scope of the paper is the proposal of our memory paradigm and, thus, we focused on the strongly convex case for this first step. However, the modularity of our line of proof based on timescale separation allows for efficiently extending the analysis to more general scenarios as it allows for only modifying the analysis of the reduced system (i.e., the gradient method) in the considered setting without modifying the other steps in the proof. If the paper is accepted, we will explicitly discuss this aspect in the revised version of the paper.\n\nWe thank you for all your useful suggestions to improve our experiments. In response, we already improved this aspect by taking three actions. First, we have added new experiments on ridge regression scenarios. Second, we have added comparisons with two additional SOTA 1-point algorithms, bringing the total number of algorithms considered in the comparisons to four. Finally, in all scenarios, we added experiments with larger problem dimensions $n = 200$ and $n = 300$, as well as problem with smaller dimensions $n = 10$ and $n = 25$.\n\nWe thank you for the observation about the tuning of $\\gamma$. However, we highlight that we manually tuned the parameters to get the best performance for each method. Using the same $\\gamma$ for all methods would not be fair as it would not allow each method to achieve its best performance.\n\n\nWe thank you for the question about the memory efficiency and computational efficiency of our method. Our method requires storing a $D$-dimensional vector compared to standard algorithms in this field (see below Algorithm 1). As for the computational complexity, we highlight that our method improves the convergence rate without increasing the oracle complexity. If the paper is accepted, we will discuss these aspects in the revised version of the paper.\n\n\nThank you again for your patience and effort in reviewing our paper!"}, {"review_id": "ZyOy9uT5QZ", "reviewer": "Reviewer_KNsy", "summary": "This paper proposes a simple memory mechanism for single-sample zeroth-order optimization. Instead of discarding past cost evaluations, the algorithm maintains a buffer that stores the most recent cost of each predefined perturbation direction. Under strong convexity, the method converges at a linear rate to a neighborhood of the optimal solution.", "questions": "1. Can the convergence analysis be generalized to weaker conditions, for example, smooth (but not strongly) convex objectives or stochastic noise in function values? \n\n2. Would the authors consider adding stronger baselines and real-world benchmarks? In addition, given that the synthetic experiment uses a modest dimension (n=50), might larger-scale experiments be included to illustrate scalability?\n\n3. Have you explored a principled or adaptive strategy for selecting the parameters (e.g., gamma and epsilon)?", "limitations": "yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n1. Appending a simple memory mechanism to a classical single-point ZO method is both novel and easy to implement. \n\n2. Convergence is rigorously established, with explicit rates provided.\n\n\n### Weaknesses\n1. Because the analysis relies on strong convexity and noiseless evaluations, the current theory may not cover many general convex problems. \n\n2. Experiments focus on a 50-dimensional synthetic logistic-regression task and compare against only two baseline methods. Additional studies on larger dimensions, real-world black-box tasks, or stronger baselines would better demonstrate practical advantages.\n\n3. The paper reports that parameters (e.g., step size gamma) were tuned manually, which may introduce bias into the comparisons.", "was_revised": true, "final_justification": "The paper presents a simple yet novel memory mechanism for single-point ZO methods with rigorous convergence analysis. The rebuttal addressed concerns such as limited baselines and the absence of larger-scale experiments. Overall, the strengths and additional evidence outweigh the remaining limitations, and I recommend a positive score.", "rebuttal": "We are very happy that you appreciated the idea characterizing our method and that you found our convergence analysis rigorous. We thank for all your useful questions and observations which allowed us to improve the quality of our paper. In the following, we provide a detailed account of the changes we would make if the paper is accepted, as well as those we have already included in the revised version.\n\n\nWe thank you for this observation about the fact that our analysis needs strongly convex cost functions. The main scope of the paper is the proposal of our memory paradigm and, thus, we focused on the strongly convex case for this first step. However, the modularity of our line of proof based on timescale separation allows for efficiently extending the analysis to more general scenarios as it allows for only modifying the analysis of the reduced system (i.e., the gradient method) in the considered setting (e.g., the one with merely convex functions) without modifying the other steps of the analysis. If the paper is accepted, we will add a remark to discuss this aspect.\n\n\nWe thank you for your suggestions on how to improve our experimental evaluation. In response, we have already added two additional state-of-the-art 1-point algorithms for comparison, bringing the total number of competitors to four. Moreover, we introduced new experiments on ridge regression problems (see Appendix A.6 in the revised version of the paper). Finally, for all scenarios, we included tests on both larger problems (with dimensions $n = 200$ and $n = 300$) and smaller ones (with $n = 10$ and $n = 25$).\n\nWe thank you for raising the point about potential bias due to manual tuning. We manually tuned the parameters to get the best performance for each method in the sense that we tested the performance of each method for a range of values of $\\gamma$ and selected the one that achieved the best performance. This tuning procedure strongly mitigates potential biases in the comparisons.\n\nThank you for your question regarding stochastic noise. We highlight that our theoretical result is grounded in a stability analysis of the proposed algorithm. As such, preliminary robustness guarantees are inherently provided by the robustness properties associated with stability. Moreover, our system-theoretic proof lays the groundwork for deriving more detailed robustness guarantees by leveraging well-established tools from the robustness theory of dynamical systems. If the paper is accepted, we will add some insights on this aspect.\n\nWe thank you for the question about adaptive parameters. Studying algorithms that adaptively update their parameters is indeed a fascinating research direction and certainly worth exploring. However, due to its complexity and the fact that it remains a relatively open field, it lies beyond the scope of the present work.\n\nThank you again for your patience and effort in reviewing our paper!"}]

Author rebuttal:

Meta-review: The paper introduces a simple yet effective memory mechanism for model‑free (zeroth‑order) optimization that reuses past cost samples to form weighted descent directions, which allows to achieve faster convergence without increasing the number of function queries per iteration. Reviewers consistently liked the novelty and elegance of the auxiliary‐memory idea, the rigorous two‑timescale dynamical‑systems analysis yielding linear convergence for strongly convex objectives, and the modest implementation overhead (only an 
m‑dimensional vector). It seems that the authors’ rebuttal addressed all major concerns: they added experiments against additional state‑of‑the‑art baselines, expanded the empirical study to ridge regression and a broader range of problem dimensions, clarified theoretical assumptions (model‑free definition, essentially‑cyclic sampling), and proposed alternative algorithmic variants that improve practicality without altering the convergence proof. While the current theory is limited to strongly convex, the noiseless setting and the experimental suite is a little small, the combination of solid theory and clear empirical gains merits acceptance of the paper.
Decision: Accept (poster)