On the Optimal Construction of Unbiased Gradient Estimators for Zeroth-Order Optimization
Authors: Shaocong Ma
Heng Huang
Venue: NeurIPS 2025 spotlight ()
Keywords: zeroth-order optimization
stochastic optimization
stochastic gradient descent
Primary area: optimization
TL;DR: We propose a class of unbiased zeroth-order gradient estimator.

Abstract: Zeroth-order optimization (ZOO) is an important framework for stochastic optimization when gradients are unavailable or expensive to compute. A potential limitation of existing ZOO methods is the bias inherent in most gradient estimators unless the perturbation stepsize vanishes. In this paper, we overcome this biasedness issue by proposing a novel family of *unbiased* gradient estimators based solely on function evaluations. By reformulating directional derivatives as a telescoping series and sampling from carefully designed distributions, we construct estimators that eliminate bias while maintaining favorable variance. We analyze their theoretical properties, derive optimal scaling distributions and perturbation stepsizes of four specific constructions, and prove that SGD using the proposed estimators achieves optimal complexity for smooth non-convex objectives. Experiments on synthetic tasks and language model fine-tuning confirm the superior accuracy and convergence of our approach compared to standard methods.

On the Optimal Construction of Unbiased Gradient Estimators for Zeroth-Order Optimization
Shaocong Ma Department of Computer Science
Heng Huang*
Department of Computer Science
University of Maryland College Park, MD 20742, USA scma0908@umd, edu
University of Maryland College Park, MD 20742, USA heng@umd edu
Abstract
Zeroth-order optimization (ZOO) is important framework for stochastic opti- mization when gradients are unavailable or expensive to compute. A potential limitation of existing ZOO methods is the bias inherent in most gradient estimators unless the perturbation stepsize vanishes. In this paper, we overcome this biased- ness issue by proposing novel family of unbiased gradient estimators based solely on function evaluations. By reformulating directional derivatives as telescoping series and sampling from carefully designed distributions, we construct estimators that climinate bias while maintaining favorable variance. We analyze their theo- retical properties, derive optimal scaling distributions and perturbation stepsizes of four specific constructions, and prove that SGD using the proposed estimators achieves optimal complexity for smooth non-convex objectives. Experiments on synthetic tasks and language model fine-tuning confirm the superior accuracy and convergence of our approach compared to standard methods.
Introduction
In this paper, we consider the problem of zeroth-order optimization (ZOO), where our goal to solve the following stochastic optimization problem:
min f(x) == (x; zel
where f(x; is a smooth loss function evaluated on data drawn from a distribution E. In many practical scenarios, gradient information is either unavailable or prohibitively expensive to compute. Due to its versatility, ZOO has been widely adopted across various domains, including black-box adversarial attacks machine learning models [Chen et al., 2017, Kurakin et al., 2016, Papernot et al., 2017, Cai et al., 2021, Zhao et al., 2020], physics-ifformed neural networks interfacing with external PDE solvers [Shen et al., 2024, Ma et al., 2025], and reinforcement learning [Choromanski et al., 2018 Lei et al., 2022, Suh et al., 2022]. Recent research on ZOO also focuses on enhancing memory efficiency [Cai al., 2022a.b, Li et al., 2024, Sugiura and Matsutani, 2025], motivated in large part by fine-tuning large language models [Malladi et al., 2023, Zhang et al., 2024, Gautam al., 2024, Tang al., 2024, Wang et al., 2024, 2025).
Unlike first-order methods that rely on stochastic gradients vf(x; ZOO uses only function evaluations, without access gradient information. To approximate gradients, several estimators have been proposed, including the one-point estimate fx;; f(=+pe:f), [Flaxman et al., 2005, Shamir, 2013, Bach and Perchet, 2016, Nesterov and Spokoiny, 2017, Berahas et al., 2022] and
"
This work was partially supported by NSF IIS 2347592, 2348169,, DBI 2405416, CCF 2348306. CNS 2347617, RISE 2536663.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

- [Ghadimi and Lan, 2013, Duchi et al., 2015, Nesterov and Spokoiny, 2017] (see Appendix A.1 for further discussions). The random direction is typically drawn from a Gaussian or uniform spherical distribution, while alternative choices have also gained increasing attention in recent years [Ghadimi and Lan, 2013, Duchi et al., 2015, Ji et al., 2019, Sahu et al., 2019, Coope and Tappenden, 2020, Kozak et al., 2023, Rando al., 2024a,b, Ma and Huang, 2025, Mi et al., 2025].
two-point estimator (ar; E)
However, despite these advancements, a critical limitation arises zeroth-order gradient estimation; that is, all widely used gradient estimators exhibit inherent bias. Specifically, unless the perturbation step size j asymptotically tends to zero, these estimators yield persistently biased approximations of the true gradient. This inherent bias motivates a central question explored in this paper:
Q1: Is it possible design unbiased zeroth-order gradient estimator using function evaluations?
Contribution 1: In this paper, we answer Q! affirmatively. Contrary to the belief that zeroth-order gradient estimators must inherently be biased due to finite-step perturbations, we demonstrate that is indeed possible to construct unbiased gradient estimators using only function evaluations. Our key idea to express (xr) (in the deterministic setting), the directional derivative along the direction as a telescoping series:
f(x + Mn - f(x)
V v((())= = lim
#n=+0
20
f(x + (1v) [(x) +
1
[(x pativ) - Pat1v)- f(x) f(x + knet) }
Σ Pu n=1
(2)
M1
Pn
put1
An
f(x + (1v) f(x)
1
f(x + Pa410) - f(x)
f(x + pnv) f(x)
[p
+
Pn
pn+1
Un
where the perturbation stepsize pmn - 0 as - 00, the sampling distribution {pr 1 form probability distribution (that is, 0 < < for all i E N and Σ
1), and the expectation representation holds under mild regularity conditions (Proposition 2.1). This formulation allows us to reinterpret the directional derivative as an expectation over n ~ {pm)m_> enabling the construction of a unbiased gradient estimator family S (Definition 2.2):
fx)() = ChamP((n.)
where P(n, v) is an unbiased estimator of
f(x + plv) f(x)
1 Pn.
f(x Pat1v) - f(x) f(x + pnv) f(x)
+
f1
px+1
Pn
Within this framework, we propose four specific estimators, denoted as Pa-estimator for k = 1, 3, 4, corresponding to the number of function evaluations required in each estimation. To the best of our knowledge, unbiased zeroth-order gradient estimators have received little attention in prior literature. The only existing work we aware of is the four-point estimator proposed by Chen [2020], which shares the same telescoping structure and can be viewed as a special case of our Py-estimator.
Contribution 2: Building on unbiased estimator construction, conduct rigorous variance analysis on our proposed Pæ-estimators. We first present a negative result for the Pı-estimator; although it requires fewer function evaluations, it may exhibit infinite variance under certain conditions (Theorem 3.1 (a)), which aligns with the one-point estimator in the randomized smoothing [Flaxman al., 2005]. Next, we characterize the relation among the variance of the Px-estimator (k: 2, 4), the perturbation stepsize sequence ) - and the sampling distribution {p,)=- (Theorem 3.1 (b)). Identifying the optimal choice of [um functional optimization problem:
100 and {pm) leads the following non-convex
Pn Mn+1
min (P.. 1
E Em-
Pn
20
subject to 0 Pn <
2
Pn 1;
pn 10
nm
n=1

---

We present an explicit analytical solution to this optimization problem (Theorem 3.2), which reveals two key insights: (1) our constructed unbiased gradient estimators can achieve the same variance as the classical two-point estimator without introducing additional bias, leading to the best-possible complexity for SGD algorithm (Corollary 3.5); (2) broad class of sampling distributions can achieve the minimum variance, extending beyond the specific choices considered in prior work [Chen, 2020]. While our theoretical results establish strong guarantees, an important practical question remains:
Q2: Given the optimal choice of {um] and (palá-1. do the proposed unbiased estimators empirically outperform existing zeroth-order methods?
To address Q2, we empirically validate our proposed approach across both synthetic and practical tasks. On estimating the gradient of mean-square and logistic losses, our method achieves significantly lower gradient estimation error compared to standard zeroth-order methods (Section 4. Furthermore, when applied to fine-tuning large language models, the proposed estimators demonstrate faster convergence and higher final accuracy under the same number of function evaluations (Section 4.2). These results confirm the practical advantages of our unbiased construction and underscore its effectiveness in modern zeroth-order optimization tasks.
Contribution 3:
The Derivation of Unbiased Zeroth-Order Estimators
We will start from the deterministic case then turn to the stochastic case in Section 3.3. In this section, we formally derive a class of unbiased estimators for approximating the gradient V fz(z) using only function evaluations. We also provide sufficient condition under which the telescoping series in Eq. (4) admits the expectation representation. All proofs are provided in the appendix.
2.1 Telescoping Series and Expectation Representation
For fixed direction e Rdd the directional derivative of differentiable function f Rd - RR at along the direction is defined
f(x + nv) - f(x)
vof(z)= lim
Then for any decreasing sequence {um) with limm Mn derivative as the limit of a convergent sequence
0, one can express this directional
(1)+))a)-1(()) {
Nn
v.f((x= = lim f(z + unu) f(x)
Un
This convergent sequence canonically induces telescoping series with the same limit:
v.f(z) f(x + (lr) f(x)
f(x + Pa41v)
Pa41v)- - f(x)
f(x + pnv) f(x)
Σ
(4)
P1
Pn+1
pn
Next, consider a probability mass function (PMF) [pm] with Pa. > 0 for all n and Em Pa When the series in Eq. (4) is absolutely convergent, we can interpret it as an expectation over a discrete random variable That is,
x Σ
fx + m1u) - f(x)
1 ( Mat1v) - f(x)
f(x + muw) f(x)
v.f(z)
Pn
11
Pn.
pn+1
pn
f(x + patn f(x)
1
f( + un+1v) - flx)) f(x + pn - (nx)-f(z)) }
E
+
P1
Pr
Pn41
Hn.
(5)
We follow the standard definition from Spivak [2008]: A series Exe a, called convergent, if the limit of finite sum lim ΣN- on exists. A series Σx- On called absolutely convergent, if the series |anl is convergent. See the formal definition in Appendix B.I.

---

On the Role of Absolute Convergence. The absolute convergence of the series in Eq. (4) plays critical role in interpreting the telescoping series an expectation. This due to the difference between the series convergence and the existence of expectation:
The series convergence: Consider the convergent series E: Peti. To evaluate its value, can calculate the finite-sum S,, := > P,zı then we have _i=1 Pazz -imm- Sn.
The existence of expectation: Consider the random variable X with IP(X = =x) Pr for i E N. Its expectation [[ is also written as 2i;=1 Pe However, the notion of expectation must be well-defined independently of any ordering of outcomes. That is, for an arbitrary permutation a N - N, all series Σ 2 (o()) should represent the same value E[X]
As result, convergent series can yield different values depending the order of summation (this result is called the Riemann series theorem [Riemann, 1868, Spivak, 2008]); however, the outcomes of random variable requires random variable's expectation to be well-defined regardless of any such ordering. While the expectation representation has been discussed in prior work (e.g., [Chen, 2020J), the lack of attention to absolute convergence has left the conditions ensuring unbiasedness underexplored.
Due to this reason, provide the following (mild) sufficient condition for ensuring the absolute convergence with adding a slightly stronger requirement on the objective function f RR _ and the sequence [ 1m-1-
Proposition 2.1. If the second-order continuously differentiable function f I2 BR has L L-Liscchitz continuous gradient and 2n Pn < 00, then the series
f(z u1v) f(x)
f(x + pn+11) f(x) f(x + pn.) - f(x)
Σ 7=
Pm
+
Pn
pn+1
Un
is absolutely convergent and its limit is V (xr)..
2.2 The Construction of Unbiased Estimators
With the expectation representation in place, are ready to define the class of unbiased estimators explicitly.
Definition 2.2. Suppose that the function f _2d IR is continuously differentiable and {palnal
a positive sequence with lim,- An = 0 such that the telescoping series
f(x + (1v) f(x)
Σ
f(x + Pa+1v) - f(x) f(x + pnu) - fx)
41
Mn+1
Hn
absolutely convergent, the sequence {pm)== forms a PME, and the distribution over Idd Then
the family of estimators T0 := (f, Lun ( J V) denote the class of random variables such that for every P(n, v) P, satisfies
f(x + prv) f(x)
1 Pn
[(x un41o) v) - f(x) f(x + unv) ((x)
([P((n,)) v) n, w]
+
A1
Mat1
Pnn
where is sampled from V. independent with {pn)2
In the following theorem, we formally prove that our proposed class S is exactly the unbiased estimator of the gradient f(x).
Theorem 2.3 (Unbiasedness). Let 90 (f. (pati 1 V) is defined as Definition 2.2. Then, for any estimator P(n, v) 9, the following hold:
(a) E[P(n, v) v] Vx)(x); that P(n,v) is an unbiased estimator of the directional derivative vxJ(x).
(b) If the random direction is chosen independently of the sampling n ~ {palm- and satisfies Elu u I, then
En-(pa) ((nn Vf(x),
that P(n, v) is unbiased estimator of the gradient f())).

---

2.3 Specific Constructions
In this subsection, we propose four concrete constructions from the estimator family (Definition 2.2)
(1, (pa) V)
based the number of function evaluations used in estimating the gradient. These constructions are designed to explore two main aspects: (1) the trade-off between the estimator variance and number of function evaluations, allowing flexibility depending on the computational budget; and (2) a fundamental question purely driven by the theoretical interest: What is the minimum number of function evaluations required to construct unbiased gradient estimator?
This estimator corresponds to the four- four-point estimator originally proposed by Chen [2020] with slightly generalizing the choice of the perturbation stepsize sequence (M,, J and the sampling distribution (pn)n-1- For a given direction v ~ V and n ~ {Pa In- the Pa-estimator is defined as
Pa-Estimator.
Pa(n, f(r + (lv) f(x)
1 Pn
f(x (n+10) f(x) f(r + unv) f(x)
(6)
pa+1
This construction requires four function evaluations 2,, x + 111 U, and x + 1n+111 exhibiting the lowest variance and the most function evaluation counts among all members of P.
P3-Estimator. We can reduce function evaluation by introducing selection random variable U2 Uniform (f0, 13)². The estimator is then defined as
Ps(n, v) f(x + mvv) [(x) U2
1 Pn.
f(x + pa+10) - f(x)) f(x + pnv) - [(x) (1 U2).
(1-U2).
A1
pn-1
Pn
(7)
This construction randomly selects one of two pathways: With probability 1/2, it uses the first term only and requires two function evaluations at fi and z;; otherwise, uses the second term and requires three function evaluations at x + Mn, + Pa41U, and æ. This estimator maintains unbiasedness as P4, with slightly higher variance.
Pi & Pæ-Estimator: The selection random variable can be naturally extended to construct P and P2-cstimators as follows:
P2(n, f(x + mar) f(x)
[(Ug=0)
(8)
111
1
f(x + un+1e) 12) - f(x)
f(x + pnv) f(x)
+
I(u,=1)
1[
Pn
Mat1
f(x + [(,))[ f(z)
Pi(n.v)
(9)
f1
1
f(x + (m++11)) (o=1e))
f(x) ( f(x + (nt) - - (v)
+ Pn
pn-1
where Ug Uniform (10,1,2)). U, Uniform (f0, (,) and IA is the indicator function, which equals the event occurs, and 0 otherwise. of the event A. Remarkably, the construction of Pr-estimator achieves unbiasedness using only single function evaluation. However, we will show that in the next section, --estimat will have infinite variance under certain condition.
Variance Analysis of Unbiased Zeroth-Ordes Estimators
In this section, we provide theoretical analysis of the variance behavior for the unbiased estimator family 9 (f. (f. (Jt,) 1 (px) V" (Definition 2.2). While the unbiasedness has been shown
Here Uniform(A) to represent the uniform distribution over the finite compact set A.

---

in Theorem 2.3, their variances can differ dramatically depending on the estimator construction. In particular, we prove that the variance becomes unbounded (i.e., infinite) for certain constructions such as Pı-estimator. We also provide iinieevariance bounds for Pu-estimators (for l = 2, 2,3, 4) with matching the optimal variance under specific choices of (Pn) and {pa).
3.1 Theoretical Analysis
In the following result, we adopt the same condition as Proposition 2.1 to ensure the expectation representation.
Theorem 3.1. Lei p= PS, (um) (pa) V) is defined Definition 2.2. Suppose that f: R E is second-order continuously differentiable and has L-Lipschitz continuous gradient, 00 2 Pn and V is the uniform distribution over the sphere with the radius vd. Define
(pa+1= - kn)²
:== P1,
2
and
:: 1=1
Pn
Pn.
==1
Then the following statements hold:
(a) If there exists a point R.d such that the Hessian V fa)) is positive definite and f(x then the variances of the for estimating ff(z) is infinite.
(b) The variance of Pa-estimator PM(n,, v)u (k = 2, for estimating V Vf(x) is given by
L
L
L
Var(P2(n, u) Var[P4(n, o]
d²
d
d=
3
12
L
L
Var(Pa(n, u) < VarlP4 o|
d3 d'p² 22+
d
L.2d3 312 d3 d²,²
Var[Px(n,oouf) < (d (d-1)\ +
2 0
Proof. Part (a) directly follows by analyzing the tail of 1 f(=+a.)) and leveraging the curvature from positive definite Hessian. For the part (b), we simply decompose the variance of P2(n, viu and P (n u'u into the variance of estimating (1 e)v using
Var u] dE P w) Var[Pa u)
for arbitrary P := P(n, v) 3 Then we apply the second-order Taylor expansions with the mean value theorem to control the finite-differenc noise. Full details and auxiliary lemmas are provided in Appendix C.
Comparison with Existing Literature. Theorem 3.1 reveals that while P1 is unbiased. its variance can be infinite under certain conditions, making them unsuitable for SGD. In contrast, Pæ-estimator (k 2, 3, 4) offer the finite variance when { and [pm) 1 are appropriately selected. We will show it later that under the optimal setting, their variances match the optimal order of classical two-point estimators [Nesterov and Spokoiny, 2017] but with zero bias:
Var[P v] Odd f(() &d
This variance will lead to the optimal function query complexity O(4) for achieving e-accuracy in the gradient norm V fz)]] [Duchi et al., 2015].
Comparison with the Noisy Oracle Setup In our work, we consider the exact function evaluation setting with noiseless values. In this case, our variance scales as d², which is worse than the d²² of some specific biased estimators, which mitigated by choosing small enough μı; the overall sample complexity remains optimal. However, in the noisy function evaluation setting, where each function evaluation may return noisy value, smaller 1 amplifies the noise, leading to degraded performance. Several recent works have provided more refined analysis under noisy setups with improved variance behavior. Notably, Akhavan et al. [2024] demonstrated that for highly smooth functions, the & -randomization can reduce the variance scaling to d2 with achieving the improved

---

performance for highly smooth objective functions, which extends the existing e--anndomzztiion proposed by Akhavan al. [2022]. Earlier work by Gasnikov et al. [2017] analyzed the variance behavior in single-point and multi-point bandit feedback settings, and more recent developments further explore the impact first-order smoothness in noisy black-box optimization [Gasnikov et al., 2022]. Notably, all of these results achieve the optimal complexity derived by Duchi et al. [2015].
3.2
On the Optimal Choices of ( and 1
In previous section, Theorem 3.1 connects the perturbation stepsize sequence {Ma), the sampling distribution {px}, and the variance upper bounds of our constructed unbiased estimators, which has received limited discussion in the existing literature. To control the variance term, one must ensure that @ Σm=1 (#n+)-pa)* (and yp Σ Hi for Pz-estimator) is sufficiently small, This Ps
Pa
observation naturally raises the question: What are the optimal sequences {Un)a=1 and {pm)== that minimize this sum? The following theorem addresses this question:
Theorem 3.2. Let Minnn-1 be a positive, decreasing sequence with En Mn < cc, and let {pm]m mn be a PMF Denote μ := fa- Then the following statements hold:
(a) The lower bound of is given by > μ², Moreover, the equality holds if and only if pp-1
Pn.
(b) The lower bound of & is given by & 3 | Un) un) )² > 122 Moreover, the equality holds if and only if Pn Σ
1 This result characterizes the choices of k== and {pa)a=1 that minimizes e (and &o for the Pz-estimator), leading to the variance upper bound of the form:
max(Var(P22ee) u)u], Var[Pa(n, v)u], Var[Pa(n, v)e]} & O(d]V)(()) + ddp²).
Here, can always choose 1o for the Pz-estimator such that μ41 Σ Pnn nearly match the lower bound (* v)))
Sampling from the Optimal Sampling Distribution (pa 1h-1 When the perturbation stepsize sequence {pn) ) is given, sampling the corresponding optimal distribution Pa, = (4.-,-1 could be difficult; in most of cases, {Px 1a0 cannot be ready-to-use distribution naively supported by existing software. Fortunately, we can do it conversely: given an arbitrary PMF ] the perturbation stepsize takes the form
Mn PPN > n), where N
(pm 3
providing practical way to implement the unbiased zeroth-order gradient estimator. To illustrate this point, we provide two concrete examples.
Example 3.3 (Geometric Px-Estimators). We consider the geometric distribution Geom(c)
(c (0, 1)). Then Pn, = (1 (1- c) c for all n E N. We define (4 by the recursion An - Pa+1 11 Pn (1 c) n-1 Summing this relation leads to the closed-form solution
Un = p4 c n-1
and the optimal value @ = nu7. This construction recovers the geometric sampling scheme used by Chen [2020]. We call the Pa-estimator constructed on the geometric distribution as the geometric
Pa-estimator. It is easy to verify that the corresponding & is given as & = - 21 n²
(I-
Example 3.4 (Zipf's Fsttmmtmatooss. We consider the Zipf's distribution Pa [(*) for all n N, where the Riemannian zeta function defined as ((s) Σ n
Zipf 9) (s > 1). Then
We detine μt, by the recursion pn closed-form solution
Pn+1 P1P
14 द(क) _ Summing this relation leads to the
11
pn = 11
((s)
7

---

This construction also leads to the optimal value 2 = =u. When estimating the upper bound of s, we additionally assume > 3. In this case, we have ;p = Σx m= (c-1))((e) pt. The detailed calculation is put in Example C.7.
((a-2)
In both examples, start with well-known easy-to-sample distribution (pa) and calculate the associated perturbation stepsize sequence [ either analytically (Geometric Px-estimators) or iteratively (Zipf's Px-estimators). While all estimators (i.e. Pa-estimator with A: 2,3, 4) achieve the optimal variance in the order with d and μ., these examples indicate a key difference between the P2-estimator and the Pæ-estimator (for k 3, 4): the variance bound of P3- and Pa Pa-estimator is parameter-agnostic: that is, once [Px) is specified, no additional tuning of distribution parameters is required attain the optimal bound p². This distinction highlight the practical advantages of P3- and P estimators.
3.3 Convergence of SGD with Unbiased Gradient Estimators
In this subsection, consider the stochastic optimization setting described in Eq. (1), where the goal to estimate the stochastic gradient V f((;;) rather than the full gradient. Under the optimal sampling distribution pr.s_- and the corresponding perturbation stepsize sequence (Pa) [unJn-1 the convergence upper bound of SGD follows directly from standard results for general unbiased stochastic gradient methods.
Corollary 3.5 (Khaled and Richtárik (2022]). Consider the stochastic optimization problem in Eq. (1),
and suppose that the individual loss fx: {) is second-order differentiable with L-Lipschitz continuous gradient in uniformly over & E. Assume the stochastic gradient is approximated using the Pk- estimator Ps(ni,0 v) o for k = 2,3, 4. Let the SGD ueeration be defined as I++ = It - x ve) ot where 1) (0, L is the stepsize. Then the iterates satisfy the following convergence guarantee:
2 nf
min 0 -1
dn
Consequently, choosing 1) ((I/VdT) and μ ((2) yields the optimal complexity T O of having 1 <
This complexity has matched the lower bound of solving a smooth non-convex optimization problem using zeroth-order gradient-bart based method [Duchi et al., 2015] and cannot be further improved without adding additional assumptions. Though we directly apply the result from Khaled and Richtárik [2022] (which is applicable for all unbiased estimators), the zeroth-order estimation can result in an additional dependence on the dimension d: this dependence has been reflected in our upper bound.
4 Experiments
To validate our theoretical results and demonstrate the effectiveness of the proposed unbiased zeroth- order gradient estimators, conduct experiments on two settings: synthetic objectives and language model optimization. Details and hyperparameter configurations are provided in Appendix E.
4.1 Synthetic Examples
We first evaluate estimators on two classic loss functions [James et al., 2013]: the quadratic loss freg: R - IR. for linear regression and the logistic loss fca. : IRd - for binary classification.
(z) = TAATAC Az, fz)(x) 1 Σ log(1 + exp(-b, (a. xx)))
where each entry of A E Rddd is independently sampled from the uniform distribution U[- 1]. each feature vector a; e Rd is sampled from the standard normal distribution Normal(O, Ia), and b, (-1, are binary labels generated based on Bernoulli distribution with the fixed sample size n. The gradient of each objective function can be explicitly evaluated; we compare the performance of different zeroth-order gradient estimator using the Mean-Square -rrror (MSE), which is defined as
MSE(V)(<)) = (x) vf(x) -f(t))
(10)
8

---

12"
Quidiatic Fundion
Logate Fundion
12'
N/ Blu
10
256 SAl
Oerroision 10
(lorma))
Figure 1: This figure presents the MSE error of four different estimators across various dimensions d ranging from 16 4096. The left panel corresponds to the quadratic loss fryg while the right panel illustrates results for the logistic loss fels. Each box plot describes the distribution of the MSE error across 100 random trials.
We compare the accuracy of estimating the gradient of two loss functions among four different gradient estimators including Zipf's P3-estimator (Example 3,4), two-point estimator with Gaussian or uniform random perturbations, and centralized two-point estimator with uniform perturbation (the batch size of two-point estimators is adjusted to exactly 3 function evaluations). For detailed hyper-parameter setting, we put in Appendix E. Several observations can be made from the results shown in Figure 1. First, comparing the same estimator across different dimensions, the MSE error for both objective functions generally increases with the dimension d, which is expected as higher- dimensional settings pose greater estimation challenges. Second, comparing different estimators, the Zipf's P3 Pa-estimator consistently achieves lower MSE compared to others. These results collectively demonstrate the effectiveness of our proposed estimator when estimating the gradient, especially in high-dimensional settings, which will be further validate in the next experiment.
4.2 Language Model Optimization
In this section, we demonstrate the practical applicability of the unbiased gradient estimators in optimizing the deep neural network. Particularly, we apply it to the task of fine-tuning a pre-trained language model. Using zeroth-order optimization to fine-tune the LLMs has been an active research field in recent years due to its effectiveness in saving memory [Malladi et al.. 2023, Zhang et al., 2024, Gautam al., 2024, Guo et al.. 2024]: it allows for fine-tuning model parameters without requiring access to the full computational graph, which can be prohibitively large for modern language models.
Loss vs. of Updates
Loss vs. # of Function Calls
Loss
0.6
Train
0.5
0
2500 5000
10000 12500 15000 17500 20000 Number of Updates
7500
10000 20000 30000 40000 50000 60000 70000 80000 Number of Function Calls
One Side Twa-PPiit Estimator (b
ZipP's P3-tttimator
Geommtricc P3-Estimator P-S-Estimalor
Two-Side Two-hoint
Figure Comparison of training loss during fine-tuning of 1.T-1 on SST-2 using different zeroth-order gradient estimators. The right panel rescales iterations by the number of function evaluations. The unbiased Zipf's P3-, Zipf's P4-, Geometric P3-, and Geometric Pa-estimators achieve faster convergence under the same number of function evaluations.

---

We conducted experiments using the OPT-1.3b model [Zhang et al., 2022] for sentiment classification on the Stanford Sentiment Treebank (SST-2) dataset [Socher al., 2013]. To ensure fair comparison. we maintained consistent parameters across experiments: the learning rate y 10-1 and the perturbation stepsize μ 10-3 (corresponding to j4r in the proposed unbiased estimators), which is taken from Malladi et al. [2023]'s Table 7 without additional tuning. For two-point estimators, we have adjusted the batch size to align function evaluations. Detailed experimental settings provided in Appendix E. As shown in Figure 2, zeroth-order optimization using the proposed unbiased zeroth-orded estimators achieved superior performance compared to other baseline methods.
Direct Comparison to the Two-Point Estimator (b = 1) Previously, we compare our proposed methods against two-point estimators under the constraint four function evaluations. It is also interesting to consider direct comparison with classical two-point estimators using batch size of b 1, which corresponding two function evaluations. Figure 3 presents this comparison.
Loss vs. of Updates
Loss vs. # of Function Calls
Train Loss
2500 5000 7500 10000 12500 15000 17500 20000 Number of Updates
10000 20000 3000 40000 50000 cooo0 70000 80000 Number of Function Calls
Sine-Side Two-Point Estimator
First-Order
Geometric P3-Estimator
Figure 3: Comparison to the two-point estimator with b under the same setting as Figure 2. We also include the performance of the first-order Adam and SGD in the left panel.
Choosing larger batch sizes gives more accurate gradient estimates, leading to lower training loss when measured by the number of updates. However, also observe that selecting the batch as b = 1 may also present its own advantage. Therefore, choosing the batch size can be non-trivial and requires balance the variance of gradient estimation against the per-step cost.
5 Conclusion
In this work, we proposed a novel class of unbiased zeroth-order gradient estimators based on a tele- scoping series expansion of directional derivatives. We established new theoretical results, including a sufficient condition for the expectation representation (Proposition 2.1). the unbiasedness of the proposed estimators (Theorem 2.3), a variance analysis for four specific constructions (Theorem 3),,, and the characterization of the optimal sampling distribution and perturbation stepsize sequence (Theorem 3.2). We further demonstrated that SGD equipped with our estimators achieves optimal sample complexity and empirically outperforms existing mini-batch two-poin estimators. These results provide principled foundation for new class of estimators in geroth-order optimization, offering both theoretical insights and practical improvements.
References
Arya Akhavan, Evgenii Chzhen, Massimiliano Pontil, and Alexandre Tsyhakov. A gradient estimator via 11-randomization for online zero-order optimization with two point feedback. Advances in Neural Information Processing Systems, 35:7685-7696, 2022.
Arya Akhavan, Evgenii Chzhen, Massimiliano Pontil, and Alexandre Tsybakov. Gradient-free optimization of highly smooth functions: improved analysis and a new algorithm. Journal of Machine Learning Research, 25 (370):1 50, 2024.
10

---

Francis Bach and Vianney Perchet. Highly-smooth zero-th order online optimization. In Conference on Learning Theory, pages 257-283. PMLR, 2016.
Albert S Berahas, Liyuan Cao, Krzys/tof Choromanski, and Katya Scheinherg. A thcoretical and empirical comparison of gradient approximations in derivative-free optimization. Foundations of Computational Mathematics, 22(2):507-560, 2022.
HanQin Cai, Yuchen Lou, Daniel McKenzie, and Wotao Yin. A zeroth-order block coordinate descent algorithm for huge-scale black-box optimization. In International Conference on Machine Learning, pages 1 1193-1200 1203. PMIRR 2021.
HanQin Cai, Daniel McKenzie, Wotao Yin, and Zhenliang Zhang. A one-bit, comparison- based gradient estimator. Applied and Computational Harmonic Analysis, 60:242 266, 2022a.
HanQin Cai, Daniel McKenzie, Wotao Yin, and Zhenliang Zhang. Zeroth-order regularized optimization (2oro);
Approximately sparse gradients and adaptive sampling. SIAM Journal on Optimization, 32(2):687-714, 2022b.
Guanting Chen, Unbiased gradient simulation for zeroth-order optimization. In 2020 Winter Simulation
Conference (WSC), pages 2947-2959, IEEE, 2020.
Lesi Chen, Jing Xu, and Luo Luo. Faster gradient-free algorithms for nonsmooth nonconvex stochastic optimization. In International Conference on Machine Learning, pages 5219-5233. PMLR, 2023.
Pin- Yu Chen, Huan Zhang, Yash Sharma, Jinfeng Yi, and Cho-Jui Hsich. Zoo: Zeroth order optimization based black-box attacks deep neural networks without training substitute models. In Proceedings of the 10th ACM workshop on artificial intelligence and security, pages 15-26, 2017.
Ziyi Chen, Shaocong Ma, and Yi Zhou. Accelerated proximal alternating gradient- seeeentddesceet-sscceee for nonconvex minimax machine learning. In 2022 IEEE international symposium on information theory (ISIT), pages 172-677, IEEE, 2022.
Krzysztof Choromanski, Mark Rowland, Vikas Sindhwani, Richard Turner, and Adrian Weller. Structured evolution with compact architectures for scalable policy optimization. In International Conference on Machine Learning, pages 970-978. PMLR, 2018.
Ian D Coope and Rachael Tappenden. Gradient and hessian approximations derivative free optimization. arXiv preprint arXiv:2001.08355 2020.
Liyi Dai. Convergence rates of finite difference stochastic approximation algorithms. arXiv preprint arXiv:1506.0921/ 2015.
Damek Davis, Dmitriy Drusvyatskiy, Yin Tat Lee, Swati Padmanabhan, and Guanghao Ye. A gradient sampling method with complexity guarantees for lipschitz functions in high and low dimensions, Advances Neural Information Processing Systems, 35:6692-6703, 2022.
Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. Saga: A fast incremental gradient method with support for non-strongly convex composite objectives. Advances neural information processing systems, 27, 2014.
John C Duchi, Michael I 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-2806, 2015.
Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. Spider:
Near-optimal non-convex optimization stochastic path-inted integrated differential estimator. Advances in neural information processing systems, 31, 2018.
Abraham D Flaxman, Adam Tauman Kalai, and H Brendan McMahan. Online convex optimization in the bandit setting: Gradient descent without gradient. In Proceedings the sixteenth annual ACM-SIAM symposium on Discrete algorithms, pages 3944, SIAM, 2005.
Gerald B. Folland. Advanced Calculus. Prentice Hall, Upper Saddle River, NJ, 2002. ISBN 0133652652.
Alexander Gasnikov. Darina Dvinskikh, Pavel Dvurechensky, Eduard Gorbunov, Aleksander Beznosikov, and Alexander Lobanov. Randomized gradient-free methods in convex optimization. arXiv preprint arXiv:2211. 2022.
Alexander Gasnikov, Ekaterina A Krymova, Anastasia A Lagunovskaya, Iluua N Usmanova, and Fedor Fedorenko, Stochastic online optimization. single-poin and multi-point non linear multi- armed bandits. convex and coroe case. Automation and remote control, 78(2):224 234, 2017.
11

---

Tanmay Gautam. Youngsuk Park. Hao Zhou, Parameswaran Raman, and Wooseok Ha. Variance-reduced zeroth-order methods for fine-tuning language models. arXiv preprint arXiv:2404. 08080, 2024.
Saced Ghadimi and Guanghui Lan. Stochastic first-and zeroth-order methods for nonconvex stochastic program- ming. SIAM journal optimization, 23(4):2341-2366 2013.
Andreas Griewank and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiatio. SIAM, 2008
Bin Gu, Xiyuan Wei, Shangqian Gao, Ziran Xiong, Cheng Deng, and Heng Huang. Black-box reductions for zeroth-order gradient algorithms to achieve lower query complexity. Journal of Machine Learning Research, 22(170):1-47, 2021.
Wentao Guo, Jikai Long, Yimeng Zeng. Zirui Liu,, Xinyu Yang. Yide Ran, Jacob R Gardner, Osbert Bastani, Christopher De Sa, Xiaodong Yu, et al. Zeroth-order fine-tuning of ulms with extreme sparsity. arXiv preprint 00x1::2406.02911, 2024.
Feihu Huang. Bin Gu, Zhouyuan Huo, Songcan Chen, and Heng Huang. Faster gradient-free proximal stochastic methods for nonconvex nonsmooth optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33-01, pages 1503-1510. 2019.
Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani, et al. An introduction to statistical learning, volume 1 12-1. Springer, 2013.
Kaiyi Ji, Zhe Wang. Yi Zhou, and Yingbin Liang. Improved zeroth-order variance reduced algorithms and analysis for nonconvex optimization. In International conference on machine learning. pages 3100-3109. PMLR, 2019.
Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. Advances in neural information processing systems, 26, 2013.
Ahmed Khaled and Peter Richtárik. Better theory for sgd in the nonconvex world. Transactions on Machine Learning Research, 2022.
Guy Kornowski and Ohad Shamir. An algorithm with optimal dimension-deeenddnne for zero-order nonsmooth nonconvex stochastic optimization. Journal of Machine Learning Research, 25(122): 14, 2024.
David Kozak, Cesare Molinari, Lorenzo Rosasco, Luis Tenorio, and Silvia Villa. Zeroth-order optimization with orthogonal random directions. Mathematical Programming,
2023.
Alexey Kurakin, Ian Goodfellow, and Samy Bengio. Adversarial machine learning scale. arXiv preprint arXiv:6611.01236. 2016.
Yuheng Lei, Jianyu Chen, Shengbo Eben Li, and Sifa Zheng. Zeroth-order actor-critic. arXiv preprint 2xXiv:2201.12511 2022.
Zeman Li, Xinwei Zhang, Peilin Zhong, Yuan Deng, Mcisam Razaviyayn, and Vahab Mirrokni. Addax: Utilizing zeroth-order gradients improve memory efficiency and performance of sgd for fine-tuning language models. arXiv preprint arXiv:2410.06441, 2024.
Tianyi Lin, Zeyu Zheng, and Michael Jordan. Gradient-free methods for deterministic and stochastic nonsmooth nonconvex optimization. Advances in Neural Information Processing Systems, 35:26160 11-2660-226155. 2022.
Sijia Liu, Bhavya Kailkhura, Pin- Yu Chen, Paishun Ting, Shiyu Chang, and Lisa Amini Zeroth-order stochastic variance reduction for nonconvex optimization. Advances in Neural Information Processing Systems, 31, 2018.
Sijia Liu, Pin- Yu Chen, Bhavya Kailkhura, Gaoyuan Zhang, Alfred O Hero III, and Pramod K Varshney, A primer zeroth-order optimizaaion in signal processing and machine learning: Principals, recent advances, and applications. IEEE Signol Processing Magazine, 37(5):33-54 2020.
Shaocong Ma and Heng Huang. Revisiting zeroth-order optimization: Minimum-variana two-point esti- mators and directionally aligned perturbations. In The Thirteenth International Conference on Learning Representations, 2025. URL
Shaocong Ma and Yi Zhou. Understanding the impact of model incoherence on convergence of incremental sgd with random reshuffle. In International Conference Machine Learning, pages 6563-6574, PMLR, 2020.
Shaocong Ma, Ziyi Chen, Yi Zhou, Kaiyi Ji, and Yingbin Liang. Data sampling affects the complexity of online sgd over dependent data. In Uncertainty in Artificial Intelligence, pages 1296 1305. PMLR, 2022.
12

---

Shaocong Ma, James Diffenderfer. Bhavya Kailkhura. and Yi Zhou. Deep learning of pde correction and mesh adaption without automatic differentiation. Machine Learning, 114(3):61, Feb 2025. ISSN 1573-0565. doi: 020040000-0900744555500 URL https: httpss///do.. http//oii.ggg/oooooeee 10070000000000 0e0---0050-0065500 9.
Sadhika Malladi, Tianyu Gao, Eshaan Nichani, Alex Damian, Jason D Lee, Danqi Chen, and Sanjeev Arora. Fine-tuning language models with just forward passes. Advances Neural Information Processing Systems, 05553088-33753 2023.
Zhendong Mi, Qitao Tan, Xiaodong Yu, Zining Zhu, Geng Yuan, and Shaoyi Huang. Kerzoo: Kernel function informed zeroth-order optimization for accurate and accelerated IIm fine-tuning. arXiv preprint arXiv:2505. 02222500.88866 18886, 2025.
Konstantin Mishchenko, Ahmed Khaled, and Peter Richtárik. Random reshuftling: Simple analysis with vast improvements. Advances in Neural Information Processing Systems, 33: 13319-19-77300, 2020.
Yurii Nesterov and Vladimir Spokoiny. Random gradient-fres minimization of convex functions. Foundations of Computational Mathematics, 17:527-566, 2017.
Lam M Nguyen, Jie Liu, Katya Scheinberg, and Martin Takáč. Sarah: A novel method for machine learning problems using stochastic recursive gradient. In International conference on machine learning, pages 2613-2621. PMLR, 2017.
Nicolas Papernot, Patrick McDaniel, Ian Goodfellow, Somesh Jha, Z Berkay Celik, and Ananthram Swami. Practical black-box attacks against machine learning. In Proceedings of the 2017 ACM on Asia conference on computer and communications security, pages 506-519 2017.
Marco Rando, Cesare Molinari, Lorenzo Rosasco, and Silvia Villa. An optimal structured zeroth-order algorithm for non-smooth optimization. Advances in Neural Information Processing Systems, 36, 2024a.
Marco Rando, Cesare Molinari, Silvia Villa, and Lorenzo Rosasco. Stochastic zeroth order descent with structured directions. Computational Optimization and Applications, pages 1-37, 2024b.
Bernhard Riemann. Über die darstellbarkeit einer function durch eine trigonometrische reihe, Abhandlungen der Königlichen Gesellschaft der Wissenschaften zu Görringen, 13:87-132, 1868. URL https: https://eudml. org/ doc/ hoc//35759.
Anit Kumar Sahu, Manzil Zaheer, and Soummya Kar. Towards gradient free and projection free stochastic optimization. In The 22nd International Conference on Artificial Intelligence and Statistics, pages PMLR, 2019.
Ohad Shamir. On the complexity of bandit and derivative-free stochastic convex optimization. In Conference on learning theory, pages 3 24. PMLR, 2013.
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.
Qianli Shen, Yezhen Wang, Zhouhao Yang, Xiang Li, Haonan Wang, Yang Zhang, Jonathan Scarlett, Zhanxing Zhu, and Kenji Kawaguchi. Memory-efficient gradient unrolling for large-scale bi-level optimization. arXiv preprint ar 2406 14095, 2024.
Richard Socher, Alex Perelygin. Jean Wu, Jason Chuang. Christopher D. Manning, Andrew Ng, and Christopher Potts. Recursive deep models for semantic compositionality over a sentiment trechank. In Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing, pages 1631-1642, Seattle, Washington, USA, Octoher 2013. Association for Computational Linguistics. URL htclab///ww......eeebb. aclweb. org/ httaaaaooo///---1111/ 1170.
Michael Spivak. Calculus. Publish or Perish, Inc., Houston, TX, USA, 4 edition, 2008. ISBN 978-0-914098--11
Keisuke Sugiura and Hiroki Matsutani. Elasticzo: A memory-efficient on-device learning with combined zeroth-and first-order optimization. arXiv preprint arXiv:2501.0 04287, 2025.
Hyung Ju Suh, Max Simehowitz, Kaiging Zhang, and Russ Tedrake. Do differentiable simulators give better policy gradients? In International Conference on Machine Learning, pages 20668-2066 PMLR, 2022.
Xinyu Tang, Ashwinee Panda, Milad Nasr, Saced Mahloujifar, and Prateek Mittal. Private fine-tuning of large language models with zeroth-order optimization. arXiv preprint arXiv:2401.0432 04343, 2024.
Fei Wang. Li Shen, Liang Ding, Chao Xue, Ye Liu, and Changxing Ding. Simultancous computation and memory efficient zeroth-ordas optimizer for fine-tuning large language models. arXiv preprint arXiv:2410.09823. 2024.
13

---

Liangyu Wang, Jie Ren, Hang Xu, Junxiao Wang, Huanyi Xie, David E Keyes, and Di Wang. Zo2: Scalable zeroth-order fine-tuning for extremely large language models with limited gpu memory. arXiv preprint arXiv:2503. 12668, 2025.
Jingzhao Zhang. Hongzhou Lin. Stefanie Jegelka, Suvrit Sra, and Ali Jadbabaie. Complexity of finding stationary points nonconvex nonsmooth functions. International Conference on Machine Learning. pages 11173-11182. PMLR, 2020.
Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li. Xi Victoria Lin, et al. Opt: Open pre pre-trained transformer language models. arXiv preprint arw.v:2205.0111.... 2022.
Yihua Zhang, Pingzhi Li. Junyuan Hong, Jiaxiang Li. Yimeng Zhang, Wenqing Zheng. Pin- Yu Chen, Jason D. Loe, Wotao Yin, Mingyi Hong, Zhangyang Wang, Sijia Liu, and Tianlong Chen. Revisiting zeroth-order optimization for memory-effeiient IIm fine-tuning: A benchmark, 2024.
Pu Zhao, Pin- Yu Chen, Siyue Wang, and Xue Lin. Towards query query-efficient black-box adversary with seroth-order natural gradient descent. In Proceedings of the AAAI Conference Artificial Intelligence, volume 34-04, pages 6909-6916, 2020.
14

---

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]
We have included the main contributions in the abstract and introduction in bold.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
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 not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We have included separate limitation section in the appendix.
Guidelines:
The answer NA means that the paper has no 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 a separate "Limitations" section 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 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 few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low 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 of 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 used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
15

---

Justification:
We have include the full assumptions in the statement of each lemma and theorem.
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 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 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 be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have included the full source codes and reproduction instructions. For synthetic experiment, we additionally include separate self-contained Jupyter notebook.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not 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 make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others 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 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) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(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).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe 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 possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
16

---

The authors should provide scripts to reproduce all experimental results for the new reproducible, they Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We have included the full source codes and reproduction instructions. For synthetic experiment, we additionally include separate self-contained Jupyter notebook.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (uttps://nips.. cc/
Sewriides/CodeSubmiise
for more details.
While we encourage the release of code and data, we understand that this might 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: / httppsscccppbccccceeeeceeecc wwwp//uiesse///ee /CoeeS/mmdSbbiiic ssi aesionPolicyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
proposed method and baselines. only subset of experiments 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 is 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:
We have included the these details in the experiment section and the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, in appendix, or as 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: In the synthetic experiment, we ran 100 independent trails to obtain the estimation. The median and the +25% percentiles reported in the boxplot.
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, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
17

---

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 one 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 not to show in tables figures symmetric error bars that would yield results that are out range (c.g. negative emror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables 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: [Yes]
Justification:
We have provided the compute resources information in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, 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 as 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 or failed experiments that didn't 't make it into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: http://eeeipppeeeeeeeeee wc/ppbliic/thiicauuiieeee
Answer: [Yes]
Justification:
We exactly follow the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is 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: [Yes]
Justification:
We have included this discussion in the appendix.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or 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.
18

---

The conference expects that many papers will be foundational research and not tied 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 an 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 a generic algorithm for optimizing neural networks could enable people to 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 being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there 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 monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
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 no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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 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. Lieenses for existing assets
Question: 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?
Answer: [Yes]
Justification: We have clearly cited the related papers of the code, data, and models.
Guidelines:
The answer NA means that the paper does 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 (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com 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.
19

---

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: [Yes]
Justification: All codes come with self-contained names and comments. The file, README.md, is included as the documentation.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl 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
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 details about compensation (if any)?
Question:
Answer: [NA]
Justification:
The paper does not involve crowdsourcing 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 in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or 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 equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing 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.
We 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, not include any information that would break anonymity (if applicable), such the institution conducting the review.
20

---

16. Declaration of LLM usage
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, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https //tp:///ee....eeeee Seeeurips.cc/ooffeeeee /LLM) for what should or should be described.
21

---

Appendix
Table of Contents
A Additional Backgrounds
22
A.1 Gradient Estimators in Zeroth-Order Optimization
22
A.2 Zeroth-Order SGD and Its Variants
23
Discussions on the Forward Auto Differentiation (AD) Approach
24
Discussions on the Difference Between Our Results and Chen [2020]
24
B Bias Analysis
25
B.I Absolute Convergence
25
B.2 Unbiasedness of Zerath-Order Estimators in P-Family
26
C Variance Analysis
27
C.I Variance of P -Estimmttrr
27
C.2 Variance of P. P.-Fistmator
28
C,3 Variance of P--Estimator
30
C.4 Variance of Pa-Estimator
31
C.5 On the Optimal Sampling Distribution and the Perturbation Stepsize Sequence
32
C.6 Discussions: Variance of Random Directional Derivative
33
Convergence Analysis
33
D.1 The Convergence Upper Bound
33
D.2 Supporting Lemmas
34
E Experiments Details
35
E.1 Synthetic Example
35
E.2 Language Model Optimization
35
F Broader Impact
36
G Limitations
36
A Additional Backgrounds
A.1 Gradient Estimators in Zeroth-Order Optimization
One-Point Zeroth-Order Estimator One-point estimators represent the simplest class, needing only a single function query per estimate. This construction makes them suitable when queries are costly or limited, like in online settings [Flaxman et al., 2005]. A common form, motivated by Gaussian smoothing [Nesterov and Spokoiny, 2017], is
( 1 --x + ww)v,
(Single-Point)
where is often drawn from Normal(O, Ia). While the expectation EI (xx ur) approximates the gradient of the smoothed function v [Eu-Norrmat0..z f +w)]. the estimator is biased regarding the true gradient V (x:). This bias diminishes as the smoothing parameter μ - 0 [Berahas et al., 2022]. However, these estimators suffer from high variance, potentially scaling with d² and exploding μ 0 [Flaxman al., 2005].
22

---

Two-Point Zeroth-Order Estimator Two-point estimators improve on one one-poin methods by using two function evaluations, often via finite difference along a random direction [Shamir, 2017]. The standard difference form is
V 2-side ((a) f(x + pu) - f(xz pu2)
(Two-Side)
U2
2/1
flæ + pv) f(x) I-side f())
(One-Side)
0,
requiring two queries [Shamir, 2017, Nesterov and Spokoiny, 2017]. This construction approximates the directional derivative [Chen, 2020]- Their expectation exactly matches the gradient of smoothed function v=Ju (z) [Nesterov and Spokoiny, 2017] and maintains (uu)-lvvll error [Ma and Huang, 2025]. Variance is significantly lower than one-point methods, often scaling linearly with dimension d [Duchi et al., 2015, Berahas et al., 2022].
Multiple-Poin Zeroth-Order Estimator Multiple-point estimators use more than two function evaluations to further enhance gradient estimate quality. Common strategies include mriie-Diffferene method [Dai, 2015] or mini-batch averaging [Duchi al., 2015]. The finite-difference method approximates the gradient using finite differences along each standard basis direction:
f(r + nE;) f(x pe;))
FFiite-Differenne
V if (2))
C2.
2p
requiring 2d queries, where is the i-th coordinate vector. Mini-batch averaging reduces variance by averaging independent two-point estimates:
= Σ f(x + pwvz) - f(x pmv)
(Mini-Batee Bani-Batch)
V batchf(2))
U/4
2p
The finiteddifference method offers low intrinsic variance but high query cost. Mini-batching reduces base estimator variance by 1/b at a cost of or 2b queries. These multi-point approaches can be combined arbitrary directional derivative estimators.
A.2 Zeroth-Order SGD and Its Variants
Vanilla Zeroth-Order SGD The convergence of SGD has been extensively studied under various settings. Ghadimi and Lan [2013] established complexity results for computing approximate solutions using first-order and zeroth-order (gradient-fired information with Gaussian smoothing. For smooth convex objective functions, Duchi et al. [2015] obtained the optimal convergence upper bound for SGD under the zeroth-order optimization (ZOO) setting. Nesterov and Spokoiny [2017] provided the optimal convergence upper bound for Gaussian smoothing. In the realm of nonconvex optimization, Ji et al. [2019] proposed two new zeroth-order seriance-reduced optimization algorithms, ZO-SVRG- Goord-Rand and
MPIDER-Coord and provided improved analysis for the existing ZO-SVRG: Coord algorithm. These methods achieved better convergence rates and function query complexities than previous approaches. Berahas et al. [2022] derived convergence analyses for finite differences, linear interpolation, Gaussian smoothing, and uniform sphere smoothing methods. Recent studies have focused on non-smooth settings. Davis et al. [2022] and Zhang et al. [2020] established the sample complexity for Lipschitz functions without assuming smoothness. Lin et al. [2022] derived the complexity upper bound of SGD while noting vd scale compared to the smooth setting. Notably, Rando et al. [2024a] and Kornowski and Shamir [2024] revealed that by applying certain techniques, the non-smooth case is not inherently more challenging than the smooth case. A potential direction for extending this line of research is to explore the intersection between zeroth-order SGD and random reshuffling [Ma and Zhou, 2020, Mishchenko et al., 2020], minimax optimization [Chen et al., 2022], or dependent data [Ma et al., 2022].
Variance- Variance-Reduced Zeroth- Seroth-Order SGD A key bottleneck in vanilla zeroth-order SGD is the high variance of gradient estimators, which arises from both stochastic data sampling and the inherent randomness the gradient estimation process. This high variance necessitates small stepsizes, leading to slow convergence [Liu et al., 2020]. To address the variance from stochastic data sampling,
23

---

variance reduction techniques originally developed for first-order methods [Fang et al., 2018, Defazio et al., 2014, Johnson and Zhang, 2013, Nguyen et al., 2017], have been adapted to the zeroth-order setting. Algorithms such as ZO-SVRG [Liu et al., 2018, Huang et al., 2019, Gu et al., 2021], SOD-SVR//SPIDRR-COO et al., 2019, Chen et al., 2023] leverage epoch- updates with pariance-reducing correction terms or recursive estimator refinements. These methods significantly improve convergence by reducing the iteration complexity.
[Ji et al., 2019], and 700-SPIDER/SRRAF [Fang et al., 2018, Ji
Memory-Efficien Zeroth-Order SGD Standard SGD typically requires storing all intermediate gradient across layers to enable chain-rule-based backpropagation, which incurs substantial memory overhead, especially when training large models. To alleviate this, MeZO [Malladi et al., 2023] introduces a memory- -fffecrrefficient approach wherein it suffices to store the random seed used to generate the perturbation vector for each layer, dramatically reducing memory consumption. This principle motivates algorithms such as Addax [Li et al., 2024], ElasticZO [Sugiura and Matsutani, 2025], and ZO2 [Wang et al., 2025], along with related benchmarking efforts [Zhang et al., 2024, Gautam et al., 2024, Wang et al., 2024, Guo et al., 2024]. Additional strategies exploit sparsity to further reduce memory usage; notable examples include ZORO [Cai et al.. 2022b], the Extreme Sparsity framework [Guo et al., 2024], and the One-Bir Bit method [Cai et al., 2022a].
A.3 Discussions on the Forward Auto-Differentiation (AD) Approach
The forward gradient [Griewank and Walther, 2008] provides the exact directional derivative (with exactly zero bias), while the zeroth-order approach offers only an approximation of the derivative gradient. As a result, the peroth-orde approximation inherently introduces additional variance (even if it can be unbiased). As pointed out by Zhang et al. [2024], this makes the Forward AD method theoretically better in terms of estimator quality. However, there still multiple scenarios where the zeroth-order method is preferable.
Implementation Difficully: The practical implementation of Forward AD heavily relies on the availability of JVP (a. the Jacobian- -ector Product). A naive implementation will not reduce the memory usage and potentially increase the computation cost.
Memory usage: Forward AD can be memory-efficient when implemented properly. However, it still presents a higher memory usage than the zeroth-order optimization. Therefore, for the edge device or other extreme cases where the memory cost is sensitive, we may still prefer the zeroth-order approach.
We also note that aaroth-order optimization is clearly advantageous in black-box settings where the forward gradient is not available. Therefore, the forward difeeeefrreniiiionnn and zeroth-order approaches are not mutually exclusive, but complementary, depending on the feasibility and the device memory.
A.4 Discussions the Difference Between Our Results and Chen [2020]
Although the telescoping structure is the same as the one used in Chen [2020] as we have commented in the introduction section, have developed more results to the unbiased zeroth-order gradient estimator beyond this telescoping structure:
1. Identify when we can have an unbiased gradient estimator: We identify the condition under which the relescoping structure admits a valid expectation representation (Proposition 2.... This condition is critical for constructing unbiased estimators, but has not been established in Chen [2020].
2. More general unbiased gradient estimator & Reduce the number of fanction evaluations from 4 to 1: The estimator in Chen [2020] is a special case of our P4-estimator. Our framework extends beyond this, answering a fundamental theoretical question: What is the minimal number of function evaluations needed to construct an unbiased gradient estimator? We improve the known answer from give by Chen [2020] 1.
3. Identify the necessary and sufficient condition of achieving the optimul variance: More- over, one of our key focuses is identifying optimal parameter sequences {um) and (pn) (Theorem 3.1), rather than proposing specific estimator.
24

---

B
Bias Analysis
B.1 Absolute Convergence
In this subsection, we derive a sufficient condition to guarantee the expectation representation of the telescoping series introduced in Eq. (5). This requires ensuring the series is absolutely convergent, a property essential for interpreting it as the expectation of well-defined random variable. The following definitions are directly taken from Folland [2002]:
Definition B.1 (Convergent series). A series Σn=1 is said to be convergent if the sequence of partial sums Sv Σn On is a convergent sequence; that is, lim Sy exists.
Definition B.2 (Absolutely convergent series). A series Σx= On is said to be absolutely convergent the series - [0,] convergent.
The following classical result, known as the Riemann series theorem [Riemann, 1868, Spivak, 2008), highlights the necessity of absolute convergence when interpreting an infinite sum as well-defined value.
Theorem B.3 (Riemann Series Theorem). Let En=1 an be a conditionally convergent series of real numbers (i.e., convergent but not absolutely convergent). Then, for any real number r E R, there exists a rearrangement N M such that Σx=1 Oz(n) r. Moreover, there exist rearrangements such that the sum diverges to +00, -C, or fails to converge at all.
This theorem underscores the critical distinction between conditional and absolute convergence: a convergent series may yield different values under different summation orders. However, the definition of expectation for random variable does not permit such ambiguity, since the outcomes have no inherent order. To guarantee a well-defined expectation, absolute convergence is required: that is, for random variable X with outcomes {xx} and probabilities {pa), must hold that [|X] -m=1 In pp
00.
We now recap Proposition 2.1 describing the condition where the telescoping series is absolutely convergent, enabling valid expectation representation.
Proposition B.4. If the second-orde continuously differentiable function f Rd - R has L- Lipschitz continuous gradient and En Pn on, then the series
f(x
u10) f(x)
[(x (n+10) f(x)
f(x pn (xx
Pm
f4
Pr
pn+1
pn
is absolutely convergent and its limit is V.)()).
Proof. First, because Em lan + bal $ Em=1 laul + Σm= [bal. it suffices to prove
Σ
f( + (ativ) - f(x) f(x + mnv) ((x))
Mn+1
Un
is absolutely convergent. To prove its absolute convergence, we estimate the magnitude of the difference term using Taylor's theorem:
f { + (n+1v) f(x)
f(x + un - f(x)
Mn+1
Un
Pn+1 Vf(x) u R(x)
pn V)(x) To R((2))
pnt1
Hp.
(10) L
Pn+1 + Hnl
where (i) applies the Taylor theorem [Folland, 2002] with setting the integral form remainder fx + ffe) dt, (ii) assumes the global Lipschitz continuous
R(x) = [ (--)))
1) > a am
gradient, which results in the uniform estimate R(x) & for all Idd. Therefore, ensure the telescoping series is absolutely continuous, it suffices require Σ Hhn cc. The limit is directly determined by the original convergent series. It concludes our proof.
D
25

---

In the following two examples, present commonly used sequences {1x} that satisfy the condition 2m=1 < o, ensuring the absolute convergence required in Proposition B.4.
Example B.5 (Exponential Decay), Let pn = a² for some constant 0 < a < Then,
Mn
00.
Example B.6 (Polynomial Decay), Let pln
for some constant 8 > Then,
m
Pn < 00,
n2"
which is well-known as the Riemann zeta function C(s).
B.2 Unbiasedness of Zeroth-Order Estimators in P-Family
We begin by recalling the definition of the unbiased zeroth-order gradient estimators, denoted by 9 (J: {1m])2111 (Pn {pm)=_1> V), as given by Definition 2.2. Under suitable conditions on the differentiable function f Rd Proposition 2.1), this definition naturally yields the desired expectation representation. Moreover, the random distributions (pm)--) and V are independent. These conditions are sufficient guarantee the unbiasedness stated in the following result.
and the sequence { (e.g.. those provided by
Theorem B.7 (Unbiasedness). Lei P := P(S (pa)m= (po 1oeea- V) is defined as Definition 2.2. Then, for any estimator P(n, v) P, the following hold:
(a) E[P(n, v) u] vx) that is, P(n, v) is an unbiased estimator of the directional derivative v.f(x).
(b) If the random direction is chosen independently of the sampling n {pa]a- 1o and satisfies E E I, then
P(n, ( Vf(x),
that P(n, v) is an unbiased estimator of the full gradient.
Proof. By Definition 2.2, the directional derivative (x) naturally has the expectation representa- tion.
(a) Denote
Xn, = f(x + +11+))100000 - f(x)
f(x f(x+ + pn+10) - f(x)
f(x + (nv) f(x)
Pr
pn-1
Then by the tower property of the conditional expectation,
v=f(x)= En-(p=)%__ [Xn v] = En (ps 12 [E| Pm,, u
w/ u[
It concludes the proof.
(b) We consider the conditional expectation E[-lu]. We have
e E/P(n, v) 0
P(n, u) u
E R Pmn, v)
(z)o
= vf(x))
Therefore, we conclude that P(n, w) is an unbiased estimator of the gradient Vz(x).
26

---

C Variance Analysis
In this section we present the proof of Theorem 3.1. Each subsection contains the proof of the corresponding item. We recap its statement below:
Theorem C.1. Let 9 (f, (uf-innn- (pa) V) is defined as Definition 2.2. Suppose that : _R - I is second-order continuously differentiable and has L-Lipschitz continuous gradient, 2 Pn < oc, and V is the uniform distribution over the sphere with the radius d. Define
(pn+1 - Mn)
and
::
Pn
Pn
Then the following statements hold:
(a) If there exists a point Rd such that the Hessian V2 is positive definite and f(x) then the variances of the P1 is infinite. That is,
Var[Pi(n, v)o] = +00.
(b) The variance of P2-estimator P2(n, w)e is given by
L2
L
La
Var[P2(n, u] Var[PP((,,))) 3
p²
d
0+
d
12
3
(c) The variance of P3-estimator Ps(n, uje given by
L2
L
VarlPa(n, VarlP Var
d=p²
d
d
(d) The variance of Pa-ssimmator P4(n, u)v is given by
3L
L'd
Var(Pa(n, & (d - (ed-1)/V// 1)[
8 +
4
Proof. For the item (a), we present the proof in Lemma C.2. For arbitrary P p(n, P, have
Var[Pu]
(1)
2 u T [=))
[P - P4(11 v) + Pa(n, (V/()))
[P Pa(n, v)) dB (n, u)] |Vf(=)]
dP Pn +)) VarlP (n, u)
where (i) applies the unbiasedness (Theorem 2.3) and (ii) applies the definition of e (P is an unbiased estimator of P4). Therefore, we start with Var[P: vj. the variance of the P4-estimator. Then it suffices evaluate (n, v) (n, "] for ls and
3. Therefore, we prove the item (d) first. The detailed proof is included in Lemma C.3. Based on this result, we obtain the variance upper bound of P2- and P3 Pa-estimators Lemma C.4 and Lemma C.5, respectively.
C.1 Variance of Esstimmtorr
Lemma C.2. Under the same setting as Theorem C.I, if there exists point E Rd such that the
Hessian (xx is positive definite, then the variances of the Pi-sstmmator at x is infinite.
Proof. Recall that for the random direction 72" Uniform(Sd-1) and the random variable U4 Uniform (f0, 1, 3)), we have the P eestimmmator defined as
(i(n, v) = f(x + (1o) f(x (U()))))))))
Mi
27

---

1 f(x + (n-1e) - -fo==f((([l))) ((4))((1400)
f(x+ pno) (()){(U)-0))
Pn
pn-1
Since this estimator is unbiased (Theorem 2.3).
Var Pi(n, b) u] dE PI(n, u) [V
so it suffices show (n, 2) ab. For brevity write :1:- P v)-
As the Hessian v² is positive definite at the point x, it has v²f(2) 2 Amia for some Amin > 0. We consider the event {N n, U3 3). one finds
Pi f(x + pnv)
with probability Pr 4
1
Patta
Hence
(nu) Pr(N = U3 ( anu)) 2
En-(pa) [Pq] >
Pnttan
[f(x + mn))]²
Path
R=1
[f(x)]² Lf(x + pnv) - - f(x)
Pol
Pnl
Without loss of generality, we assume f(x) O. the case f(x) < O, we apply the L-smoothness to obtain lower bound instead. By the second-order Taylor expansion [Spivak, 2008] (or the strong convexity near the point x),
f(r + pv) 3 f(x) + Vf(x) p=² Amin T.
so we have
[f(x + pnv) - f(x)]f(x)
vf(x) ddminf(x)
Pnf22
Pntfi
2pn
As the result, we have
1
f(x))
dd ((x)
En (pn) t Uniform
|P?]
Σ
+
3
Pr
2pn
71
As {Pa] is PMF, it must diverge to infinite.
C.2 Variance of Py-Estimator
Lemma C.3. Under the same setting as Theorem C.I, the variance of P4-estimator Pa(n, v) upper bounded by
3T228
L.d
|M Mml
x
Σ
Pr
n=1
Proof. Recall that
PΛ(n, v) f(x + (1v) f(x)
1 f(r + pn pn-1v) f(x)
f(x + un f(x))
+
11
Pn
Pn41
Our goal to bound Var Var[P4(1 (4 v) wJ. For each n, define the remainder term
f(x (nv) - f(x) - pn Vf(x)
(v)
thn
f(x + 11,11 - [(x) 5m(()) + T. A,
pn
28

---

Then
Pa(n, u) = 01(0) + 0nt1(v) - 6n(v)
Pr.
Hence our vector estimator is
Pa(n, u vv Vf(x)) T
+ (81(0) fn (0)- v.
Pn
Since we have shown that Pa(ni, v) u unbiased (Theorem 2.3) and [[(v))]
Ia we have
E (1)
dn (x111())_8.11))
v 0.
(11)
Pr
The variance is given by
v) -
-(d-1)]V/(=) dE 61(0) + 6x+1(0) - 5n(r)) 2
Pr
+ 22EV f(x) (0((()) + Sn+1(u) 5n(1)
@ (d - 1)) de (1((1) + Sn+1lv) - 6n(v)
pn
(d - 1)) dES{(0)
(v)
d
+ 2dE
(6m+1(1) 6,(( + E
(8n+1(1) - (8n+1(0)-8((t)) =
Pn
ph
( (d- 1) 1) 3dEd (v) + (6n+1(u) - 6n(()))
pz
where (i) applies the unbiasedness Eq. (11) and (ii) applies the telescoping series Eq. (2). Now remains upper bound the remainder term (v) and the remainder difference term 6n+1(v) 6, (v).
Bound the remainder term (v): By L-smoothness of ʃ Rd R, we have
L 1012
f(x p.)) - f(x) < M4, f(x) v ++
2 ||ll
dLp
Pa
where (i) applies Lvv. La (this condition ensures that |o||2 Tr (1|ul13) Tr (ov") d). As the result,
Ld (v) & 2 pn.
(12)
or we may use the following almost-sure upper bound
L.d P
[b](1]] <
Bound the remainder difference term 5n+1 (v) 5 (u): We define the remainder term as a function in μe; that is,
6(H)= = f(x + pv) - f(x) - V
It automatically gives
8 8p+1(0) -v,(v) = @(un+1) p(pn)-
29

---

As Φ 11n-1.Pul Pnl - R is a continuous differentiable function (by our assumption that f Rdd - R is second-order continuously differentiable), we can apply the mean-value theorem [Folland, 2002]: There exists [ [n] such that
o(un+1) 6( Ø' (Mn+) -M))
We again applying the L-smoothness of f : R² - R (essentially the bounded Hessian assumption) to o' (s), which leads
18n+1(v) - -x((0)] & Ld
Mn+1 Mxl-
(13)
As the result, we obtain the upper bound as
3L2
2d
Ld
Var(Pr(n, < (d 1)[
Mn+1
pml
Pn
n-1
T,2
|a+11 - pml 12
< (d - 1)V +
4
2
Pn
R=1
where (i) applies the expectation over ~ (pa)a=i. which cancels out one
C.3
Variance of Pa-Estimator
Lemma C.A. Under the same setting as Theorem C.I, the variance of P2-estimator Pz(n, v) is upper bounded by
La
(Un - Mn+1)²
L
L
Var(Pz(n,,) < (n.
d
Σ
d
+
d
12
Pr
3
Pr.
Proof. Recall that P2(n, v) is defined as
P2(n, f(x + piv) f(x)
L(Ug-0)
P1
f(x + pn+10) - f(x)
I(Us-1) f(x + pn f(x)
+
[[Us-2)
Pn
pn+1
pn.
where Ug- ~ Uniform ((0, 1,2]) is a selection variable. Then
([((2(n, - P4(nn =)²² w]
f(x un f(x) 1
1
f(x 1n+11) f(x)
=P(U; 0)
Pr
Pn+1
Un
f(x + plu) - f(x)
1
f(x + Mn - f(x)
+ +PU(s 1)
+
141
Pn
Hhu
f(x (ut) f(x)
f(x Ma+10) f(x)
+ PUU 2)
+
141
Pa
14--1
001 (6)+1(0) (v))
[8)(1) +
[8, (v))²
+ [0)((0]]2 D
[6n+
3 3
P
1 (5 (v) (n)) 2
2 1
(v))2
18 (v))²
3 p2 [snti
3
3 3
p
where (i) applies (a + 6)2 & 20² + 26 and 8n(v) = [attn ((x) V bounded this term in Lemma C.3, have
As we have
14
[[(P2(n, v) Pa(n, ) m]
30

---

Σ Pn
(P(((, v)
Pa(n, v)) v)) m,
1=
ddd
L (pon Mn+1)²
dd
Σ
Pn
11 +
3
12
Ph
3p
L
(enn
)
L
L
H
d
+
dd
+
d
12
Pn.
3
3
Pn
n=1
n=1
Therefore, we finally obtain
Var(P2(n, v) [R[[(P2(,,)) v) 2a(n, + Var Var(Pa(n,,) v]
I
T
(Un Ma+1)2
I
d
Σ
d
P1
+
d
12
Pn.
3
3
Pn
n=1
completes the proof.
C.4
Variance of P3-Estimator
Lemma C.5. Under the same setting as Theorem C.I, the variance of Pa-estimator Py(n, v) is upper bounded by
I,2
(Pn 1n+1)
122
Var(Ps(n, & Var o))
+
d
+ 8
dd μ²
8
Pn
n=1
Proof. Recall that Ps(n, v) defined as
Pa(n, v)
f(x p1v) f(z) U2 1
f(x un+10) f(x)
f(x + pnv) - f(x)
(1 U2).
f1
Pr.
Pnt1
Pn
where U2 Uniform ((0, 1}) is selection variable. Then [[(()(n, P1(n, 0))
n, u]
f(x + 1n+1v) [(x) J(x + (nv) - [(x)
-P(U2- = 0)
Pa 21)
Pr
Pnt1
Un
f(x (uv)) - f((x) - -Pa(n, v) ?
+ IF(22 = 1)
11
f(x + (n+10) f(x)
f(x + un f(x) _1
f(x + (lv) f(x) 2
+
2 Pn
Mut1
Hh
1/1
As we have bounded this term in Lemma C.3, we have
E[(Ps Pa(n, v]
pn (P(((( P,(n,v)) [
L2d
L'dd
|n+1 - 2 +
8
8
7
I
(phn 1n+112
I
d
+
de
d
Pn
8
Therefore, we finally obtain
Var(Pa(n. -dB[((P - (P(n,, v)) + Var Var(Pa(n,1)) u]
I
(pn Mn+1) )
I
<)
Σ
+
d3
Pn
8
It concludes the proof.
31

---

C.5 On the Optimal Sampling Distribution and the Perturbation Stepsize Sequence We recap the full statement of Theorem 3.2.
Theorem C.6. Let (umln= be positive, decreasing sequence with Σn= Mn 00, and let {po]a=1 be a PMF Then the following statements hold:
(a) Define On
[an+1 -aa
The lower bound of o is given by
Pa
8
On
n
Moreover, equality holds if and only if Pm.
Un Ma--aan Ph+1
(b) Define Pn = B., The lower bound of yp is given by
m1.
Pn
Mn
n=1
Moreover, equality holds if and only Pm
In
Pr
Proof.
Write an pn+11 Un, SO
On
l41 and En Pn
By Cauchy-Schwarz inequality,
Σ |anl
Pr
On.
Pn
which yields the claimed lower bound in (a). Equality occurs exactly when pa c/ i.c. (un Pp411)/pp p/p.. Similarly, By Cauchy-Schwarz inequality,
Σ
Pn
Pni
Pr
1=1
which yields the claimed lower bound in (b).
In the following example, we include the omitted details of Example 3.4. Example C.7. We consider the Zipf distribution n Zipf(s) (s 1). Then
1 ((3) n
((s) Σ
Pr.
n-1
We define (1x) by the recursion
1
Hnn Pn+1
Pn 11 ((s) n**
so that summing gives the closed-form
Un = μi (1 _ = 11
Σgen
((8)
((8)
direct check shows this choice attains the lower bound = 1³ on [(un - 12pp) Now we turn to bound yp:
((s)
Σ
n² (1-
p
Σ
n
(())
(Ejngga )
((3)
((s)
Pa
For > 3, use the integral bound
8
nn
1-s
I
dx
32

---

n² to get
11
&
Pn (s 1)2 ((s)
Since s > 3 the series
-n n 5
<(s 2) converges, giving the clean bound
((s 2) (8- 1)2 (()) 3
C.6 Discussions: Variance of Random Directional Derivative
In this subsection, analyze the variance of gradient estimate based random directional derivative. Let e be a random vector uniformly sampled from the sphere with the dimension d. We approximate the gradient (z) using the random directional derivative defined as
v.f(z) = vv V fz)).
Assuming that the expectation satisfies Elve" '1 I, it is essential to evaluate the variance of this estimator. Specifically, we compute:
[[V)(x) Gu ue Vf(=)] = /
d]V/(>))]
Thus, the variance is given by
Var V vff(f)))]
|V (d- (d-1-dddll((())) 1)IV
This result indicates that even when the exact directional derivative is available, the variance still scales with (dd relative to the gradient norm. Consequently, it is not avoidable to remove the dependence on the dimension
Convergence Analysis
In this section, we present the proof of Corollary 3.5. Recall that goal is to solve the stochastic optimization problem Eq. (1):
min f(=)====:::= = Bcee:(()) z
where the second-order continuously differentiable function f(=:E) Rd - R has L-Lipschitz gradient for every E. We consider the convergence upper bound of the classical stochastic gradient descent (SGD) algorithm with the constant learning rate given the initialization Zn.
Ti+1 = (z1),
(SGD)
where g(xx) is an unbiased estimator of the full gradient V (x)).
D.1 The Convergence Upper Bound
The full statement of Corollary 3.5 is given as follow:
Corollary D.1. Consider the stochastic optimization problem in Eq. (1), and suppose that the individual loss f(r; is second-order differentiable with L-Lipschitz continuous gradient in Assume the stochastic gradient approximated using the Px-estimator (R v) for k = 2, 3,3,4 4. Let the SGD iteration be defined as
uniformly over E
F1+1 = X, - - vz) 1z
where 7) e (0, I is the stepsize. Then the iterates satisfy the following convergence guarantee:
min [(()()) < Ln [Cdd 2dL [E==)(²]] T/T c
0
33

---

2812"
5. k 2,
13
where å f(xo)= f and Ck
= 3, (k -2,344 is the estimator-dependent
1
512
k == 4,
error term. Consequently, choosing 1j = (11/VdT) and = - yields the optimal complexity T = (2) of having B c.
Proof. Let g(x) = Pa(n, u) for = 2,3, 4. By Lemma D.4, we have
Blg(E)) < d[]V[(()|| Czd u² + 2d1 Be-=fş)
Then we = dL and C C2d p.2 2dL Be- in Lemma D.3. As the result, when & D
we have
2
min
E Vf
& Ln Czd /² 2dL(f* Ee-=f()] n/T 5.
OEtET
The complexity of making minogreT E (z:) & given by
125L
Cnd 2 2dL(s Eg=ES()
T
max
dI,
2
2²
Setting 1 (2), the optimal complexity is given by T 0(4).
D.2 Supporting Lemmas
Lemma D.2. Let f* == min, f(x) and f == min, f(= for each E. f(x;() has L-Lipschitz gradient, then
E V < LIV -
Proof. This lemma is directly taken from Proposition 2, Khaled and Richtárik [2022].
Lemma D.3. Let the second-order continuously differentiable function f(: : Rd - IR be lower bounded by f* := ming fxx and have L-Lipschitz gradient for every E, and g(x) be an unbiased estimator of fzx). Suppose that 9(x) and fxx satisfy the expected smoothness condition
[] < B [[V( +
If the learning rate LB and define 6 == f(xo) f then the T-th iteration of SGD satisfies
min 0stET
fx LCn+ + nT
Proof. This lemma is directly taken from Theorem 2, Khaled and Richtárik [2022].
Lemma D.4. Let g(x) be the Ps-estimator (for 2, 4), f min, f(x), and f ming fxx E). Then the expected smoothness condition is given by
Ellg(E)] vL Cxds 1² + 2d Bx-=fE)
where
2,
12
C
31.2
k 3, 51
k == 4,
T
the error tern introduced by the zeroth-order gradient estimation. For P2-estimator we choose {px} and {pn) such that @ 212 and up & 2μ².
34

---

Proof. Let (x; be the Pæ-estimator for k = 2,3,4. First, we notice the following variance- decomposition holds;
v - - vf(=;))
Vf(=()) V/ar)) () Iaaar E]
By Theorem 3.1 and Theorem 3.2, we set
Var [P& < (d [>; pxd
where Ck is determined by the estimator; we also assume that the optimal distribution and perturba- tions taken obeying Theorem 3.2. As the result,
(== B (xx + +El 1)I C2d
(2;; 01 Czdl
dI vl)) + Cgd 122 2dI[f* Re-=/(),
where (i) applies Lemma D.2. It completes the proof.
E Experiments Details
In this section, we describe the detailed experiment setting and the hyperparameter configurationss
E.1 Synthetic Example
In the synthetic example, we compared gradient estimators across varying dimensions (d = 16, 64, 256, 1024, 4096) using both quadratic and logistic functions. For fair comparison. all es- timators used consistent number of function evaluations of 3 and the perturbation stepsize 10-5 (for Pæ-estimator, we set 141 (t). The P3-estimator was configured with parameter s 2.0 and followed the same perturbation stepsize scheduling as Example 3.4. We evaluated four gradient estimators: Zipf's Pa-eestimtoor, one-side two-point estimator with the Gaussian smoothing, one-side two-point estimator with the uniform smoothing, and two-side two-point estimator with the Gaussian smoothing. Each configuration was tested over 100 independent trials with a fixed random seed for reproducibility.
Code Availability and System Requirements All source code is included in the supplementary materials. No specific hardware is required; any machine supporting Python 3.10.10 10 should suffice. Jupyter notebook version is also provided for convenient execution on Google Colab.
E.2 Language Model Optimization
In the language model optimization experiment, we compare the performance of different gradient estimators within vanilla SGD framework for fine-tuning language model sentiment clas- sification task. The learning rate of SGD is fixed at r = 10-4, the batch size SGD fixed at 16 (this batch size corresponds the number of stochastic samples and different from the batch size multiple-pon proposed unbiased estimators, do not sample the extreme tail of the distribution. Instead, we truncate the sampling distribution by enforcing Pa 10-3 to remaining hyperparameters are summarized in Table 1.
zeroth-order estimator), and the perturbation stepsize is set to μ = 10-3 (for
μ41 μ)). Due limitations in numerical precision,
numerical stability and avoid excessively small probabilities. All
Code Availability and System Requirements All source code and reproduction instructions are provided in the supplementary materials. Experiments were conducted on a cluster running RHEL 8, equipped with dual AMD EPYC 7763 processors, 512 GB of memory, and seven NVIDIA RTX 5000 GPUs. reproduce the language model optimization experiments, we recommend using CUDA-compn cUDA-compatible GPU with at least 24 GB of VRAM.
35

---

Table Summary of gradient estimators used in the language model optimization experiment.
Estimation Method
Estimator Formula
Batch Size
# Function Calls
Notes
fle prs)) - f(z)
One-Side Two-Puint Estimator
12 Normal(), [g)
Two-Side Two-Point Estimator
fi pmi) [i 102)))
2
44 Normal Lz)
2μ
Zipf's Py-Estimator
Eq. (7)
Example 3.4
Geometric Pa-Estimator
Eq. (7)
0.5, Example
Zipl's P4 Pa-Estimator
Eq. (6)
1.5, Example 3.4
Geometric g-Estimator
Eq. (6)
4
0.5, Example 3.3
F Broader Impact
This work introduces a new class of unbiased zeroth-order gradient estimators that offer both theoretical guarantees and practical advantages. By eliminating bias without increasing variance, our method enhances the reliability of optimization in settings where gradient information is unavailable or costly obtain. These include fine-tuning large language models under memory constraints, conducting black-box adversarial robustness evaluations, and solving scientific computing tasks such as physics-informed neural networks. On the theoretical side, our estimators advance the understanding of zeroth-order optimization and provide new tools for the zeroth-order gradient estimation. This work opens a promising direction for future research gradient-fre optimization and its broad applications in machine learning and beyond.
G Limitations
Despite the theoretical guarantees and empirical improvements demonstrated by our proposed unbiased zeroth-order gradient estimators, several limitations remain. First, the estimators rely on sampling from an infinite sequence of perturbation steps, but practical implementations must truncate this sequence, which may reintroduce bias or affect variance control. Second, the proposed estimators are based on directional derivatives, which inherently exhibit a dependence on the problem dimension d; this dimensional dependence generally unavoidable for this class of methods. Lastly, while we validate the approach on synthetic tasks and language model fine-tuning, we have not extensively evaluated its performance across broader range of optimization problems, and the observed empirical gains may not fully generalize to settings involving non-smooth objectives or high levels of evaluation noise.
36


Reviews:
[{"review_id": "lje2PxYGKH", "reviewer": "Reviewer_sGPT", "summary": "The authors considered the problem of zeroth-order optimization (ZOO) in a stochastic optimization setup. They argued that most existing ZOO methods have a biased estimate of the gradient. Therefore, they proposed a family of unbiased estimators of the gradient based on merely function evaluations. In particular, they constructed four unbiased estimators that depend on 1 to 4 function evaluations in estimating the gradient. The main idea in devising these estimators is to write $\\nabla_v f(x)$ in a direction of vector $v$ as a telescoping series of gradient estimators with diminishing perturbation stepsize $\\mu_n$. This formulation provides an unbiased estimator given in (6) (which is called $\\text{P}_4$ estimator). Moreover, they suggested three other estimators with a reduced number of function evaluations. Under some regularity assumption on the objective function, they showed that the proposed family of the estimator is unbiased. Moreover, they gave upper bounds on the variance of gradient estimators in terms of perturbation step sizes $\\{\\mu_n\\}$ and sampling distribution $\\{p_n\\}$. In particular, for the estimator with one evaluation function, the variance might be unbounded, but for the rest of the estimators, by tuning $\\mu_n$ and $p_n$, the variances match the optimal order of classical two-point estimators. Moreover, for ease of implementation, they proposed a method that, given a sampling distribution, derived the perturbation stepsize, and they gave a formulation for the perturbation step size for two distributions of Geometric and Zipf distributions. The experimental results showed improvement over some previous work on synthetic functions and training of a language model.", "questions": "- What is the motivation behind having an unbiased estimator of the gradient? For instance, in the literature of variance reduction methods in stochastic optimization, most proposed methods are biased estimators while they reduce the variance. Overall, they can achieve the best known sample complexity while the estimator is biased.\n\n- Compared to biased estimators in the ZOO setup, what are the advantages of the current work in terms of sample complexity? It is good to clarify any advantages in the paper.\n\n- The explanation in lines 89-97 is not so clear to me. In particular, would you please explain more about ``the outcomes of a random variable have no naturally given order (which makes it different from a series), which requires a random variable’s expectation 95 to be well-defined regardless of any such ordering''. It seems that the assumption in Proposition 2.1 is needed because of this argument, which limits the applicability of the current work.\n\n- Would you please compare the current method with (Chen 2020) in Figures 1 and 2? Moreover, it is also good to add optimization algorithms such as SGD and Adam in these figures as a reference.\n\n- Could you elaborate on the impact of sampling distribution on the performance of the estimator or the optimization algorithm in general? In particular, in what cases, which one of the two considered sampling distributions should be used?", "limitations": "The assumptions are mentioned in the paper, but it is good to explain what the impacts of these assumptions are and for what applications the proposed estimators can be used. It seems that there was no discussion on the potential negative societal impact of the work. It is good to add any societal impacts or otherwise, mention that there is no immediate impact.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n- The authors proposed a family of unbiased estimators of the gradient in the ZOO setting. It seems that there is only one unbiased estimator in the literature by (Chen 2020), which was a special case of $\\text{P}_4$ estimator.\n\n- The authors provided theoretical guarantees on the unbiasedness of the proposed family of estimators and also gave upper bounds on the variance for the three estimators mentioned in Section 2.3.\n\n- They derived optimal sampling distributions for given perturbation step sizes that minimize the upper bound on the variance (Theorem 3.2). They also provided an efficient implementation of the estimator by fixing first the sampling distribution and then adjusting the perturbation step sizes accordingly. \n\n### Weaknesses:\n- The motivation behind the unbiased estimator is not well explained. In fact, it is not clear what the advantages of having an unbiased estimator in the ZOO setup are.\n\n- The assumption on the objective function being a second-order continuously differentiable function is not well motivated, and it is not clear in which theorems this assumption is needed. Apparently, it is needed in all theorems and propositions, but it is not explicitly mentioned in the statements. I think this assumption is also restrictive, as neural networks with ReLU activation may not satisfy this assumption.\n\n- In the experiments, it seems that the authors did not compare with the closest related work of (Chen 2020).", "was_revised": true, "final_justification": "Overall, I consider this a good work that can inspire further improvements by other researchers in the area of zeroth-order methods. However, comparisons with (Chen, 2020), SGD, and Adam—particularly with (Chen, 2020)—should be included in the revised version. Regarding the assumption of second-order continuous differentiability, I still believe this is a strong condition that may not hold in many practical applications. Therefore, I would like to keep my score unchanged.", "rebuttal": "We deeply appreciate the reviewer’s positive evaluation and constructive feedback. Below, we provide our point-by-point responses to each of the concerns raised.\n\n## Response to questions\n* (**Motivation behind having an unbiased estimator of the gradient**) What is the motivation behind having an unbiased estimator of the gradient? For instance, in the literature of variance reduction methods in stochastic optimization, most proposed methods are biased estimators while they reduce the variance. Overall, they can achieve the best known sample complexity while the estimator is biased.\n\n    **Response:**  We thank the reviewer for raising this important question. We highlight two key advantages and one potential advantage of using an unbiased estimator: \n\n    (i) *Practical benefit:* It removes the need to balance bias and variance as a coupled trade-off. This simplifies the hyperparameter tuning process, allowing us to focus solely on controlling variance.  As indicated by our Theorem 3.2, the variance upper bounds of the $P_3$ and $P_4$ estimators are parameter-agnostic; that is, they do not introduce any additional hyperparameters compared to classical biased estimators.\n\n    (ii) *Theoretical benefit:* The absence of bias enables the direct application of standard optimization theory, which leads to tighter convergence guarantees without needing to account for the complications introduced by estimator bias. While our current analysis focuses on vanilla SGD (as in Corollary 3.5), the same theoretical framework readily extends to a broad range of optimization methods.\n\n    (iii) *Potential benefit:* The unbiased estimator opens a novel and underexplored direction in the design of zeroth-order gradient estimators. While existing (typically biased) estimators have been extensively studied and refined, the space of unbiased estimators remains largely untapped. We believe our work creates new opportunities for developing more effective and theoretically grounded zeroth-order estimators. For example, in our Proposition 2.1, we identify what kinds of assumptions are needed to represent a *Taylor series* as an expectation; this result could be extended to many other series and may motivate new zeroth-order estimators that can be tailored to the specific structure of the objective function.\n\n* (**Advantage in terms of sample complexity**) Compared to biased estimators in the ZOO setup, what are the advantages of the current work in terms of sample complexity? It is good to clarify any advantages in the paper. \n\n    **Response:** Corollary 3.5 shows that SGD with $P_3$/$P_4$-estimator achieves $\\mathcal{O}(d/\\epsilon^4)$ oracle complexity. This complexity is information‑theoretic optimal (that is, it doesn't depend on any specific form of the gradient estimator) and cannot be further improved (see Proposition 1 [Duchi et al. 2015]). From this perspective, all existing zeroth-order gradient estimators cannot outperform our proposed method in terms of sample complexity. \n\n    * [Duchi et al. 2015] Duchi, John C., et al. \"Optimal rates for zero-order convex optimization: The power of two function evaluations.\" IEEE Transactions on Information Theory 61.5 (2015): 2788-2806.\n\n* (**Explain  lines 89-97**) The explanation in lines 89-97 is not so clear to me. In particular, would you please explain more about ``the outcomes of a random variable have no naturally given order (which makes it different from a series), which requires a random variable’s expectation 95 to be well-defined regardless of any such ordering''. It seems that the assumption in Proposition 2.1 is needed because of this argument, which limits the applicability of the current work.\n\n    **Response:** Thank you for the insightful comment. We will revise lines 89–97 to clarify the underlying idea. The goal of this paragraph is to explain when an infinite series can be interpreted as the expectation of a random variable. Consider the convergent series $\\sum_{i=1}^\\infty p_i x_i$. \n\n    1. What is the value of this infinite series?\n\n        We can calculate the finite-sum $S_n := \\sum_{i=1}^n p_i x_i$; then we let $n$ tend to infinite. \n\n    2. What is the value of an expectation of the random variable? \n\n        We consider the random variable $X$ with $P(X=x_i) = p_i$. Then the expectation $E[X]$ can of course be written as  $\\sum_{i=1}^\\infty p_i x_i$. However, the notion of expectation must be well-defined independently of any ordering of outcomes (that is, **\"regardless of any such ordering\"**). That is, for an arbitrary permutation $\\sigma: \\mathbb{N}  \\to \\mathbb{N} $, all series   $\\sum_{i=1}^\\infty p_{\\sigma(i)} x_{\\sigma(i)}$ should represent the same expectation $E[X]$. \n\n    This invariance under reordering is guaranteed only when the series is absolutely convergent. Proposition 2.1 provides a sufficient condition to ensure this property. While the condition itself may appear stronger than absolute convergence, it has the advantage of being easier to verify in practice.\n\n* (**Compare (Chen 2020), SGD, and Adam**) Would you please compare the current method with (Chen 2020) in Figures 1 and 2? Moreover, it is also good to add optimization algorithms such as SGD and Adam in these figures as a reference.\n\n    **Response:** Thank you for the suggestion. Yes, we will include all three methods in Figure 2 to facilitate a direct comparison, as requested. Regarding Figure 1, it is designed specifically to evaluate the gradient estimation error and does not depend on any particular optimization algorithm. Therefore, in Figure 1, we will include a comparison with (Chen, 2020), but not with SGD or Adam, as they are not directly relevant in this context. \n\n* (**Impact of sampling distribution on the performance**) Could you elaborate on the impact of sampling distribution on the performance of the estimator or the optimization algorithm in general? In particular, in what cases, which one of the two considered sampling distributions should be used?\n\n    **Response:** Our theoretical analysis (Theorem 3.2) shows that both sampling distributions lead to estimators with optimal variance. This result is, in our view, quite surprising at the first glance. Because, intuitively, using the $P_3$ or $P_4$ estimator seems challenging in practice due to the need to tune an infinite-dimensional hyperparameter $\\{ p_n \\}$ and $\\{ \\mu_n \\}$. However, Theorem 3.2 reveals a key insight: optimal complexity can still be achieved as long as these weights satisfy a simple algebraic relation.\n\n    To demonstrate this, we have provided an example based on the Zipf distribution, which illustrates how to construct $P_3$ or $P_4$ estimators when a closed-form solution is not available. Our empirical comparison between the Zipf and Geometric estimators in Figure 2 further supports Theorem 3.2, showing that their performances are indeed comparable, with no significant difference observed in practice.\n\n    Therefore, Zipf, Geometric, or otherwise, does not fundamentally affect the convergence rate or the variance of the estimator. All of them will achieve the optimal sample complexity.\n\n\n## Response to weaknesses\n\n* (**The second-order continuously differentiable sssumption**) The assumption on the objective function being a second-order continuously differentiable function is not well motivated, and it is not clear in which theorems this assumption is needed. Apparently, it is needed in all theorems and propositions, but it is not explicitly mentioned in the statements. I think this assumption is also restrictive, as neural networks with ReLU activation may not satisfy this assumption.\n\n    **Response:** We thank the reviewer for this insightful comment. We would like to clarify the motivation for assuming second-order continuous differentiability. This assumption is introduced because our theoretical analysis relies on applying Taylor’s theorem to the second order, which requires the existence of Hessian $\\nabla^2 f$. In contrast, assuming only Lipschitz continuity of the gradient (i.e., $L$-smoothness) is insufficient to guarantee this existence.\n\n    We acknowledge that this assumption excludes certain widely used neural network architectures, such as those employing ReLU activations (we also note that the $L$-smoothness assumption will also exclude this example). However, we would like to emphasize that this assumption is in line with prior work in zeroth-order optimization.  For example, the MeZO paper [MeZO] relies on the existence of the Hessian in its proof, even though they do not state this assumption explicitly.\n\n    * [MeZO] Malladi, Sadhika, et al. \"Fine-tuning language models with just forward passes.\" *Advances in Neural Information Processing Systems* 36 (2023): 53038-53075."}, {"review_id": "i8mqEk5Uxa", "reviewer": "Reviewer_5co9", "summary": "Authors propose new way to build unbiased estimate of gradient via gradient-free setup. The idea is interesting and according to authors claim propose for free a way to eliminate bias (by conserving the variance estimate).", "questions": "Why do authors use such a finite difference? It seems that symmetric one will be better\nShamir O. An optimal algorithm for bandit and zero-order convex optimization with two-point feedback //Journal of Machine Learning Research. – 2017. – V. 18. – №. 52. – P. 1-11.", "limitations": "Main limitation-- if $f$ is available with noise like in Gasnikov A. et al. The power of first-order smooth optimization for black-box non-smooth problems //International Conference on Machine Learning. – PMLR, 2022. – V. 7241-7265. In this case the construction seems to be quite sensitive and may work worse than classical one approach.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "So the paper is well written and propose a useful way to reduce bias. But I'm not sure that this way leads to the same variance at the end as it will be for biased case, see e.g. Tsybakov A. et al. http://www.jmlr.org/papers/volume25/23-0733/23-0733.pdf or Gasnikov A. V. et al. Stochastic online optimization. Single-point and multi-point non-linear multi-armed bandits. Convex and strongly-convex case // Automation and remote control. – 2017. – V. 78. – P. 224-234. Roughly speaking, for me it seems that $d^3\\mu^2$ terms could be improved as $d^2\\mu^2$. May be it's worth to include more accurate discussion about this.", "was_revised": true, "final_justification": "I consider the paper will be useful for gradient-free community, though in initial variant in some places the paper was not accurate in interpretation of details. Now authors fix this drawbacks that is why my score is 5.", "rebuttal": "We deeply appreciate the reviewer’s positive evaluation and constructive feedback. Below, we provide our point-by-point responses to each of the concerns raised.\n\n## Response to weakness\n\n- (**Include more accurate discussions on the variance term $d^3 \\mu^2$ vs. $d^2 \\mu^2$**) So the paper is well written and propose a useful way to reduce bias. But I'm not sure that this way leads to the same variance at the end as it will be for biased case, see e.g. \n\n    - Akhavan, Arya, et al. \"Gradient-free optimization of highly smooth functions: improved analysis and a new algorithm.\" *Journal of Machine Learning Research* 25.370 (2024): 1-50.\n    - Gasnikov, Alexander V., et al. \"Stochastic online optimization. Single-point and multi-point non-linear multi-armed bandits. Convex and strongly-convex case.\" *Automation and remote control* 78.2 (2017): 224-234.\n\n    Roughly speaking, for me it seems that $d^3 \\mu^2$ terms could be improved as $d^2 \\mu^2$. May be it's worth to include more accurate discussion about this.\n\n    (**Function value noise**) If is available with noise like in:\n\n    - Gasnikov, Alexander, et al. \"The power of first-order smooth optimization for black-box non-smooth problems.\" *International Conference on Machine Learning*. PMLR, 2022.\n\n    In this case the construction seems to be quite sensitive and may work worse than classical one approach.\n\n- **Response:** We deeply appreciate the reviewer's insightful comment; it indeed identifies an important setup of zeroth-order optimization which we lack sufficient discussion. When considering the scenario where the function evaluation is exact (considered in our setting), the variance error term contains $d^3 \\mu^2$ is worse than $d^2 \\mu^2$ but it could be resolved by setting a sufficiently small $\\mu$. So, it still leads to the optimal complexity. However, in the noisy oracle setting (as studied in above references), this term will have negative impact on final result. These works spend additional effort to delivery a tighter characterization on the upper bound. To better clarify the difference between these settings, we will add the following discussion right after our current Comparison with Existing Literature.\n\n    > (Comparison with the Noisy Oracle Setup)  In our work, we consider the exact function evaluation setting with noiseless values. In this case, our variance scales as $d^3 \\mu^2$, which is worse than the $d^2 \\mu^2$ of some specific biased estimators. However, this can be mitigated by choosing a small enough $\\mu$, so the overall sample complexity remains optimal. However, in the noisy function evaluation setting, where each function evaluation may return a noisy value, a smaller $\\mu$ amplifies the noise, leading to degraded performance. Several recent works have provided more refined analysis under noisy setups with improved variance behavior. Notably, Akhavan et al. (2024) demonstrated that for highly smooth functions, the $\\ell_1$-randomization can reduce the variance scaling to $d^2 \\mu^2$ with achieving the improved performance for highly smooth objective functions, which extends the existing $\\ell_1$-randomization proposed by Akhavan et al. (2022). Earlier foundational work by Gasnikov et al. (2017) analyzed the variance behavior in single-point and multi-point bandit feedback settings, and more recent developments further explore the impact of first-order smoothness in noisy black-box optimization [Gasnikov et al., 2022]. Notably, all of these results achieve the optimal complexity derived by Duchi et al., 2015. \n\n## Response to questions\n\n- Why do authors use such a finite difference? It seems that symmetric one will be better according to [Shamir2017].\n\n    **Response:**  We appreciate the reviewer’s insightful comment. We agree that the symmetric one (also called the two-side two-point estimator) is commonly considered better than the naive one-side two-point estimator. However, there are two reasons for still considering the one-side estimator: \n    (1) On the theoretical side, they still achieve the same sample complexity for smooth objective functions. This result is implied by the lower bound results in [Duchi2015]. \n    (2) For the practical consideration, while the symmetric estimator uses $f(x+\\mu v)$ and $f(x-\\mu v)$ to compute the gradient, it commonly requires evaluating $f(x)$ for logging purposes. This operation will introduce an additional function evaluation.\n\n    Moreover,  our choice is motivated by a deeper theoretical interest in a relatively underexplored class of estimators characterized by structured hyperparameters (i.e. $\\{\\mu_n\\}$ and $\\{p_n\\}$). These hyper-parameters do not appear in other types of gradient estimators.  Remarkably, our analysis (Theorem 3.2) identifies the optimal structure for these parameters, showing that a simple relation between them suffices to achieve optimal sample complexity. Therefore, we believe there still could be a large space to improve this approach.\n\n    * [Shamir2017] Shamir, Ohad. \"An optimal algorithm for bandit and zero-order convex optimization with two-point feedback.\" *Journal of Machine Learning Research* 18.52 (2017): 1-11.\n    * [Duchi2015] Duchi, John C., et al. \"Optimal rates for zero-order convex optimization: The power of two function evaluations.\" *IEEE Transactions on Information Theory* 61.5 (2015): 2788-2806."}, {"review_id": "1evp75Zbci", "reviewer": "Reviewer_nG3r", "summary": "The paper proposes new  unbiased estimates of the directional derivative, based on a telescopic series representation of the directional derivative. The paper is very well written, clear, and the proofs are easy to follow. My only real concern is how this work fits within Neurips. I understand there is still plenty of legacy code (often written in Fortran) for simulations that requires zero order training. But here the connection to Neurips + ML is through the proposed application of model fine-tuning, and to save of memory when using autodiff. But here I have a question/issue. See [Application to LLM fine-tuning] below. Based on the response of the authors on this question, and the other reviewers, I am open to reconsidering my score.", "questions": "1. The reference lower bound in [Duchi et al. 2015] is for zero-order optimization, and not first-order. The fact that you have an unbiased estimator, means you can leverage the proof structure of SGD, such as the one in [Khaled and Richtarik]. Because of this, your statements in the main text and abstract such as \"... prove that SGD using the proposed estimators achieves optimal complexity for smooth non-convex objectives\" can be misleading. A reader could interpret that you have matched the optimal complexity of \"SGD for smooth non-convex objectives\" which is \n\n\n$$\\frac{ \\sqrt{2} L (f(x^0) - \\inf f) + (\\inf f - \\mathbb{E}[\\inf f\\_{\\xi}])}{\\sqrt{T}} $$\n\nsee Theorem 5.12 in \nGarrigos et. al, Handbook of Convergence Theorems for (Stochastic) Gradient Methods, arXiv, 2024\n\nI think you should re-word this contribution, to make it clear you are talking about the optimal zero-order complexity.\n\n2. How do your methods compare to other zero order methods used for fine-tuning, such as Mezo in [Malladi et al 2023]? Or is Mezo equivalent here to what you call Two-side point estimator? On this note, I see a small mistake in your Two-side point estimator in Table 1. You are off by a factor of 2, since it should be\n$$ \\frac{f(x+\\mu v_i)- f(x+\\mu v_i)}{2\\mu} v_i$$\nThis same mistake is also repeated on page 14 in three locations. \nA factor of two would make no difference if you tuned the learning rate. But since you have a fixed $\\eta = 10^{-4}$, this factor could make a difference.\n\n\n\n[Application to LLM fine-tuning]. To justify this application, you have cited some prior work, and stated \"storing the full computational graph ... can be prohibitively large for modern language models\". I checked your references, and only one [Malladi et al 2023] has been peer reviewed. There, their claims is the using zero-order uses a fraction of the memory of backpropagtion. This I agree with, since backpropagation requires storing all weights, and adjoints of your computational graph. But here is where I see an issue. These zero order methods should instead be compared to the cost/memory of computing the directional derivative $\\nabla_v f(x)$ when using forward propagation, which can be done with a very small memory overhead. To justify the case of using zero-order for fine-tuning, you would need to compare to using forward prop for computing $\\nabla_v f(x)$, and compare both the time, memory usage, and the resulting performance. Otherwise I don't think it is a compelling argument to apply these methods to fine-tuning.\n\n### Forward-Mode AD: Flow and Memory Usage\n\nIn forward-mode automatic differentiation, to compute the directional derivative $\\nabla f(x)^T d $, we propagate both the primal values $ x_\\ell $ and tangent vectors $\\dot{x}_\\ell $ through each layer of the function composition $f = f_L \\circ \\cdots \\circ f_1 $:\n\n$x_{\\ell} = f_{\\ell}(x_{\\ell-1}),$\n\n$\\dot{x}\\_\\ell = J\\_\\ell(x\\_{\\ell-1}) \\dot{x}\\_{\\ell-1}$\n\nwhere $J\\_{\\ell}$ is the Jacobian of the $\\ell$th layer.\n**Memory usage** is minimal: only the current pair $(x_{\\ell-1}, \\dot{x}_{\\ell-1})$ must be stored to compute the next. Both values and tangents can be discarded once the next layer is evaluated, making forward-mode memory-efficient for directional derivatives.\n\n\n**Minor points**\n\n1. Line 625, if $a_n = \\mu_{n+1}-\\mu_n$ then $\\sum_{n=1}^{\\infty} a_n = -\\mu_1$ instead of $\\mu_1.$ The proof here is nonetheless correct", "limitations": "No explicit mention of limitations.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper is very clear, the quality of the writing, and proofs is high. The unbiased zero order estimators are mostly new, excluding a 2020 paper by Chen.", "was_revised": true, "final_justification": "The authors engaged with my main issue, which was including forward mode differentiation as a baseline. They found limitations when using forward mode (I.e flash attention), were transparent on the potential benefits of using forward mode differentiation and promised to include forward mode in their background material. I was very pleased with open and scientific attitude of the authors, and also as mentioned in my initial review, see many merits in the paper. So I’m increasing my score to an accept.", "rebuttal": "We deeply appreciate the reviewer’s careful reading and constructive feedback. We have fixed the minor issues and typos raised by the reviewer. Below, we provide our point-by-point responses to each of the questions raised.\n\n1. **Misleading statement: “optimal complexity for smooth non‑convex objectives” could be read as first‑order optimality.**\n\n    **Response:** We appreciate the reviewer's suggestion. We agree that our sentence here is indeed misleading. We will revise it by correcting the misuse of citing  [Duchi et al. 2015] from the first-order to the zeroth-order optimization lower bound and re-wording the statement of applying [Khaled and Richtarik]. The revised paragraph is now given as the follow:\n\n    > (Revised paragraph) This complexity has matched the lower bound of solving a smooth non-convex optimization problem using **zeroth-order** gradient-based method [Duchi et al., 2015] and cannot be further improved without adding additional assumptions. **Though we directly apply the result from Khaled and Richtárik [2022] (which is applicable for all unbiased estimators), the zeroth-order estimation can result in an additional dependence on the dimension $d$; this dependence has been reflected in our upper bound.**  \n\n2. **(1) Comparison to other zero-order methods used for fine-tuning, such as Mezo in [Malladi et al 2023].** \n\n    **(2) Off a factor $2$ in the two-point estimator.**\n\n    **Response:** We appreciate the reviewer's careful reading and insightful comments. \n\n    - **Regarding the MeZO algorithm**, MeZO is not a new gradient estimator; instead, it is a memory-efficient implementation of the zeroth-order SGD to make existing gradient estimator (e.g. the two-side two-point estimator) a feasible approach to fine-tune LLMs. Before MeZO, people may enumerate each layer and assign the estimated gradient to `param.grad`. After the assignment is done, people take `optimizer.step()` to update the parameter exactly as the first-order method does. The MeZO paper has a smart observation: it is not necessary to save the estimated gradient to `param.grad`; instead, it suffices to save the random seed and the function value difference $[f(x+\\mu v)-f(x-\\mu v)]/(2\\mu)$; when we need to recover the gradient (for the gradient update of each layer), we just need to use the random seed to re-generate the random vector $v$. This approach will significantly reduce the memory cost of the zeroth-order optimization approach.\n\n        As the result, our two-side two-point estimator is the same as the one used in the MeZO paper (with using the zeroth-order batch size $b=2$; we also note that in the original MeZO paper, $b$ is set to be $1$).  \n\n    - **Regarding the omitted factor $2$**, we appreciate the reviewer's careful reading; it is fortunate that it is just a typo, and in our experimental implementation (located in the class given below), we have correctly included this factor. \n\n        - The synthetic experiment: `synthetic-experiments/exp_dim_comparison.ipynb`: the `CentralizedUniformEstimator` class. \n        - The LLM training experiment: `llm-fine-tuning/zo-bench/trainer.py`: the `zo_step` method in the `OurTrainer` class. \n\n        We will fix this typo by adding the factor $2$ in the table and on page 14. \n\n3. **Forward‑mode AD may achieve directional derivatives with low memory and better performance.**\n\n    **Response:** We appreciate the reviewer's thoughtful input. We enjoy reading this comment and learn a lot from it. We agree with the reviewer that there could be other memory-efficient approach (e.g. Forward-mode AD) that could be better than the existing zeroth-order optimization approach. We would be glad to add a discussion of this topic in the background. \n\n    As the goal of our paper is to present a deeper understanding of the unbiased zeroth-order gradient estimator, it is indeed out of the scoop of our paper to include experiments for comparing the performance of Forward-mode AD with the zeroth-order optimization approach on LLM fine-tuning applications. Here we would like to recap our theoretical contributions: (1) we identify the conditions under which constructing an unbiased estimator is possible; (2) we derive the optimal structure that achieves the best possible complexity; and (3) we determine the minimal number of function evaluations required to obtain an unbiased gradient estimator, improving the best-known result from $4$ evaluations [Chen2020] to $1$. We sincerely hope the reviewer will take these theoretical insights into account when assessing our submission."}, {"review_id": "S3ebf76FCE", "reviewer": "Reviewer_Ea9x", "summary": "The central contribution of this paper is the introduction of a novel family of provably unbiased zeroth-order gradient estimators, which the authors term Pk-estimators (where k = 1, 2, 3, 4 denotes the number of function queries required). A detailed variance analysis reveals that while the P1 estimator may exhibit infinite variance, the P2, P3, and P4 estimators are shown to possess finite variance under certain conditions.", "questions": "Clarification is requested regarding the plots in Figure 2. The x-axis of the right panel is labeled \"function calls,\" while the left panel shows the \"number of updates.\" For the one-sided and two-sided two-point estimators, one would expect the total number of function calls to be twice the number of gradient updates. However, the plots for these estimators conclude at approximately 80,000 function calls respectively, rather than the expected 40,000 corresponding to 20,000 updates. Could you please explain this apparent discrepancy? If I am wrong, please correct me.", "limitations": "N/A", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "- The paper's objective is clearly articulated, and the manuscript is well-written, making the core concepts and contributions accessible.\n\n- The novelty of the proposed technique is a point of concern, as it appears to share similarities with methods presented in Chen (2023). Furthermore, a more robust justification for the pursuit of an unbiased estimator is needed. Given that the resulting optimization algorithm demonstrates a convergence rate comparable to methods using biased estimators, what is the fundamental advantage of this approach? Why should one prefer this unbiased estimator over a simpler, biased alternative?", "was_revised": true, "final_justification": "The author presents an unbiased estimate of the gradient of the original function, which is mathematically meaningful; however, its contribution to the machine learning community remains subject to debate. Moreover, the rebuttal addresses most of my concerns. Therefore, I have decided to raise my score.", "rebuttal": "We deeply appreciate the reviewer’s constructive feedback. Below, we provide our point-by-point responses to each of the concerns raised.\n\n## Responses to weakness\n\n- (**Similarity to Chen (2020)**) The novelty of the proposed technique is a point of concern, as it appears to share similarities with methods presented in Chen (2023).\n\n    **Response:** The telescoping structure is the same as the one used in Chen (2020); however, our contribution goes much beyond this telescoping structure and we provide a deeper understanding to the unbiased zeroth-order gradient estimator:\n\n    1.  **Identify when we can use unbiased gradient estimator:** We identify the condition under which the telescoping structure admits a valid expectation representation (Proposition 2.1). This condition is critical for constructing unbiased estimators, but has never been established in Chen (2020) or other prior work. \n\n    2. **More general unbiased gradient estimator \\& Reduce the number of function evaluations from 4 to 1:** The estimator in Chen (2020) is a special case of our $P_4$-estimator. Our framework extends beyond this, answering a fundamental theoretical question: *What is the minimal number of function evaluations needed to construct an unbiased gradient estimator?* We improve the known answer from 4 give by Chen (2020) to 1. Although our one-point estimator comes with a negative result (this is also a new result), we believe it remains of independent theoretical interest.\n\n    3. **Identify the necessary and sufficient condition of achieving the optimal variance:** Moreover,  one of our key focuses is identifying optimal parameter sequences $\\\\\\{\\mu_n\\\\}$ and $\\\\{p_n\\\\\\}$ (Theorem 3.2), rather than proposing a specific estimator. This result indicates that we cannot simply take an arbitrary sequence as being done in Chen (2020), instead, we need to require them to follow a specific algebraic structure. \n\n    These substantial contributions are beyond the shared telescoping structure between Chen (2020) and our work, and we sincerely hope the reviewer will take these theoretical insights into account when assessing our submission.\n\n    \n\n- (**Advantage compared to biased estimator**) Furthermore, a more robust justification for the pursuit of an unbiased estimator is needed. Given that the resulting optimization algorithm demonstrates a convergence rate comparable to methods using biased estimators, what is the fundamental advantage of this approach? Why should one prefer this unbiased estimator over a simpler, biased alternative?\n\n    **Response:** We thank the reviewer for raising this important question. We highlight two key advantages and one potential advantage of using an unbiased estimator: \n\n    (i) **Practical benefit:** *First*, using the unbiased gradient estimator removes the need to balance bias and variance as a coupled trade-off. This simplifies the hyperparameter tuning process, allowing us to focus solely on controlling variance.  As indicated by our Theorem 3.2, the variance upper bounds of the $P_3$ and $P_4$ estimators are parameter-agnostic; that is, they do not introduce any additional hyperparameters compared to classical biased estimators. *Second*, our empirical results indicate that the reduced biased could benefit the gradient estimation, which could result in better performance in training neural network or other applications of zeroth-order optimization. \n\n    (ii) **Theoretical benefit:** The zero bias enables the direct application of standard optimization theory without needing to account for the complications introduced by estimator bias. While our current analysis focuses on vanilla SGD (as in Corollary 3.5), the same theoretical framework readily extends to a broad range of optimization methods or learning theory results, in which the theoretical analysis only relies on the unbiased estimation of the true gradient. \n\n    (iii) **Potential benefit:** The unbiased estimator opens a novel and underexplored direction in the design of zeroth-order gradient estimators. While existing (typically biased) estimators have been extensively studied and refined, the space of unbiased estimators remains largely under-explored. We believe our work creates new opportunities for developing more effective and theoretically grounded zeroth-order estimators.\n\n## Responses to questions\n\n- (**Figure 2 discrepancy: function calls vs. updates seem mismatched**) Clarification is requested regarding the plots in Figure 2. The x-axis of the right panel is labeled \"function calls,\" while the left panel shows the \"number of updates.\" For the one-sided and two-sided two-point estimators, one would expect the total number of function calls to be twice the number of gradient updates. However, the plots for these estimators conclude at approximately 80,000 function calls respectively, rather than the expected 40,000 corresponding to 20,000 updates. Could you please explain this apparent discrepancy? If I am wrong, please correct me.\n\n    **Response:** We appreciate the reviewer's careful reading. We would like to clarify that the x-axis is correctly labeled here, as we are using the mini-batch two-point gradient estimator to align the number of function evaluations to $4$; that is, we increase the batch size to ensure that each gradient update takes the same number of function evaluations as the $P_4$-estimator.  More explicitly, the exact number is described as the following: \n    - For the two-side two-point estimator, we evaluate four points: $f(x+\\mu v_1)$, $f(x-\\mu v_1)$. $f(x+\\mu v_2)$, $f(x-\\mu v_2$). Then we estimate the gradient using $\\frac{1}{2}[\\frac{f(x+\\mu v_1) - f(x-\\mu v_1)}{2 \\mu}v_1+ \\frac{f(x+\\mu v_2) - f(x-\\mu v_2)}{2 \\mu}v_2]$. It corresponds to the batch size $b=2$. \n    - For the one-side two-point estimator, we evaluate four points: $f(x)$, $f(x+\\mu v_1)$, $f(x+\\mu v_2)$, $f(x+\\mu v_3)$ Then we estimate the gradient using $\\frac{1}{3}[\\frac{f(x+\\mu v_1) - f(x)}{ \\mu}v_1+ \\frac{f(x+\\mu v_2) - f(x)}{ \\mu}v_2+ \\frac{f(x+\\mu v_3) - f(x)}{ \\mu}v_3]$. It corresponds to the batch size $b=3$. \n\n    We will also update the figure to include the case of batch size $b=1$ in the revision for a direct comparison."}, {"review_id": "PdgAHhEgrS", "reviewer": "Reviewer_Yvd5", "summary": "This paper introduces a framework for constructing unbiased estimators of the true directional derivative of a smooth function using only zeroth-order queries. The authors reformulate the directional derivative as an infinite series and, leveraging the absolute convergence of this series, derive an expectation-based formulation. Within this framework, they propose and analyze specific estimators, providing theoretical guarantees on the expected stationarity gap in the context of non-convex stochastic optimization.", "questions": "1. The result that unbiased estimators can attain the same level of variance in the smooth setting as standard single- or two-point estimators is certainly interesting. However, it remains unclear how this observation contributes either theoretically or practically. Specifically, does exact unbiasedness offer any concrete advantages over the existing analyses of zeroth-order optimization that rely on the gradient of a smoothed function?\n\n2. The construction of these estimators involves selecting a countably infinite set of distributions $p_n$ and corresponding smoothing parameters $\\mu_n$. It is unclear how this construction guarantees exact unbiasedness in practice. Could this be interpreted instead as an approximation or a limiting form of unbiasedness? If so, why would this approach be preferable to, or more suitable than, using standard smoothed gradient estimators in algorithmic applications?", "limitations": "Yes. The authors have adequately addressed the limitations of this work.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\n\n1. The paper introduces a general framework for analyzing a class of zeroth-order estimators that estimate directional derivatives in an unbiased manner. This provides a unifying perspective on various estimation strategies.\n2. The observation that an unbiased single-point estimator of the true gradient can have infinite variance—contrasting with the bounded variance of zeroth-order estimators for smoothed functions—is both surprising and insightful.\n\n**Weaknesses:**\n\n1. The contributions and practical implications of the work are not clearly articulated. The motivation for introducing this framework, as well as its advantages over existing methods, remains vague. Furthermore, the writing lacks a clear sense of direction regarding potential applications or future developments based on the proposed theory.\n2. The proposed technique relies on the construction of a countably infinite set $p_n$ and corresponding smoothing parameters $\\mu_n$ for $n = 1, \\cdots, \\infty$. It is not immediately clear how this formulation translates into a practical implementation. Further clarification on how this can be approximated or truncated in practice would be valuable. If truncation is necessary, then using standard smoothed gradient estimators—which are better understood and widely adopted—might be more practical and theoretically grounded.\n3. The experimental setup is limited in scope. To better demonstrate the practical utility of the proposed method, it would be beneficial to include experiments on larger-scale models or more diverse benchmark tasks that highlight tangible improvements over existing approaches.", "was_revised": true, "final_justification": "My initial concern was centered on the novelty and significance of the work. I had mistakenly assumed that the authors were directly encoding the directional derivative as an infinite sum, which would offer limited improvements and result in a biased estimator—an approach already explored in prior literature. However, the authors clarified this point using the example of the $P_4$ estimator, demonstrating how *randomness in $n$* allows the sum to represent the directional derivative *in expectation*. This leads to an unbiased estimator for zeroth-order oracles using a finite number of function evaluations.\nGiven the involvement of random variables with complex dependencies, the authors appropriately employed more advanced analytical tools to rigorously justify their results. This careful treatment adds to the work’s credibility and technical depth. I now recognize this contribution as novel within the context of zeroth-order optimization. In particular, it opens the door to deriving high-probability bounds for zeroth-order estimators, which justifies my revised evaluation from a 4 (weak accept) to a 5 (accept).\n\nMoreover, despite introducing randomness in $n$, the authors still provide performance guarantees that match those of existing zeroth-order methods. This highlights the practical viability of their approach, while also expanding the scope of zeroth-order optimization into the realm of unbiased stochastic optimization.", "rebuttal": "We appreciate the reviewer’s constructive feedback. Below, we provide our point-by-point responses to each of questions raised. \n\n1. (**how this observation contributes either theoretically or practically?  any concrete advantages?**) The result that unbiased estimators can attain the same level of variance in the smooth setting as standard single- or two-point estimators is certainly interesting. However, it remains unclear how this observation contributes either theoretically or practically. Specifically, does exact unbiasedness offer any concrete advantages over the existing analyses of zeroth-order optimization that rely on the gradient of a smoothed function?\n\n    **Response:**   We thank the reviewer for raising this important question. We highlight two key advantages and one potential advantage of using an unbiased estimator: \n\n    (i) **Practical benefit:** *First*, using the unbiased gradient estimator removes the need to balance bias and variance as a coupled trade-off. This simplifies the hyperparameter tuning process, allowing us to focus solely on controlling variance.  As indicated by our Theorem 3.2, the variance upper bounds of the $P_3$ and $P_4$ estimators are parameter-agnostic; that is, they do not introduce any additional hyperparameters compared to classical biased estimators. *Second*, our empirical results indicate that the reduced biased could benefit the gradient estimation, which could result in better performance in training neural network or other applications of zeroth-order optimization. \n\n    (ii) **Theoretical benefit:** We especially highlight the theoretical benefit of considering the unbiased gradient estimator. The zero bias enables the direct application of standard optimization theory without needing to account for the complications introduced by estimator bias. While our current analysis focuses on vanilla SGD (as in Corollary 3.5), the same theoretical framework readily extends to a broad range of optimization methods or learning theory results, in which the theoretical analysis only relies on the unbiased estimation of the true gradient. \n\n    (iii) **Potential benefit:** The unbiased estimator opens a novel and underexplored direction in the design of zeroth-order gradient estimators. While existing (typically biased) estimators have been extensively studied and refined, the space of unbiased estimators remains largely under-explored. We believe our work creates new opportunities for developing more effective and theoretically grounded zeroth-order estimators.  \n\n    To illustrate the potential benefit (and the *'potential applications or future developments based on the proposed theory'* noted in the **Weakness 1**), we would like to clarify that in recent years, most zeroth-order gradient estimators have been constructed using the random smoothing technique. This approach approximates the objective function by a smoothed version of the form $\\mathbb{E}_{v \\sim V}[f(x + \\mu v)]$. In contrast, our work is based on a fundamentally different principle: we express the gradient $\\nabla f(x)$ as the expectation of a random variable that depends solely on function evaluations. Here we list two **potential future developments** based on this new expectation representation perspective (developed in our Section 2.1): \n    * > (Designing a new gradient estimator using frequency information): If we replace the random variable solely depending on function evaluations with the random variable depending on frequency of gradient (i.e. use Fourier series to replace the Taylor series), we can access the gradient using the information from the frequency domain. Our work can motivate a branch of further developments in using other side signals to recover the gradient.  \n    * > (Reducing the memory usage of Adam-like optimizer): In Adam, the first- or second- order momentum is represented as a (truncated) series. If we replace it with the expectation representation (described in our Section 2.1), we can recover the whole trajectory using a single item of this series, avoding to save the momentum information on GPU. \n\n    We will include these future directions in our discussion and hope it can resolve your concern described in *Weakness 1*. \n\n2. (**how selecting a countably infinite set of hyperparameters?**) The construction of these estimators involves selecting a countably infinite set of distributions and corresponding smoothing parameters . It is unclear how this construction guarantees exact unbiasedness in practice. Could this be interpreted instead as an approximation or a limiting form of unbiasedness? If so, why would this approach be preferable to, or more suitable than, using standard smoothed gradient estimators in algorithmic applications?\n\n    **Response:**   We appreciate the reviewer's insightful comment. These two parameters  $\\\\\\{p_n\\\\\\}$  and  $\\\\\\{ \\mu_n\\\\\\}$ are indeed the main obstacle of implementing the unbiased zeroth-order gradient estimator, and our work exactly contributes to solve this problem.  These two sequences form an infinite-dimensional hyper-parameter space; by intuition, it must be hard to tune these parameters. However, we make a key observation derived in Theorem 3.2: if we require them to satisfy a simple algebraic relation; e.g. $p_n = \\frac{\\mu_n - \\mu_{n+1}}{\\mu}$ (for the $P_3$- and $P_4$-estimator), then it is guaranteed to achieve the optimal complexity. Moreover, the condition presented in our Theorem 3.2 is the necessary and sufficient condition. That is, if the presented relation is not satisfied, the variance must be sub-optimal (even though they may still achieve the optimal complexity in the asymptotic sense).  As the result, introducing these two parameters do not add additional difficulty to tune the hyper-parameter. \n\n    We also have included two examples (Example 3.3 and Example 3.4) to illustrate how to obtain $\\\\\\{p_n\\\\\\}$  and  $\\\\\\{ \\mu_n\\\\\\}$  from this relation. To better illustrate its practical implementation, we further explain it here: Typically, we can select an existing distribution over $\\mathbb{N}$ (e.g. the Geometric distribution or the Zipf's distribution) that satisfies (1) we can directly sample from using `numpy`, `scipy`, or `pytorch`, and (2) we can access the density function $p_n$ through  these packages. Then the estimation procedure is decribed as follows:\n    1. We sample $i \\sim  \\\\\\{p_n\\\\\\}$ (the distribution we just selected) and calculate $\\mu_i$ using our algebraic relation derived from our Theorem 3.2. \n    2. We plugging it back to our gradient estimator formula. \n\n    Then we directly use this estimator as the true gradient. The truncation occurs only if the numerical stability is absent by a specific $\\mu_i$ (i.e. $\\mu_i$ is too small). This could also happen for a classical zeroth-order estimation; though the classical method uses a constant $\\mu$, at some point $x$ or for some data $\\xi$, this specific $\\mu$ could be too small. Therefore, the practical implementation doesn't present significant difference between the unbiased estimator and other gradient estimators. We also hope this clarification could resolve your concern presented in **Weakness 2**. \n\nWe thank the reviewer again for these insightful comments. In summary, the goal of our work is to present a deeper understanding of the unbiased zeroth-order gradient estimator, which is a topic that has been rarely explored. From our perspective, we have successfully achieved our goal: (1) we identify the conditions under which constructing an unbiased estimator is possible; (2) we derive the optimal structure that achieves the best possible complexity; and (3) we determine the minimal number of function evaluations required to obtain an unbiased gradient estimator, improving the best-known result from 4 evaluations to 1. We sincerely hope the reviewer will take these theoretical insights into account when assessing our submission."}]

Author rebuttal:

Meta-review: This paper proposes a family of unbiased zeroth-order gradient estimators based on a telescoping series representation of directional derivatives. The authors establish the conditions under which unbiasedness is achievable, derive optimal perturbation stepsizes and sampling distributions, and prove that their P3/P4 estimators achieve the optimal  complexity for smooth non-convex optimization. They show the minimal number of function evaluations for unbiased estimation can be reduced from 4 to 1. The rebuttal effectively addresses concerns: unbiasedness is achieved via randomization over an infinite sequence but w/o requiring infinite evaluations, and the method remains practical due to constraints. Most reviewers have raised their scores after discussion, confirming its merit. I recommend acceptance.
Decision: Accept (spotlight)