Provable Faster Zeroth-order Method for Bilevel Optimization with Optimal Dependency on Error and Dimension

We sincerely appreciate the reviewers' valuable suggestions regarding the writing of our paper, which have been very helpful in improving our work. Following your recommendations, we have revised the structure of the paper, included comparisons with more algorithms in the tables, and added additional experiment.

Q1: Can the authors explain the motivation of studying zeroth-order bilevel algorithms? The experiments can not reflect if VRZSBO can still be a better choice in large-scale settings

1.Bilevel is a highly general framework that can model many real-world problems. Min-max problems and single-level optimization problems are both special cases of bilevel problems. Zeroth-order optimization has already been extensively studied in single-level and min-max problems, attracting significant interest. Therefore, studying zeroth-order algorithms under the bilevel framework to obtain more general conclusions compared to the aforementioned cases is highly appealing.
2.Moreover, reformulating some zeroth-order optimization tasks using the bilevel framework can potentially achieve better results. For example, as shown in [1] and [2], combining adversarial attacks (a common application of zeroth-order algorithms) with meta-learning (a common bilevel application) led to improved performance over existing methods. Given that [3] demonstrates the use of zeroth-order algorithms for fine-tuning large language models (LLMs), achieving up to 12× memory reduction and up to 2× GPU-hour reduction, and [4] highlights the close connection between the bilevel framework and LLM , zeroth-order bilevel algorithms hold significant potential for large-scale applications. Our VRZSBO algorithm, with its strong theoretical guarantees and low complexity, can serve as an important baseline for future research in zeroth-order bilevel optimization.

Q2 Can the authors explain why $(d _1 + d _2)$ dependency is optimal? For example, is there a reference for the lower bound of the dimension dependency.

Thank you for raising this question. We answer it in the General Response .

Q3: the cost in (Aghasi and Ghadimi, 2024) is not in expectation   while our oracle cost is in expectation

This problem is actually results in the difference between PAGE[5] and Spider[6].To achieve a single-loop algorithm, PAGE selects a large batch size with a probability of $p$, while SPIDER uses a large batch size once every $1/p$ iterations. This results in both methods having the same complexity, but PAGE's results are presented in an expected form. This reflects a trade-off between the single-loop structure and randomness. We do not consider this a significant issue, as in single-level problems, comparing the complexities of SPIDER and PAGE is not regarded as unfair.

Reference：

[1]Boosting Black-Box Adversarial Attacks with Meta Learning

[2]Simulating Unknown Target Models for Query-Efficient Black-box Attacks

[3]Fine-Tuning Language Models with Just Forward Passes

[4]LLMS ARE HIGHLY-CONSTRAINED BIOPHYSICAL SE QUENCE OPTIMIZERS

[5]PAGE: A Simple and Optimal Probabilistic Gradient Estimator for Nonconvex Optimization

[6]SPIDER: Near-Optimal Non-Convex Optimization via Stochastic Path Integrated Differential Estimator

Q: Based on the proofs of Theorems 1 and 2, the iteration count T in these theorems depends on $\psi _{0}$, which may in turn depend on n (in the finite-sum case) and ϵ (in the expectation case). Therefore, to derive the complexity analysis results (i.e., the total oracle cost), does this require certain conditions on the initial points? Additionally, does it require an assumption of a lower bound for Φ(x)? A bit more discussion on these points would be helpful.

Good question. recall that

For the variance terms $\mathbb{E}[||v _{x, 0}-\hat{h} _{x, 0}|| ^{2}],\mathbb{E}[|| v _{y, 0} -\hat{h} _{y, 0} || ^{2}],\mathbb{E}[|| v _{z, 0} -\hat{h} _{z, 0} || ^{2}]$, by computeing $v _{z, 0}=\frac{1}{B}\sum _{i=1} ^{B} D\hat{H} _{v}(t+1;\xi _{i})-\hat{\nabla} _{y}f(x _{t+1}, y _{t+1};\zeta _{i})$, $v _{x,0}=\frac{1}{B}\sum _{i=1} ^{B}\hat{\nabla} _{x}f(x _{t+1}, y _{t+1};\zeta _{i})-D\hat{J} _{v}(t+1;\xi _{i})$, $v _{y,0}=\frac{1}{B}\sum _{i=1} ^{B}\hat{\nabla} _{y}g(x _{t+1}, y _{t+1};\xi _{i})$ with batchsize $B=n$ (finite-sum case) and $\mathcal{O}(\epsilon ^{-2})$ (the expectation case) when $t=0$, the variance terms $\mathbb{E}[||v _{x, 0}-\hat{h} _{x, 0}|| ^{2}],\mathbb{E}[|| v _{y, 0} -\hat{h} _{y, 0} || ^{2}],\mathbb{E}[|| v _{z, 0} -\hat{h} _{z, 0} || ^{2}]$ will decrease to 0（finite-sum case）and $\mathcal{O}(1)$ （expectation case）, which makes $\mathbb{E}[\psi _0]$ independt of $n$ and $\epsilon$, we will add this initialization step to our algorithm. The remaining terms $\Phi(x _{0})$ , $|| y _{0}- y ^{\star} _{0}|| ^{2}$ and $|| z _{0}- z ^{\star} _{0}|| ^{2}$ , to the best of our knowledge, appear in almost all bilevel optimization literature. Since $y _{0} ^*$ and $z _{0} ^{\star}$ are solutions to strongly convex problems, we can run algorithm $\log(\epsilon)$ times to decrease $|| y _{0}- y ^{\star} _{0}|| ^{2}$ and $|| z _{0}- z ^{\star} _{0}|| ^{2}$ to nearly 0 .This trick is called "warm start", which is effective in certain experiments. And of course we require an assumption of a lower bound for $\Phi(x _{0})$, which is standard in bilevel optimization, we will add this assumption to improve clarity.

Q: In Theorems 1 and 2, the zeroth-order tuning parameter h does not appear—why is this the case? Is there a relationship between v and h, potentially stemming from Corollary 1 and Lemma 3 through an intermediate parameter δ? Additionally, in Figure 1, ablation studies examine the effects of different choices of h and v on the zeroth-order estimator. Were these experiments conducted with consideration of the relationship between v and h?

Your understanding is correct. In Theorems 1 and 2, $h$ does not explicitly appear because its effect is incorporated into $v$ through the intermediate parameter $\delta$, as shown in Corollary 1 and Lemma 3, we will include the value of $h$ in the statements of Theorems 1 and 2 to improve clarity. In the experiments, we fix v to the best-performing value, then test h , and then swap by fixing h to the best-performing value and testing v .

Q: The theoretical results in this work suggest that smaller zeroth-order tuning parameters v and h yield better performance, yet the experimental results in Figure 1 do not support this conclusion. Could the authors explain this discrepancy?

Good question. We would like to clarify that, theoretically, the value of v indeed improves with smaller values. However, for h, there is an optimal value: it should not be too large or too small, as shown in the proof of Lemma D.6. In our experiments, the optimal value of h corresponds to a middle-range value, which is consistent with our theoretical analysis. As for v , the experimental trend does not fully align with our theoretical expectations, fortunately, we observed that the best performance in the experiments still corresponds to the smallest value of v. Based on our experience, selecting a sufficiently small v (e.g., 0.001 or 0.00001) remains effective.

Reference：

[1]On the Complexity of First-Order Methods in Stochastic Bilevel Optimization

[2]Using Complex Variables to Estimate Derivatives of Real Functions

We sincerely thank you for taking the time to provide such detailed and thoughtful feedback. Your suggestions are helpful for improving our paper. We have supplemented the experiments and corrected the typos based on your comments in the revised version, and we will now address the another questions you raised.

W: finite-difference approaches are often impractical due to potential inaccuracies in function evaluations (e.g., rounding errors), which prevent the distance of finite differences from being arbitrarily small.

1.To the best of our knowledge, achieving the best theoretical upper bound for a Hessian-free method appears challenging without employing finite differences. First-order methods based on penalty functions can at best achieve a complexity of $\mathcal{O}(\epsilon ^{-4})$ ([1], ICML 2024), whereas finite-difference-based method achieves a significantly improved complexity of $\mathcal{O}(\epsilon ^{-3})$. 2.Although we did not employ this approach in our paper, we would like to introduce an intersting method that can effectively address the issue of inaccuracies in function evaluations, provided that complex number computations are acceptable[1]. Using a Taylor’s series expansion ( $e _{i}$ denote the $i$ th standard basis vector，v is a smoothing parameter, z is a vector):

and

we can estimate the gradient and hessian-vector product as : $\nabla f(x)=[Im(f(x+\mathrm{i} ve _{1})),Im(f(x+\mathrm{i} ve _{2})),\dots,Im(f(x+\mathrm{i} ve _{d}))] ^T$ and $\nabla f(x)z=Im(\nabla f(x+\mathrm{i} vz))$. This estimation method is not subject to rounding errors because it doesnt involve subtraction,serves as a good complement to our approach, particularly when the estimation might be affected by rounding errors.

 Q: Why does the proposed algorithm achieve an optimal result with respect to dependence on dimension? A bit more discussion is needed.

Thank you for raising this question. We answer it in the General Response.

Q: In Assumption 2, why is there no stochastic assumption on $\nabla _{xy} ^2 g(x,y)$?

Sorry about this typo. We have corrected it.

Thank you for taking the time to review our paper. We appreciate the opportunity to clarify our contributions and provide a detailed comparison between our method and ZDSBA.

Q: Gaussian Smoothing vs. Coordinate-wise Smoothing

To simplify notation, let $d = d _1 + d _2$. The improved efficiency of our method with respect to dimensionality $d$ stems from two key factors: (1) Our single-loop structure eliminates the inner loop in [1], saving $\mathcal{O}(d)$ complexity from solving inner problems. (2) We use a finite-difference approximation for Hessian-vector or Jacobian-vector products, requiring only two gradient evaluations, with $\mathcal{O}(d)$ oracle cost. This is significantly more efficient than the Gaussian smoothing method in [1, Proposition 2.5], where approximating the Hessian-vector product incurs a variance of $\mathcal{O}(d ^3)$, leading to $\mathcal{O}(d ^3)$ complexity for hypergradient estimation and $\mathcal{O}(d ^3)$ outer loop iterations. Our method reduces this cost by $\mathcal{O}(d ^2)$. Combining the two above, our complexity is $d ^3$ faster than that of [1].

1. Gradient Estimation: Coordinate-wise smoothing requires $\mathcal{O}(d)$ oracle calls for accurate gradient estimation, while Gaussian smoothing uses $\mathcal{O}(1)$ calls but introduces a variance of $\mathcal{O}(d || \nabla f(x) || ^2)$, ultimately requiring $\mathcal{O}(d)$ oracle calls to cancel this variance. Hence, the overall complexity for both remains $\mathcal{O}(d)$ times that of standard first-order methods.
2. Hessian-Vector/Jacobian-Vector Estimation: Finite-difference approximation via coordinate-wise smoothing requires $\mathcal{O}(d)$ oracle calls for accurate estimation, while [1, Proposition 2.5] requires $\mathcal{O}(1)$ calls but introduces a variance of $\mathcal{O}(d ^3)$ due to term $(d+4)(d+2)\left||\nabla _{x x} ^2 q\right|| _F ^2=\mathcal{O}(d ^3\left||\nabla _{x x} ^2 q\right|| _2 ^2)$, leading to $\mathcal{O}(d ^3)$ complexity in the outer loop. Our method reduces this by $\mathcal{O}(d ^2)$. Conclusion: Combining the $\mathcal{O}(d ^2)$ improvement in Hessian-vector/Jacobian-vector estimation and the $\mathcal{O}(d)$ speedup from the single-loop structure, our method achieves an overall $\mathcal{O}(d ^3)$ improvement compared to ZDSBA [1]. Additionally, our coordinate-wise results suggest the potential for Gaussian smoothing to achieve $\mathcal{O}(1)$ samples with $\mathcal{O}(d)$ variance, as both smoothing approaches perform similarly in gradient estimation, we have incorporated this comparison into the latest version of the paper.

W: Lack of Novelty Thank you for your comment. We have already addressed a similar question raised by Reviewer L5fV, and our response can be found in their section. We kindly refer you to that response for a detailed explanation, moreover, we would like to emphasize that our differences with SPABA also extend to comparisons with other methods that achieve near-optimal complexity [2],[3]. We believe that our answer should address your concern as well.

Reference： [1]Fully Zeroth-Order Bilevel Programming via Gaussian Smoothing

[2]Achieving O (ϵ−1.5) Complexity in Hessian/Jacobian-free Stochastic Bilevel Optimization

[3]A Lower Bound and a Near-Optimal Algorithm for Bilevel Empirical Risk Minimization

We appreciate the reviewer’s thoughtful comments and feedback. In response to the concerns raised, we address the following points:

W: Concerning the use of the coordinate-wise zeroth-order estimator in Definition 1.

We understand the reviewer’s concern regarding the practical applicability of the coordinate-wise zeroth-order estimator. In response to this, we offer the following points:

1.We chose the coordinate-wise zeroth-order estimator in order to achieve stronger theoretical results.As discussed in our response to Reviewer Nm7C, the use of the coordinate-wise zeroth-order estimator enables us to achieve a complexity of $\mathcal{O}(d _{1}+d _{2})$. This result serves as an important baseline, as it matches the complexity of the most well-known zeroth-order min-max methods, which are often used in practice.

2. Early Stage of Zeroth-Order Bilevel Optimization Research We remind the reviewers that the first fully zeroth-order bilevel optimization algorithm was only recently proposed, so research in this area is still in its early stages. Given this, the improvement in our work from $(d _{1}+d _{2}) ^4$ to $d _{1}+d _{2}$ is a significant advancement. Our theoretical analysis provides valuable insights that future researchers can use to design methods tailored to specific problems. For example, a hypergradient estimator requiring only $\mathcal{O}(1)$ function evaluations or potentially even achieving faster speeds than first-order methods in practical tasks.

3. Moreover, reformulating some zeroth-order optimization tasks using the bilevel framework can potentially achieve better results. For example, as shown in [10] and [11], combining adversarial attacks (a common application of zeroth-order algorithms) with meta-learning (a common bilevel application) led to improved performance over existing methods.

I appreciate the authors' detailed response. However, I still have reservations regarding the use of the coordinate-wise zeroth-order estimator. This coordinate-wise finite-difference approach makes it not much different from a standard gradient oracle, in terms of both the computational cost and convergence analysis. While I understand that the authors aim to establish a baseline for zeroth-order methods in bilevel optimization, it is well-known that gradients can be well approximated using $O(d)$ function value queries, and thus a first-order method can be simulated using a zeroth-order oracle. In some sense, the contribution of this paper appears to be formalizing this observation for bilevel optimization problems.

In addition, while there are some technical differences in the convergence analysis between this paper and SPABA, ultimately they appear to achieve similar convergence rates. It remains unclear whether any of these differences are essential for the algorithm to work properly with a coordinate-wise zeroth-order estimator.

As a side note, I noticed that the two plots in Figure 3 appear to be identical. Intuitively, one would expect a larger dimension to result in a longer runtime. Could the authors double-check if this is correct?