First-Order Methods for Linearly Constrained Bilevel Optimization

Authors: guy  kornowski
swati  padmanabhan
kai  wang
jimmy  zhang
suvrit  sra
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 08:29:53

Abstract: Algorithms for bilevel optimization often encounter Hessian computations, which are prohibitive in high dimensions. While recent works offer first-order methods for unconstrained bilevel problems, the constrained setting remains relatively underexplored. We present  first-order linearly constrained optimization methods with finite-time hypergradient stationarity guarantees. For linear equality constraints, we attain $\epsilon$-stationarity in $\widetilde{O}(\epsilon^{-2})$ gradient oracle calls, which is nearly-optimal. 
For linear inequality constraints, we attain $(\delta,\epsilon)$-Goldstein stationarity in $\widetilde{O}(d{\delta^{-1} \epsilon^{-3}})$ gradient oracle calls, where $d$ is the upper-level dimension. 
Finally, we obtain for the linear inequality setting dimension-free rates of $\widetilde{O}({\delta^{-1} \epsilon^{-4}})$ oracle complexity under the additional assumption of oracle access to the optimal dual variable. Along the way, we develop new nonsmooth nonconvex optimization methods with inexact oracles. Our numerical experiments verify these guarantees.

Main paper: 
 introduction: 
 Bilevel optimization [1] [2] [3] [4] , an important problem in optimization, is defined as follows: minimize x∈X F (x) := f (x, y * (x)) subject to y * (x) ∈ arg min y∈S(x) g(x, y).(1.1)Here, the value of the upper-level problem at any point x depends on the solution of the lowerlevel problem. This framework has recently found numerous applications in meta-learning [5] [6] [7] [8] , hyperparameter optimization [9] [10] [11] , and reinforcement learning [12] [13] [14] [15] . Its growing importance has spurred increasing efforts towards designing computationally efficient algorithms for it.As demonstrated by [16] , a key computational step in algorithms for bilevel optimization is estimating dy * (x)/dx, the gradient of the lower-level solution. This gradient estimation problem has been extensively studied in differentiable optimization [17, 18] by applying the implicit function theorem to the KKT system of the given problem [19] [20] [21] [22] [23] [24] . However, this technique typically entails computing (or estimating) second-order derivatives, which can be prohibitive in high dimensions [25] [26] [27] .Besides being natural from the viewpoint of complexity theory, this question is well-grounded in applications such as mechanism design [29, 30] , resource allocation [31] [32] [33] [34] , and decision-making under uncertainty [20, 35, 36] . Our primary contribution is an affirmative answer to the highlighted question for bilevel programs with linear constraints, an important problem class often arising in adversarial training, decentralized meta learning, and sensor networks (see [37] ). While there have been some other recent works [37] [38] [39] on this problem, our work is first-order (as opposed to [37] ) and offers, in our view, a stronger guarantee on stationarity (compared to [38, 39] )-cf. Section 1.2.We provide first-order algorithms (with associated finite-time convergence guarantees) for linearly constrained bilevel programs (Problem 1.1). By "first-order", we mean that we use only zeroth and first-order oracle access to f and g. Our assumptions for each of our contributions are in Section 2.1.(1) Linear equality constraints. As our first contribution, we design first-order algorithms for solving Problem 1.1 where the lower-level constraint set S(x) := {y : Ax -By -b = 0} comprises linear equality constraints, and X a convex compact set. With appropriate regularity assumptions on f and g, we show in this case smoothness in x of the hyperobjective F . Inspired by ideas from Kwon et al. [40] , we use implicit differentiation of the KKT matrix of a slightly perturbed version of the lower-level problem to design a first-order approximation to ∇F . Constructing our first-order approximation entails solving a strongly convex optimization problem on affine constraints, which can be done efficiently. With this inexact gradient oracle in hand, we then run projected gradient descent, which converges in O(ϵ -2 ) iterations for smooth functions. Theorem 1.1 (Informal; cf. Theorem 3.1). Given Problem 1.1 with linear equality constraints S(x) := {y : Ax -By -b = 0} and X a convex compact set, under regularity assumptions on f and g (Assumptions 2.2 and 2.3), there exists an algorithm, which in O(ϵ -2 ) oracle calls to f and g, converges to an ϵ-stationary point of F .For linear equality constrained bilevel optimization, this is the first first-order result attaining ϵstationarity of F with assumptions solely on the constituent functions f and g and none on F -cf. Section 1.2 for a discussion of the results of Khanduri et al. [37] for this setting.(2) Linear inequality constraints. Next, we provide first-order algorithms for solving Problem 1.1 where the lower-level constraint set S(x) := {y : Ax -By -b ≤ 0} comprises linear inequality constraints, and the upper-level variable is unconstrained.Our measure of convergence of algorithms in this case is that of (δ, ϵ)-stationarity [41] : for a Lipschitz function, we say that a point x is (δ, ϵ)-stationary if within a δ-ball around x there exists a convex combination of subgradients of the function with norm at most ϵ (cf. Definition 2.1).To motivate this notion of convergence, we note that the hyperobjective F (in Problem 1.1) as a function of x could be nonsmooth and nonconvex (and Lipschitz, as we later prove). Minimizing such a function in general is known to be intractable [42] , necessitating local notions of stationarity. Indeed, not only is it impossible to attain ϵ-stationarity in finite time [43] , even getting near an approximate stationary point of an arbitrary Lipschitz function is impossible unless the number of queries is exponential in the dimension [44] . Consequently, for this function class, (δ, ϵ)-stationarity has recently emerged [43] to be a natural and algorithmically tractable notion of stationarity. We give the following guarantee under regularity assumptions on f and g. Theorem 1.2 (Informal; Theorem 4.1). Consider Problem 1.1 with linear inequality constraints S(x) := {y : Ax -By -b ≤ 0}. Under mild assumptions on f and g (Assumption 2.2) and the lower-level primal solution y * (Assumption 2.4), there exists an algorithm, which converges to a (δ, ϵ)-stationary point of F in O(dδ -1 ϵ -3 ) oracle calls to f and g, where d is the upper-level variable dimension.To the best of our knowledge, this is the first result to offer a first-order finite-time stationarity guarantee on the hyperobjective for linear inequality constrained bilevel optimization (cf. Section 1.2 for a discussion of related work [37] [38] [39] ). We obtain our guarantee in Theorem 1.2 by first invoking a result by Zhang and Lan [45] to obtain inexact hyperobjective values of F using only O(1) oracle calls to f and g. We also show (Lemma 4.3) that this hyperobjective F is Lipschitz. We then employ our inexact zeroth-order oracle for F in Algorithm 2 designed to minimize Lipschitz nonsmooth nonconvex functions (in particular, F ), with the following convergence guarantee. While such algorithms using exact zeroth-order access already exist [46] , extending them to the inexact gradient setting is non-trivial; we leverage recent ideas connecting online learning to nonsmooth nonconvex optimization by Cutkosky, Mehta, and Orabona [47] (cf. Section 4).(3) Linear inequality under assumptions on dual variable access. For the inequality setting (i.e., Problem 1.1 with the lower-level constraint set S(x) := {y : Ax -By -b ≤ 0}), we obtain dimension-free rates under an additional assumption (Assumption 2.5) on oracle access to the optimal dual variable λ * of the lower-level problem. We are not aware of a method to obtain this dual variable in a first-order fashion (though in practice, highly accurate approximations to λ * are readily available), hence the need for imposing this assumption. We believe that removing this assumption and obtaining dimension-free first-order rates in this setting would be an important direction for future work. Our guarantee for this setting is summarized below. Theorem 1.4 (Informal; Theorem 4.4 combined with Theorem 5.3). Consider Problem 1.1 with linear inequality constraints S(x) := {y : Ax -By -b ≤ 0} and unconstrained upper-level variable. Under mild regularity assumptions on f and g (Assumption 2.2), on y * (Assumption 2.4), and assuming oracle access to the optimal dual variable λ * (Assumption 2.5), there exists an algorithm, which in O(δ -1 ϵ -4 ) oracle calls to f and g converges to a (δ, ϵ)-stationary point for F .We obtain this result by first reformulating Problem 4.1 via the penalty method and constructing an inexact gradient oracle for the hyperobjective F (cf. Section 5). We then employ this inexact gradient oracle within an algorithm (Algorithm 3) designed to minimize Lipschitz nonsmooth nonconvex functions (in particular, F ), with the following convergence guarantee.Theorem 1.5 (Informal; Theorem 4.4). Given Lipschitz F : R d → R and ∥ ∇F (•) -∇F (•)∥ ≤ ϵ, there exists an algorithm that, in T = O(δ -1 ϵ -3 ) calls to ∇F , outputs a (δ, 2ϵ)-stationary point of F .Our Algorithm 3 is essentially a "first-order" version of Algorithm 2. Similar to Algorithm 2, despite the existence of algorithms with these guarantees with access to exact gradients [48] , their extensions to the inexact gradient setting are not trivial and also make use of the new framework of Cutkosky, Mehta, and Orabona [47] . We believe our analysis for this general task can be of independent interest to the broader optimization community. Lastly, we also use a more implementation-friendly variant of Algorithm 3 (with slightly worse theoretical guarantees) in numerical experiments.The vast body of work on asymptotic results for bilevel programming, starting with classical works such as Anandalingam and White [49] , Ishizuka and Aiyoshi [50] , White and Anandalingam [51] , Vicente, Savard, and Júdice [52] , Zhu [53] , and Ye and Zhu [54] , typically fall into two categories: those based on approximate implicit differentiation: Amos and Kolter [17] , Agrawal et al. [18] , Domke [55] , Pedregosa [56] , Gould et al. [57] , Liao et al. [58] , Grazzi et al. [59] , and Lorraine, Vicol, and Duvenaud [60] and those via iterative differentiation: Franceschi et al. [9] , Shaban et al. [10] , Domke [55] , Grazzi et al. [59] , Maclaurin, Duvenaud, and Adams [61] , and Franceschi et al. [62] .Another recent line of work in this category includes Khanduri et al. [37] , Liu et al. [63] , Ye et al. [64] , and Gao et al. [65] , which use various smoothing techniques.The first non-asymptotic result for bilevel programming was provided by Ghadimi and Wang [16] , which was followed by a flurry of work: for example, algorithms that are single-loop stochastic: Chen, Sun, and Yin [66] , Chen et al. [67] , and Hong et al. [68] , projection-free: Akhtar et al. [69] , Jiang et al. [70] , Abolfazli et al. [71] , and Cao et al. [72] , use variance-reduction and momentum: Khanduri et al. [73] , Guo et al. [74] , Yang, Ji, and Liang [75] , and Dagréou et al. [76] , those for single-variable bilevel programs: Jiang et al. [70] , Sabach and Shtern [77] , Amini and Yousefian [78, 79] , and Merchav and Sabach [80] , and for bilevel programs with special constraints: Khanduri et al. [37] , Abolfazli et al. [71] , Tsaknakis, Khanduri, and Hong [81] , and Xu and Zhu [82] .The most direct predecessors of our work are those by Khanduri et al. [37] , Yao et al. [38] , Lu and Mei [39] , Kwon et al. [40] , and Liu et al. [63] . As alluded to earlier, Liu et al. [28] recently made a significant contribution by providing for bilevel programming a fully first-order algorithm with finite-time stationarity guarantees. This was extended to the stochastic setting by Kwon et al. [40] (which we build upon), simplified and improved by Chen, Ma, and Zhang [83] , and extended to the constrained setting by Khanduri et al. [37] , Yao et al. [38] , and Lu and Mei [39] .The works of Yao et al. [38] and Lu and Mei [39] study the more general problem of bilevel programming with general convex constraints. However, they use KKT stationarity as a proxy to the hypergradient stationarity. Our Theorem 1.4 is restricted to linear inequality constraints, we provide stationarity guarantees directly in terms of the objective of interest. Moreover, Yao et al. [38] assumes joint convexity of the lower-level constraints in upper and lower variables to allow for efficient projections, while we require convexity only in the lower-level variable.The current best result for the linearly constrained setting is that of Khanduri et al. [37] . However, this work requires Hessian computations (and is therefore not fully first-order). Moreover, Khanduri et al. [37] imposes strong regularity assumptions on the hyperobjective F , which are, in general, impossible to verify. In contrast, Theorem 1.1 imposes assumptions solely on the constituent functions f and g, none directly on F , thus making substantial progress on these two fronts. 
preliminaries: 
 We follow standard notation (see Appendix A), with only the following crucial definition stated here. Definition 2.1. Consider a locally Lipschitz function f : R d → R, a point x ∈ R d , and a parameter δ > 0. The Goldstein subdifferential [41] off at x is the set ∂ δ f (x) := conv(∪ y∈B δ (x) ∂f (y)), where ∂f (x) = conv {lim n→∞ ∇f (x n ) : x n → x, x n ∈ dom(∇f )} is the Clarke subdifferential [84] of f and B δ (x) denotes the Euclidean ball of radius δ around x. A point x is called (δ, ϵ)-stationary if dist(0, ∂ δ f (x)) ≤ ϵ, where dist(x, S) := inf y∈S ∥x -y∥.We consider Problem (i) Upper-level: The objective f is C f -smooth and L f -Lipschitz continuous in (x, y).(ii) Lower-level: The objective g is C g -smooth. Fixing any x ∈ X , g(x, •) is µ g -strongly convex.(iii) We assume that the linear independence constraint qualification (LICQ) condition holds for the LL problem at every x and y, i.e., the constraint h(x, y) := Ax -By -b has a full row rank B. Assumption 2.3. For Problem 3.1 (with linear equality constraints), we additionally assume that the set X is convex and compact, and that the objective g is S g -Hessian smooth, that is, ∇ 2 g(x, y) -∇ 2 g(x, ȳ) ≤ S g ∥(x, y) -(x, ȳ)∥ ∀x, x ∈ X , and y, ȳ ∈ R dy . Assumption 2.4. For Problem 4.1 (with linear inequality constraints), we additionally assume that y * is L y -Lipschitz in x, where y * is the LL primal solution y * (x), λ * (x) = arg max y min β≥0 g(x, y) + β ⊤ h(x, y), where h(x, y) := Ax -By -b. Assumption 2.5. We provide additional results for Problem 4.1 under additional stronger assumptions stated here: Denote the LL primal and dual solution y * (x), λ * (x) = arg max y min β≥0 g(x, y) + β ⊤ h(x, y), where h(x, y) := Ax -By -b; then, we assume exact access to λ * and that ∥λ * (x)∥ ≤ R. Assumptions 2.2(i) and 2.2(ii) are standard in bilevel optimization. Assumption 2.2(iii) is the same as the complete recourse assumption in stochastic programming [85] , that is, the LL problem is feasible y for every x ∈ R dx . Assumption 2.3 is used only in the equality case and guarantees smoothness of F . Assumption 2.4 is used in the inequality case and implies Lipschitzness of F . We need the stronger assumption in Assumption 2.5 for our dimension-free result for the linear inequality case.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: This is discussed in Section 7Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
lower-level problem with linear equality constraint: 
 We first obtain improved results for the setting of bilevel programs with linear equality constraints in the lower-level problem. Our formal problem statement is: minimize x∈X F (x) := f (x, y * (x)) subject to y * (x) ∈ arg min y:h(x,y)=0 g(x, y),where f , g, h(x, y) := Ax -By -b, and X satisfy Assumptions 2.2 and 2.3. The previous best result on Problem 3.1 providing finite-time ϵ-stationarity guarantees, by Khanduri et al. [37] , required certain regularity assumptions on F as well as Hessian computations. In contrast, our finite-time guarantees require assumptions only on f and g, not on F ; indeed, in our work, these desirable properties of F are naturally implied by our analysis. Specifically, our key insight is that the hypergradient∇F (x) := ∇ x f (x, y * ) + dy * (x) dx ⊤ ∇ y f (x, y * ) for Problem 3.1 is Lipschitz-continuous and admits an easily computable -yet highly accurate -finite-difference approximation. Therefore, O(ϵ -2 ) iterations of gradient descent on F with this finite-difference gradient proxy yield an ϵ-stationary point.Specifically, for any fixed x ∈ X , our proposed finite-difference gradient proxy approximating the non-trivial-to-compute component dy * (x) dx ⊤ ∇ y f (x, y * ) of the hypergradient is given byEQUATIONwhere (y Output: ∇F = v x + ∇ x f (x, y * )Notice that the finite-difference term in (3.2) avoids differentiating through the implicit function y * . Instead, all we need to evaluate it are the values of (y * , λ * , y * δ , λ * δ ) (and gradients of g and h). Since (y * , λ * ) are solutions to a smooth strongly convex linearly constrained problem, they can be approximated at a linear rate. Similarly, since the approximation error in (3.2) is proportional to δ (cf. Lemma 3.2), a small enough δ in the perturbed objective g + δf in (3.3) ensures that it is dominated by the strongly convex and smooth g, whereby accurate approximates to (y * δ , λ * δ ) can also be readily obtained. Putting it all together, the proposed finite-difference hypergradient proxy in (3.2) is efficiently computable, yielding the following guarantee. Theorem 3.1. Consider Problem 3.1 under Assumption 2.2, and let κ = C g /µ g be the condition number of g. Then Algorithm 5 finds an ϵ-stationary point (in terms of gradient mapping, see (B.14)) after T = O(C F (F (x 0 ) -inf F ) √ κϵ -2) oracle calls to f and g, whereC F := 2(L f + C f + C g )C 3H S g (L g + ∥A∥) 2 is the smoothness constant of the hyperobjective F .We now sketch the proof of Theorem 3.1. The complete proofs may be found in Appendix B.We briefly outline the two key technical building blocks alluded to above, that together give us Theorem 3.1: the approximation guarantee of our finite-difference gradient proxy ((3.2)) and the smoothness of hyperobjective F (for Problem 3.1). The starting point for both these results is the following simple observation obtained by implicitly differentiating, with respect to x, the KKT system associated with y * = arg min h(x,y)=0 g(x, y) and optimal dual variable λ * :dy * (x) dx dλ * (x) dx = ∇ 2 yy g(x, y * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 -1 -∇ 2 yx g(x, y * ) -∇ x h(x, y * ) (3.4)The invertibility of the matrix in the preceding equation is proved in Corollary B.3: essentially, this invertibility is implied by strong convexity of g and ∇ y h(x, y * ) = B having full row rank. This in conjunction with the compactness of X implies that the inverse of the matrix is bounded by some constant C H (cf. Corollary B.3 for details). Our hypergradient approximation guarantee follows: Lemma 3.2. For Problem 3.1 under Assumption 2.2, with v x as in (3.2), the following holds:v x -dy * (x) dx ⊤ ∇ y f (x, y * ) ≤ O(C F δ).Proof sketch; see Appendix B. The main idea is that the two terms being compared are essentially the same by the implicit function theorem. First, we use the expression for dy * (x) dx from (3.4):dy * (x) dx ⊤ ∇ y f (x, y * ) = -∇ 2 yx g(x, y * ) -∇ x h(x, y * ) ⊤ ∇ 2 yy g(x, y * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 -1 ∇ y f (x, y * ) 0 .We now examine v x .For simplicity of exposition, we instead considerlim δ→0 ∇x[g(x,y * δ )+⟨λ * δ ,h(x,y * )⟩]-∇x[g(x,y * )+⟨λ * ,h(x,y * )⟩] δ, which, by the fundamental theorem of calculus and Assumption 2.3 , equals v x up to an O(δ)-additive error. Note that this expression is:EQUATIONSince y * δ is minimizes a strongly convex function over a linear equality constraint (3.3), the reasoning that yields (3.4) also gives the following, which, when combined with (3.5), finishes the proof: EQUATION:= O{(L f + C f + C g )C 3 H S g (L g + ∥A∥) 2 }.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: The assumptions, theorem statements, and proof sketches are included in the main paper. The full proofs are included in the appendix.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
nonsmooth nonconvex optimization with inexact oracle: 
 We now shift gears from the case of linear equality constraints to that of linear inequality constraints. Specifically, defining h(x, y) = Ax -By -b, the problem we now consider is minimize x F (x) := f (x, y * (x)) subject to y * (x) ∈ arg min y:h(x,y)≤0 g(x, y).As noted earlier, for this larger problem class, the hyperobjective F can be nonsmooth nonconvex, necessitating our measure of convergence to be the now popular notion of Goldstein stationarity [43] .Our first algorithm for solving Problem 4.1 is presented in Algorithm 2, with its convergence guarantee in Theorem 4.1. At a high level, this algorithm first obtains access to an inexact zeroth-order oracle to F (we shortly explain how this is done) and uses this oracle to construct a (biased) gradient estimate of F . It then uses this gradient estimate to update the iterates with a rule motivated by recent works reducing nonconvex optimization to online optimization [47] . We explain this in Section 4. = min δ 2 , F (x0)-inf F L f Ly , ν = δ-ρ, D = Θ νϵ 2 ρ 2 dxρ 2 L 2 f L 2 y +α 2 d 2 x, and η = Θνϵ 3 ρ 4 (dxρ 2 L 2 f L 2 y +α 2 d 2 x ) 2outputs x out such that E[dist(0, ∂ δ F (x out ))] ≤ ϵ + α with T oracle calls to f and g, where:T = O √ κd x (F (x 0 ) -inf F ) δϵ 3 • L 2 f L 2 y + α 2 d x δ 2 + d x L 2 f L 2 y (F (x 0 ) -inf F ) 2 • log(L f /α) .Algorithm 2 Nonsmooth Nonconvex Algorithm with Inexact Zero-Order oracle 1: Input: Initialization x 0 ∈ R d , clipping parameter D > 0, step size η > 0, smoothing parameter ρ > 0, accuracy parameter ν > 0, iteration budget T ∈ N, inexact zero-order oracle F : R d → R. 2: Initialize: ∆ 1 = 0 3: for t = 1, . . . , T do 4: Sample s t ∼ Unif[0, 1], w t ∼ Unif(S d-1 ) 5:x t = x t-1 + ∆ t , z t = x t-1 + s t ∆ t 6: g t = d 2ρ ( F (z t + ρw t ) -F (z t -ρw t ))w t 7: ∆ t+1 = clip D (∆ t -η g t ) ▷ clip D (z) := min{1, D ∥z∥ } • z 8: M = ⌊ ν D ⌋, K = ⌊ T M ⌋ 9: for k = 1, . . . , K do 10: x k = 1 M M m=1 z (k-1)M +m 11: Sample x out ∼ Unif{x 1 , . . . , x K } 12: Output: x out .Algorithm 2 is a variant of gradient descent with momentum and clipping, with g t the inexact gradient, ∆ t a clipped accumulated gradient (hence accounts for past gradients, which serve as a momentum), and the clipping ensuring that consecutive iterates of the algorithm reside within a δ-ball of each other. While similar algorithms have appeared in prior work on nonsmooth nonconvex optimization (e.g. [47] ), none of them account for inexactness in the gradient, crucial in our setting.We can obtain inexact zeroth-order oracle access to F because (as formalized in Lemma 4.2) despite potential nonsmoothness and nonconvexity of F in Problem 4.1, estimating its value F (x) at any point x amounts to solving a single smooth and strongly convex optimization problem, which can be done can be done in O(1) oracle calls to f and g by appealing to a result by Zhang and Lan [45] . Having computed an inexact value of the hyperobjective F , we now show how to use it to develop an algorithm for solving Problem 4.1. To this end, we first note that F , despite being possibly nonsmooth and nonconvex, is Lipschitz and hence amenable to the use of recent algorithmic developments in nonsmooth nonconvex optimization pertaining to Goldstein stationarity. Lemma 4.3. Under Assumption 2.2 and 2.5,F in Problem 4.1 is O(L f L y )-Lipschitz in x.With this guarantee on the Lipschitzness of F , we prove Theorem C.1 for attaining Goldstein stationarity using the inexact zeroth-order oracle of a Lipschitz function. Our proof of Theorem C.1 crucially uses the recent online-to-nonconvex framework of Cutkosky, Mehta, and Orabona [47] . Combining Lemma 4.2 and Theorem C.1 then immediately implies Theorem 4.1.In Section 5, we provide a way to generate approximate gradients of F . Here, we present an algorithm that attains Goldstein stationarity of Problem 4.1 using this inexact gradient oracle. While there has been a long line of recent work on algorithms for nonsmooth nonconvex optimization with convergence to Goldstein stationarity [43, 48, [86] [87] [88] , these results necessarily require exact gradients. This brittleness to any error in gradients renders them ineffective in our setting, where our computed (hyper)gradient necessarily suffers from an additive error. While inexact oracles are known to be effective for smooth or convex objectives [89] , utilizing inexact gradients in the nonsmooth nonconvex regime presents a nontrivial challenge. Indeed, without any local bound on gradient variation due to smoothness, or convexity that ensures that gradients are everywhere correlated with the direction pointing at the minimum, it is not clear a priori how to control the accumulating price of inexactness throughout the run of an algorithm. To derive such results, we use the recently proposed connection between online learning and nonsmooth nonconvex optimization by Cutkosky, Mehta, and Orabona [47] . By controlling the accumulated error suffered by online gradient descent for linear losses (cf. Lemma C.3), we derive guarantees for our setting of interest, providing Lipschitz optimization algorithms that converge to Goldstein stationary points even with inexact gradients. This algorithm matches the best known complexity in first-order nonsmooth nonconvex optimization [43, 47, 48] , merely replacing the convergence to a (δ, ϵ)-stationary point by (δ, ϵ + α)-stationarity, where α is the inexactness of the gradient oracle.Algorithm 3 Nonsmooth Nonconvex Algorithm with Inexact Gradient Oracle 1: Input: Initialization x 0 ∈ R d , clipping parameter D > 0, step size η > 0, accuracy parameter δ > 0, iteration budget T ∈ N, inexact gradient oracle ∇F : R d → R d . 2: Initialize: ∆ 1 = 0 3: for t = 1, . . . , T do 4: Sample s t ∼ Unif[0, 1] 5: x t = x t-1 + ∆ t , z t = x t-1 + s t ∆ t 6: g t = ∇F (z t ) 7: ∆ t+1 = clip D (∆ t -η g t ) ▷ clip D (z) := min{1, D ∥z∥ } • z 8: M = ⌊ δ D ⌋, K = ⌊ T M ⌋ 9: for k = 1, . . . , K do 10: x k = 1 M M m=1 z (k-1)M +m 11: Sample x out ∼ Unif{x 1 , . . . , x K } 12: Output: x out . Theorem 4.4. Suppose F : R d → R is L-Lipschitz and that ∥ ∇F (•) -∇F (•)∥ ≤ α. Then running Algorithm 3 with D = Θ( δϵ 2 L 2 ), η = Θ( δϵ 3 L 4 ), outputs a point x out such that E[dist(0, ∂ δ F (x out ))] ≤ ϵ + α, with T = O (F (x0)-inf F )L 2 δϵ 3 calls to ∇F (•).We defer the proof of Theorem 4.4 to Appendix C.1. Plugging the complexity of computing inexact gradients, as given by Theorem 5.3, into the result above, we immediately obtain convergence to a (δ, ϵ)-stationary point of Problem 1.1 with O(δ -1 ϵ -4 ) gradient calls overall.Implementation-friendly algorithm. While Algorithm 3 matches the best-known results in nonsmooth nonconvex optimization, it could be impractical due to several hyperparameters which need tuning. Arguably, a more natural application of the hypergradient estimates would be simply plugging them into gradient descent, which requires tuning only the stepsize. Since F is neither smooth nor convex, perturbations are required to guarantee differentiability along the trajectory. We therefore complement Theorem 4.4 by analyzing perturbed (inexact) gradient descent in the nonsmooth nonconvex setting (Algorithm 7) and state its theoretical guarantee in Theorem C.4. Despite its suboptimal worst-case theoretical guarantees, we find this algorithm easier to implement in practice.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: We provide our full code in the supplemental material, and it can be used to reproduce the experimental results.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
inequality constraints: constructing the inexact gradient oracle: 
 Computing a stationary point of F of Problem 4.1 via any first-order method would require:EQUATIONfor which the key challenge lies in computing dy * (x)/dx. This requires differentiating through an argmin operator, which typically requires second-order derivatives. Instead, here we differentiate (using the implicit function theorem) through the KKT conditions describing y * (x) and get:∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h ∇ y h ⊤ I diag(λ * I )∇ y h I 0 dy * (x) dx dλ * I (x) dx = - ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h diag(λ * I )∇ x h I (5.2)where given x, we assume efficient access to the optimal dual solution λ * (x) of the LL problem in Problem 4.1. In (5.2), we use Observe that as currently stated, (5.2) leads to a second-order computation of dy * (x)/dx. In the rest of the section, we provide a fully first-order approximate hypergradient oracle by constructing an equivalent reformulation of Problem 4.1 using a penalty function.I := {i ∈ [d h ] : h i (x, y) = 0, λ * i > 0}We begin by reformulating Problem 4.1 into a single level constrained optimization problem:EQUATIONwhere g * (x) := min y:h(x,y)≤0 g(x, y) = g(x, y * (x)) and λ * (x) is the optimal dual solution. The equivalence of this reformulation to Problem 4.1 is spelled out in Appendix D. From (5.3), we define the following penalty function, crucial to our analysis:EQUATIONwhere α = [α 1 , α 2 ] ≥ 0 are the penalty parameters. Notably, we can compute its derivative with respect to x of (5.4) in a fully first-order fashion by the following expression:∇ x L λ * ,α (x, y) = ∇ x f (x, y)+α 1 (∇ x g(x, y)+∇ x h(x, y) ⊤ λ * -∇ x g * (x))+α 2 ∇ x h I (x, y) ⊤ h I (x, y).To give some intuition for our choice of penalties in (5.4), we note that the two constraints in (5.3) behave quite differently. The first constraint g(x, y) + λ * (x) ⊤ h(x, y) ≤ g * (x) is one-sided, i.e., can only be violated or met, and hence just needs a penalty parameter α 1 to weight the "violation". As to the second constraint h(x, y) ≤ 0, it can be arbitrary. To allow for such a "two-sided" constraint, we penalize only the active constraints I, i.e., we use ∥h I (x, y)∥ 2 to penalize deviation.Next, we define the optimal solutions to the penalty function optimization by: y * λ * ,α (x) := arg min y L λ * ,α (x, y).(5.5)We now show that this minimizer is close to the optimal solution of the LL problem, while suffering only a small constraint violation. Lemma 5.1. Given any x, the corresponding dual solution λ * (x), primal solution y * (x) of the lower optimization problem in Problem 4.1, and y * λ * ,α (x) as in (5.5), satisfy:y * λ * ,α (x) -y * (x) ≤ O(α -1 1 ) and h I (x, y * λ * ,α (x)) ≤ O(α -1/2 1 α -1/2 2).(5.6)The proof of Lemma 5.1 is based on the Lipschitzness of f and strong convexity of g for sufficiently large α 1 . The aforementioned constraint violation bound on h I (x, y) is later used in Lemma 5.2 to bound the inexactness of our proposed gradient oracle.The main export of this section is the following bound on the approximation of the hypergradient. This, together with the bounds in Lemma 5.1, validate our use of the penalty function in (5.4). Lemma 5.2. Consider F as in Problem 4.1, L as in (5.4), a fixed x, and y * λ * ,α as in (5.5). Then under Assumptions 2.2 and 2.5, we have:EQUATION2 ).The proof can be found in Appendix G. With this hypergradient approximation guarantee, we design Algorithm 4 to compute an inexact gradient oracle for the hyperobjective F .Theorem 5.3. Given any accuracy parameter α > 0, Algorithm 4 outputs∇ x F (x) such that ∥ ∇F (x) -∇F (x)∥ ≤ α within O(α -1) gradient oracle evaluations.The full proof of this result may be found in Appendix H. Inner level dimension  
experiments: 
 We generate instances of the following constrained bilevel optimization problem:EQUATIONwhere h(x, y) = Ay -b is a d h -dim linear constraint. The PSD matrix Q ∈ R dy×dy , c ∈ R dy , P ∈ R dx×dy , and constraints A ∈ R d h ×dy , b ∈ R d h are randomly generated from normal distributions (cf. Appendix K). We compare Algorithm 3 with a non-fully first-order method using cvxpyLayer [18] . Both algorithms use Adam [90] to control the learning rate, and are averaged over 10 random seeds.Figure 2a shows that both the algorithms converge to the same optimal solution at the same rate. Simultaneously, the colorful bars represent the gradient differences between two methods, showing the inexactness of our gradients. Figure 2b additionally varies this inexactness to demonstrate its impact on convergence with standard deviation plotted. Figure 2c compares the computation costs for different lower-level problem sizes. Our fully first-order method significantly outperforms, in computation cost, the non-fully first-order method implemented using differentiable optimization method. The implementation can be found in https://github.com/guaguakai/constrained-bilevel-optimization. 
limitations and future directions: 
 One limitation to our approach is that the inexact gradient oracle we constructed in Section 5 requires access to the exact dual multiplier λ * . For a first-order algorithm, the closest proxy one could get to this would be a highly accurate approximation (which could be computed up to ϵ error in O(log(1/ϵ)) gradient oracle evaluations). Removing this "exact dual access" assumption (Assumption 2.5) would be an important result.Another important direction for future work would be to extend the hypergradient stationarity guarantee of Theorem 1.4 to bilevel programs with general convex constraints. To this end, we conjecture that the use of a primal-only gradient approximation oracle could be potentially effective.Finally, our current rate of O(δ -1 ϵ -4 ) oracle calls for reaching (δ, ϵ)-Goldstein stationarity is not necessarily inherent to the problem; indeed, it might be the case that an alternate approach could improve it to the best known rate of O(δ -1 ϵ -3 ) for nonsmooth nonconvex optimization.We use ⟨•, •⟩ to denote inner products and ∥ • ∥ for the Euclidean norm. Unless transposed, all vectors are column vectors. For f : R d2 → R d1 its Jacobian with respect to x ∈ R d2 is ∇f ∈ R d1×d2 . For f : R d → R, we overload ∇f to refer to its gradient (the transposed Jacobian), a column vector. We use ∇ x to denote partial derivatives with respect to x.A function f : R n → R m is L-Lipschitz if for any x, y, we have ∥f (x) -f (y)∥ ≤ L∥x -y∥. A differentiable function f : R n → R is convex if for any x, y ∈ R n we have f (y) ≥ f (x) + ∇f (x) ⊤ (y -x); it is µ-strongly convex if f -µ 2 ∥ • ∥ 2 is convex; it is β-smooth if ∇f is β-Lipschitz.For a Lipschitz function f , a point x is (δ, ϵ)-stationary if within a δ-ball around x, there exists a convex combination of subgradients of f with norm at most ϵ. For a differentiable function f , we say that x is ϵ-stationary if ∥∇f (x)∥ ≤ ϵ.In this section, we provide the full proofs of claims for bilevel programs with linear equality constraints, as stated in Section 3. We first state a few technical results using the implicit function theorem that we repeatedly invoke in our results for this setting. Lemma B.1. Fix a point x. Given y * = arg min y:h(x,y)=0 g(x, y) where g is strongly convex in y and λ * is the dual optimal variable for this problem, define L eq (x, y, λ) = g(x, y) + ⟨λ, h(x, y)⟩. Then, we have∇ 2 yy L eq (x, y * , λ * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0H for linear equality constraintsdy * dx dλ * dx = -∇ 2 yx g(x, y * ) -∇ 2 yx ⟨λ * , h(x, y * )⟩ -∇ x h(x, y * ) .Proof. Since g is strongly convex, by linear constraint qualification, the KKT condition is both sufficient and necessary condition for optimality. Hence, consider the following KKT system obtained via first order optimality of y * , with dual optimal variable λ * : ∇ y g(x, y * ) + ∇ y ⟨λ * , h(x, y * )⟩ = 0, and h(x, y * ) = 0. (B.1)Differentiating the system of equations in (B.1) with respect to x and rearranging terms in a matrixvector format yields: Proof. Let u = [y, λ]. We show that Hu = 0 implies u = 0, which in turn implies invertibility of H. If ∇ y h(x, y * )y ̸ = 0, then by construction of u and H, we must also have Hu ̸ = 0. Otherwise if ∇ y h(x, y * )y = 0 and y ̸ = 0, the quadratic form u ⊤ Hu is positive, as seen by∇ 2 yy g(x, y * ) + ∇ 2 yy ⟨λ * , h(x, y * )⟩ ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 dy * dx dλ * dx = -∇ 2 yx g(x, y * ) -∇ 2 yx ⟨λ * , h(x, y * )⟩ -∇ x h(x, y * ) (B.2) Noting that ∇ 2 yy L eq (x, y, λ) = ∇ 2 yy g(x, y) + ∇ 2 yy ⟨λ, h(x,u ⊤ Hu = y ⊤ ∇ 2yy L eq (x, y * , λ * )y > 0, where the final step is by the assumption of L eq being positive definite over the defined tangent plane T = {y : ∇ y h(x, y * )y = 0}. If y = 0 while Hu = 0, then ∇ y h having full rank implies λ = 0. Combined with y = 0, this means u = 0, as required when Hu = 0. This concludes the proof.Corollary B.3. For Problem 3.1 under Assumption 2.2 and Assumption 2.3, the matrix H (as defined in (3.4)) is non-singular. Further, there exists a finite C H such that ∥H -1 ∥ ≤ C H .Proof. Since we are assuming strong convexity of g, Lemma B.2 applies, yielding the claimed invertibility of H. Combined with the boundedness of variables x (per Assumption 2.3) and continuity of the inverse implies a bound on ∥H -1 ∥.We now show how to construct the inexact gradient oracle for the objective F in Problem 3.1. As sketched in Section 3, we then use this oracle in a projected gradient descent algorithm to get the claimed guarantee. Lemma B.4. Consider Problem 3.1 under Assumption 2.2 and Assumption 2.3. Let y * δ be as defined in (3.3) . Then, for any δ ∈ [0, ∆] with ∆ ≤ µ g /2C f , the following relation is valid:∥y * δ -y * ∥ ≤ M (x)δ, with M (x) := 2 µ g ∥∇ y f (x, y * )∥ ≤ 2L f µ g .Proof. The first-order optimality condition applied to g(x, y) + δf (x, y) at y * and y * δ gives ⟨∇ y g(x, y * δ ) + δ∇ y f (x, y * δ ), y * -y * δ ⟩ ≥ 0, which upon adding and subtracting ∇ y f (x, y * ) transforms into⟨∇ y g(x, y * δ ) + δ[∇ y f (x, y * δ ) -∇ y f (x, y * )] + δ∇ y f (x, y * ), y * -y * δ ⟩ ≥ 0. (B.3)Similarly, the first-order optimality condition applied to g at y * and y * δ gives ⟨∇ y g(x, y * ),y * δ -y * ⟩ ≥ 0. (B.4)Adding Inequality (B.3) and Inequality (B.4) and rearranging yields⟨∇ y g(x, y * δ ) -∇ y g(x, y * ) + δ[∇ y f (x, y * δ ) -∇ y f (x, y * )], y * δ -y * ⟩ ≤ ⟨δ∇ y f (x, y * ), y * -y * δ ⟩.Applying to the left side above a lower bound via strong convexity of g + δf and to the right hand side an upper bound via Cauchy-Schwarz inequality, we haves∥y * δ -y * ∥ ≤ δ∥∇ y f (x, y * )∥, (B.5)where s is the strong convexity of g + δf . Since f is C f -smooth, the worst case value of this is s = µ g -δC f = µ g - µg 2C f C f = µ g /2,lim δ→0 ∇ x [g(x, y * δ (x)) + λ * δ h(x, y * )] -∇ x [g(x, y * (x)) + λ * h(x, y * )] δ = dy * (x) dx ⊤ ∇ y f (x, y * (x)).Proof. Recall that by definition, g is strongly convex and y * = arg min y:h(x,y)=0 g(x, y). Hence, we can apply Lemma B.1. Combining this with Lemma B.2 and further applying that linearity of h implies ∇ 2 yy h = 0 and ∇ 2 xy h = 0, we obtain the following:dy * dx dλ * dx = ∇ 2 yy g(x, y * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 -1 -∇ 2 yx g(x, y * ) -∇ x h(x, y * ).So we can express the right-hand side of the claimed equation in the lemma statement bydy * (x) dx ⊤ ∇ y f (x, y * (x)) = dy * dx ⊤ dλ * dx ⊤ ∇ y f (x, y * (x)) 0 ,which can be further simplified to-∇ 2 yx g(x, y * ) ⊤ -∇ x h(x, y * ) ⊤ ∇ 2 yy g(x, y * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 -1 ∇ y f (x, y * (x)) 0 . (B.6)We now apply Lemma B.1 to the perturbed problem defined in (3.3) . We know from Lemma B.4 that lim δ→0 y * δ = y * . The associated KKT system is given by δ∇ y f (x, y * δ ) + ∇ y g(x, y * δ ) + ∇ y ⟨λ * δ , h(x, y * δ )⟩ = 0 and h(x, y * δ ) = 0. (B.7)Taking the derivative with respect of (B.7) gives the following implicit system, where we used the fact that h is linear and hence ∇ 2 yy h = 0:δ∇ 2 yy f (x, y * δ ) + ∇ 2 yy g(x, y * δ ) ∇ y h(x, y * δ ) ⊤ ∇ y h(x, y * δ ) 0 H δ dy * δ dδ dλ * δ dδ = -∇ y f (x, y * δ ) ⊤ 0 . (B.8)For a sufficiently small δ, we have ∇ | δ=0 = ∇ 2 yy g(x, y * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 -1 -∇ y f (x, y * ) 0 .So we can express the left-hand side of the expression in the lemma statement by lim δ→0 ∇ x [g(x, y * δ (x)) + ⟨λ * δ , h(x, y * )⟩] -∇ x [g(x, y * (x)) + ⟨λ * , h(x, y * )⟩] δ = ∇ 2 xy g(x, y * ) dy * δ dδ + ∇ x h(x, y * ) ⊤ dλ * δ dδ = ∇ 2 xy g(x, y * ) ∇ x h(x, y * ) ⊤ ∇ 2 yy g(x, y * ) ∇ y h(x, y * ) ⊤ ∇ y h(x, y * ) 0 -1 -∇ y f (x, y * ) 0 , which matches (B.6) (since (∇ 2 yx g) ⊤ = ∇ 2 xy g),F := O{(L f + C f + C g )C 3 H S g (L g + ∥A∥) 2 }.Proof. Rearranging (3.4) and applying Corollary B.3, we havedy * dx dλ * dx = ∇ 2 yy g(x, y * ) B ⊤ B 0 -1 -∇ 2 yx g(x, y * ) -∇ x h(x, y * ) .This implies a Lipschitz bound of C H • (C g + ∥A∥). Next, note that in the case with linear equality constraints, the terms in (B.2) involving second-order derivatives of h are all zero; differentiating (B.2) with respect to x, we notice that the linear system we get again has the same matrix H from before. We can therefore again perform the same inversion and apply the bound on ∥H -1 ∥ and on the third-order derivatives of g (Assumption 2.3) to observe that∥ d 2 y * dx 2 ∥ ≤ O(C H • S g ∥ dy * dx ∥ 2 ) = O(C 3 H • S g • (C g + ∥A∥) 2 ), where we are hiding numerical constants in the Big-Oh notation. As a result, we can calculate the Lipschitz smoothness constant associated with the hyper-objective F by∥∇F (x) -∇F (x)∥ ≤ ∥ dy * (x) dx ∇ y f (x, y * (x)) - dy * (x) dx ∇ y f (x, y * (x))∥ + ∥∇ x f (x, y * (x)) -∇ x f (x, y * (x))∥ ≤ [C f C H (L g + ∥A∥) + C f C 2 H (L g + ∥A∥) 2 + L f C 3 H S g (L g + ∥A∥) 2 ]∥x -x∥ + [C f + C f C H (L g + ∥A∥)]∥x -x∥ ≤ 2(L f + C f + C g )C 3 H S g (L g + ∥A∥) 2 C F ∥x -x∥.Lemma 3.2. For Problem 3.1 under Assumption 2.2, with v x as in (3.2), the following holds:v x -dy * (x) dx ⊤ ∇ y f (x, y * ) ≤ O(C F δ).Proof. For simplicity, we adopt the following notation throughout this proof: g xy (x, y) = ∇ 2 xy g, and g xyy denotes the tensor such that its ijk entry is given by ∂ 3 g ∂xi∂yj ∂y k . We first consider the terms involving g. By the fundamental theorem of calculus, we have∇ x g(x, y * δ (x)) -∇ x g(x, y * (x)) = δ t=0 g xy (x, y * t (x)) dy * t (x) dt dt.As a result, we have∇ x g(x, y * δ (x)) -∇ x g(x, y * (x)) δ -g xy (x, y * (x)) dy * t (x) dt | t=0 = 1 δ δ t=0 g xy (x, y * t (x)) dy * t (x) dt -g xy (x, y * (x)) dy * t (x) dt | t=0 dt = 1 δ δ t=0 g xy (x, y * t (x)) dy * t (x) dt -g xy (x, y * (x)) dy * t (x) dt | t=0 dt = 1 δ δ t=0 (g xy (x, y * t (x)) -g xy (x, y * (x))) dy * t (x) dt dt + 1 δ δ t=0 g xy (x, y * (x)) • dy * t (x) dt - dy * t (x) dt | t=0 dt. (B.9)We now bound each of the terms on the right-hand side of (B.9). For the first term, we have∥ 1 δ δ t=0(g xy (x, y * t (x)) -g xy (x, y * (x))dt)dy * t (x) dt ∥ ≤ 1 δ δ t=0 ∥ dy * t (x) dt ∥ • t s=0 ∥g xyy (x, y * s (x))∥∥ dy * s (x) ds ∥ds • dt ≤ 1 δ δ t=0 ∥ dy * t (x) dt ∥ • max s∈[0,δ] ∥g xyy (x, y * s (x))∥ • ∥ dy * s (x) ds ∥tdt ≤ 1 δ • max u∈[0,δ] ∥g xyy (x, y * u (x))∥ • δ 2 • max t∈[0,δ] ∥ dy * t (x) dt ∥ 2 ≤ δ • max u∈[0,δ] ∥g xyy (x, y * u (x))∥ • max t∈[0,δ] ∥ dy * t (x) dt ∥ 2 = δ • S g • M 2 y , (B.10)where M y is the Lipschitz bound on y * as shown in Lemma 3.3, and S g is the smoothness of g from Assumption 2.3. For the second term on the right-hand side of (B.9), we have∥ 1 δ δ t=0 g xy (x, y * (x)) • dy * t (x) dt - dy * (x) dt ∥ ≤ 1 δ • ∥g xy (x, y * (x))∥ • δ t=0 t s=0 ∥ d 2 ds 2 y * s (x)∥ds dt ≤ 1 δ • ∥g xy (x, y * (x))∥ • max s∈[0,δ] ∥ d 2 ds 2 y * s (x)∥ • δ 2 ≤ δ • ∥g xy (x, y * (x))∥ • max s∈[0,δ] ∥ d 2 ds 2 y * s (x)∥ = δ • C g • C y , (B.11)where C g is the bound on smoothness of g as in Assumption 2.3, and C y is the bound on ∥ d 2 y * dx 2 ∥ from Lemma 3.3. For the terms involving the function h, we have∥ λ * δ -λ * δ - dλ * δ dδ | δ=0 ∥ = 1 δ δ t=0 ∥ dλ * t dt - dλ * δ dδ | δ=0 ∥dt = 1 δ δ t=0 t s=0 ∥ d 2 ds 2 λ * s ∥ds • dt ≤ 1 δ max s∈[0,δ] ∥ d 2 ds 2 λ * s ∥ • δ 2 ≤ δ • max s∈[0,δ] ∥ d 2 ds 2 λ * s ∥ = δ • C ℓ , (B.12)where C ℓ is the bound on ∥ δ • (S g M 2 y + C g C y + C ℓ ) ≤ O(δ • (S g • C 3 H • (C g + ∥A∥) 2 • (C g + C f + L f ))).Algorithm 5 The Fully First-Order Method for Bilevel Equality Constrained Problem 1: Input: Current x 0 , accuracy ϵ, perturbation δ = ϵ 2 /8C 2 F R X with C F = 2(L f + C f + C g )C 3 H S g (L g + ∥A∥)T = O(C F (F (x 0 ) -inf F ) √ κϵ -2) oracle calls to f and g, whereC F := 2(L f + C f + C g )C 3H S g (L g + ∥A∥) 2 is the smoothness constant of the hyperobjective F .Proof. We first show the inexact gradient ∇F (x t ) generated in Algorithm 5 is an δ-accurate approximation to the hyper-gradient ∇F (x t ). Consider the inexact gradient defined in (3.2)∥v t -vt ∥ ≤ 1 δ {∥[∇ x g(x t , ŷ * δ ) -∇ x [g(x t , ŷ * )] -[∇ x g(x t , y * δ ) -∇ x [g(x t , y * )∥ + ∥ λ * δ -λ * -[λ * δ -λ * ∥∥A∥} ≤ 2 δ [C g + ∥A∥] δ.Thus we get∥ ∇F (x t ) -∇F (x t )∥ ≤ ∥∇ x f (x t , y * ) -∇ x f (x t , ŷ * )∥ + vt -v t + ∥v t - dy * (x t ) dx ∇ y f (x t , y * (x t ))∥ ≤ C f δ + 2 δ [C g + ∥A∥] δ + C F δ ≤ 2 δ δ [C f + C g + ∥A∥] + C F δ ≤ ϵ 2 4C F R X .Applied to the C F -smooth hyper-objective F , such an inexact gradient oracle satisfies the requirement for Proposition B.6. Thus an ϵ-stationary point with ∥G F (x t )∥ ≤ ϵ (see Eq. (B.14) for the definition of gradient mapping) must be found inN = O( C F (F (x 0 )-F * ) ϵ 2) iterations. Noting the evaluation of inexact solutions (ŷ * , λ * , ŷ * δ , λ * δ ) requires Õ( C g /µ g ) first order oracle evaluations, we arrive at the total oracle complexity of Õ(C g /µ g C F (F (x 0 )-F * ) ϵ 2) for finding an ϵ-stationary point.In this subsection, we state the number of iterations required by projected gradient descent method to find an ϵ-stationary point using inexact gradient oracles. Specifically, we consider the following non-convex smooth problem where the objective F is assumed to be C F -Lipschitz smooth:minimize x∈X F (x). (B.13)Since the feasible region X is compact, we use the norm of the following gradient mapping G F (x) as the stationarity criterionG F (x) := C F (x -x + ) where x + = arg min z∈X z -x - 1 C F ∇F (x) 2 . (B.14)Initialized to some x 0 and the inexact gradient oracle ∇F , the updates of the inexact projected gradient descent method is given by For t=1,2,..., N do:Set x t ← arg min z∈X z -x t-1 - 1 C F ∇F (x t-1 ) 2 . (B.15)The next proposition calculates the complexity result. Proposition B.6. Consider the constrained optimization problem in (B.13) with F being C F -Lipschitz smooth and X having a radius of R. When supplied with a δ = ϵ 2 /4C F R -inexact gradient oracle ∇F , that is, ∥∇F (x) -∇F (x)∥ ≤ δ, the solution generated by the projected gradient descent method (B.15) satisfiesmin t∈[N ] ∥G F (x t )∥ 2 ≤ C F (F (x 0 ) -F * ) N + δC F R, that is, it takes at most O( C F (F (x 0 )-F * ) ϵ 2) iterations to generate some x with ∥G F (x)∥ ≤ ϵ.Proof. By C F -smoothness of F , we havef (x t+1 ) = f (x t - 1 C F G F (x t )) ≤ f (x t ) - 1 C F G F (x t ) ⊤ ∇f (x t ) + 1 2C F ∥ G F (x t )∥ 2 = f (x t ) - 1 2C F ∥ G F (x t )(x t )∥ 2 + 1 C F G F (x t ) ⊤ ( G F (x t ) -∇f (x t )). (B.16) We now show that 1 β G F (x t ) ⊤ ( G F (x t ) -∇f (x t )) ≤ 0. Let y t = x t -1 C F ∇F (x t ), and let y t = x t -1 C F ∇f (x t ). Then have that 1 C F G F (x t ) ⊤ ( 1 C F G F (x t ) -∇f (x t )) = C F (x t -proj X ( y t )) ⊤ (y t -proj X ( y t )) = C F (x t -proj X ( y t )) ⊤ ( y t -proj X ( y t )) + C F (x t -proj X ( y t )) ⊤ (y t -y t ) ≤ C F (x t -proj X ( y t )) ⊤ (y t -y t ) ≤ δC F R,where the penultimate inequality uses the fact that X is a convex set, and R is the diameter of the set X. Combining this with Inequality (B.16), we have that the function decrease per iteration isF (x t+1 ) ≤ F (x t ) - 1 2C F ∥ G F (x t )∥ 2 + δC F R.Summing over N iterations telescopes the terms, we getmin t∈[N ] ∥ G F (x t )∥ 2 ≤ 1 N C F (F (x 0 ) -F * ) + δC F R. Substituting in N = 4 ϵ 2 C F (F (x 0 ) -F * ) and the choice of δ = ϵ 2 /4C F R, we get min t∈[N ] ∥ G F (x t )∥ 2 ≤ ϵ 2 2 .Taking into account the fact that ∥ G F (x t ) -G F (x t )∥ ≤ ∥∇F (x t ) -∇F (x t )∥ ≤ δ, we obtain the desired result.In this subsection, we address the issue of generating approximations to the primal and dual solutions (y * , λ * ) associated with the lower-level problem in Problem 3.1. These approximations are required for computing the approximate hypergradient in Algorithm 1. For notational simplicity, we are going to consider the following constrained strongly convex problem:minimize y∈R d g(y) subject to By = b. (B.17)We propose the following simple scheme to generate approximate solutions to Problem B.17.Compute a feasible ŷ such that ∥ŷ -y * ∥ ≤ δ. Then solveλ = arg min λ∈R m ∥∇ y g(ŷ) -B ⊤ λ∥ 2 . (B.18)The following lemma tells us that λ is close to λ * if B has full row rank. Lemma B.7. Suppose g in Problem B.17 is a C g -Lipschitz smooth, and the matrix B has full row rank such that the following matrix M B is invertibleM B = I B ⊤ B 0 .Then the approximate solution ( λ, ŷ) from (B.18) satisfies ∥ λ -λ * ∥ ≤ ∥M -1 B ∥(1 + C g )δ.Proof. Since (λ * , y * ) satisfy the KKT conditions, they are the solution to the following linear systemI B ⊤ B 0 =M B y * λ * = -∇ y g(y * ) + Iy * b . (B.19)That isy * λ * = M -1 B -∇ y g(y * ) + Iy * b .On the other hand, the approximate solutions (ŷ, λ) in (B.18) satisfiesI B ⊤ B 0 ŷ λ = B ⊤ λ + I ŷ b .We show the right hand side (r.h.s) of the above equation to be close to the r.h.s of (B. 19 ). Let S := {B ⊤ λ : λ ∈ R m } denote the subspace spanned by the rows of B. We can rewrite B ⊤ λ as the projection of ∇g(ŷ) onto S, that is,B ⊤ λ = arg min s∈S ∥∇ y g(ŷ) -s∥ 2 -∇ y g(y * ) =B ⊤ λ * = arg min s∈S ∥∇ y g(y * ) -s∥ 2 ,where the second relation follows from the KKT conditon associated with (λ * , y * ). Since the projection is an non-expansive operation, we have∥B ⊤ λ -(-∇ y g(y * ))∥ = ∥B ⊤ λ -B ⊤ λ * ∥ ≤ ∥∇ y g(ŷ) -∇g(y * )∥ ≤ C g ∥ŷ -y * ∥ ≤ C g δ.We can rewrite (ŷ, λ) as solutions to the following linear system with some ∥τ ∥ ≤ (1+ C g )δ, ŷ λ = M -1 B -∇ y g(y * ) + Iy * + τ b .Thus we get∥ ŷ λ - y * λ * ∥ = ∥M -1 B ∥∥ τ 0 ≤ ∥M -1 B ∥(1 + C g )δ.Now we can just use the AGD method to generate a close enough approximate solution ŷ and call up the Subroutine in (B.18) to generate the approximate dual solution λ.The Projected Gradient Method to Generate Primal and Dual Solutions for a Linearly Constrained Problem 1: Input: accuracy requirement ϵ > 0 and linearly constrained problem min y:By=b g(y). Proof. With N := ⌈4 C g /µ g log∥y * ∥∥M -1 B ∥(Lg+1) µgϵ⌉, Theorem 3.7 in [91] shows that ∥y N -ŷ∥ ≤ ϵ/∥M -1 B ∥(1 + L g ). Then we can apply Lemma B.7 to obtain the desired bound.Our algorithms are based on the Lipschitzness of F , which we prove below. Lemma 4.3. Under Assumption 2.2 and 2.5,F in Problem 4.1 is O(L f L y )-Lipschitz in x.Proof. By Lemma 2.1 of [16] , the hypergradient of F computed with respect to the variable x may be expressed as∇ x F (x) = ∇ x f (x, y * (x)) + dy * (x) dx ⊤ • ∇ y f (x, y * (x)). Since we impose Lipschitzness on f and y * , we can bound each of the terms of ∇ x F (x) by the claimed bound.In this section we analyze Algorithm 2, which as stated in Section 4, requires evaluating only the hyperobjective F (as opposed to estimating the hypergradient in Algorithm 3).The motivation for designing such an algorithm, is that while evaluating ∇F up to α accuracy requires O(α -1 ) gradient evaluations, the hyperobjective value can be estimated at a linear rate:Lemma 4.2 (Proof in Appendix C.1). Given any x, we can return F (x) such that |F (x) -F (x)| ≤ α using O( C g /µ g log(L f /α))first-order oracle calls to f and g.Proof of Lemma 4.2. We note that it suffices to find ỹ * such that ∥ỹ * -y * (x)∥ ≤ α/L f , since setting F (x) := f (x, ỹ * ) will then satisfy| F (x) -F (x)| = |f (x, ỹ * ) -f (x, y * (x))| ≤ L f • α L f =α by Lispchitzness of f , as required. Noting that y * (x) = arg min h(x,y)≤0 g(x, y) is the solution to a constrained smooth, strongly-convex problem with condition number C g /µ g , it is possible to approximate it up to α/L f with O( C g /µ g log(L f /α)) first-order oracle calls using the result of Zhang and Lan [45] .Accordingly, we consider Algorithm 2, which is a zero-order variant of Algorithm 3, whose guarantee is summarized is the theorem below. = min δ 2 , F (x0)-inf F L , ν = δ -ρ, D = Θ νϵ 2 ρ 2 dρ 2 L 2 +α 2 d 2 , η = Θ νϵ 3 ρ 4 (dρ 2 L 2 +α 2 d 2 ) 2 , outputs a point x out such that E[dist(0, ∂ δ F (x out ))] ≤ ϵ + α with T = O d(F (x 0 ) -inf F ) δϵ 3 • L 2 + α 2 ( d δ 2 + dL 2 (F (x 0 ) -inf F ) 2 ) calls to F (•).Combining the result of Theorem C.1 with the complexity of hyperobjective estimation, as given by Lemma 4.2, we obtain convergence to a (δ, ϵ)-stationary point of Problem 4.1 with O(d x δ -1 ϵ -3 ) gradient calls overall.Denoting the uniform randomized smoothingF ρ (x) := E ∥z∥≤1 [F (x + ρ • z)]where the expectation, here and in what follows, is taken with respect to the uniform measure, it is well known [92, Lemma 10 ] thatE ∥w∥=1 d 2ρ (F (x + ρw) -F (x -ρw))w = ∇F ρ (x) , E ∥w∥=1 ∇F ρ (x) -d 2ρ (F (x + ρw) -F (x -ρw))w 2 ≲ dL 2 . (C.1)We first show that replacing the gradient estimator with the inexact evaluations F (•) leads to a biased gradient estimator of F .EQUATIONE ∥w∥=1 ∥ g x ∥ 2 ≲ α 2 d 2 ρ 2 + dL 2 .Proof. For the first bound, we haveE ∥w∥=1 ∥g x -g x ∥ ≤ d 2ρ (2α)E ∥w∥=1 ∥w∥ = αd ρ ,while for the second boundE ∥w∥=1 ∥ g x ∥ 2 = E ∥w∥=1 ∥ g x -g x + g x ∥ 2 ≤ 2E ∥w∥=1 ∥ g x -g x ∥ 2 +2E ∥w∥=1 ∥g x ∥ 2 ≲ d 2 ρ 2 •α 2 +dL 2, where the last step invoked (C.1).We are now ready to analyze Algorithm 2. We denoteα ′ = αd ρ , G = α 2 d 2 ρ 2 + dL 2 . Since x t = x t-1 + ∆ t , we have F ρ (x t ) -F ρ (x t-1 ) = 1 0 ⟨∇F ρ (x t-1 + s∆ t ), ∆ t ⟩ ds = E st∼Unif[0,1] [∇F ρ (x t-1 + s t ∆ t ), ∆ t ] = E [⟨∇F ρ (z t ), ∆ t ⟩] . By summing over t ∈ [T ] = [K × M ], we get for any fixed sequence u 1 , . . . , u K ∈ R d : inf F ρ ≤ F ρ (x T ) ≤ F ρ (x 0 ) + T t=1 E [⟨∇F ρ (z t ), ∆ t ⟩] = F ρ (x 0 ) + K k=1 M m=1 E ∇F ρ (z (k-1)M +m ), ∆ (k-1)M +m -u k + K k=1 M m=1 E ∇F ρ (z (k-1)M +m ), u k ≤ F ρ (x 0 ) + K k=1 Reg M (u k ) + K k=1 M m=1 E ∇F ρ (z (k-1)M +m ), u k ≤ F ρ (x 0 ) + KD G √ M + Kα ′ DM + K k=1 M m=1 E ∇F ρ (z (k-1)M +m ), u kwhere the last inequality follows by combining Lemma C.2 and Lemma C.3. By settingu k := -D M m=1 ∇Fρ(z (k-1)M +m ) ∥ M m=1 ∇Fρ(z (k-1)M +m )∥, rearranging and dividing by DT = DKM we obtain1 K K k=1 E 1 M M m=1 ∇F ρ (z (k-1)M +m ) ≤ F ρ (x 0 ) -inf F ρ DT + G √ M + α ′ = F ρ (x 0 ) -inf F ρ Kν + α 2 d 2 ρ 2 + L 2 d √ M + αd ρ ≤ F ρ (x 0 ) -inf F ρ Kν + αd ρ √ M + L √ d √ M + αd ρ . (C.2)Finally, note that for all m ∈ [M ] :z (k-1)M +m -x k ≤ M D ≤ ν, therefore ∇F ρ (z (k-1)M +m ) ∈ ∂ ν F ρ (x k ) ⊂ ∂ δ F (x k ),where the last containment is due to [46, Lemma 4 ] by using our assignment ρ + ν = δ. Invoking the convexity of the Goldstein subdifferential, this implies that1 M M m=1 ∇F ρ (z (k-1)M +m ) ∈ ∂ δ F (x k ) ,thus it suffices to bound the first three summands in (C.2) by ϵ in order to finish the proof. This happens as long asFρ(x0)-inf Fρ Kν ≤ ϵ 3 , αd ρ √ M ≤ ϵ 3 , and L √ d √ M ≤ ϵ 3 , which imply K ≳ Fρ(x0)-inf Fρ νϵ , M ≳ α 2 d 2 ρ 2 ϵ 2 , and M ≳ L 2 d ϵ 2. By our assignments of ρ and ν, these result inT = KM = O F ρ (x 0 ) -inf F ρ νϵ • α 2 d 2 ρ 2 ϵ 2 + L 2 d ϵ 2 = O (F (x 0 ) -inf F )d δϵ 3 • α 2 d ρ 2 + L 2 = O (F (x 0 ) -inf F )d δϵ 3 • α 2 d • max 1 δ 2 , L 2 (F (x 0 ) -inf F ) 2 + L 2 ,completing the proof.We recall Theorem 4.4 below to keep this section self-contained. Theorem 4.4.Suppose F : R d → R is L-Lipschitz and that ∥ ∇F (•) -∇F (•)∥ ≤ α. Then running Algorithm 3 with D = Θ( δϵ 2 L 2 ), η = Θ( δϵ 3 L 4 ), outputs a point x out such that E[dist(0, ∂ δ F (x out ))] ≤ ϵ + α, with T = O (F (x0)-inf F )L 2 δϵ 3 calls to ∇F (•).Our analysis is inspired by the reduction from online learning to nonconvex optimization given by [47] . To that end, we start by proving a seemingly unrelated result, asserting that online gradient descent minimizes the regret with respect to inexact evaluations. Recalling standard definitions from online learning, given a sequence of linear losses ℓ m (•) = ⟨g m , •⟩, if an algorithm chooses ∆ 1 , . . . , ∆ M we denote the regret with respect to u as Reg M (u) := M m=1 ⟨g m , ∆ m -u⟩ .E [Reg M (u)] ≤ D 2 η M + G 2 M m=1 η m + αDM ,where the expectation is with respect to the (possible) randomness of ( g m ) M m=1 . In particular, settingη m ≡ D G √ M yields E [Reg M (u)] ≲ D G √ M + αDM .Proof. For any m ∈ [M ] :∥∆ m+1 -u∥ 2 = ∥clip D (∆ m -η m g m ) -u∥ 2 ≤ ∥∆ m -η m g m -u∥ 2 = ∥∆ m -u∥ 2 + η 2 m ∥ g m ∥ 2 -2η m ⟨∆ m -u, g m ⟩ , thus ⟨ g m , ∆ m -u⟩ ≤ ∥∆ m -u∥ 2 -∥∆ m+1 -u∥ 2 2η m + η m 2 ∥ g m ∥ 2 ,from which we get thatE ⟨g m , ∆ m -u⟩ = E ⟨ g m , ∆ m -u⟩ + E ⟨g m -g m , ∆ m -u⟩ ≤ ∥∆ m -u∥ 2 -∥∆ m+1 -u∥ 2 2η m + η m 2 E ∥ g m ∥ 2 + E ∥g m -g m ∥ • ∥∆ m -u∥ ≤ ∥∆ m -u∥ 2 -∥∆ m+1 -u∥ 2 2η m + η m 2 G 2 + αD .Summing over m ∈ [M ], we see thatE [Reg M (u)] ≤ M m=1 ∥∆ m -u∥ 2 1 η m - 1 η m-1 + G 2 2 M m=1 η m + M αD ≤ D 2 η M + G 2 M m=1 η m + αDM .The simplification for η m ≡ D G √ M readily follows.We are now ready to analyze Algorithm 3 in the inexact gradient setting.Proof of Theorem 4.4. Since Algorithm 3 has x t = x t-1 + ∆ t , we haveF (x t ) -F (x t-1 ) = 1 0 ⟨∇F (x t-1 + s∆ t ), ∆ t ⟩ ds = E st∼Unif[0,1] [⟨∇F (x t-1 + s t ∆ t ), ∆ t ⟩] = E [⟨∇F (z t ), ∆ t ⟩] .By summing over t∈ [T ] = [K × M ],we get for any fixed sequence u 1 , . . . , uK ∈ R d : inf F ≤ F (x T ) ≤ F (x 0 ) + T t=1 E [⟨∇F (z t ), ∆ t ⟩] = F (x 0 ) + K k=1 M m=1 E ∇F (z (k-1)M +m ), ∆ (k-1)M +m -u k + K k=1 M m=1 E ∇F (z (k-1)M +m ), u k ≤ F (x 0 ) + K k=1 Reg M (u k ) + K k=1 M m=1 E ∇F (z (k-1)M +m ), u k ≤ F (x 0 ) + KD G √ M + KαDM + K k=1 M m=1 E ∇F (z (k-1)M +m ), u kwhere the last inequality follows from Lemma C.3 forG = √ L 2 + α 2 , η = D G √ M , since ∥ g t -∇F (z t )∥ ≤ α (deterministically) for all t ∈ [T ] by assumption. Letting u k := -D M m=1 ∇F (z (k-1)M +m ) ∥ M m=1 ∇F (z (k-1)M +m )∥, rearranging and dividing by DT = DKM , we obtain1 K K k=1 E 1 M M m=1 ∇F (z (k-1)M +m ) ≤ F (x 0 ) -inf F DT + G √ M + α = F (x 0 ) -inf F Kδ + G √ M + α . (C.3) Finally, note that for all k ∈ [K], m ∈ [M ] : z (k-1)M +m -x k ≤ M D ≤ δ, therefore ∇F (z (k-1)M +m ) ∈ ∂ δ F (x k ).Invoking the convexity of the Goldstein subdifferential, we see that1 M M m=1 ∇F (z (k-1)M +m ) ∈ ∂ δ F (x k ) ,thus it suffices to bound the first two summands on the right-hand side in (C.3) by ϵ in order to finish the proof. This happens as long asF (x0)-inf F Kδ ≤ ϵ 2 and G √ M ≤ ϵ 2 . These are equivalent to K ≥ 2(F (x0)-inf F ) δϵ and M ≥ 4 G 2 ϵ 2 , which results in T = KM = O F (x 0 ) -inf F δϵ • L 2 + α 2 ϵ 2 = O (F (x 0 ) -inf F )L 2 δϵ 3 ,completing the proof.Algorithm 7 Perturbed Inexact GD1: Input: Inexact gradient oracle ∇F : R d → R d , initialization x 0 ∈ R d , spatial parameter δ > 0, step size η > 0, iteration budget T ∈ N. 2: for t = 0, . . . , T -1 do 3: Sample w t ∼ Unif(S d-1 ) 4: g t = ∇F (x t + δ • w t ) 5:x t+1 = x t -η g t 6: Output: 4 calls to ∇F (•).x out ∼ Unif{x 0 , . . . , x T -1 }. Theorem C.4. Suppose F : R d → R is L-Lipschitz, and that ∥ ∇F (•) -∇F (•)∥ ≤ α. Then running Algorithm 7 with η = Θ ((F (x0)-inf F )+δL) 1/2 δ 1/2 T 1/2 L 1/2 d 1/4 (α+L) outputs a point x out such that E[dist(0, ∂ δ F (x out ))] ≤ ϵ + √ αL, with T = O (F (x 0 ) -inf F + δL)L 3 √ d δϵProof. Throughout the proof we denotez t = x t + δ • w t . Since F is L-Lipschitz, F δ (x) := E w∼Unif(S d-1 ) [F (x + δ • w)] is L-Lipschitz and O(L √ d/δ)-smooth.By smoothness we getF δ (x t+1 ) -F δ (x t ) ≤ ⟨∇F δ (x t ), x t+1 -x t ⟩ + O L √ d δ • ∥x t+1 -x t ∥ 2 = -η ⟨∇F δ (x t ), g t ⟩ + O η 2 L √ d δ • ∥ g t ∥ 2 = -η ⟨∇F δ (x t ), ∇F (z t )⟩ -η ⟨∇F δ (x t ), g t -∇F (z t )⟩ + O η 2 L √ d δ • ∥ g t ∥ 2 .Noting that E[∇F (z t )] = ∇F δ (x t ) and that ∥ g t ∥ ≤ ∥ g t -∇F (z t )∥ + ∥∇F (z t )∥ ≤ α + L, we see thatE[F δ (x t+1 ) -F δ (x t )] ≤ -ηE ∥∇F δ (x t )∥ 2 + ηLα + O η 2 L √ d δ (α + L) 2 ,which impliesE ∥∇F δ (x t )∥ 2 ≤ E[F δ (x t )] -E[F δ (x t+1 )] η + Lα + O ηL √ d(α + L) 2 δ .Averaging over t = 0, . . . , T -1 and noting thatF δ (x 0 ) -inf F δ ≤ (F (x 0 ) -inf F ) + δL results in E ∇F δ (x out ) 2 = 1 T T -1 t=0 E ∥∇F δ (x t )∥ 2 ≤ (F (x 0 ) -inf F ) + δL ηT +Lα+O ηL √ d(α + L) 2 δ .By Jensen's inequality and the sub-additivity of the square root,E ∇F δ (x out ) ≤ (F (x 0 ) -inf F ) + δL ηT + √ Lα + O   ηL √ d(α + L) 2 δ   . Setting η = √ ((F (x0)-inf F )+δL)δ √ T L √ d(α+L) 2yields the final boundE ∇F δ (x out ) ≲ ((F (x 0 ) -inf F ) + δL) 1/4 L 1/4 d 1/8 (α + L) 1/2 δ 1/4 T 1/4 + √ Lα ,and the first summand is bounded byϵ for T = O ((F (x0)-inf F )+δL)L √ d(L+α) 2 δϵ 4 .Theorem D.1 (Reformulation equivalence). When λ * matches to an optimal dual solution to the lower level problem y * = arg min y g(x, y) s.t. h(x, y) ≤ 0, we show that for each x, the reformulation has the same feasible region of y.Proof. We first show that lower-level feasibility implies feasibility of the reformulated problem. Let y * , λ * = min y max β≥0 g(x, y) + β ⊤ h(x, y) be the primal and the dual solution to the lower level problem with parameter x. We can verify that y * satisfies all the constraints in the reformulation problem. The feasibility condition h(x, y * ) is automatically satisfied. We just need to check:g * (x) := min θ g(x, θ) + (λ * ) ⊤ h(x, θ) = g(x, y * ) + (λ * ) ⊤ h(x, y * ). (D.1)Therefore, x, y * is a feasible point to the reformulation problem.We now show the other direction, i.e., that feasibility of the reformulaed problem implies that of the lower-level problem. Given λ * , let us assume y satisfies g(x, y) ≤ g * λ * (x) and h(x, y) ≤ 0. On the other hand, assume y * , λ * = min y max β≥0 g(x, y) + β ⊤ h(x, y) be the primal and the dual solution. We can show that:g(x, y) + (λ * ) ⊤ h(x, y) ≤ g * (x) := min θ g(x, θ) + (λ * ) ⊤ h(x, θ). (D.2)By the strong convexity of g + (λ * ) ⊤ h, we know that y matches to the unique minimum y * , which implies that y = y * is also a feasible point to the original bilevel problem.By computing the derivative of the KKT conditions in Section 2.1, we get:(∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h) + (∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h) dy * dx + (∇ y h) ⊤ dλ * dx = 0 (E.1) diag(λ * )∇ x h + diag(λ * )∇ y h dy * dx + diag(h) dλ * dx = 0. (E.2) Let I = {i ∈ [d h ]|h(x, y * ) i = 0, λ * i > 0}be the set of active constraints with positive dual solution, and I 1 = {i|h(x, y * ) i ̸ = 0} be the set of inactive constraints and I 2 = {i|h(x, y * ) i = 0, λ * i = 0}. We know that Ī = I 1 ∪ I 2 . For each i ∈ I 1 , due to complementary slackness, we know that λ * i = 0.For i ∈ I 1 in (E.1), we have λ * i ∇ x h(x, y * ) i +λ * i ∇ y h(x, y * ) i dy * dx +h(x, y * ) i dλ * i dx = 0, which implies h(x, y * ) i dλ * i dx = 0 because λ * i = 0. This in turn implies dλ * i dx = 0 because h(x, y * ) i < 0.That means the dual variable λ * i = 0 and has zero gradient dλ * i dx = 0 for any index i ∈ I 1 . Therefore, we can remove row i ∈ I 1 in (E.2) and obtain λ * i = 0 anddλ * i dx = 0.For i ∈ I 2 , the KKT condition in (E.2) is degenerate. Therefore, dλ * i dx can be arbitrary, i.e., nondifferentiable. As a subgradient choice, we can set dλ * i dx = 0 for such i. This choice will also eliminate its impact on the KKT condition in (E.1) because dλ * i dx is set to be 0. By this choice of subgradient, we can also remove row i ∈ I 2 (E.2).I = [h i ] i∈I and λ * I = [λ * i ] i∈I : diag(λ * )∇ x h I + diag(λ * I )∇ y h I dy * dx + diag(h I ) dλ * I dx = 0 =⇒ diag(λ * )∇ x h I + diag(λ * I )∇ y h I dy * dx = 0 (due to h I (x, y * ) = 0). (E.3)In (E.1), due todλ * i dx = 0 for all i ∈ Ī, we can remove dλ * i dx ∀i ∈ Ī in (E.1) by:0 = (∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h) + (∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h) dy * dx + (∇ y h) ⊤ dλ * dx = (∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h) + (∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h) dy * dx + (∇ y h I ) ⊤ dλ * I dx . (E.4)Combining (E.4) and (E.3), we get:(∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h) + (∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h) dy * dx + (∇ y h I ) ⊤ dλ * I dx = 0 diag(λ * )∇ x h I + diag(λ * I )∇ y h I dy * dx = 0,which can be written in its matrix form:∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h ∇ y h ⊤ I diag(λ * I )∇ y h I 0 dy * dx dλ * I dx = - ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h diag(λ * I )∇ x h I (E.5)This concludes the derivation of the derivative of constrained optimization in (5.2).F Inequality case: bounds on primal solution error and constraint violation Lemma 5.1. Given any x, the corresponding dual solution λ * (x), primal solution y * (x) of the lower optimization problem in Problem 4.1, and y * λ * ,α (x) as in (5.5), satisfy:EQUATION).(5.6)Proof. We first provide the claimed bound on ∥y * α1,α2 -y * (x)∥. Part 1: Bound on the convergence of y.Since y * λ * ,α minimizes L α,λ * (x, y), the first-order condition gives us:0 = ∇ y L α,λ * (x, y * λ * ,α). Similarly, we can compute the gradient of L α,λ * (x, y) at y * :∇ y L α (x, y * ) = ∇ y f (x, y * ) + α 1 (∇ y g(x, y * ) + (λ * ) ⊤ ∇ y h(x, y * )) + α 2 ∇ y h I (x, y * ) ⊤ h I (x, y * ) = ∇ y f (x, y * ),where the second step is due to the property of the primal and dual solution: ∇ y g(x, y * ) + (λ * ) ⊤ ∇ y h(x, y * ) = 0 by the stationarity condition in the KKT conditions, and by definition of the active constraints h I where the optimal y * must have h I (x, y * ) = 0.Since, for a sufficiently large α 1 , the penalty function isα 1 µ g -L f ≥ α1µg 2strongly convex in y, we have:α 1 µ g 2 y * -y * λ * ,α ≤ ∇ y L α,λ * (x, y * ) -∇ y L α,λ * (x, y * λ * ,α ) = ∥∇ y f (x, y * )∥ ≤ L f .Therefore, upon rearranging the terms, we obtain the claimed bound:y * -y * α,λ * ≤ 2L f α 1 µ g .Part 2: bound on the constraint violation.When we plug y * into (5.4), we get:L α,λ * (x, y * ) = f (x, y * ) + α 1 (g(x, y * ) + (λ * ) ⊤ h(x, y * ) -g * λ * (x)) + α 2 2 ∥h I (x, y * )∥ 2 = f (x, y * ).Plugging in y * α,λ * , we may obtain:L α,λ * (x, y * λ * ,α ) = f (x, y * λ * ,α ) + α 1 (g(x, y * λ * ,α ) + (λ * ) ⊤ h(x, y * λ * ,α ) -g * (x)) + α 2 2 h I (x, y * λ * ,α ) 2 = f (x, y * λ * ,α ) + α 1 (g(x, y * λ * ,α ) + (λ * ) ⊤ h(x, y * λ * ,α ) -g(x, y * ) -(λ * ) ⊤ h(x, y * )) + α 2 2 h I (x, y * λ * ,α ) 2 ≥ f (x, y * λ * ,α ) + α 1 µ g 2 y * -y * λ * ,α 2 + α 2 2 h I (x, y * λ * ,α ) 2 ,where we used the strong convexity (with respect to y) of g(x, y) + (λ * ) ⊤ h(x, y) and the optimality of y * for g(x, y) + (λ * ) ⊤ h(x, y). By the optimality of y * λ * ,α for L α,λ * , we know thatf (x, y * ) = L α,λ * (x, y * ) ≥ L α,λ * (x, y * λ * ,α ) ≥ f (x, y * λ * ,α ) + α 1 µ g 2 y * -y * λ * ,α 2 + α 2 2 h I (x, y * λ * ,α ) 2 .Therefore, by the Lipschitzness of the function f in terms of y, and the bound ∥y * -y * λ * ,α ∥ ≤ 2L f α1µg , we know that:α 2 2 h I (x, y * λ * ,α ) 2 ≤ f (x, y * ) -f (x, y * λ * ,α ) -α 1 µ g 2 y * -y * λ * ,α 2 ≤ L f y * -y * λ * ,α -α 1 µ g 2 y * -y * λ * ,α 2 ≤ L f y * -y * λ * ,α = O(α -1 1 ).Rearranging terms then gives the claimed bound.The bound on the constraint violation in Lemma 5.1 is an important step in the following theorem.G Proof of Lemma 5.2: gradient approximation for inequality constraints Lemma 5.2. Consider F as in Problem 4.1, L as in (5.4), a fixed x, and y * λ * ,α as in (5.5) . Then under Assumptions 2.2 and 2.5, we have:EQUATION2 ).Proof. First, we recall (5.4) here:L λ * ,α (x, y) = f (x, y) + α 1 g(x, y) + (λ * ) ⊤ h(x, y) -g * (x) + α 2 2 ∥h I (x, y)∥ 2 .Next, recall from Equation D.1, we can express g * (x) = g(x, y * ) + (λ * ) ⊤ h(x, y * ), which we use in the first step below:∇ x F (x) - d dx L λ * ,α (x, y * λ * ,α ) = ∇ x f (x, y * ) + dy * dx ⊤ ∇ y f (x, y * ) -∇ x f (x, y * λ * ,α ) + α 1 (∇ x g(x, y * λ * ,α ) + ∇ x h(x, y * λ * ,α ) ⊤ λ * -α 1 (∇ x g(x, y * ) + ∇ x h(x, y * ) ⊤ λ * ) + α 2 ∇ x h I (x, y * λ * ,α ) ⊤ h I (x, y * λ * ,α ) =∇ x f (x, y * ) -∇ x f (x, y * λ * ,α ) (G.1) + dy * dx ⊤ ∇ y f (x, y * ) - dy * dx ⊤ ∇ y f (x, y * λ * ,α ) (G.2) + dy * dx ⊤ ∇ y f (x, y * λ * ,α ) -α 1 ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h diag(λ * I )∇ x h I ⊤ y * λ * ,α -y * 0 added term 1 -α 2 ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx hλ * diag(λ * I )∇ x h I ⊤ 0 diag(1/λ * I )h I (x, y * λ * ,α ) added term 2 (G.3) + α 1 ∇ x g(x, y * ) -∇ x g(x, y * λ * ,α ) + ∇ x h(x, y * ) ⊤ λ * -∇ x h(x, y * λ * ,α ) ⊤ λ * + ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h diag(λ * I )∇ x h I ⊤ y * λ * ,α -y * 0 added term 1 (G.4) -α 2 ∇ x h I (x, y * λ * ,α ) ⊤ h I (x, y * λ * ,α ) + α 2 ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx hλ * diag(λ * I )∇ x h I ⊤ 0 diag(1/λ * I )h I (x, y * λ * ,α ) added term 2.(G.5)According to (5.2) and (E.5), we letH = ∇ 2 yy g + (λ * ) ⊤ ∇ 2 yy h ∇ y h ⊤ I diag((λ * ) * I )∇ y h I 0 ,which is invertible by Assumption 2.2(iii) and by the fact that we remove all the inactive constraints. We now bound the terms in (G.1), (G.2), (G.3), (G.4), and (G.5).Bounding (G.1) and (G.2): (G.1) can be easily bounded by the smoothness of f in terms of x and y, and the bound on ∥y * -y * (λ * ) * ,α ∥ ≤ O(α -1 1 ) from Lemma 5.1. Therefore, we know:∇ x f (x, y * ) -∇ x f (x, y * (λ * ) * ,α ) ≤ C f y * -y * (λ * ) * ,α ≤ C f • O(α -1 1 ).Similarly, given Assumption 2.5 by which y * (x) is L y -Lipschitz in x, we have the bound dy * dx ≤ L y . Therefore, (G.2) can be bounded by:dy * dx ⊤ ∇ y f (x, y * ) - dy * dx ⊤ ∇ y f (x, y * (λ * ) * ,α ) ≤ C f dy * dx y * -y * (λ * ) * ,α ≤ C f L y • O(α -1 1 ).Using (5.2) to solvedy * dx d(λ * ) * dx = -H -1 ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h diag((λ * ) * I )∇ x h I, we can write:dy * dx ⊤ ∇ y f (x, y * (λ * ) * ,α ) = ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h diag((λ * ) * I )∇ x h I ⊤ (H -1 ) ⊤ -∇ y f (x, y * (λ * ) * ,α ) 0 = - dy * dx ⊤ α 1 ∇ y g(x, y * (λ * ) * ,α ) + ∇ y h(x, y * (λ * ) * ,α ) ⊤ (λ * ) * 0 + α 2 ∇ y h I (x, y * (λ * ) * ,α ) ⊤ h I (x, y * (λ * ) * ,α ) 0 , (G.6)where we use the optimality of y * (λ * ) * ,α from (5.5):∇ y f (x, y * (λ * ) * ,α ) + α 1 ∇ y g(x, y * (λ * ) * ,α ) + ∇ y h(x, y * (λ * ) * ,α ) ⊤ (λ * ) * (G.7) + α 2 ∇ y h I (x, y * (λ * ) * ,α ) ⊤ h I (x, y * (λ * ) * ,α ) = 0.Further, recall that H is non-degenerate by Assumption 2.2, as a result of which, the added term 1 in (G.3) can be modified as follows:∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h diag((λ * ) * I )∇ x h I ⊤ α 1 (y * (λ * ) * ,α -y * ) 0 = ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h diag((λ * ) * I )∇ x h I ⊤ (H -1 ) ⊤ H ⊤ α 1 (y * (λ * ) * ,α -y * ) 0 =α 1 ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h diag((λ * ) * I )∇ x h I ⊤ (H -1 ) ⊤ (∇ 2 yy g + ((λ * ) * ) ⊤ ∇ 2 yy h) ⊤ (y * (λ * ) * ,α -y * ) ∇ y h I (x, y * )(y * (λ * ) * ,α -y * ) . (G.8)The added term 2 in (G.3) can be expanded to:α 2 ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h(λ * ) * diag((λ * ) * I )∇ x h I ⊤ 0 diag(1/(λ * ) * I )h I (x, y * (λ * ) * ,α ) =α 2 ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h(λ * ) * diag((λ * ) * I )∇ x h I ⊤ (H -1 ) ⊤ H ⊤ 0 diag(1/(λ * ) * I )h I (x, y * (λ * ) * ,α ) =α 2 ∇ 2 yx g + ((λ * ) * ) ⊤ ∇ 2 yx h(λ * ) * diag((λ * ) * I )∇ x h I ⊤ (H -1 ) ⊤ ∇ y h I (x, y * ) ⊤ h I (x, y * (λ * ) * ,α ) 0 (G.9)Therefore, we can compute the difference between (G.6), (G.8), and (G.9) to bound (G.3), and use the fact that ∇ y g(x, y * ) + (λ * ) ⊤ ∇ y h(x, y * ) = 0:dy * dx ⊤ ∇ y f (x, y * λ * ,α ) -added term 1 -added term 2 = ∇ 2 yx g + (λ * ) ⊤ ∇ 2 yx h diag(λ * I )∇ x h I ⊤ (H -1 ) ⊤ α 1 ∇ y g(x, y * λ * ,α ) -∇ y g(x, y * ) -∇ 2 yy g(x, y * )(y * λ * ,α -y * ) 0 (G.10) + α 1 ∇ y h(x, y * λ * ,α ) ⊤ λ * -∇ y h(x, y * ) ⊤ λ * -∇ 2 yy h(x, y * ) ⊤ λ * (y * λ * ,α -y * ) 0 (G.11) -α 1 0 ∇ y h I (x, y * )(y * λ * ,α -y * ) (G.12) + α 2 ∇ y h I (x, y * λ * ,α ) ⊤ h I (x, y * λ * ,α ) 0 - ∇ y h I (x, y * ) ⊤ h I (x, y * λ * ,α ) 0 (G.13)The terms in (G.10) and (G.11) can both be bounded by α 1 C g L y ∥y * λ * ,α -y * ∥ 2 and α 1 RC h L y ∥y * λ * ,α -y * ∥ 2 by the smoothness of g and h ⊤ λ * . Further, plugging in ∥y * -y * λ * ,α ∥ ≤ O(α -1 1 ) from Lemma 5.1 bounds both these terms by O(α -1 1 ).To bound the term in (G.12), we use:h I (x, y * λ * ,α ) -h I (x, y * ) -∇ y h I (x, y * )(y * λ * ,α -y * ) ≤ C h y * λ * ,α -y * 2. Therefore, we have:∇ y h I (x, y * )(y * λ * ,α -y * ) ≤ h I (x, y * λ * ,α ) + ∥h I (x, y * )∥ + C h O( y * λ * ,α -y * 2 ) ≤ O(α -1/2 1 α -1/2 2 ) + 0 + O(α -2 1 ) = O(α -1/2 1 α -1/2 2 + α -21 ), which upon scaling by α 1 gives us the following bound on the term in (G.12):α 1 ∇ y h I (x, y * )(y * λ * ,α -y * ) ≤ O(α 1/2 1 α -1/2 2 + α -1 1 ) .The term in (G.13) can be bounded by:α 2 ∇ x h I (x, y * λ * ,α ) ⊤ h I (x, y * λ * ,α ) -∇ x h I (x, y * ) ⊤ h I (x, y * λ * ,α ) =α 2 ∇ x h I (x, y * λ * ,α ) -∇ x h I (x, y * ) O( h I (x, y * α,λ * ) ) =α 2 • O(α -1 1 )O(α -1/2 1 α -1/2 2 ) =O(α -3/2 1 α 1/2 2 ) (G.14)Bounding (G.4): This can be easily bounded by the smoothness of g and h, and the bound on the dual solution ∥λ * ∥ ≤ R. Thus (G.4) can be bounded by REQUATION). Bounding (G.5): By the same argument in (G.14), we get:EQUATION2 ) .Combining all upper bounds gives the claimed bound.H Proof of the main result (Theorem 5.3): convergence and computation cost Theorem 5.3. Given any accuracy parameter α > 0, Algorithm 4 outputs ∇ x F (x) such that ∥ ∇F (x) -∇F (x)∥ ≤ α within O(α -1 ) gradient oracle evaluations.Proof. First, given the bound in Lemma 5.2, we choose α 1 = α -2 and α 2 = α -4 to ensure the inexactness of the gradient oracle is bounded by α. In the later analysis, we will still use α 1 and α 2 in the penalty function for clarity.Lower-level problem. Given the oracle access to the optimal dual solution λ * (x), we can recover the primal solution y * (x) efficiently (e.g, by [45] ). Therefore, we can use the primal and dual solutions to construct the penalty function L λ * ,α (x, y) in (5.4).Penalty function minimization problem. The second main optimization problem is the penalty minimization problem in Line 4 of Algorithm 4. Recall from (5.4) thatL λ,α (x, y) = f (x, y) + α 1 g(x, y) + λ ⊤ h(x, y) -g * (x) + α 2 2 ∥h I (x, y)∥ 2 , (H.1)where we use the approximate dual solution λ as opposed to the optimal dual solution λ * .Given (H.1), we solve the penalty minimization problem:y ′ λ,α (x) := arg min y L λ,α (x, y).The penalty minimization is a unconstrained strongly convex optimization problem, which is known to have linear convergence rate. We further analyze its convexity and smoothness below to precisely estimate the computation cost:Moreover, there is an additional error because we need λ * as opposed to a near-optimal λ to make the analysis in Appendix G work. The error between using λ and λ * in (H.3) can be bounded by:∇ y h(x, y ′ λ,α ) ⊤ (λ -λ * ) ≤ L h α, (H.4)where we use the local Lipschitzness of the function h in an open ball near y * . Therefore, the additional error is also O(α), which matches the inexactness of the inexact gradient oracle.Therefore, we conclude that in order to bound the inexactness of the gradient oracle, we just need an efficient inexact dual solution with α accuracy.I Practical oracle to optimal (approximate) dual solutionHere we discuss how practical the assumption on the oracle access to the optimal dual solution is.For linear inequality constraint h(x, y) = Ax -By -b, the LL problem is a constrained strongly convex smooth optimization problem. To show that we can compute an approximate solution to the optimal dual solution for linear inequality constraints, we apply the result from [93] : Corollary I.1 (Application of Corollary 3.1 in [93] ). When h(x, y) = Ax -By + b is linear in y, the primal and dual solutions can be written as:y * , λ * = arg min y max λ g(x, y) + (λ * ) ⊤ h(x, y) = g(x, y) -(λ * ) ⊤ By + R(x) ⇐⇒ y * , λ * = arg min y max λ g(x, y) -(λ * ) ⊤ By (I.1)where g is strongly convex in y and B is of full rank by Assumption 2.2. According to Corollary 3.1 from [93] , the primal-dual gradient method guarantees a linear convergence. More precisely, in t = O(log 1 α ) iterations, we get:y t -y * ≤ α and λ t -λ * ≤ α. (I.2)Given Corollary I.1, we can efficiently approximate the primal and dual solutions up to high accuracy with O(log 1 α ) oracle calls when the inequality constraints are linear. This gives us an efficient approximate oracle access to the dual solution. Remark I.2. Under the assumption of an optimal dual solution oracle, all the analyses mentioned in Section 5 hold for the general convex inequality constraints. However, the main technical challenge is that the dual solution oracle for general convex inequality cannot be guaranteed in practice. In fact, to the best of our knowledge, there is no iterate convergence in the dual solution λ for general convex inequality constraints. Most of the literature in strongly-convex-concave saddle point convergence only guarantees dual solution convergence in terms of its duality gap or some other merit functions. We are not aware of any successful bound on the dual solution iterate convergence, which is an important research question to answer by itself. This is the main technical bottleneck for general convex inequality constraints as well. Remark I.3. On the other hand, we need the dual solution iterate convergence with rate O(1/α) to ensure the error to be bounded. But this is not a necessary condition. To ensure a bound on the error propagation, we just need to bound some forms of merit functions ((H.2) and (H.4)) of the dual solutions, which we believe that this is much more tractable than the actual iterate dual solution convergence. We leave this as a future direction and this will generalize the analysis from linear inequality constraints to general convex inequality constraints.J The role of λ * (x) in the derivative of Equation (5.4) Notice that Equation (5.4), we treat the dual solution λ * (x) as a constant to define the penalty function derivative. Yet, the dual solution λ * (x) is in fact also a function of x. Therefore, in theory, we should also compute its derivative with respect to x.However, notice that the following:∇ x (λ * (x)) ⊤ h(x, y) = ∇ x h(x, y) ⊤ λ * + dλ(x) dx ⊤ h(x, y) (J.1)The later term in Equation (J.1) can be divided into two cases:• For active constraint i ∈ I with h(x, y * ) = 0, we know that y * λ,α is close to y * by Lemma 5.1. Therefore, the derivative dλ(x) dx ⊤ h(x, y * λ,α ) ≤ L h L λ α 1 = O(α 1 ) = O(α 2) by the local smoothness of h near y * and the Lipschitzness assumption of λ * in Assumption 2.5.• For inactive constraint i ∈ Ī and λ * i > 0, we can solve the KKT conditions and get dλ(x) dx = 0. Therefore, the second term becomes 0. • For inactive constraint i ∈ Ī and λ * i = 0, the KKT system degenerates and we need to use subgradient. By solving the KKT system, we find that dλ(x) dx = 0 is a valid subgradient. Therefore, by choosing this subgradient, the second term also vanishes. Therefore, we do not need to compute the derivative of λ * as the terms involved its derivative is negligible compared to other major terms.All experiments were run on a computing cluster with Dual Intel Xeon Gold 6226 CPUs @ 2.7 GHz and DDR4-2933 MHz DRAM. No GPU was used, and we used 1 core with 8GB RAM per instance of the experiment. The cutoff time for running the algorithms is set to be 6 hours. All experiments were run and averaged over 10 different random seeds. All parameters in the constrained bilevel optimization in Section 6, including the objective parameters and the constrain parameters, are randomly generated from a normal distribution with 0-mean and standard deviation 1.For our fully first-order algorithm, we implement Algorithm 3, where the inexact gradient oracle subroutine is provided by implementing Algorithm 4. All algorithms are implemented in PyTorch [94] to compute gradients, and using Cvxpy [95] to solve the LL problem and the penalty minimization problem. We implement our fully first-order method based on the solutions returned by Cvxpy with certain accuracy requirement, and use PyTorch to compute the inexact gradient discussed in Section 5. We implement the non-fully first-order method using the CvxpyLayer [18] , which is a Cvxpy compatible library that can differentiate through the LL convex optimization problem. where h i (x, y) = x ⊤ A i y -b ⊤ i x ∀i ∈ [d h ] is a d h -dim bilinear constraint, where the constraint bilinear matrix A i ∈ R dx×dy , b i ∈ R dx for all i ∈ [d h ] are randomly generated from normal distributions. The bilinear (nonlinear) constraint of the lower-level problem is the major difference compared to the experiment in Section 6. We are interested in whether our algorithms work beyond the linear constraints where our theory guarantees.The rest of the parameters are the same as in Section 6. The PSD matrix Q ∈ R dy×dy , c ∈ R dy , P ∈ R dx×dy . We compare our Algorithm 3 with a non-fully first-order method implemented using cvxpyLayer [18] . Both algorithms use Adam [90] to control the learning rate in gradient descent. All the experiments are averaged over ten random seeds. Justification: We provide algorithms and corresponding theoretical guarantees for all our claims in the abstract. We provide experiments (and relevant code) as claimed.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
 
 Appendix:  


Decision: This work studies constrained bilevel optimization problems, for which it presents fully first-order methods (thereby avoiding Hessian computations). There is overall consensus among the reviewers that this is an interesting work that is worthy of acceptance, and I am inclined to agree. It would be strongly suggested that, for the camera-ready version, the authors clarify the presentation/conventions of the convergence bounds, in order to better facilitate comparison with existing upper and lower bounds.

Reviews: review: paper_summary: paper_summary: This paper studies the linearly constrained bilevel optimization problem and provides a fully first-order method with solid theoretical analysis. To approximate hypergradient, the penalty method seems novel to me and it is applied in two settings where the LL problem is linear inequality or equality constraints.
strength_weakness: strength_weakness: 1. The experiment does not seem sufficient.
2. The checklist should be behind the appendix as I memorize.
questions: questions: 1. Could the authors explain the design of Algorithm 2? Why is there a clip? 
2. Why do authors not include the experiment of Algorithm 4? 
3. Could authors compare the method in [a] though it studies the unconstrained bilevel optimization problem?
[a] Yang Y, Xiao P, Ji K. Achieving ${O}(\epsilon^{-1.5}) $ Complexity in Hessian/Jacobian-free Stochastic Bilevel Optimization.
limitations: limitations: Please check weaknesses and problems.
paper_summary: paper_summary: The paper provides algorithms for bilevel optimization with linear equality and inequality constraints. The main contribution is that the algorithms are fully first order and do not require Hessian computations. This is achieved by reformulating the linearly constrained bilevel optimization problem using the penalty method and assuming access to the upper-level variable $x$ and approximate dual optimal solution. In this setting an inexact gradient oracle for the problem is constructed. The authors provide stationarity guarantees for their algorithms under certain standard assumptions. They also describe how to use their inexact gradient oracle for non-convex non smooth optimization. The algorithms and theoretical results are backed by small set of proof-of-concept experiments.
strength_weakness: strength_weakness: 1) There are no experiments highlighting the effectiveness of Algorithm 2.
2) The experiments given here are also only of a proof -of -concept nature and are not comprehensive. Some more comparisons with existing algorithms (based on time) and ablation studies may provide more insights regarding the practical applicability of the methods.
questions: questions: There is a lot of literature available on approximating Hessian computations faster for 2nd order methods (for example using iterative sketching etc.) though not necessarily for bilevel optimization. I would appreciate it if you can comment on such methods when compared to yours in terms of practicality.
limitations: limitations: Yes, the limitations have been clarified in the paper.
paper_summary: paper_summary: This works deals with constrained bilevel optimization problems where the lower-level has linear inequality or equality constraints. A set of algorithms is developed that does not require access to the Hessian, but only to zeroth and first-order information. In the case of inequality constraints, convergence is established to an $(\delta,\epsilon)$-Goldstein stationary point, due to the non-smoothness of the problem, using $O(\delta^{-1}\epsilon^{-4})$ or $O(d\delta^{-1}\epsilon^{-3})$ oracle calls (d: dimension of the upper-level).For equality-constrained problems, the convergence rate is nearly optimal with a rate $O(\epsilon^{-2})$.
strength_weakness: strength_weakness: * A major weakness of this work is the small number of experiments. There is only a single example bilevel problem over which experiments are performed. In addition, there are no experiments on real applications (e.g., in machine learning).
* In the experiments section the proposed method is compared only with a single baseline, which does not correspond to any published bilevel method. Why aren’t there any comparisons with other bilevel algorithms that deal with constrained bilevel problems? There are both value function and implicit gradient methods that can, at least in theory, deal with the example problem used here, e.g. [37,38]. 
* The algorithms require access to exact solutions of certain problems, such as $y^{\ast}(x)$ of the lower-level or $y_{\lambda^{\ast},\alpha}^{\ast}(x)$ of problem (3.6). This is not the case in practice.
questions: questions: * In the gradient of the penalty function (3.5.), there is the term $\nabla_x g^{\ast}(x)$. How do you compute this gradient given that $g^{\ast}(x)$ is a min function? The authors should explain this. 
* The authors should explain in more detail Algorithm 2. It is not clear what the utility of each step is.
* The authors derive the equations for the gradient of $F$ and $y^{\ast}(x)$ (eq. 3.3). However, as they note in section 1.1, the hyperfunction $F$ is non-smooth. Then, how can we derive a formula for the gradient of $y^{\ast}(x)$? How do we determine if a given point is differentiable?
* Perhaps the least common assumption in literature is Assumption 2.3. The authors should justify whether this assumption is easy to be satisfied. For instance, is it typically the case that the dual solution is bounded? How about the Lipschitz property of $y^{\ast}(x)$? Are there other relevant (bilevel) works that use the same or similar assumption?
* How does the cvxpyLayer baseline work? This is not explained clearly in the text, besides a minor reference in the Appendix. As this is the only baseline currently used, I believe that some more details are required.
limitations: limitations: * See weaknesses above.
* The authors discuss limitations of their work in section 7.
* No negative societal impact
paper_summary: paper_summary: This paper studies fully first-order methods for bilevel optimization with strongly convex lower-level objectives and linear lower-level constraints, including both inequality constraints and equality constraints. For linear inequality constraints, the paper used a penalty methods to construct hypergradient estimators and demonstrate the approximation error. Then combined with a gradient method for nonsmooth nonconvex optimization to obtain a Goldstein stationarity point. While for linear inequality constraints, the paper proposed to use zeroth-order approximation of the hyper-gradient and demonstrate complexity.
strength_weakness: strength_weakness: 1. How to solve linear inequality constrained case with simple gradient methods to achieve optimal complexity bounds still remain unknown. The current algorithm is a bit complicated and hard to tune in practice. 
2. It is unclear if the complexity bounds are optimal for both linearly inequality and equality constrained settings. Following the conventional notation for stationarity such that $||\nabla F(x)||\leq \epsilon^2$, the complexity bound for Alg 1 is $O(\epsilon^{-8})$, for Alg 2 is $O(\epsilon^{-10})$ if further taking into account the complexity to get $tilde \nabla F$, and for Alg 3 is $O(\epsilon^{-4})$. Note that the paper consider deterministic bilevel optimization, it implies that there are a lot of room for improvements. The paper should mention comparison to unconstrained lower-level bilevel optimization so that readers know the results can be improved. Despite the complexity bounds, I still believe that this paper makes good contribution to the community of bilevel optimization.
3. The methods proposed in the paper are not aligning with each other well. Despite the constrained lower-level bilevel optimization setting, the paper lacks a coherent narrative in terms of algorithmic design.
questions: questions: 1. In the Alg 3 perturbed GD, what is the intuition to use perturbed GD instead of just inexact GD $\tilde \nabla F(x_t)$. In other words, why one needs to perturb using samples from sphere distribution?
2. In Sec 3, although the problem is nonsmooth and nonconvex, given Thm 3.3, I would assume that the $\nabla F(x)$ exists but it is not Lipschitz continuous? 
3. Please be specific about Lipschitz, i.e., does it mean Lipschitz continuous or Lipschitz smooth.
4. What is the reason to consider zeroth-order methods in the linear equality constrained setting? Is it possible to derive fully first-order gradient estimator for the linear equality constrained setting? If possible, why not do it? 
5. In Sec 3.1, assuming knowing $\lambda^*$ means that one also knows $y^*$ and $g^*$. 
6. The definition of complexity is not very clear. In line 229, it says solving a single smooth and strongly convex optimization problems amounts to $\tilde O(1)$ oracle calls. But line 226 implies that one of such oracle call requires $\epsilon^{-1}$ gradient evaluations. Then in line 322, the complexity seems to become $O(log(1/\epsilon))$. Please specify how expensive to obtain $y^*$ and $\lambda^*$ in an $\epsilon$ error and clean oracle calls and complexity bounds in a consistent manner. 
7. What is the motivation of (5.3) the perturbed lower-level problem? 
8. Alg 1 requires computing $y^*_{\lambda^*,\alpha}$, this requires either $f$ to be convex in $y$ or requires $\alpha$ to be large enough to ensure $\mathcal{L}$ is convex in $y$. Please specify explicitly the requirements. Is $\mathcal{L}$ strongly convex in $y$?

I am happy to raise the score if these issues are well-addressed.
limitations: limitations: The paper includes several settings. Yet why adapting a specific strategy to solve each setting is not sufficiently motivated. Such as why one needs to use Alg 2 and Alg 3 instead of directly using the hypergradient estimator with gradient descent. Why one needs to use zeroth-order oracle for the linearly equality constrained setting instead of first-order oracle?

Comments: title: Thank you!
Comment by Reviewer
Replying to further discussions
Discussions
Replying to further discussion
Discussions
comment: The responses addressed my concerns. Please consider improve the presentation for the final version based on the discussions. I raise my score from 5 to 6.
We want to again thank the reviewer for their insightful review and great questions. We hope we have satisfactorily answered them. 

We would like to respectfully ask if, per their original message, the reviewer would consider increasing their score? Thanks again!
Thank you for the clarification. Just to be sure, [3] studied the stochastic case and their complexity bound is $O(\epsilon^{-4})$ with a corresponding lower bounds for stochastic cases.
My main concern during the review was the limited number of experiments. The authors have addressed this issue by providing additional experiments. Given that this is primarily a theoretical work, I also understand the current absence of large-scale experiments.

I am raising my score to 6.
Thank you for answering my questions and they are resolved. I would like to keep my score.
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.
For the linear equality case, when the hyperobjective function is smooth, as the exponent of the stationary notion (under the smooth setting) is 1, the complexity bounds $O(\epsilon^{-2})$ is optimal if considering stochastic first order algorithms. However, the paper considers deterministic first-order algorithms and thus this complexity bound $O(\epsilon^{-2})$ is no longer optimal. Could you clarify what do you mean by "Nonetheless, we wish to point out that our rate for the linear equality case is in fact optimal" in your initial responses?
I have read the rebuttal and the reviews. As of now I will keep my score.
Thank you very much for your questions! We answer them below. 

---
1. ## The exponent of stationarity.

   Both stationarity notions, namely $\epsilon$-stationarity for the equality case and $(\delta, \epsilon)$-stationarity for the inequality case, are defined and addressed throughout the paper with respect to the norm with exponent 1.

---
2. ## Need for different notions of stationarity. 

   The need for different stationarity notions for the two settings (and hence, the different approaches) stems from the crucial difference that **the hyperobjective is not smooth in the inequality case, whereas in the equality case we prove that it is** — as a result, in the equality setting, we can resort to the simpler (more classically understood) notion of $\epsilon$-stationarity [1, 2]. Moreover, for smooth objectives, these two stationarity notions are equivalent, see e.g. Proposition 6 (ii) in [3]. 

   In our paper, we motivate $(\delta, \epsilon)$-stationarity in the (nonsmooth) inequality setting in Lines 45-52. To improve our presentation, we will add a more thorough explanation for these two notions of stationarity in the appendix. 

---
3. ## Does the linear equality algorithm incur a dimension dependence? 

   No, our algorithm for the linear equality setting is dimension-free (cf. Theorem 5.1). This is because, as we show, our functions in this setting are already smooth. 

   The purpose served by the perturbation in this setting is to help us generate a sufficiently good hypergradient approximation (as explained in our previous response). We will clarify this better in the text.  

---
Please feel free to let us know if we can provide any further clarifications! Thank you once again for all your questions and feedback!

---
# References 

[1] Lan, G. (2020). First-order and stochastic optimization methods for machine learning. Springer.

[2] Nesterov, Y., and Polyak, B. Cubic regularization of Newton method and its global performance. Mathematical programming (2006).

[3] Zhang, J., Lin, H., Jegelka, S., Sra, S., & Jadbabaie, A. (2020, November). Complexity of finding stationary points of nonconvex nonsmooth functions. In International Conference on Machine Learning (pp. 11173-11182). PMLR.
Thanks for the detailed responses. Most of my concerns are addressed. I realize that the definition for the stationary notion in the linear equality case is with exponent 2 while in the linear inequality cases is with exponent 1, making it very confusing to comprehend. The complexity bounds in the linear equality case seems indeed optimal in terms of accuracy if the definition is with exponent 2. In addition, with the perturbation methods, does the complexity bound in the linear equality constrained case depend on the dimension of the problem? 

The presentation of the work requires serious improvements.