Bilevel Optimization under Unbounded Smoothness: A New Algorithm and Convergence Analysis

Weakness 1. Too many loops in Algorithm 1,2 which may lead to time-consuming.

A: There is a minor issue in the pseudo-code for Algorithm 2 in the original version of the paper and we are sorry for the ambiguity and confusion brought to the readers. The outer loop in Algorithm 2 with the iteration variable $k$ spanning from $0$ to $K-1$ is indeed repeatedly defined, and we already removed it in the new version of our paper. Algorithm 1 is actually a double-loop algorithm instead of a triple-loop one. Notably, we only need to execute the inner loop occasionally since updating for the lower-level variable is periodic rather than in each iteration, which helps reduce the running time.

Weakness 2. Some symbols are confusing. Both Algorithms 1, and 2 use $K$. In lemmas, what is $K_0$?

A: Please check the algorithm framework in new version of our paper, we fix a typo in the pseudo-code for Algorithm 2 (remove the unnecessary loop in line 1). As illustrated before, we only need $K$ in Algorithm 1, and we do not need it in Algorithm 2. As for the notation $K_0$ in lemmas, please check "Remark" paragraph below Assumption 2 in the main text, we mentioned that $K_0$ and $K_1$ are relaxed smoothness parameters for the function $\Phi(x)$. For the definition of $K_0$ and $K_1$, please check Lemma 9 in Appendix C for details, where we define $K_0$ and $K_1$ in equation (16) and we prove that function $\Phi(x)$ satisfies the standard relaxed smoothness condition.

Weakness 3. Too many hyperparameters in algorithms. More experimental results of different hyperparameters are needed.

A: Nearly all papers about bilevel optimization have many hyperparameters, and our problem setting is more difficult to solve due to Assumption 1, i.e., $(L_{x,0},L_{x,1}, L_{y,0}, L_{y,1})$-smoothness assumption. We conduct careful ablation studies to explore the impact of hyperparameter $I$ (the update period of the lower-level variable) and $N$ (the number of iterations for updating the lower-level variable during each period). The experimental results in Figure 7 (a) (in Appendix J) show that the performance of BO-REP algorithm decreases slightly when increasing the update period $I$ from $2$ to $8$ while fixing the number of iterations in each update period $N$. That demonstrates empirically our algorithm is not sensitive to the update period hyperparameter $I$. When the value of $I$ is too large ($I \geq 16$), we observe a significant performance degradation. In our experiments in the main text, we choose $I=2$ for all experiments and get good performance universally. The ablation result for $N$ is shown in Figure 6 (b), where the update period $I$ is fixed. The figure shows that the performance of BO-REP algorithm would be better as $N$ increases. The algorithm can achieve the best performance when $N\geq 5$. In our experiments, it is good enough to choose $N=3$ to achieve good performance universally for all tasks.

Due to these ablation studies, it means that the algorithm does not need lots of tuning efforts despite these hyperparameters (e.g., $I$ and $N$ ): some default values of $I$ and $N$ (e.g., $I = 2, N = 3$) work very well for a wide range of tasks in practice.

Dear Reviewer gGYZ,

Thank you for reviewing our paper. We have carefully answered your concerns on the nested loops, single-loop, hyperparameters, the new baseline SABA, and verification of unbounded smoothness.

Please let us know if our answers accurately address your concerns. If our response resolves your concerns, we kindly ask you to consider raising the rating of our work. Thank you very much for your time and efforts! We would like to discuss any additional questions you may have.

Best, Authors

Weakness 1. The theoretical complexity result of the proposed algorithm, Theorem 1, lacks persuasiveness. Specifically, the choice of the iteration number $N$ in Algorithm 2 for lower-level variable updates is left unspecified. Furthermore, the selection of the parameter $I$ appears improper, particularly in comparison to the selection of the parameter $K$. This peculiarity arises because, according to Theorem 1, $I=\mathcal{O}(1/\epsilon^2)$, while $K=\mathcal{O}(1/\epsilon)$. Consequently, for sufficiently small values of $\epsilon$, $I$ is significantly larger than $K$. However, both the implementation of the proposed algorithm and the theoretical analysis, as exemplified by Lemma 2, necessitate $I$ to be smaller than $K$.

A: Thank you for pointing it out. We actually state the choice of the iteration number $N$ in the proof of Lemma 2, but without explicitly pointing it out in statement of Lemma 2. We have fixed it in the new version of the paper. For the exact value of $N$, you can check equation (45) and (50) in proof of Lemma 2 in Appendix D.3. The iteration complexity $N$ in Algorithm 2 for lower-level variable updates (within a single subroutine) should be (by eq. (45) and (50)) $N = 64^2\sigma_{g,1}^2K_0^2\left(\lambda+\sqrt{\lambda+1}\right)^2/\mu^2\epsilon^2$, where $\lambda$ is some logarithmic factors. Thus we conclude $N=\widetilde{\mathcal{O}}(\sigma_{g,1}^2K_0^2/\mu^2\epsilon^2)$ in Algorithm 2.

For Theorem 1, $I$ is indeed much smaller than $K$ for small enough $\epsilon$. Note that we choose update period $I$ to be $I = \sigma_{g,1}^2K_0^2/\mu^2\epsilon^2 = \mathcal{O}(1/\epsilon^2)$, and there is a learning rate dependency (i.e. $1/\eta$) in $K$: $K=\frac{4\Delta}{\eta\epsilon}=\mathcal{O}(1/\epsilon^4)$ because $\eta=\mathcal{O}(\epsilon^3)$. In this case, $I$ is much smaller than $K$ for small enough $\epsilon$.

Weakness 2. Despite the authors' assertion that they do not require the boundedness of the norm of the gradient of the upper-level objective function, $\\|\nabla_yf(x,y)\\|$, as stipulated in Assumption 2 (i), it is worth noting that they still demand the boundedness of $\\|\nabla_y f(x, y^\*(x))\\|$. The rationale behind the relative ease of achieving this condition in practical scenarios, as opposed to the former requirement, remains unclear. Moreover, it remains unaddressed whether the problems examined in the numerical experiments exhibit unbounded smoothness yet satisfy Assumption 2 (i).

A: Thank you for this insightful question. First of all, to the best of our knowledge, this is an assumption (i.e., $\\|\nabla_yf(x,y^*(x))\\|\leq M$) which is required by all bilevel optimization papers [1, 2, 3]. Also, the boundedness of $\\|\nabla_y f(x, y^*(x))\\|$ does not contradict with relaxed smoothness assumption, i.e., Assumption 1. We will illustrate it in the following. As can be seen from Figure 4 in Appendix G, we show that the smoothness parameter scales linearly in terms of the gradient norm of $x$ (parameter of the recurrent layers) and $y$ (parameter of the last linear layer).

In addition, we show that if we denote $x$ as previous layers of recurrent neural networks and $y$ as the last layer (i.e., classification layer), the condition $\\|\nabla_yf(x,y^*(x))\\|\leq M$ indeed holds both theoretically and experimentally. Figure 4 (b) in Appendix G shows that the gradient norm in terms of $y$ is bounded, i.e., the maximal value of the horizontal axis can be bounded by $M=10$. Besides, when considering the setting of hyper-representation, we have included a proof in Appendix H, indicating that this condition holds true for any $(x, y)$ and thus also holds for $(x, y^*(x))$.

[1] Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. arXiv preprint arXiv:1802.02246, 2018.

[2] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International conference on machine learning, pp. 4882–4892. PMLR, 2021.

[3] Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert D Nowak. A fully first-order method for stochastic bilevel optimization. In International Conference on Machine Learning, pp. 18083-18113. PMLR, 2023.

Question 1. It is important for the authors to acknowledge that not all bilevel optimization papers require the assumption $\\|\nabla\_yf(x,y^*(x)\\|\leq M$. For instance, the work Shen & Chen (2023).

A: The work Shen & Chen (2023) indeed makes this assumption, but with a different notation. It's worth noting that for Assumption 1 in Shen & Chen (2023), the author assumes that '' There exists $L>0$ that given any $x\in\mathcal{C}$, $f(x,\cdot)$ is $L$-Lipschitz continuous on $\mathcal{U}(x)$", which actually implies $\\|\nabla\_yf(x,y)\\|\leq L$ for any given $x\in\mathcal{C}$ under their notation. The difference is that Shen & Chen (2023) does not assume strong convexity of the lower-level problem, so the optimal solution set $\mathcal{S}(x)$ for the lower-level problem may not be unique. But since for any given $x\in\mathcal{C}$, $\\|\nabla\_yf(x,y)\\|\leq L$ holds for all $y\in\mathcal{U}(x)$, then this condition also holds for any $y\in\mathcal{S}(x)$.

In fact, this condition is quite important for the analysis in Shen & Chen (2023). For example,

• Formula (A.2) for proof of Theorem 1 in Appendix A.1, on page 21
• Formula (A.8) and (A.10) for proof of Theorem 6 in Appendix A.2, on pages 23 and 24
• $L$-dependency for parameter choice $\epsilon=L^2\mu/\gamma^2$, $\gamma\geq L\sqrt{3\mu\delta^{-1}}$ in part (a), and $\epsilon=L^2/(4\gamma^2\sigma^2)$, $\gamma\geq\max\\{L\sqrt{2\mu\delta^{-1}}, L\sqrt{\delta^{-1}}/\sigma\\}$ in part (b) for proof of Proposition 2 in Appendix B.2, on page 25
• Formula (B.5) for proof of Theorem 7 in Appendix B.3, on page 26
• Formula (B.33) for proof of Theorem 8 in Appendix B.4, on page 32
• Formula (C.2) for proof of Proposition 3 in Appendix C.1, on page 33
• $L$-dependency for parameter choice $\gamma^*=L\rho$ in formula (D.3) for proof of Theorem 9 in Appendix D.1, on page 37
• Formula (D.6) for proof of Theorem 10 in Appendix D.1, on page 38

Question 2. The validation of Assumption 2(i) in Appendix H does not substantiate one of the motivation underlying this work, which asserts that Assumption 2(i) is more practical to satisfy than the requirement that the norm of the gradient of the upper-level objective function, $\\|\nabla\_yf(x,y)\\|$, is bounded. Appendix H demonstrates that the latter condition holds for the specific hyper-representation problem.

A: Please note that the assumption "$\\|\nabla\_yf(x,y^*(x))\\|\leq M$ for any $x\in\mathbb{R}^{d\_x}$" is weaker than the assumption "$\\|\nabla\_yf(x,y)\\|\leq M$ for any $x\in\mathbb{R}^{d\_x}$ and $y\in\mathbb{R}^{d\_y}$", that is, the latter assumption implies the former assumption. In theory, we only need the former (i.e., weaker) assumption to obtain the convergence result, and in Appendix H we demonstrate that the latter (i.e., stronger) condition holds for the hyper-representation problem, thus the former condition in the theoretical part also holds naturally.

In addition, we do not assert Assumption 2(i) (i.e., "$\\|\nabla\_yf(x,y^*(x))\\|\leq M$ for any $x\in\mathbb{R}^{d\_x}$") is one of the motivation in our paper. Our motivation focuses on addressing unbounded smoothness phenomenon for specific problem settings, and then we introduce the new $(L\_{x,0}, L\_{x,1}, L\_{y,0}, L\_{y,1})$-smoothness assumption under bilevel setting, but existing bilevel optimization algorithms do not have any convergence guarantees in the presence of unbounded smoothness (please read our second paragraph in Section 1 for details). We also do not claim that the assumption "$\\|\nabla\_yf(x,y^*(x))\\|\leq M$ for any $x\in\mathbb{R}^{d\_x}$" is more practical to satisfy than the assumption "$\\|\nabla\_yf(x,y)\\|\leq M$ for any $x\in\mathbb{R}^{d\_x}$ and $y\in\mathbb{R}^{d\_y}$". We would like to emphasize that we only need a weaker assumption for theoretical analysis although we can show that stronger assumption indeed holds in practice.

Weakness 1. This paper relaxes the bounded assumptions but is not based unbounded assumption as stated in the abstract. The example pf RNN and LSTMs are misleading unless the authors are able to prove that these networks satisfy the Assumption 1.

A: Thank you for your suggestion. First, we want to emphasize that the unbounded smoothness assumptions were empirically observed in many prior works ([1, 2]) for RNNs, LSTMs and Transformers. Second, we have empirically verified that the Recurrent Neural Network model satisfies Assumption 1 in the revised version. In Figure 4 in Appendix G, we plot the estimated smoothness at different iterations during training neural networks. In particular, we conduct the hyper-representation experiments to verify Assumption 1. We adopt the same recurrent neural network as that in Section 4.1. The lower-level variable is the parameter of the last linear layer, and the upper-level variable is the parameter of the previous layers (2 hidden layers). In each training iteration, we calculate the gradient norm w.r.t. the upper variable $\\|\nabla_x f(x,y)\\|$ and the lower variable $\\|\nabla_y f(x,y)\\|$, and use the same method as described in Appendix H.3 in [1] to estimate the smoothness of $x$ and $y$. From Figure 4 in Appendix G, we can find that the smoothness parameter scales linearly in terms of gradient norm for both layers. This verifies the Assumption 1 empirically. Note that these results are consistent with the results in the literature, e.g., Figure 1 in [1] and Figure 1 in [2].

Weakness 2. The second stated contribution in introduction section is not valid. As the author also mentions in the footnote 3, the SOTA complexity results under bounded smoothness is $\widetilde{\mathcal{O}}(1/\epsilon^3)$. The proposed method does not achieve the SOTA complexity rate, even after considering the logarithmic factors.

A: As stated in footnote 3, there are a few works which achieve $\mathcal{O}(\epsilon^{-3})$ oracle complexity when the stochastic function is mean-squared smooth [3, 4, 5], but our paper does not have such an assumption. Without mean-squared smoothness assumption, our proposed method is optimal up to logarithmic factors as proved in [6]. We have made it more clear in the abstract as well.

Weakness 3. Momentum technique has been incorporated in the bilevel optimization area before. For example, MRBO [3] achieves the $\widetilde{\mathcal{O}}(1/\epsilon^3)$ after momentum introduction. The authors are expected to improve the convergence rate of the proposed algorithm.

A: Please note that these two momentum techniques are completely different. Our momentum is simply doing moving-average estimate for history stochastic gradients, but the momentum used in [3] is a type of recursive momentum which is based on the STORM estimator [7]. The STORM estimator requires the mean-squared smoothness assumption while our momentum does not. Without this assumption, our complexity result is already optimal up to logarithmic factors and it is impossible to achieve the $\widetilde{\mathcal{O}}(1/\epsilon^3)$ rate as proved in [6]. We could also improve our complexity result to $\widetilde{\mathcal{O}}(1/\epsilon^3)$ by imposing the mean-squared smoothness assumption and adopting the STORM technique, but it is out of the scope of our current paper.

[1] Jingzhao Zhang, Tianxing He, Suvrit Sra, and Ali Jadbabaie. Why gradient clipping accelerates training: A theoretical justification for adaptivity. International Conference on Learning Representations, 2020b.

[2] Michael Crawshaw, Mingrui Liu, Francesco Orabona, Wei Zhang, and Zhenxun Zhuang. Robustness to unbounded smoothness of generalized signsgd. Advances in neural information processing systems, 2022.

[3] Junjie Yang, Kaiyi Ji, and Yingbin Liang. Provably faster algorithms for bilevel optimization. Advances in Neural Information Processing Systems, 34:13670-13682, 2021.

[4] Zhishuai Guo, Quanqi Hu, Lijun Zhang, and Tianbao Yang. Randomized stochastic variance-reduced methods for multi-task stochastic bilevel optimization. arXiv preprint arXiv:2105.02266, 2021.

[5] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A near-optimal algorithm for stochastic bilevel optimization via double-momentum. Advances in neural information processing systems, 34:30271–30283, 2021.

[6] Yossi Arjevani, Yair Carmon, John C Duchi, Dylan J Foster, Nathan Srebro, and Blake Woodworth. Lower bounds for non-convex stochastic optimization. Mathematical Programming, 199(1-2): 165–214, 2023.

[7] Ashok Cutkosky and Francesco Orabona. Momentum-based variance reduction in non-convex sgd. Advances in neural information processing systems, 32, 2019.

Question 2. Is the achievement of $\widetilde{\mathcal{O}}(1/\epsilon^3)$ convergence rate is always based on such mean-squared smooth or have to be completed with STORM based technique? If not, the contribution is still not valid.

A: Yes. All the bilevel optimization literature (Yang et al., 2021; Guo et al., 2021; Khanduri et al., 2021) with $\widetilde{\mathcal{O}}(1/\epsilon^3)$ complexity use individual smoothness assumption, which is even stronger than the mean-squared smoothness assumption as we mentioned in answer to Question 1. Also, their algorithm design depends on the STORM technique (Cutkosky & Orabona, 2019). In addition, it is shown in Arjevani et al. (2023) that without this assumption, $\Omega(1/\epsilon^4)$ complexity is necessary for stochastic first-order optimization algorithms for single level problems. This indicates that our complexity is optimal up to logarithmic factors for the bilevel problems we considered in this paper.

Let us give more details to explain this. First we consider potentially non-convex single level optimization problems as considered in Arjevani et al. (2023) . Given differentiable objective function $F: \mathbb{R}^d \rightarrow \mathbb{R}$, our goal is to find an $\epsilon$-stationary point using stochastic first-order algorithms. Assume the algorithms access the function $F$ through a stochastic first-order oracle consisting of a noisy gradient estimator $\nabla F: \mathbb{R}^d\times\Xi \rightarrow \mathbb{R}^d$ and distribution $P_\xi$ on $\Xi$ satisfying

For potentially non-convex single level optimization problems, Theorem 3 in Arjevani et al. (2023) states that

• Under Assumption (5), any stochastic first-order methods requires $\Omega(1/\epsilon^4)$ queries to find an $\epsilon$-stationary point in the worst case (please check Contribution 1 in Section 1.1 on page 2, and formula (23) in Theorem 3 on page 14 in Arjevani et al. (2023) for details);
• Under Assumptions (5), plus mean-squared smoothness assumption, any stochastic first-order methods require $\Omega(1/\epsilon^3)$ queries to find an $\epsilon$-stationary point in the worst case (please check Contribution 2 in Section 1.1 on page 2, and formula (24) in Theorem 3 on page 14 in Arjevani et al. (2023) for details).

Note that our problem class is more expressive than the function class considered in Arjevani et al. (2023) and hence our problem is harder. For example, if we consider an easy function where the upper-level function does not depend on the lower-level problem and does not have mean-squared smoothness, the $\Omega(1/\epsilon^4)$ lower bound can be applied in our setting. Therefore the oracle complexity $\widetilde{\mathcal{O}}(1/\epsilon^4)$ achieved in our paper is already optimal up to logarithmic factors.

Question 1. What is the definition of mean-squared smooth you are referring? Is it related to outer function or inner function? Is it with regard to $x$ or $y$?

A: For single level problems, given differentiable objective function $F: \mathbb{R}^d \rightarrow \mathbb{R}$, we say the function $F$ satisfies mean-squared smoothness property (formula (4) on page 4 in Arjevani et al. (2023)) if for any $x,y\in\mathbb{R}^d$ and any $\xi\sim P_{\xi}$,

where we use $\xi$, $\nabla F(x;\xi)$ and $L$ to denote the random data sample, the stochastic gradient estimator and the mean-squared smoothness constant (the corresponding notations in Arjevani et al. (2023) are $z$, $g(x,z)$ and $\bar{L}$, please check formula (4) on page 4 in Arjevani et al. (2023) for details).

A slightly stronger condition than mean-squared smoothness is the individual smoothness property (please check the statement "We assume that $f(x,\xi_t)$ is differentiable, and $L$-smooth as a function of $x$ with probability 1." in Section 3 on page 3 in Cutkosky & Orabona (2019) for details). We say function $F$ has individual smoothness property if for any $x,y\in\mathbb{R}^d$ and any $\xi\sim P_{\xi}$,

For bilevel problems, in order to obtain $\widetilde{\mathcal{O}}(1/\epsilon^3)$ oracle complexity bound, Yang et al. (2021), Guo et al. (2021) and Khanduri et al. (2021) actually require individual smoothness assumption jointly in $(x,y)\in\mathbb{R}^{d_x}\times\mathbb{R}^{d_y}$ for both upper-level (i.e., outer function $f(x,y)$) and lower-level problems (i.e., inner function $g(x,y)$). To be more specific, for upper level problem they require for any $u=(x,y), u'=(x',y')\in\mathbb{R}^{d_x}\times\mathbb{R}^{d_y}$ and any $\xi$,

and for lower level problem they require for any $u=(x,y), u'=(x',y')\in\mathbb{R}^{d_x}\times\mathbb{R}^{d_y}$ and any $\zeta$,

These three inequalities above, namely formulas (3) and (4), can be viewed as definition of individual smoothness property under bilevel optimization setting.

For more details, please check Assumption 2 in Section 3.1 on page 6 in Yang et al. (2021), Assumption 2 in Section 2 on page 5 and Assumption 4 in Section 4 on page 8 in Guo et al. (2021), Assumption 2 in Section 2 on page 4 in Khanduri et al. (2021). Please note that our paper does not assume any of these assumptions.

Thanks for your further responses. It addressed my concern and I increase my score to 6.

Note that all we discussed are not included in the paper. I strongly encourage you to write an individual paragraph to explain it.

Besides, I read your mentioned works and it looks like variance reduction technique can also help to improve the convergence rate. I suggest you to discuss both STORM and variance reduction techniques in your revision.

Weakness. The main drawback of the proposed algorithm is that the structure is complex, and each subroutine contains additional hyperparameters that need to be tuned.

A: Our problem setting is much more challenging to solve than the standard bilevel optimization for smooth objectives. The reason is that, under Assumption 1, i.e., $(L_{x,0},L_{x,1}, L_{y,0}, L_{y,1})$-smoothness assumption, we have to carefully update both upper-level and lower-level variables. Therefore our proposed algorithm has more complex structure compared to other bilevel algorithms. However, we want to emphasize that our algorithm is not sensitive to additional hyperparameters such as $I$ and $N$ (please see our ablation study in Appendix J). In fact, for each subroutine of the periodic updates, we do not need to carefully tune these values in practice: some default values of $I$ and $N$ (e.g., $I=2, N=3$) do not seem to hurt the performance for a wide range of bilevel optimization tasks.

Question 1. In initialization refinement step, what's the purpose to invoke a multi-epoch algorithm to obtain an accurate estimate of lower-level variables and why periodic update doesn't use such a multi-epoch strategy?

A: For the initialization refinement step, to achieve $\\|y_0-y_0^*\\|\leq\epsilon/8K_0$ with high probability guarantee, the iteration (and also oracle) complexity for multi-epoch SGD (by detailed statement of Lemma 1 in Appendix D.2) and SGD (by Lemma 12 in Appendix D.2) are $\widetilde{\mathcal{O}}(1/\epsilon^2)$ and $\widetilde{\mathcal{O}}(1/\epsilon^4)$, respectively. Thus multi-epoch SGD algorithm has a better convergence rate compared to SGD for initialization refinement.

When $k$ is a multiple of $I$, the periodic update step performs actual updates for the lower-level variable $y$. Note that the initial distance between $y_k$ and $y_k^*$ satisfies $\\|y_k-y_k^*\\|\leq\epsilon/4K_0$ with high probability. By the detailed statements of Lemma 12 (SGD) and Lemma 1 (multi-epoch SGD) in Appendix D.2, and Lemma 2 in Appendix D.3 (periodic updates), we can see that the high probability bound for SGD benefits from the small initial distance $\\|y_k-y_k^*\\|$ (e.g., the polynomial dependency in terms of $D$ in Lemma 12 in formula (20)), while the high probability bound for multi-epoch SGD only marginally benefits from the small initial distance (e.g., the logarithmic dependency of $\mathcal{V}\_0$ in Lemma 1 in formula (32)). Therefore, with $\\|y_k-y_k^*\\|\leq\epsilon/4K_0$, both multi-epoch SGD and SGD guarantee $\\|y_{k+1}-y_k^*\\|\leq\epsilon/8K_0$ with high probability within $\widetilde{\mathcal{O}}(1/\epsilon^2)$ iterations. So we choose SGD instead of multi-epoch SGD for the periodic update step due to its easy implementation in practice.

I would like to thank the authors' detailed responses and explanations which have addressed my previous concerns. So I will increase my score to 8.

Weakness 1. Lack of lower bounds does not allow to understand optimality of the proposed method under given assumption.

A: Thank you for your insightful question. We want to emphasize that our oracle complexity of finding an $\epsilon$-stationary point is indeed optimal up to logarithmic factors by directly applying the lower bound results in [1]. Let us give more details here. For the upper-level problem, we assume the stochastic gradient estimator is unbiased and has bounded variance. Under this stochastic gradient oracle, the lower bound for finding $\epsilon$-stationary points in non-convex smooth single-level stochastic optimization with first-order oracle is at least $\Omega(1/\epsilon^4)$ (Theorem 1 of [1]). Since the single-level optimization is a special case of bilevel optimization (e.g., bilevel optimization becomes a single-level optimization when the upper-level function does not depend on the optimal solution of the lower-level problem) and standard smoothness is a special case of relaxed smoothness (e.g., relaxed smoothness with $L_1=0$ correspond to standard smoothness), our $\widetilde{\mathcal{O}}(1/\epsilon^4)$ complexity result is optimal up to logarithmic factors. We have added descriptions in the main paper to make it more clear.

Weakness 2. Proposed guarantees holds only for finding an $\epsilon$-stationary point in expectation.

A: Please note that we do not assume light tail assumption for upper-level problems in our paper, so our convergence guarantees hold only in expectation. Note that the high probability result even in the single-level stochastic nonconvex optimization setting requires the light tail assumption (e.g., Corollary 2.5 in [2]).

Question 1. Is it possible to relax the strong convexity assumption of the lower function to log-concavity?

A: Strong convexity guarantees that the optimal solution $y^*(x)$ is unique for the lower-level problem: this is very important for our analysis. However, this property does not always hold for log-concave functions. For example, let $g(x,y)=(\sin{x}+2)\exp(-y^2/2)$, then $g(x,y)$ is log-concave in $y$ for any $x$, and $g(x,y)$ also satisfies other assumptions in Assumption 2. For this specific function, $y^*(x)$ does not even exist, since for any fixed $x$, we could not find a minimum for $y$ if $y\in\mathbb{R}^{d_y}$.

Question 2. What are challenges to provide high-probability guarantees for $\epsilon$-stationary point for the given algorithm?

A: For Corollary 2.5 in [2], it requires stochastic estimators to be unbiased and have bounded variance, and also the light-tail assumption for stochastic gradient noise to provide high probability bound for single-level non-convex problems. The main challenge for bilevel optimization is the following: even if we have all of these assumptions in bilevel problems, the hypergradient estimator is still biased (check Lemma 4 in our paper for details) and therefore the techniques in the single-level stochastic optimization may not apply.

[1] Yossi Arjevani, Yair Carmon, John C Duchi, Dylan J Foster, Nathan Srebro, and Blake Woodworth. Lower bounds for non-convex stochastic optimization. Mathematical Programming, 199(1-2): 165–214, 2023.

[2] Saeed Ghadimi and Guanghui Lan. Stochastic first and zeroth-order methods for nonconvex stochastic programming. SIAM Journal on Optimization, 23(4):2341–2368, 2013.