We are extremely grateful to the reviewer for their thoughtful questions and are encouraged that they found numerous aspects of our contribution interesting and useful for follow-up papers.

Response to Weaknesses

1. Complexity of Algorithms.

   We agree with the reviewer that some of our results are probably suboptimal in terms of convergence rates. Nonetheless, we wish to point out that our rate for the linear equality case is in fact optimal. Moreover, despite this suboptimality, ours is the first result with dimension-free rates for these settings for our choice of stationarity and oracle. We will add an explicit comparison with the unconstrained setting.

2. Coherence of narrative.

   We thank the reviewer for raising this valid point. Our algorithm for the linear equality case crucially relies on our novel observation that the hyperobjective $F$ is smooth; in contrast, this does not necessarily hold in the linear inequality case. Therefore, the two cases require different approaches, and unifying them would be an interesting future direction. We will clarify this better.

Response to Questions

1. Need for perturbation.

   Because we do not assume a bound on the smoothness of the hyperobjective (a strong assumption by several prior works), it is known [1] to be impossible to optimize $F$ directly in general – in particular by inexact GD, as suggested.

   Since the hyperobjective is Lipschitz (cf. Lemma 4.1), it is differentiable almost everywhere by Rademacher's theorem; we can therefore implicitly optimize its uniform smoothing via random perturbations. This method of using perturbations to avoid non-differentiable points is widely used in the literature on nonsmooth optimization. We will clarify this.

2. Meaning of ${\nabla} F$.

   We thank the reviewer for pointing out this subtlety, which we will better explain: As correctly noted, $F$ can be nonsmooth nonconvex; in this setting, algorithms generally require a subgradient, i.e., an element of the Clarke subdifferential (see Definition 2.1), which is what we mean by $\nabla F$ (which, as the reviewer notes, is not Lipschitz continuous). Furthermore, the $\tilde{\nabla} F$ we construct is such that its distance to some subgradient of $\nabla F$ is at most $\alpha$.

3. Lipschitz notation.

   We acknowledge the confusion caused by our notation and will switch all notations to either $L$-Lipschitz or $S$-smooth.

4. Zeroth-order approach for linear equality setting.

   We apologize for the confusion: Please note that for the linear equality case, our algorithm accesses $\nabla f$ and $\nabla g$, which makes it first-order. So, even though it looks like a zeroth-order approach (due to the finite differences), it is actually a first-order method. We will clarify this.

5. Assumption on $\lambda^\*$.

   Access to $\lambda^\*$ indeed implies access to $y^\*$. That said, the core difficulty in bilevel programs is access to $\frac{d}{dx}y^\*$. Our technical novelty is precisely to approximate this first-order oracle. Moreover, in practice, we can easily approximate $\lambda^\*$ to high accuracy due to the linear convergence of strongly convex problems with linear constraints --- our experiments demonstrate our algorithm's robustness to the use of such a proxy for $\lambda^\*$.

6. Clarifying complexity.

   The notation $\tilde{O}(1)$ in Line 229 refers to $\log(1/\epsilon)$, matching the cost in Line 322 — We apologize for our oversight in not stating this and will fix it. Next, the complexities in Lines 226 and 229 are for different operations: Line 226 is the cost of estimating $\nabla F(x)$ and Line 229 that of estimating $F(x)$. For the linear inequality case, the ACGD method in [2] implies linear rate to obtain $y^\*$ and $\lambda^\*$. We will explicitly state these.

7. Motivating the perturbed lower-level problem.

   For intuition on the perturbation, we demonstrate our method in the unconstrained setting. Let $y_{\delta}^\*:=\text{argmin} (g(x,y)+\delta f(x,y)).$ Optimality of $y_{\delta}^\*$ and differentiating w.r.t. $\delta$ yields $\frac{dy_{\delta}^\*}{d\delta}=-\left(\delta\cdot\nabla_{yy}^{2}f(x,y_{\delta}^{\ast})+\nabla_{yy}^{2}g(x,y_{\delta}^{\ast})\right)^{-1}\nabla_{y}f(x,y_{\delta}^{\ast}), \implies \frac{dy_{\delta}^{\ast}}{d\delta}\vert_{\delta=0} = -\left(\nabla_{yy}^{2}g(x,y^\*)\right)^{-1}\nabla_{y}f(x,y^\*).$

   Therefore, $v_{x}:=\frac{\nabla_{x}g(x,y_{\delta}^{\ast})-\nabla_{x}g(x,y^{\ast})}{\delta}\approx \lim_{\delta\rightarrow0}\frac{\nabla_{x}g(x,y_{\delta}^{\ast})-\nabla_{x}g(x,y^{\ast})}{\delta}=\nabla_{xy}^{2}g(x,y^{\ast})\cdot\frac{dy_{\delta}^{\ast}}{d\delta}\vert_{\delta=0} = \frac{dy^\*}{dx}^{\top} \nabla_y f(x, y^\*)$ is a valid gradient proxy. Thus, the motivation behind the perturbation is to allow for the construction of a hypergradient approximation by an application of the implicit function theorem.

8. Strong convexity of $\mathcal{L}$.

   $\mathcal{L}$ is strongly convex when $\alpha_1 > C_f / \mu$. This may be deduced from the smoothness and strong convexity of the components of $\mathcal{L}$. We will clarify this.

Conclusion

We again thank the reviewer for their effort in raising pertinent questions; we will incorporate these answers in our paper.

Please let us know if any questions remain. If we adequately answered all the questions, we would like to respectfully ask the reviewer to re-consider their score. We are happy to answer any further questions.

References

[1] Kornowski, G., and Shamir, O. Oracle complexity in nonsmooth nonconvex optimization. JMLR (2022).

[2] Zhang, Z., & Lan, G. (2022). Solving Convex Smooth Function Constrained Optimization Is Almost As Easy As Unconstrained Optimization.

We greatly appreciate the prompt response from the reviewer and clarify our point about optimality below.

1. Optimality in the equality setting

   The work of [1] showed the optimal oracle complexity (in the deterministic setting) to be $O(\epsilon^{-2})$ for unconstrained bilevel programs; our result matches this under the more complicated setting with equality constraints.

   More specifically, for our bilevel problem with equality constraints, the upper level objective function is smooth and possibly non-convex, so it covers the more traditional single-level deterministic smooth non-convex optimization as a special case. Under this setting, it has been shown in [2] that the number of oracle evaluations required to find an $\epsilon$-stationary point (i.e., a point $x$ satisfying $\|\nabla f(x)\|\leq \epsilon$) is lower bounded by $O(\epsilon^{-2})$. Thus, the $O(\epsilon^{-2})$ oracle complexity we achieved is unimprovable. This point was also recently stated in the work of [3].

   We will clarify this point better in the paper. Thank you again for your feedback, and please let us know if we can answer any further questions!

References

[1] Chen, L., Ma, Y., & Zhang, J. (2023). Near-optimal fully first-order algorithms for finding stationary points in bilevel optimization. arXiv preprint arXiv:2306.14853.

[2] Carmon, Y., Duchi, J. C., Hinder, O., & Sidford, A. (2020). Lower bounds for finding stationary points I. Mathematical Programming, 184(1), 71-120.

[3] Kwon, J., Kwon, D., & Lyu, H. On The Complexity of First-Order Methods in Stochastic Bilevel Optimization. In Forty-first International Conference on Machine Learning.

We are extremely grateful to the reviewer for their thoughtful questions and are encouraged by their positive assessment of our contributions, soundness, and presentation.

Response to Weaknesses

1. Limited experiments.

We acknowledge the importance of testing our algorithms with large-scale experiments. However, our focus has primarily been theoretical since prior work had far weaker theoretical guarantees for our choice of stationarity and oracle.

In the code, we also have results on bilevel problems with bilinear constraints (in figures/exp_bilinear). Per reviewer pjir's suggestion, we have now implemented additional experiments for Algorithm 2 (please see the top-level response), which we will add to our paper. In follow-up work, we hope to expand our implementation to a wider range of problems and large-scale experiments.

2. Access to $y^\*$

Our algorithm for the inequality setting indeed requires access to $\lambda^\*$ (the lower-level dual solution). In practice, though, we can easily approximate $\lambda^\*$ to high accuracy due to the linear convergence of strongly convex problems with linear constraints --- our experiments demonstrate our algorithm's robustness to the use of such approximations as proxy for the exact solution. That said, we believe that developing an algorithm for this problem without this assumption is an important direction for future work.

Response to Questions

1. Computing $\nabla_x g^*(x)$

We compute this gradient as $\nabla_x g(x,y^\*) + \nabla_x h(x,y^\*)^\top \lambda^\*$.

For some intuition, note that $g^\*(x) := \min_y \max_{\lambda \geq 0} g(x,y) + \lambda^\top h(x,y),$ to which applying Theorem 4.13 from [1] immediately yields the claim. We will add this explanation in the text.

2. Explaining Algorithm 2

Algorithm 2 is a variant of gradient descent with momentum and clipping, where $\tilde{g}\_{t}$ is the inexact gradient, $\Delta\_{t}$ is a clipped accumulated gradient (hence accounts for past gradients, which serve as a momentum), and the main iterate is updated as $x_{t+1}=x_t+\Delta_t$. Similar algorithms have appeared in prior work on nonsmooth nonconvex optimization (e.g. [2]). However, none of them accounted for inexactness in the gradient, crucial in the bilevel setting. We will add this in the final version.

3. Computing $\nabla y^*(x)$

Our algorithm does not require actually computing $\nabla y^*(x)$. The only place where $\nabla y^*(x)$ appears is in our analysis when bounding the difference between our inexact and the actual hypergradients. We use $\nabla y^*(x)$ to mean a Clarke subgradient of $y^{\ast}$. We will clarify this better.

4. Assumption 2.3

We acknowledge that this assumption is somewhat technical — however, prior work also imposes similar assumptions. For instance, Khanduri et al [3] assume (weak/strong) convexity of the hyperobjective $F$ to obtain finite-time guarantees --- our assumption is strictly weaker. Further, Yao et al [4] also assume boundedness of the optimal dual variable.

As to $y^\*(x)$ being Lipschitz, there are some useful settings with this property: e.g., if the lower level feasible region $Y(x):=\\{y: h(x,y) \leq 0\\}$ is independent of $x$, an argument in [5] shows it to be $\frac{C_g}{\mu_g}$ Lipschitz. The situation is more complicated when $Y(x)$ changes with $x$. If $g(x, y)$ is quadratic in $(x,y)$ (e.g., in model predictive control), $y^\*(x)$ is piecewise affine and Lipschitz [6]. For more general settings, this is related to sensitivity analysis, and some constraint qualifications can establish the property [7].

5. cvxpyLayer

cvxpylayer is a Python library [8, 9] that uses a second-order method to compute the derivative of the optimal solution $y^*(x)$ of a parametric convex program with respect to the parameter $x$. We invoke cvxpylayer on the lower-level problem (which is parametrized by $x$) to compute $\frac{dy^*(x)}{dx}$. We will add more details in the text.

Conclusion

We again thank the reviewer for their time and effort in bringing up pertinent questions; we will incorporate these answers in our paper.

Please let us know if any questions remain. If we adequately answered all the questions, we would like to respectfully ask the reviewer to re-consider their score. We are happy to answer any further questions.

References

[1] Bonnans, J. F., & Shapiro, A. (2013). Perturbation analysis of optimization problems.

[2] Cutkosky, A., Mehta, H., & Orabona, F. (2023). Optimal stochastic non-smooth non-convex optimization through online-to-non-convex conversion. ICML

[3] Khanduri, P., Tsaknakis, I., Zhang, Y., Liu, J., Liu, S., Zhang, J., & Hong, M. (2023). Linearly constrained bilevel optimization: A smoothed implicit gradient approach. ICML

[4] Yao, W., Yu, C., Zeng, S., & Zhang, J. (2024) Constrained Bi-Level Optimization: Proximal Lagrangian Value Function Approach and Hessian-free Algorithm. ICLR.

[5] Nesterov, Yu. Smooth minimization of non-smooth functions. Mathematical programming (2005).

[6] Borrelli, F., Bemporad, A., & Morari, M. (2017). Predictive control for linear and hybrid systems.

[7] Minchenko, L. I., and P. P. Sakolchik. Hölder behavior of optimal solutions and directional differentiability of marginal functions in nonlinear programming. Journal of optimization theory and applications (1996).

[8] Amos, Brandon, and J. Zico Kolter. Optnet: Differentiable optimization as a layer in neural networks. ICML, 2017.

[9] Agrawal, A., Amos, B., Barratt, S., Boyd, S., Diamond, S., & Kolter, J. Z. (2019). Differentiable convex optimization layers. Advances in neural information processing systems.

We are very grateful to the reviewer for their effort in reviewing our submission and are encouraged by their positive assessment of the theory, presentation, and significance of our work.

Response to Weaknesses

1. Implementing Algorithm 2

In our paper, we implemented the simpler algorithm (Algorithm 3). Per the reviewer's suggestion, we additionally implemented Algorithm 2 and repeated our experiments --- we have added the results in the pdf with our top-level response. These experiments show that Algorithm 2 outperforms both Algorithm 3 and the baseline cvxpylayer in terms of convergence, as suggested by our theory (see Theorem 4.2 and Theorem 4.3). We appreciate the reviewer’s suggestion and will update the paper with these experiments.

2. Large-scale implementation

We acknowledge the importance of testing our algorithms with large-scale experiments. However, our focus has primarily been theoretical since prior work had far weaker theoretical guarantees for our choice of stationarity and oracle. As the reviewer notes, our experiments indeed mainly provide a proof of concept for our algorithms.

In the code, we also have results on bilevel problems with bilinear constraints (in figures/exp_bilinear); in follow-up work, we hope to expand our implementation to a wider range of problems and large-scale experiments.

Response to Questions

1. Hessian approximations

The key step in algorithms for bilevel programming is that of computing $\frac{d y^\*}{dx}$, which in turn is composed of the product of a matrix, an inverse Hessian, and vector --- this therefore is the primary computational bottleneck in these algorithms. As the reviewer suggests, there have been many approaches for Hessian approximations, including in the context of bilevel programming.

For instance, [1] uses conjugate gradient to compute the Hessian-inverse-vector product, [2] uses Neumann approximation (essentially a geometric series for approximating matrix inverse) to approximate the inverse Hessian, and [3] uses a method inspired by Gauss-Newton to approximate the Hessian as the outer product of the corresponding gradients. However, all these works deal with unconstrained bilevel programming. In constrained settings (e.g., the ones we consider), the Hessian becomes significantly more complicated, and it is unclear if these techniques would apply to them (but would, nevertheless, be an interesting question to consider).

Another line of Hessian-free approaches by [4-6] uses value-function reformulation to circumvent the use of Hessians. Finally, we also mention the approach of [7], which approximates the Hessian-inverse-vector product by a linear system solver in conjunction with a finite difference. As we detail in our response to reviewer wH2u, the runtime of our algorithm (when applied in the unconstrained setting) improves upon that of this algorithm by [7]. It would be interesting to extend these approaches to the constrained setting we consider with our choice of stationarity and oracle access.

Conclusion

We again thank the reviewer for their time and effort in raising thoughtful questions. Please let us know if there are any further questions that we could help clarify!

References

[1] Pedregosa, F. (2016). Hyperparameter optimization with approximate gradient. In International conference on machine learning. PMLR.

[2] Lorraine, J., Vicol, P., & Duvenaud, D. (2020). Optimizing millions of hyperparameters by implicit differentiation. In International conference on artificial intelligence and statistics. PMLR.

[3] Giovannelli, T., Kent, G., & Vicente, L. N. (2021). Bilevel stochastic methods for optimization and machine learning: Bilevel stochastic descent and darts. arXiv preprint arXiv:2110.00604.

[4] Liu, B., Ye, M., Wright, S., Stone, P., & Liu, Q. (2022). Bome! bilevel optimization made easy: A simple first-order approach. Advances in neural information processing systems.

[5] Kwon, J., Kwon, D., Wright, S., & Nowak, R. D. (2023). A fully first-order method for stochastic bilevel optimization. In International Conference on Machine Learning. PMLR.

[6] Yao, W., Yu, C., Zeng, S., & Zhang, J. Constrained Bi-Level Optimization: Proximal Lagrangian Value Function Approach and Hessian-free Algorithm. In The Twelfth International Conference on Learning Representations.

[7] Yang, Y., Xiao, P., & Ji, K. (2023). Achieving O (ε-1.5) complexity in hessian/jacobian-free stochastic bilevel optimization. In Proceedings of the 37th International Conference on Neural Information Processing Systems.

We are very grateful to the reviewer for their effort in reviewing our submission and are encouraged by their positive assessment of our theory and presentation.

Response to Questions

1. Intuition behind clipping

Intuitively, the clipping ensures that consecutive iterates of the algorithm reside within a $\delta$-ball of each other — this in turn allows for a guarantee in terms of $(\delta, \epsilon)$-Goldstein stationarity (which, as we recall, is defined in terms of gradients in a $\delta$-ball). We remark that following this intuition, all previous papers on Goldstein stationarity, as discussed in the paper, also utilize normalized and/or clipped steps (a notable exception is [1], where the clipping was dropped by relaxing the stationarity notion accordingly).

2. Implementation

We acknowledge that the implementation and testing of our algorithm for the setting in question is an extremely important task. That said, our focus in this section was primarily theoretical since before our result, there did not exist any work at all with finite-time theoretical guarantees using a first-order oracle for the problem in question. We hope to implement our algorithm and perform large-scale experiments in follow-up works.

3. Comparison with Yang-Xiao-Ji [2]

The core idea in both our paper and that by Yang-Xiao-Ji (and indeed, by all papers on bilevel optimization) is an efficient approximation of $\frac{dy^\ast}{dx}^{\top} \nabla_y f(x, y^\ast)$. Recall that this term simplifies to $\frac{dy^\ast}{dx}^{\top} \nabla_y f(x, y^\ast)=-\nabla_{xy}^{2}g(x,y^{\ast})\cdot\nabla_{yy}^{2}g(x,y^{\ast})^{-1}\nabla_{y}f(x,y^{\ast}).$

The way the paper of Yang-Xiao-Ji handles this term is by separately approximating parts of this product: in particular, it approximates $\nabla_{xy}^{2}g(x,y^{\ast})$ via a finite-difference method and the term $\nabla_{yy}^{2}g(x,y^{\ast})^{-1}\nabla_{y}f(x,y^{\ast})$ by a linear system solver.

In contrast to this method, our paper approximates the entire term via a novel application of the implicit function theorem. We now elaborate on this point by illustrating our method for the unconstrained setting.

As stated in our Equation (5.2), our approximation (for this unconstrained setting) of the term $\frac{dy^\ast}{dx}^{\top} \nabla_y f(x, y^\ast)$ is $v_{x}:=\frac{\nabla_{x}g(x,y_{\delta}^{\ast})-\nabla_{x}g(x,y^{\ast})}{\delta}$. To see the validity of this approximation, we first note that $\lim_{\delta\rightarrow0}\frac{\nabla_{x}g(x,y_{\delta}^{\ast})-\nabla_{x}g(x,y^{\ast})}{\delta}=\nabla_{xy}^{2}g(x,y^{\ast})\cdot\frac{dy_{\delta}^{\ast}}{d\delta}\vert_{\delta=0}$. Next, to approximate $\frac{dy_{\delta}^{\ast}}{d\delta}\vert_{\delta=0},$ we observe that $y_{\delta}^{\ast}:=\text{argmin} (g(x,y)+\delta f(x,y)).$ Applying first-order optimality of $y_{\delta}^{\ast}$ and taking the derivative with respect to $\delta$ yields $\frac{dy_{\delta}^{\ast}}{d\delta}\vert_{\delta=0}=-\left(\nabla_{yy}^{2}g(x,y^{\ast})\right)^{-1}\nabla_{y}f(x,y^{\ast}).$ Combining this with the first step proves the claimed approximation.

Finally, when measured according to the stationarity criterion of Yang-Xiao-Ji (in the unconstrained setting), our work's oracle complexity is $O(\epsilon^{-1})$ (improving upon their result of $O(\epsilon^{-1.5})$).

Conclusion

We again thank the reviewer for their time and effort. Please let us know if there are any further questions that we could help clarify!

References

[1] Zhang, Qinzi, and Ashok Cutkosky. "Random scaling and momentum for non-smooth non-convex optimization." arXiv preprint arXiv:2405.09742 (2024).

[2] Yang, Y., Xiao, P., & Ji, K. (2023). Achieving O (ε-1.5) complexity in hessian/jacobian-free stochastic bilevel optimization. In Proceedings of the 37th International Conference on Neural Information Processing Systems.