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

Authors: Yuxing Peng, Yuanyuan Liu, Fanhua Shang, Hongying Liu, Zhouchen Lin
Primary area: optimization
Keywords: Stochastic bilevel optimization, Hessian-free algorithms, near-optimal complexity

Abstract:
In this paper, we study and analyze zeroth-order stochastic approximation algorithms for solving black-box bilevel optimization problems, where only the upper and lower function values can be obtained. \citep{Saeed2024} proposed the first full zeroth-order bilevel method that utilizes Gaussian smoothing to estimate the first- and second-order partial derivatives of functions with two independent blocks of variables. However, this method suffers from a high dimensional dependency of $\mathcal{O}((d_{1}+d_{2})^{4})$, where $d_{1}$ and $d_{2}$ are the dimensions of the outer and inner problems, respectively. They left an open question: can this dimension dependency be improved? To answer this question, we propose a single-loop accelerated zeroth-order bilevel algorithm, which achieves a dimension dependency of $\mathcal{O}(d_{1}+d_{2})$ by incorporating coordinate-wise smoothing gradient estimators (coord). 
	We develop a new theoretical analysis for the proposed algorithm, which converges to a stationary point of $\Phi(x)$ with a complexity of $\mathcal{O}((d_{1}+d_{2})\epsilon^{-3})$ in expectation settings and $\mathcal{O}((d_{1}+d_{2})\sqrt{n}\epsilon^{-2})$ in finite sum settings. These complexities are both best-known with respect to dimension and error $\epsilon$. We also provide  experiment to validate the effectiveness of the proposed algorithm.

Decision: reject
Reviewer scores: 5, 3, 6, 6

Reviews and discussion:
<thread>
creation date:
1730002019683

thread_type:
official_review

<thread_object>
creator: reviewer_syfa

<content>
summary:
This paper considers bilevel optimization. In particular, the authors aim at solving the stochastic bilevel optimization via zeroth-order methods. Motivated by the fact that a recent work has high dimensional dependency in the complexity bound, the authors propose VRZSBO, and prove that it requires less oracle complexity than the existing work, in the context of zeroth-order bilevel optimization algorithms. Numerical simulations are provided to validate the performance of their algorithm under different settings and their better convergence as compared to existing algorithms.

strengths:
1. The authors provide a better upper bound for the oracle complexity of zeroth-order bilevel optimization, without the need of additional assumptions.

2. The technical claims are in general sound.

weaknesses:
1. Incomplete related work. Some important bilevel optimization algorithms are missing, even if they are not based on zeroth-order methods. It would be good if the authors could include a more comprehensive table containing existing algorithms with comparisons of assumptions, complexities, and types of oracles needed. See, for example, the related works in table 1 in [1].

2. Presentation. 
    1. the authors claim one contribution is that they do not have additional assumptions like 3rd Lipschitz. However, this is not well reflected in Table 1. 
    2. It seems that Section 1.1 can be moved to other sections for better presentation — for those who are familiar with bilevel optimization this section seems unnecessary and for those who are unfamiliar with bilevel optimization this section is full of technical details and hard to follow.
    3. It would be better if the authors could provide a well-organized Section 3 — the details of analysis described in 3.5 - 3.7 can be moved to the appendix to make the main context more readable. For example, proof sketch in section 3.6 is fairly standard and unnecessary to be included in the main body of the paper.

3. There are many important baselines missing in the experiments. 

4. Overfull box in lines: 257 - 262, 409, 453, etc.

Reference

[1] BILEVEL OPTIMIZATION UNDER UNBOUNDED SMOOTHNESS: A NEW ALGORITHM AND CONVERGENCE ANALYSIS

questions:
1. 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.

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

3. The oracle cost mentioned in Remark 2 is in expectation (see (14) - (16)). However, the cost in (Aghasi and Ghadimi, 2024) is not in expectation, as the estimator therein does not require an additional sampling step (i.e., (14) - (16)). I am wondering if the authors could provide more context on this unfair comparison.
</content>
</thread_object>

</thread>
<thread>
creation date:
1730571224205

thread_type:
official_review

<thread_object>
creator: reviewer_h91m

<content>
summary:
Using finite-difference approximations and the PAGE method for nonconvex optimization, this work introduces a novel zeroth-order algorithm, VRZSBO, for solving bilevel optimization problems with a strongly convex lower-level problem. The authors show that VRZSBO converges to an $\epsilon$-stationary point with a function query complexity of $\mathcal{O}((d_1 + d_2)\sqrt{n} \epsilon^{-2})$ in the finite-sum setting and $\mathcal{O}((d_1 + d_2) \epsilon^{-3})$ in the expectation setting, where $d_1$ and $d_2$ denote the dimensions of the upper- and lower-level variables, respectively. Experimental results on hyper-representation tasks using linear and two-layer network embedding models further validate the algorithm’s effectiveness.

strengths:
$\textbf{S1:}$ The paper is easy to follow. Zeroth-order algorithms for bilevel optimization remain underexplored, even in the nonconvex-strongly-convex setting, making the introduction of the novel VRZSBO algorithm a valuable algorithmic contribution to this area. 

$\textbf{S2:}$ The proposed VRZSBO algorithm achieves a significant  improvement over the best-known complexity of ZDSBA (Aghasi and Ghadimi, 2024), reducing the query complexity from $\mathcal{O}((d_1 + d_2)^4 \epsilon^{-3})$ to $\mathcal{O}((d_1 + d_2) \epsilon^{-3})$ in the expectation setting. Additionally, this work provides a complexity analysis for the finite-sum setting.

weaknesses:
$\textbf{W1:}$ It is well-known that the Hessian-vector product can be implemented through finite-difference approximations of the gradient along a given direction. Given this, the novelty of the proposed results seems limited compared to existing methods. 

Moreover, as noted in Chapter 9 of Numerical Optimization by Jorge Nocedal and Stephen J. Wright, 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.

$\textbf{W2:}$ The experiments could be more comprehensive. Incorporating updated first-order algorithms, such as F2SA (Kwon et al., ICML 2023), would provide a more thorough comparison. Furthermore, including commonly used bilevel optimization tasks, such as hyperparameter selection and data hyper-cleaning, would further validate VRZSBO’s effectiveness, showcasing its reliability across various bilevel optimization problems. 

I would be willing to increase the score based on the rebuttal.

questions:
Apart from the questions raised in the Weaknesses section, some additional questions are as follows:

$\textbf{Q1:}$ Why does the proposed algorithm achieve an optimal result with respect to dependence on dimension? A bit more discussion is needed.

$\textbf{Q2:}$ In Assumption 2, why is there no stochastic assumption on $\nabla_{xy}^2 g$? 

$\textbf{Q3:}$ 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 $\epsilon$ (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 $\Phi(x)$? A bit more discussion on these points would be helpful.

$\textbf{Q4:}$ 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 $\delta$? 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$?

$\textbf{Q5:}$ 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?
 

$\textbf{Q6:}$ In the experiments, why was there no comparison with ZDSBA from Aghasi and Ghadimi (2024)?

$\textbf{Suggestions for improvement that did not affect the score:}$

In the convergence analysis, for example, in Theorem 1, it is stated that $\psi_{t+1}-\psi_{t} \leq -\frac{\eta_x}{2} ( | \nabla \Phi(x_t) |^2 + \epsilon^2)$. This implies that the Lyapunov function $\psi_t$ is decreasing. However, how is this possible when $\Phi(x)$ is nonconvex? It seems that the sign before $\epsilon^2$ may be incorrect.

In Theorem 1, why is there no $m$ in the finite-sum case? Is there an assumption that $m=\mathcal{O}(n)$?

Page 1, Line 045: In $g(x,y)$, $n$ should be $m$. 

Page 2, Lines 102-103: Verify the form of $\bar{\nabla}\Phi(x)$.  

Page 3, Line 145: Replace “$y(x)$, $y_\lambda(x)$, $\nabla \mathcal{L}\lambda(x)$” with “$y^*(x)$, $y\lambda^(x)$, $\nabla \mathcal{L}_\lambda^(x)$”. 

Page 5: Line 228: The referenced equation (2) may not be appropriate. 

Page 5: Lines 230-231: $R(x,y,z)$ does not include an inverse matrix operator.

Page 5, Lemma 1: It would be clearer to replace $y^$ with $y^(x_t)$ since $y^*$ is a function.

Page 5, Definition 1: In $\hat{\nabla}_x g(x,y)$, replace $d_2$ with $d_1$, and in $\hat{\nabla}_y f(x,y)$, replace $d_1$ with $d_2$.  

Page 6, Lemma 3: Clarify the meaning of $r_z$ in the definition of $h$. 

Page 7, Line 327: “Section 3.1” should be changed to “Sections 3.2–3.3”.  

Page 7, Equations (14), (16): Keep the notations consistent with those in (10)–(11).  

Page 7, Line 369: Should $-p$ replace $-\frac{1}{p}$ according to Lemma 9? 

Page 7, Line 372:  Clarify what $\eta$ represents in the definition of $\Psi_t$. 

Page 8, Lines 383-387: Should $-p$ replace $-\frac{1}{p}$ based on Lemmas 6–7? 

Page 8, Line 398: Replace $\eta$ with $\eta_x$.  

Page 8, Lemmas 4-5: Define $\ell_{Z^}$ and $\ell_{Y^}$.   

Page 10, Figure 1: Correct typos, particularly in the caption.
</content>
</thread_object>

</thread>
<thread>
creation date:
1730697415410

thread_type:
official_review

<thread_object>
creator: reviewer_l5fv

<content>
summary:
In this paper, the authors propose a zeroth-order method for solving stochastic bilevel optimization problems, where the lower-level problem is strongly convex and only function values are accessible for both levels. Their method can be regarded as a zeroth-order approximation of the SPABA method proposed by (Chu et al., 2024), where the gradients and Hessian-vector products are replaced by their zeroth-order estimators. Specifically, the authors employ a coordinate-wise zeroth-order estimator, approximating each coordinate of the gradient by finite differences. They analyze their proposed algorithm under both general stochastic and finite-sum settings, demonstrating that the function value query complexity scales linearly with the problem’s dimension and has optimal dependence on the target accuracy.

strengths:
In the context of zeroth-order algorithms for bilevel optimization, the proposed algorithm improves dimensional dependence compared to prior work and introduces a simpler single-loop update format.

weaknesses:
- My major concern is the use of the coordinate-wise zeroth-order estimator in Definition 1. As noted by the authors, this estimator requires $\mathcal{O}(d_1+d_2)$ function value queries in total, which may be as computationally expensive as obtaining the exact gradient by backpropagation. While it could be argued that only function values are accessible in certain applications, this does not hold for most common bilevel optimization applicationss, including the hyper-representation problem used in the experiments. 
- In addition, the novelty of the paper appears somewhat limited. If my understanding is correct, the algorithmic framework closely follows the SPABA method by (Chu et al., 2024) such as the use of the PAGE variance reduction technique, and the main modidfication is to substitute the exact gradients and Hessian-vector products by their zeroth-order estimators. However, since the coordinate-wise zeroth-order estimator is employed, the gradient estimator error can theoretically be made arbitrarily small by selecting a sufficiently small $v$ (see Lemma 2). Indeed, the authors choose $v = \mathcal{O}(\epsilon^2)$ in Theorems 1 and 2, which means that the gradient error is also on the order of $\mathcal{O}(\epsilon^2)$. Consequently, the convergence analysis may not significantly differ from that of the original SPABA method.

----
Tianshu Chu, Dachuan Xu, Wei Yao, and Jin Zhang. SPABA: A single-loop and probabilistic stochastic bilevel algorithm achieving optimal sample complexity. ICML 2024

questions:
- The author claims that the dependence on the dimension is optimal. While this appears reasonable, a rigorous lower bound is necessary to support this claim. 
- The experiments should include the SPABA method as a baseline, as it is the most relevant algorithm for comparison with the proposed method. Additionally, comparisons across different algorithms should be based on runtime, given the potentially significant differences in per-iteration costs.
</content>
</thread_object>

</thread>
<thread>
creation date:
1730702688806

thread_type:
official_review

<thread_object>
creator: reviewer_nm7c

<content>
summary:
This paper explores zeroth-order algorithms for addressing stochastic bilevel optimization problems. Expanding on earlier work that introduced a zeroth-order bilevel method using Gaussian smoothing, the authors enhance the dimensional dependency in sample complexity from $(d_1 + d_2)^4$ to $(d_1 + d_2)$. They also utilize variance reduction techniques to optimize the epsilon dependency in sample complexity from $\epsilon^{-4}$ to $\epsilon^{-3}$, achieving optimal bounds. The new algorithm undergoes theoretical analysis, demonstrating convergence to a stationary point with optimal complexity in both expectation and finite sum setting. Additionally, experiments are conducted to validate the algorithm's effectiveness.

strengths:
The paper presents several key strengths that enhance its contribution to the field of stochastic bilevel optimization. Notably, it significantly improves the dimensional dependency in sample complexity from $(d_1 + d_2)^4$ to $(d_1 + d_2)$ and employs variance reduction techniques to optimize epsilon dependency, achieving optimal bounds. The thorough literature review effectively situates the work within existing research, underscoring its relevance. Additionally, the clear and structured presentation of the methodology supports reproducibility, while the empirical validation through experiments demonstrates the algorithm's practical effectiveness.

weaknesses:
Lack of Novelty: The contributions of this paper build upon established concepts, integrating ideas from quadratic auxiliary functions, PAGE-type variance reduction, and zeroth-order gradient estimation. While these improvements are valuable, they may not be unexpected, in my humble opinion. For instance, ZDSBA built on previous work as a double-looped algorithm utilizing Hessian inverse approximation, so it is understandable that the incorporation of auxiliary quadratic functions would enhance dimension dependency. Additionally, the introduction of various variance reduction techniques is likely to improve epsilon dependency. Although I recognize the meaningful nature of these contributions, I feel that the paper does not fully meet the acceptance criteria for ICLR. Therefore, I will reserve my opinion on acceptance until after the rebuttal stage.

questions:
Gaussian Smoothing vs. Coordinate-wise Smoothing: Based on my previous experiences, I have not observed a significant difference in the theoretical bounds between these two smoothing techniques. Consequently, I believe the improvements presented in the paper may not be attributed to coordinate-wise smoothing. The paper does not address the rationale behind choosing coordinate-wise smoothing over Gaussian smoothing. Could the authors please clarify the reasoning for this choice?
</content>
</thread_object>

</thread>
<thread>
creation date:
1732215862530

thread_type:
official_comment

<thread_object>
creator: authors

<content>
title:
General Response

comment:
We would like to thank all the reviewers for their thoughtful feedback and for raising important questions about our work. Several of these questions share common themes, and we have provided a general response below to address these shared concerns comprehensively.

**G: the optimality of $\epsilon$ and $d _{1},d _{2}$, and lowerbound.**
After we carefully check the lower bounds under our assumptions,we indeed find that our claim on "optimal" is not rigorous. We have decided to clarify the statements related to "optimal" to "best-known", and add a discussion on the dependencies on $\epsilon$ and $d$  in our revised version, based on the following facts:

1. Our dimension dependency $\mathcal{O}(d _{1} + d _{2})$ matches the best-known result in simpler  minmax problems [1, 2], and the $\epsilon$ dependency aligns with the best-known  complexity for first-order NC-SC bilevel problems [3, 4, 5], making our results "best-known."
2. Since [6] establishes a $\Omega(d)$ lower bound for single-level smooth convex problems, which are simpler than NC-SC bilevel problems (e.g., by setting $f(x, y) = f(x)$), a $\Omega(d _{1})$ dependency is inevitable.
3. Our assumptions align with those used to derive the lower bound in [4] for NC-SC  finite-sum bilevel problems, confirming that the $\mathcal{O}(\sqrt{n}\epsilon ^{-2})$ dependency is inevitable in the finite-sum case .
4. The $\mathcal{O}(\epsilon ^{-3})$ lower bound for single-level problems[8] in the expectation case is constructed under the mean-square-smooth assumption, which is slightly stronger than  the typical assumption of smoothness about smoothness of $f(x, y, \xi)$ for all $\xi$. Thus, strictly speaking, $O(\epsilon ^{-3})$ is not necessarily optimal in $\epsilon$ for the expectation case.

Based on existing lower bounds, the dependency on $d _1$ is inevitable. However, the optimality of the $d _2$ dependency remains unclear. It would be particularly interesting if the exponent of $d _2$ could be reduced to less than 1 for lower level strongly convex optimization, as this would imply a theoretical speedup of zeroth-order algorithms compared to first-order algorithms (assuming the gradient computation cost is $\mathcal{O}(d _1 + d _2)$).



[1]Zeroth-order algorithms for nonconvex–strongly-concave minimax problems with improved complexities

[2]Zeroth-Order Alternating Gradient Descent Ascent Algorithms for a Class of Nonconvex-Nonconcave Minimax Problems

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

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

[5]SPABA: A Single-Loop and Probabilistic Stochastic Bilevel Algorithm Achieving Optimal Sample Complexity

[6]Optimal Rates for Zero-Order Convex Optimization: The Power of Two Function Evaluations

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

[8]Lower bounds for non-convex stochastic optimization
</content>
</thread_object>

</thread>


Metareview:
metareview:
This paper studies zeroth-order method for bilevel optimization. In particular, the authors proposed a single-loop accelerated zeroth-order algorithm for bilevel optimization. The reviewers believe that the proposed algorithm closely follows some existing works, and the novelty is limited. Moreover, it requires significant amount of work to improve the presentation and readability of the paper.