On the Optimal Construction of Unbiased Gradient Estimators for Zeroth-Order Optimization

We appreciate the reviewer’s constructive feedback. Below, we provide our point-by-point responses to each of questions raised.

1. (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?

   Response: We thank the reviewer for raising this important question. We highlight two key advantages and one potential advantage of using an unbiased estimator:

   (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.

   (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.

   (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.

   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):

   • (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.

   • (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.

   We will include these future directions in our discussion and hope it can resolve your concern described in Weakness 1.

2. (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?

   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.

   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:

   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.
   2. We plugging it back to our gradient estimator formula.

   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.

We 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.

I would sincerely like to thank the authors for clearing out my concerns. This completely clears out all doubts on my end. I apologize for the mistake on my end. I did not realize that $n$ itself is a random variable, which in case would imply that $n$ has infinite support not an infinite number of parameters, as mistakenly stated in my earlier comment. I think that this work is novel and very interesting now and would like to raise my score from a 3 to a 5.

We deeply appreciate the reviewer’s constructive feedback. Below, we provide our point-by-point responses to each of the concerns raised.

Responses to weakness

• (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).

  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:

  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.

  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.

  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.

  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.

• (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?

  Response: We thank the reviewer for raising this important question. We highlight two key advantages and one potential advantage of using an unbiased estimator:

  (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.

  (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.

  (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.

Responses to questions

• (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.

  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:

  • 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$.
  • 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$.

  We will also update the figure to include the case of batch size $b=1$ in the revision for a direct comparison.

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.

1. Misleading statement: “optimal complexity for smooth non‑convex objectives” could be read as first‑order optimality.

   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:

   (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.

2. (1) Comparison to other zero-order methods used for fine-tuning, such as Mezo in [Malladi et al 2023].

   (2) Off a factor $2$ in the two-point estimator.

   Response: We appreciate the reviewer's careful reading and insightful comments.

   • 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.

     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$).

   • 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.

     • The synthetic experiment: synthetic-experiments/exp_dim_comparison.ipynb: the CentralizedUniformEstimator class.
     • The LLM training experiment: llm-fine-tuning/zo-bench/trainer.py: the zo_step method in the OurTrainer class.

     We will fix this typo by adding the factor $2$ in the table and on page 14.

3. Forward‑mode AD may achieve directional derivatives with low memory and better performance.

   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.

   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.

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.

Response to weakness

• (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.

  • 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.
  • 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.

  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.

  (Function value noise) If is available with noise like in:

  • Gasnikov, Alexander, et al. "The power of first-order smooth optimization for black-box non-smooth problems." International Conference on Machine Learning. PMLR, 2022.

  In this case the construction seems to be quite sensitive and may work worse than classical one approach.

• 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.

  (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.

Response to questions

• Why do authors use such a finite difference? It seems that symmetric one will be better according to [Shamir2017].

  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: (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]. (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.

  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.

  • [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.
  • [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.

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.

Response to questions

• (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.

  Response: We thank the reviewer for raising this important question. We highlight two key advantages and one potential advantage of using an unbiased estimator:

  (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.

  (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.

  (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.

• (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.

  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.

  • [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.

• (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.

  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$.

  1. What is the value of this infinite series?

     We can calculate the finite-sum $S_n := \sum_{i=1}^n p_i x_i$; then we let $n$ tend to infinite.

  2. What is the value of an expectation of the random variable?

     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]$.

  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.

• (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.

  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.

• (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?

  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.

  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.

  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.

Response to weaknesses

• (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.

  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.

  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.

  • [MeZO] Malladi, Sadhika, et al. "Fine-tuning language models with just forward passes." Advances in Neural Information Processing Systems 36 (2023): 53038-53075.