An Accelerated Algorithm for Stochastic Bilevel Optimization under Unbounded Smoothness

Authors: xiaochuan  gong
jie  hao
mingrui  liu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 07:08:40

Abstract: This paper investigates a class of stochastic bilevel optimization problems where the upper-level function is nonconvex with potentially unbounded smoothness and the lower-level problem is strongly convex. These problems have significant applications in sequential data learning, such as text classification using recurrent neural networks. The unbounded smoothness is characterized by the smoothness constant of the upper-level function scaling linearly with the gradient norm, lacking a uniform upper bound. Existing state-of-the-art algorithms require $\widetilde{O}(\epsilon^{-4})$ oracle calls of stochastic gradient or Hessian/Jacobian-vector product to find an $\epsilon$-stationary point. However, it remains unclear if we can further improve the convergence rate when the assumptions for the function in the population level also hold for each random realization almost surely (e.g., Lipschitzness of each realization of the stochastic gradient). To address this issue, we propose a new Accelerated Bilevel Optimization algorithm named AccBO. The algorithm updates the upper-level variable by normalized stochastic gradient descent with recursive momentum and the lower-level variable by the stochastic Nesterov accelerated gradient descent algorithm with averaging. We prove that our algorithm achieves an oracle complexity of $\widetilde{O}(\epsilon^{-3})$ to find an $\epsilon$-stationary point. Our proof relies on a novel lemma characterizing the dynamics of stochastic Nesterov accelerated gradient descent algorithm under distribution drift with high probability for the lower-level variable, which is of independent interest and also plays a crucial role in analyzing the hypergradient estimation error over time. Experimental results on various tasks confirm that our proposed algorithm achieves the predicted theoretical acceleration and significantly outperforms baselines in bilevel optimization.

Main paper: 
 introduction: 
 Bilevel optimization receives tremendous attention recently in the machine learning community, due to its applications in meta-learning [25, 54] , hyperparameter optimization [25, 23] , data hypercleaning [39] , continual learning [6, 34] , and reinforcement learning [42] . The bilevel optimization problem has the following formulation: where f and g are upper-level and lower-level functions respectively. For example, in metalearning [24, 25] , x denotes the layers of neural networks for shared representation learning, y denotes the task-specific head encoded in the last layer, and the formulation (1) aims to learn the a common representation learning encoder x such that it can be quickly adapted to downstream tasks by only updating the task-specific head y. In machine learning, people typically consider stochastic optimization setting such that f (x, y) = E ξ∼D f [F (x, y; ξ)] and g(x, y) = E ζ∼Dg [G(x, y; ζ)],where D f and D g are the underlying unknown data distributions for f and g respectively, and one can access noisy observations of f and g based on sampling from D f and D g .There emerges a wave of studies for algorithmic design and analysis for solving the bilevel optimization problem (1) under different assumptions of f and g. Most theoretical work assumes the upper-level function is smooth (i.e., gradient is Lipschitz) and nonconvex, and the lower-level function is strongly convex [28, 39, 38, 31, 43] . However, as pointed out by [70, 14] , certain neural networks such as recurrent neural networks [20] , long-short term memory networks [37] and transformers [60] have smoothness constants that scale with gradient norm, potentially leading to unbounded smoothness constants (i.e., gradient Lipschitz constant can be infinity). Motivated by this, Hao et al. [35] designed the first bilevel optimization algorithm to handle the cases where f is nonconvex with potentially unbounded smoothness and g is strongly convex. The algorithm in [35] achieves O(ϵ -4 ) oracle complexity for finding an ϵ-stationary point (i.e., a point x such that ∥∇Φ(x)∥ ≤ ϵ). Gong et al. [30] proposed an single-loop algorithm under the same setting as in [35] and also achieved O(ϵ -4 ) oracle complexity. This complexity result is worse than the O(ϵ -3 ) oracle complexity under the relatively easier setting where f has a Lipschitz gradient, and each realization of the stochastic oracle calls is Lipschitz with respect to its argument (e.g., almost-sure Lipschitz oracle) [66, 41, 17] . This naturally motivates us to study the following question:Is it possible to improve the O(ϵ -4 ) oracle complexity for bilevel optimization problems where the upper-level function is nonconvex with unbounded smoothness and the lower-level function is strongly convex, by assuming that the properties of the function at the population level also hold almost surely for each random realization?In this paper, we give a positive answer to this question by designing a new algorithm named AccBO with an improved oracle complexity of O(ϵ -3 ). Our algorithm is inspired by momentum-based variance reduction techniques used in nonconvex smooth optimization [17] under the almost-sure Lipschitz stochastic gradient oracle framework. The innovation of AccBO lies in its update rules: it employs normalized stochastic gradient descent with recursive momentum for the upper-level variable and stochastic Nesterov accelerated gradient descent with averaging for the lower-level variable. Our approach differs from existing accelerated bilevel optimization algorithms, such as those proposed by [66, 41] in two key ways: (i) while these algorithms use recursive momentum for the upperlevel variable update, AccBO utilizes normalized recursive momentum to address the unbounded smoothness of the upper-level function; (ii) for the lower-level variable update, we use stochastic Nesterov accelerated gradient descent with averaging, in contrast to the recursive momentum method used by the other algorithms. The primary challenge in analyzing the convergence rate of AccBO arises from the need to simultaneously control errors from both upper-level and lower-level variables, given the unbounded smoothness, large learning rate, and recursive momentum in the upper-level problem. Our main contributions are summarized as follows.• We design a new algorithm named AccBO for solving bilevel optimization problems where the upper-level function is nonconvex with unbounded smoothness and the lower-level function is strongly convex. AccBO leverages normalized recursive momentum for the upperlevel variable and Nesterov momentum for the lower-level variable under the stochastic setting to achieve acceleration. To the best of our knowledge, the simultaneous usage of these two techniques in stochastic bilevel optimization is novel and has not been previously explored in the bilevel optimization literature.• We prove that the AccBO algorithm requires O(ϵ -3 ) oracle calls for finding an ϵ-stationary point. This complexity strictly improves the state-of-the-art oracle complexity for unbounded smooth nonconvex upper-level problem and strongly-convex lower-level problem as described in [35, 30] . To achieve this result, we introduce novel proof techniques for analyzing the dynamics of stochastic Nesterov accelerated gradient descent under distribution drift with high probability for the lower-level variable, which are crucial for analyzing the hypergradient error and also of independent interest.• We empirically verify the effectiveness of our proposed algorithm on various tasks, including deep AUC maximization and data hypercleaning. Our algorithm indeed achieves the predicted theoretical acceleration and significantly outperforms baselines in bilevel optimization. 
related work: 
 Relaxed Smoothness. The concept of relaxed smoothness was initially introduced by [70] , inspired by the loss landscapes observed in recurrent neural networks and long-short term memory networks. They show that techniques such as gradient clipping and normalization could improve the performance compared with gradient descent in these scenarios. It inspired further investigations that concentrate on various aspects, including improved analysis on gradient clipping and normalization [69, 40] , adaptive algorithms [14, 46, 63, 22] , federated algorithms [49, 15, 16] , generalized assumptions [14, 13] , and recursive momentum based methods with faster rates [55, 51] . The work of [35, 30] considered a relaxed smoothness condition for the upper-level problem in the bilevel optimization setting.Bilevel Optimization. Bilevel optimization refers to a special kind of optimization where one problem is embedded within another. It was first introduced by [8] . Early works developed specific bilevel optimization algorithms with asymptotic convergence analysis [62, 1, 64] . Ghadimi and Wang [28] initiated the study of non-asymptotic convergence for gradient-based methods in bilevel optimization where the upper-level problem is smooth and the lower-level problem is strongly convex. This field saw further advancements with improved complexity results [38, 39, 10, 19, 11] and fully first-order algorithms [43, 47] . There is a line of work which leverages almost-sure Lipschitz oracle (e.g., stochastic gradient) to obtain improved convergence rates of bilevel optimization algorithms [66, 41] . When the lower-level function is not strongly convex, several algorithmic framework and approximation schemes were proposed [56, 58, 44, 57, 50, 9] . The setting considered in this paper is very close to [35, 30] , where the upper-level function is nonconvex and unbounded smooth, and the lower-level function is strongly convex. However, the work of [35, 30] do not have an accelerated rate O(ϵ -3 ) for finding an ϵ-stationary point as established in this paper.Nesterov Accelerated Gradient and Variants. Nesterov Accelerated Gradient (NAG) method was introduced by [53] for deterministic convex optimization problems. The stochastic version of NAG (SNAG) was extensively studied in the literature [3, 4, 61, 12] . To the best of our knowledge, none of them provide a high probability analysis for SNAG. In the online learning setting, there is a line of work focusing on the perspective of sequential stochastic/online optimization with distributional drift [5, 65, 52, 18] . While these studies provide valuable insights into adaptive techniques and performance bounds under distributional drift, they do not explore the potential integration of such methods with bilevel optimization problems, nor do they consider the application of SNAG within this framework. 
problem setup and preliminaries: 
 Define ⟨•, •⟩ and ∥ • ∥ as the inner product and the Euclidean norm. Throughout the paper, we use asymptotic notation O(•), Θ(•), Ω(•) to hide polylog factors in ϵ -1 and 1/δ. Denote f : R dx × R dy → R as the upper-level function, and g: R dx × R dy → R as the lower-level function. The hypergradient ∇Φ(x) has the following form [28] :EQUATIONTo avoid the Hessian inverse computation, we typically use the following Neumann series to approximate the hypergradient [28, 41, 38] . In particular, for the stochastic setting, define∇f (x, y; ξ) = ∇ x F (x, y; ξ) - Q l g,1 ∇ 2 xy G(x, y; ζ (0) ) q(Q) i=1 I - ∇ 2 yy G(x, y; ζ (i) ) l g,1 ∇ y F (x, y; ξ),where q(Q) ∼ Uniform{0, . . . , Q -1}, ξ := {ξ, ζ (0) , . . . , ζ (q(Q)) } and we use the convention thatj i=1 A i = I if j = 0.Then Eξ[ ∇f (x, y; ξ)] is a good approximation of ∇Φ(x) if y and y * (x) are close [28] .Throughout the paper, we make the following assumptions. Assumption 3.1 ((L x,0 , L x,1 , L y,0 , L y,1 )-smoothness [35] ). Let z = (x, y) and z ′ = (x ′ , y ′ ), there exists L x,0 , L x,1 , L y,0 , L y,1 > 0 such that for all z, z ′ , if ∥z -z ′ ∥ ≤ 1/ 2(L 2 x,1 + L 2 y,1 ), then ∥∇ x f (z) -∇ x f (z ′ )∥ ≤ (L x,0 + L x,1 ∥∇ x f (z)∥)∥z -z ′ ∥ and ∥∇ y f (z) -∇ y f (z ′ )∥ ≤ (L y,0 + L y,1 ∥∇ y f (z)∥)∥z -z ′ ∥.Remark: Assumption 3.1 is introduced by [35] for describing the bilevel optimization problems with recurrent neural networks. This assumption can be regarded as a block-wise relaxed smoothness assumptions for two blocks x and y, which is a variant of the relaxed smoothness assumption [70] and the coordinate-wise relaxed smooth assumption [14] . Assumption 3.2. Suppose the followings hold for objective functions f and g: (i) f is continuously differentiable and (L x,0 , L x,1 , L y,0 , L y,1 )-smooth in (x, y); (ii) For every x, ∥∇ y f (x, y)∥ ≤ l f,0 for all y; (iii) For every x, g(x, y) is µ-strongly-convex in y for µ > 0; (iv) g is l g,1 -smooth jointly in (x, y); (v) g is twice continuously differentiable, and ∇ 2 xy g, ∇ 2 yy g are l g,2 -Lipschitz jointly in (x, y).Remark: Assumption 3.2 is standard in the bilevel optimization literature [43, 35, 28] . Assumption 3.2 (i) characterizes the unbounded smoothness of the upper-level function and is empirically observed in recurrent neural networks [35] . Assumption 3.3. The following stochastic estimators are unbiased and have the following properties:E ξ∼D f [∥∇xF (x, y; ξ) -∇xf (x, y)∥ 2 ] ≤ σ 2 f,1 , E ξ∼D f [∥∇yF (x, y; ξ) -∇yf (x, y)∥ 2 ] ≤ σ 2 f,1 , Pr(∥∇yG(x, y; ξ) -∇yg(x, y)∥ ≥ λ) ≤ 2 exp(-2λ 2 /σ 2 g,1 ) ∀λ > 0, E ζ∼Dg [∥∇ 2 xy G(x, y; ζ) -∇ 2 xy g(x, y)∥ 2 ] ≤ σ 2 g,2 , E ζ∼Dg [∥∇ 2 yy G(x, y; ζ) -∇ 2 yy g(x, y)∥ 2 ] ≤ σ 2 g,2 .Remark: Assumption 3.3 assumes the stochastic oracle for the upper-level problem has bounded variance, which is standard in nonconvex stochastic optimization [26] [27] [28] . It also assumes the stochastic oracle for the lower-level problem is light-tailed, which is common for the high probability analysis for the lower-level problem [45, 36] . Note that the same assumption is also made in [35, 30] for the bilevel problems with a unbounded smooth upper-level function. Remark: Assumption 3.4 assumes that each random realization of the upper-and lower-level functions satisfies the same property as in the population level. Note that this condition is the key to achieve an improved O(ϵ -3 ) oracle complexity under various settings, including both single-level nonconvex smooth problems [21, 17, 59 ] and bilevel problems with nonconvex smooth upper-level objectives [66, 41] . Furthermore, this assumption is shown to be necessary for achieving improved oracle complexity in single-level problems [2] .4 Algorithm and Analysis 
main challenges and algorithm design: 
 Main Challenges. We begin by explaining why existing bilevel optimization algorithms and their corresponding analysis techniques are insufficient in our setting. First, most algorithms developed for bilevel optimization require the upper-level function is smooth (i.e., the gradient of the upperlevel function is Lipschitz) [28, 39, 38, 66, 41, 19, 43] . They characterize the estimation error of the optimal solution for the lower-level problem, utilize an approximate hypergradient descent approach and the descent lemma for L-smooth functions to prove the convergence. In particular, they demonstrate that a potential function, incorporating both the function value and the bilevel error from the lower-level problem, progressively decreases in expectation. However, when the upper-level function is (L x,0 , L x,1 , L y,0 , L y,1 )-smooth as illustrated in Assumption 3.1, the previous algorithms and analyses relying on L-smoothness do not work. The reason is that the hypergradient bias depends on the approximation error of the lower-level variable as well as the hypergradient itself: these elements are statistically dependent and the standard potential function argument with an expectation-based analysis would not work. Second, the recent work of Hao et al. [35] and Gong et al. [30] considered that the upper-level function is unbounded smooth and addressed this issue by performing normalized stochastic gradient with momentum for the upper-level variable and periodic updates or stochastic gradient descent for the lower-level variable, but their oracle complexity is notAlgorithm 1 STOCHASTIC NESTEROV ACCELERATED GRADIENT METHOD (SNAG) 1: Input: x, ỹ-1 , ỹ0 , α, T 0 # SNAG(x, ỹ0 , α, T 0 ) 2: for t = 0, 1, . . . , T 0 -1 do 3:Sample πt from distribution D g 4: zt = ỹt + γ(ỹ t -ỹt-1 )ỹt+1 = zt -α∇ y G(x, zt ; πt ) 6: end for Algorithm 2 ACCELERATED BILEVEL OPTIMIZATION ALGORITHM (ACCBO) 1: Input: α init , α, β, γ, η, τ, I, S, T 0 , T , set x 0 , y init 0 = 0 2: y 0 = SNAG(x 0 , y init 0 , α init , T 0 ), and set y -1 = ŷ0 = y 0 # Warm-start 3: for t = 0, 1, . . . , T -1 do for j = 0, 1, . . . , N -1 do 14:EQUATIONz j t = y j t + γ(y j t -y j-1 t )15: x t+1 = x t -η mt ∥mt∥ 25: end for better than O(ϵ -4 ). These facts indicate that we need new algorithm design and analysis techniques to get potential acceleration.y j+1 t = z j t -α∇ y G(x t ,Algorithm Design. To obtain potential acceleration, our key idea is to update the upper-level variable by normalized stochastic gradient descent with recursive momentum and update the lower-level variable by the stochastic Nesterov accelerated gradient (SNAG) method. Different from [35, 30] , the key innovation of our algorithm design is that we achieve acceleration for both upper-level and lower-level problems simultaneously but without affecting each other. The upper-level update rule can be regarded as a generalization of the acceleration technique (e.g., the momentum-based variance reduction technique) [17, 51] from single-level to bilevel problems. The main challenge is that we need to deal with the accumulated error of the recursive momentum over time due to the hypergradient bias, which is caused by the inaccurate estimation of the optimal lower-level variable. Therefore we require a very small tracking error between the iterate of the lower-level variable and the optimal lower-level solution defined by the upper-level variable (i.e., y * (x)) at every iteration. This requirement is satisfied by executing SNAG method under the distribution drift for the lower-level problem, where the drift is caused by the change of the upper-level variable over time.The detailed description of our algorithm is illustrated in Algorithm We first introduce some useful notations. Let σ(•) be the σ-algebra generated by the random variables in the argument. We define the following filtrations for t ≥ 1: F init = σ(π 0 , . . . , πT0-1 ), F t = σ( ξ0 , . . . , ξt-1 ), F 1 t = σ(π 0 , . . . , π t-1 ), and we also defineF 2 t = σ(π 0 t , . . . , π N -1 t) when t is a multiple of I. We use E t , E Ft and E to denote the conditional expectation E[• | F t ], the expectation over F t and the total expectation over F T respectively. Theorem 4.1. Suppose Assumptions 3.1 to 3.4 hold. Let {x t } be the iterates produced by Algorithm 2. For any given δ ∈ (0, 1) and small enough ϵ (see exact choice in (53)), we set parameters α init , α, β, γ, η, τ, I, N, S, Q, T 0 (see exact choices in (54), ( 55), ( 56), (57) , and (58)) asα init = Θ(ϵ 4 ), α = Θ(ϵ 2 ), 1 -β = Θ(ϵ 2 ), η = Θ(ϵ 2 ), τ = Θ(ϵ), γ = O(1), T0 = O(ϵ -2 ), I = O(ϵ -1 ), N = O(ϵ -1 ), Q = O(1), S = O(1).Then with probability at least 1 -2δ over the randomness in σ(F init ∪ F 1 T ) (for Option I) or σ(F init ∪ (∪ t≤T F 2 t )) (for Option II), Algorithm 2 guarantees 1 T T t=1 E∥∇Φ(x t )∥ ≤ 20ϵ within T = 4d0 ηϵ = O(ϵ -3) iterations, where d 0 := Φ(x 0 ) -inf x Φ(x) and the expectation is taken over the randomness over F T . For Option I, it requires T 0 + SQT = O(ϵ -3 ) oracle calls of stochastic gradient or Hessian/Jacobian vector product. For Option II, it requires T 0 + N T I + SQT = O(ϵ -3 ) oracle calls of stochastic gradient or Hessian/Jacobian vector product.Remark: Theorem 4.1 established an improved O(ϵ -3 ) oracle complexity for finding an ϵ-stationary point. This complexity result strictly improves the O(ϵ -4 ) obtained by [35, 30] when the upper-level function is nonconvex and unbounded smooth. This complexity result also matches that in the single-level unbounded smooth setting [51] and is nearly optimal in terms of the dependency on ϵ [2] . The full statement of Theorem 4.1 is included in Theorem E.2.In this section, we provide a roadmap of proving Theorem 4.1 and the main steps. The detailed proofs can be found in Appendix D and E. The key idea is to prove two things: (1) the lower-level iterate is very close to the optimal lower-level variable at every iteration; (2) two consecutive iterates of the lower-level iterates are close to each other. In particular, define y * t = y * (x t ), and we aim to prove that ∥ŷ t -y * t ∥ ≤ O(ϵ) and ∥ŷ t+1 -ŷt ∥ ≤ O(ϵ 2 ) for every t. These two requirements are essential to control the hypergradient estimation error (i.e., ∥m t -∇Φ(x t )∥) caused by inaccurate estimate of the lower-level problem. Lemma 4.7 provides the guarantee for the lower-level problem, and Lemma 4.8 characterizes the hypergradient estimation error. Equipped with these two lemmas, we can adapt the momentum-based variance reduction techniques [17, 51] to the upper-level problem and prove the main theorem.The main technical contribution of this paper is to provide a general framework for proving the convergence of SNAG under distributional drift in Section 4.3.1, which can be leveraged as a tool to control the lower-level error in bilevel optimization and derive the Lemma 4.7, as illustrated in Section 4.3.2. In particular, we can regard the change of the upper-level variable x at each iteration as the distributional drift for the lower-level problem: the drift is small due to the normalization operator of the upper-level update rule and also the Lipschitzness of y * (x). Once we have the general lemma for tracking the minimizer for any fixed distributional drift over time, this lemma can be applied to our algorithm analysis and establish guarantees for the bilevel problem.In this section, we study the sequences of stochastic optimization problems min w∈R d ϕ t (w) indexed by time t ∈ N. We denote the minimizer and the minimal value of ϕ t as w * t and ϕ * t , and we define the minimizer drift at time t to be ∆ t := ∥w * t -w * t+1 ∥. With a slight abuse of notation2 , we consider the SNAG algorithm applied to the sequence {ϕ t } T t=1 , where T is the total number of iterations:EQUATIONwhere g t = ∇ϕ t (z t ; ξ t ) is the stochastic gradient evaluated at z t with random sample ξ t . Define ε t = g t -∇ϕ t (z t ) as the stochastic gradient noise at t-th iteration. Define H t = σ(ξ 1 , . . . , ξ t-1 ) as the filtration, which is the σ-algebra generated by all random variables until t-th iteration. We make the following assumption, which is the same as Assumption 3 in [18] for high probability analysis. Assumption 4.2. Function ϕ t : R d → R is µ-strongly convex and L-smooth for constants µ, L > 0. Also, there exists constants ∆, σ > 0 such that the drift ∆ 2 t is sub-exponential conditioned on H t with parameter ∆ 2 and the noise ε t is norm sub-Gaussian conditioned on H t with parameter σ/2. Lemma 4.3. Suppose Assumption 4.2 holds and let {w t } be the iterates produced by the update rule (3) with constant learning rate α ≤ min{1/2L, 1/10 3 µ, µ/10 3 L 2 }, and set γ =1- √ µα 1+ √ µα . Define θ t = [(w t -w * t ) ⊤ , (w t-1 -w * t )] ⊤ ∈ R 2d, and the potential function V t asV t = θ ⊤ t Pθ t + ϕ t (w t ) -ϕ t (w * t ), where P = 1 2α 1 √ µα -1 √ µα -1 (1 - √ µα) 2 ⊗ I d .Then for any given δ ∈ (0, 1) and all t ≥ 0, the following holds with probability at least 1 -δ over the randomness in H t (here e denotes the base of natural logarithms):(i) (With drift) Let ϕ t (w) := µ 2 ∥w -w * t ∥ 2 and w ∈ R, then V t ≤ 1 - √ µα 4 t V 0 + 2ασ 2 + 80∆ 2 α ln eT δ .(ii) (Without drift) Let ϕ t (w) ≡ ϕ(w) be any general functions and w ∈ R d with ∆ = 0, thenV t ≤ 1 - √ µα 4 t V 0 + 2ασ 2 ln eT δ .Remark: When {ϕ t } T t=1 is a sequence of one-dimensional quadratic functions with moving minimizers, Lemma 4.3 provides a high probability tracking guarantee for SNAG with distributional drift, which is useful to provide guarantees for Option I in Algorithm 2. Note that this guarantee strictly improves the guarantee of stochastic gradient descent with distributional drift (e.g., Theorem 6 in [18] ) and therefore is of independent interest. In particular, for small α, the decaying factor in the first term is improved from 1 -µα 2 to 1 -√ µα 4 , and the drift term is improved from∆ 2 α 2 to ∆ 2 α .To the best of our knowledge, this is the first high probability guarantee with improved rate for SNAG under distributional drift. When there is no drift, Lemma 4.3 also provides a high probability guarantee for SNAG. It holds for any smooth and strongly convex function ϕ with any dimension d, and it is useful to provide guarantees for Option II of Algorithm 2.Inspired by [35, 30] , we can regard ϕ t (•) as ϕ t (•) := g(x t , •) in the bilevel setting, and then we have ∆ t = ηl g,1 /µ for every t due to the upper-level update rule and the Lipschitzness of y * (x). Therefore we can focus on the high probability analysis on the lower-level variable without worrying about the randomness from the upper-level. Throughout, we assume Assumption 3.1, 3.2, 3.3 and 3.4 hold. In addition, the failure probability δ ∈ (0, 1) and ϵ > 0 are chosen in the same way as in Theorem 4.1.t } be the iterates produced by line 2 of Algorithm 2. Set α init = Θ(ϵ 4 ) and ϕ t (y) ≡ g(x 0 , y). Then ∥y init T0 -y * 0 ∥ ≤ µα 32 ϵ L0 holds with probability at least 1 -δ over the randomness in F init (we denote this event as E init ) in T 0 = O(ϵ -2 ) iterations.Remark: Lemma 4.4 shows that for fixed initialization x 0 , running SNAG for at most T 0 = O(ϵ -2 ) iterations can guarantee that the Euclidean distance between the lower-level variable y init T0 and the optimal solution y * (x 0 ) is at most O(ϵ), with high probability. Lemma 4.5 (Option I). Under event E init , let {y t } be the iterates produced by Option I. Set α = Θ(ϵ 2 ) and ϕ t (y) = g(x t , y) = µ 2 ∥y -y * t ∥ 2 with y ∈ R. Then for any t ∈ [T ], Algorithm 2 guarantees with probability at least 1 -δ over the randomness in F 1 T (we denote this event asE 1 y ) that ∥y t -y * t ∥ ≤ ϵ/2L 0 . Lemma 4.6 (Option II). Under event E init , let {y t } be the iterates produced by Option II. Set α = Θ(ϵ 2 ), N = O(ϵ -1 ), I = O(ϵ -1) and ϕ t (y) = g(x t , y) when t is a multiple of I (i.e., x t is fixed for each update round of Option II so g can be general functions). Then for any t ∈ [T ], Algorithm 2 guarantees with probability at least 1 -δ over the randomness in σ(∪ t≤T F 2 t ) (we denote this event as E 2 y ) that ∥y t -y * t ∥ ≤ ϵ/L 0 .Remark: Lemma 4.5 and Lemma 4.6 show that, under event E init and both option I and option II, the algorithm guarantees that each iterate y t is O(ϵ)-close to the the optimal lower-level variable y * t at every iteration t with high probability. Lemma 4.7 (Averaging). Under eventE init ∩ E 1 y (Option I) or E init ∩ E 2 y (Option II), set τ = √ µαin the averaging step (line 21 of Algorithm 2). Then for any t ≥ 0 we have∥ŷ t -y * t ∥ ≤ 2ϵ L0 and ∥ŷ t+1 -ŷt ∥ ≤ µϵ 2 24L 2 0 σg,1 =: ϑ.Remark: Lemma 4.7 shows that after performing averaging operations over the sequence {y t } T t=1 , the averaged sequence enjoys stronger guarantees. First, each averaged iterate ŷt is still O(ϵ)-close to the optimal lower-level variable y * t ; Second, two consecutive averaged iterates (i.e., ŷt and ŷt+1 ) is O(ϵ 2 )-close to each other. The stronger guarantees are crucial to control the hypergradient estimation error as described in Lemma 4.8. Lemma 4.8. Under event E init ∩ E 1 y (Option I) or E init ∩ E 2 y (Option II), define ϵ t = m t -E t [ ∇f (x t , ŷt ; ξt )], then we have the following averaged cumulative error bound:1 T T -1 t=0 E∥ϵt∥ ≤ σ T (1 -β) + 1 -β σ + L0 √ 1 -β 2(η 2 + ϑ 2 ) S + L1 2(η 2 + ϑ 2 ) S(1 -β) 1 T T -1 t=0 E∥∇Φ(xt)∥,where S denotes the batch size, and σ, L0 , L1 are defined in Lemmas B.4 and B.6.Remark: Lemma 4.8 characterizes the upper-level hypergradient estimation error under the good event that the lower-level error can be controlled. One can choose hyperparameters appropriately such that the cumulative error (i.e., LHS) grows only sublinearly in terms of T , which is important for establishing the fast convergence of our algorithm. 
experiments: 
 Deep AUC Maximization with Recurrent Neural Networks. AUC (Area Under the ROC Curve) [33] is a critical metric in evaluating the performance of binary classification models. It measures the ability of the model to distinguish between positive and negative classes, and it is defined as the probability that the prediction score of a positive example is higher than that is a negative example [32] . Deep AUC maximization [48, 67] can be formulated as a minmax optimization problem [48] :min w∈R d ,(a,b)∈R 2 max α∈R f (w, a, b, α) := E z [F (w, a, b, α; z)],where samples in the population. This min-max formulation is an special case of the bilevel problem with g = -f in (1), which can be reformulated as the following:F (w, a, b, α; z) = (1 -r)(h(w; x) -a) 2 I [c=1] + r(h(w; x) -b) 2 I [c=-1] + 2(1 + α)(rh(w; x)I [c=-1] -(1-r)h(w; x)I [c=1] )-r(1-r)α 2 ,EQUATION) where (w, a, b) denotes the upper-level variable, and α denotes the lower-level variable. In this case, the lower-level is a one-dimensional quadratic function in terms of α and is strongly convex, and the upper-level function is non-convex function with potential unbounded smoothness when using a recurrent neural network as the predictive model.We aim to perform imbalanced text classification task and maximize the AUC metric. The Deep AUC maximization experiment is performed on imbalanced Sentiment140 [29] dataset (under the license of CC BY 4.0), which is a binary text classification task. Specifically, we follow [68] to make training set imbalanced with a pre-defined imbalanced ratio (r), and leave the test set unchanged. Given r, we randomly discard the positive samples (with label 1) in original training set until the portion of positive samples equals to r. The imbalance ratio r is set to 0.2 in our experiment, which means only 20% data is positive in the training set. We use a two-layer recurrent neural network with input dimension=300, hidden dimension=4096, and output dimension=2 for the model prediction.We compare with some bilevel optimization baselines, including StocBio [39] , TTSA [38] , SABA [19] , MA-SOBA [11] , SUSTAIN [41] , VRBO [66] and BO-REP [35] . We show the training and test AUC result with 25 epochs in (a) (b) of Figure 1 EQUATION)where w is the model parameter of a neural network, and σ(x) = 1 1+e -x is the sigmoid function. We perform bilevel optimization algorithms on the noisy text classification dataset Stanford Natural Language Inference (SNLI) [7] (under the license of CC BY 4.0) with a three-layer recurrent neural network with input dimension=300, hidden dimension=4096, and output dimension=3 for the label prediction. Each of sentence-pairs manually labeled as entailment, contradiction, and neutral. Specifically, the label of each training data is randomly flipped to one of the other two labels with probability p. We set p = 0.1 and p = 0.2 in the experiments, respectively. We compare all the baselines used in the deep AUC maximization experiment. Different from the formulation (4) for the deep AUC maximization, the lower-level function in (5) is not quadratic function of the lower-level variable. Therefore we choose Option II in Algorithm 2, i.e., periodic updates for the lower-level variable. The results are presented in Figure 2 (p = 0.1 and p = 0.2). Our algorithm AccBO exhibits the highest classification accuracy on training and test set among all the bilevel baselines, and also shows a high runtime efficiency. More detailed parameter tuning and selection can be found in Appendix G. All the experiments are run on the device of NVIDIA A6000 (48GB memory) GPU and AMD EPYC 7513 32-Core CPU. 
conclusion: 
 In this paper, we propose a new algorithm named AccBO for solving bilevel optimization problems where the upper-level is nonconvex and unbounded smooth and the lower-level problem is strongly convex. The algorithm achieved O(ϵ -3 ) oracle complexity for finding an ϵ-stationary point, which matches the rate of the state-of-the-art single-level relaxed smooth optimization [51] and is nearly optimal in terms of dependency on ϵ [2] .Limitations. One limitation of our work is that the convergence analysis for the Option I of our algorithm relies on the lower-level problem being a one-dimensional quadratic function: only under this case the algorithm becomes a single-loop procedure. However, it remains unclear how to design single-loop algorithms for more general lower-level strongly convex functions in high dimension. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our paper does not involve the release of any data or models.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: Our paper uses existing text classification datasets and are cited in Section 5 and their licenses are mentioned.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix: In this section, we will introduce a few useful lemmas. The following technical lemma on recursive control is crucial for providing high probability guarantee of the lower-level variables y t and ŷt in Algorithm 2 at anytime. We follow a similar argument as in [18, Proposition 29] with a slight generalization.Lemma A.1 (Recursive control on MGF). Consider scalar stochastic processes (V t ), (V ′ t,1 ), (V ′ t,2 ), (D t,1 ), (D t,2 ) and (X t ) on a probability space with filtration (H t ), which are linked by the inequalityV t+1 ≤ α t V t + D t,1 V ′ t,1 + D t,2 V ′ t,2 + X t + κ tfor some deterministic constants α t ∈ (-∞, 1] and κ t ∈ R. Suppose the following properties hold.• V t , V ′ t,1 and V ′ t,2 are non-negative and H t -measurable. • D t,i is mean-zero sub-Gaussian conditioned on H t with deterministic parameter σ i , andV ′ t,i ≤ V t for i = 1, 2: E[exp(λD t,i ) | H t ] ≤ exp(λ 2 σ 2 i /2) for all λ ∈ R.• X t is non-negative and sub-exponential conditioned on H t with deterministic parameter ν t :E[exp(λX t ) | H t ] ≤ exp(λν t ) for all 0 ≤ λ ≤ 1/ν t .Then the estimateE[exp(λV t+1 )] ≤ exp(λ(ν t + κ t ))E[exp(λ(1 + α t )V t /2)]holds for any λ satisfying 0 ≤ λ ≤ min1-αt 2(σ 2 1 +σ 2 2 ) , 1 2νt .Proof of Lemma A.1. For any index t ≥ 0 and any scalar λ ≥ 0, the law of total expectation impliesE[exp(λV t+1 )] ≤ E exp(λ(α t V t + D t,1 V ′ t,1 + D t,2 V ′ t,2 + X t + κ t )) = exp(λκ t )E exp(λα t V t )E exp λ D t,1 V ′ t,1 + D t,2 V ′ t,2 exp(λX t ) | H t .Hölder's inequality in turn yieldsE exp(λα t V t )E exp λ D t,1 V ′ t,1 + D t,2 V ′ t,2 exp(λX t ) | H t ≤ E exp 2λ D t,1 V ′ t,1 + D t,2 V ′ t,2 | H t • E [exp(2λX t ) | H t ] ≤ exp 2λ 2 (σ 2 i V ′ t,1 + σ 2 i V ′ t,2 ) exp(2λν t ) ≤ exp λ 2 (σ 2 1 + σ 2 2 )V t exp(λν t ) provided 0 ≤ λ ≤ 1/2ν t , where we use V ′ t,i ≤ V t for i = 1, 2 in the last inequality. Therefore, under the condition that0 ≤ λ ≤ min 1 -α t 2(σ 2 1 + σ 2 2 ) , 1 2ν t ,the following estimate hold for all t ≥ 0:E[exp(λV t+1 )] ≤ exp(λκ t )E exp(λα t V t ) exp λ 2 (σ 2 1 + σ 2 2 )V t exp(λν t ) = exp(λ(ν t + κ t ))E exp λ(α t + λ(σ 2 1 + σ 2 2 ))V t ≤ exp(λ(ν t + κ t ))E[exp(λ(1 + α t )V t /2)],where the last inequality follows by the given range of λ. Thus the proof is completed.Next, we introduce the following Young's inequality beyond Euclidean norm cases. This lemma serves as an important role when dealing with distributional drift for high probability SNAG analysis.Lemma A.2 (Young's inequality). For any vectors v 1 , v 2 ∈ R d , positive semidefinite (PSD) matrix Q ∈ R d×d , and scalar c > 0, it holds that 3∥v 1 + v 2 ∥ 2 Q ≤ (1 + c)∥v 1 ∥ 2 Q + 1 + 1 c ∥v 2 ∥ 2 Q .Proof of Lemma A.2. By definition of ∥ • ∥ Q , we haveEQUATIONSince Q ∈ R d×d is PSD, let Q = UU ⊤ be the Cholesky decomposition, thenEQUATIONwhere we use Young's inequality and definition of ∥ • ∥ Q for the second and third lines, respectively. Combing ( 6) and (7) gives the result as claimed.In this section, we provide important properties of the objective function Φ in bilevel optimization problems, as well as characterizations (such as variance and bias) for stochastic hypergradient estimator ∇f (x, y; ξ) based on Neumann series. For readers' convenience, we only list the results here and defer the detailed proofs to Appendix F. ∥∇Φ(x)-∇Φ(x ′ )∥ ≤ (L 0 +L 1 ∥∇Φ(x ′ )∥)∥x-x ′ ∥ if ∥x-x ′ ∥ ≤ 1 2(1 + l 2 g,1 /µ 2 )(L 2 x,1 + L 2 y,1 ), where (L 0 , L 1 )-smoothness constant L 0 and L 1 are defined asL 0 = 1 + l 2 g,1 µ 2 L x,0 + L x,1 l g,1 l f,0 µ + l g,1 µ (L y,0 + L y,1 l f,0 ) + l f,0 l g,1 l g,2 + l g,2 µ µ 2 and L 1 = 1 + l 2 g,1 µ 2 L x,1 .Lemma B.3 (Descent inequality, [35, Lemma 10] ). Suppose Assumptions 3.1 and 3.2 and 3.2 hold. Then for any x, x ′ we haveΦ(x) ≤ Φ(x ′ )+⟨∇Φ(x ′ ), x-x ′ ⟩+ L 0 + L 1 ∥∇Φ(x ′ )∥ 2 ∥x-x ′ ∥ 2 if ∥x-x ′ ∥ ≤ 1 2(1 + l 2 g,1 /µ 2 )(L 2 x,1 + L 2 y,1 ).Lemma B.4 ([41, Lemma B.1]). Under Assumptions 3.1 to 3.4, the bias of the stochastic hypergradient estimate of the upper-level objective satisfies∥ ∇f (x, y) -Eξ[ ∇f (x, y; ξ)]∥ ≤ l g,1 l f,0 µ 1 - µ l g,1 Q, where Q is the number of samples chosen to approximate the Hessian inverse. Moreover, we haveEξ[∥ ∇f (x, y; ξ) -Eξ[ ∇f (x, y; ξ)]∥ 2 ] ≤ σ 2 f,1 +Lemma B.5. Under Assumptions 3.1 to 3.4, we have∥ ∇f (x, y) -∇Φ(x)∥ ≤ ( L + L x,1 ∥∇Φ(x)∥)∥y -y * (x)∥,where constant L is defined asL := L x,0 + L x,1 l g,1 l f,0 µ + l g,1 µ (L y,0 + L y,1 l f,0 ) + l f,0 µl g,2 + l g,1 l g,2 µ 2 ≤ L 0 .Lemma B.6. Under Assumptions 3.1 to 3.4, we have (i) For any fixed y ∈ R dy and any x, x ′ ∈ R dx ,Eξ∥ ∇f (x, y; ξ) -∇f (x ′ , y; ξ)∥ 2 ≤ ( L2 0 + L2 1 ∥∇Φ(x)∥ 2 )∥x -x ′ ∥ 2 .(ii) For any fixed x ∈ R dx and any y, y ′ ∈ R dy ,Eξ∥ ∇f (x, y; ξ) -∇f (x, y ′ ; ξ)∥ 2 ≤ ( L2 0 + L2 1 ∥∇Φ(x)∥ 2 )∥y -y ′ ∥ 2 .In the above expressions, we define L0 and L1 asL0 = 4 L x,0 + L x,1 l g,1 l f,0 µ + L x,0 + L x,1 l g,1 l f,0 µ ∥y -y * (x)∥ 2 + 6Q 2µl g,1 -µ 2 l 2 g,1 (L y,0 + L y,1 l f,0 ) 2 + l 2 f,0 l 2 g,2 + l 2 f,0 l 2 g,1 l 2 g,2 Q 2 (l g,1 -µ) 2 1/2 , L1 = 2L x,1 (1 + L x,1 ∥y -y * (x)∥).Note that in Lemma B.6, constant L0 depends on the value of ∥y -y * (x)∥. When we consider this term in Algorithm 2, it turns into ∥y t -y * t ∥ or ∥ŷ t -y * t ∥, which are both as small as O(ϵ) (and thus bounded) with high probability by Lemmas 4.5 to 4.7. In other words, we can treat this term as another constant for our algorithm and analysis.For convenience, we will restate a few concepts included in Section 4.3.1 here. We consider the sequences of stochastic optimization problemsEQUATIONindexed by time t ∈ N. We denote the minimizer and the minimal value of ϕ t as w * t and ϕ * t , and we define the minimizer drift at time t to be ∆ t := ∥w * t -w * t+1 ∥. With a slight abuse of notation, we consider the SNAG algorithm applied to the sequence {ϕ t } T t=1 , where T is the total number of iterations:EQUATIONwhere g t = ∇ϕ t (z t ; ξ t ) is the stochastic gradient evaluated at z t with random sample ξ t . Define ε t = g t -∇ϕ t (z t ) as the stochastic gradient noise at t-th iteration. Define H t = σ(ξ 1 , . . . , ξ t-1 ) as the filtration, which is the σ-algebra generated by all random variables until t-th iteration. We will make the following standard assumption, as illustrated below 4 .Assumption C.1. The sequences of time-varying functions satisfy that, each function ϕ t : R d → R is µ-strongly convex and L-smooth for some constants µ, L > 0. Assumption C.2 (Sub-Gaussian drift and noise). There exists constants ∆, σ > 0 such that the following holds for all t ≥ 0:(i) (Drift) The drift ∆ 2 t is sub-exponential conditioned on H t with parameter ∆ 2 : E exp(λ∆ 2 t ) | H t ≤ exp(λ∆ 2 ) for all 0 ≤ λ ≤ ∆ -2 .(ii) (Noise) The noise ε t is norm sub-Gaussian conditioned on H t with parameter σ/2:Pr {∥ε t ∥ ≥ ϱ | H t } ≤ 2 exp(-2ϱ 2 /σ 2 ) for all ϱ > 0.The following lemma characterize the one-step improvement for stochastic Nesterov accelerated gradient method. Although part of our analysis is similar to [12, 4] , our final goal is quite different: we aim to derive a careful formulation (see (13) ) such that we can apply Lemma A.1 to recursively control the moment generating function of V t with distributional drift, thus leading to a high probability bound for V t at anytime (see Lemma C.5), while [12, 4] only show the convergence in expectation without distributional drift.Lemma C.3 (Distance recursion, with drift). Suppose that Assumptions C.1 and C.2 hold. Let {w t } be the iterates produced by update rule (9) with constant learning rate α ≤ 1/2L and set constants γ, ρ > 0, and matrix P ∈ R 2d×2d asEQUATION)Define θ t = [(w t -w * t ) ⊤ , (w t-1 -w * t )] ⊤ ∈ R 2d, also define the potential function and u t,1 , u t,2 asEQUATION)Then for all t ≥ 0, it holds thatEQUATIONSpecifically, if ϕ t (w) := µ 2 ∥w -w * t ∥ 2 with d = 1, then we haveEQUATION)Proof of Lemma C.3. We first apply Lemma A.2 withv 1 + v 2 = θ t+1 = w t+1 -w * t+1 w t -w * t+1 , v 1 = w t+1 -w * t w t -w * t , v 2 = w * t -w * t+1 w * t -w * t+1 , Q = P to obtain Vt+1 = θ ⊤ t+1 Pθt+1 + ϕt+1(wt) -ϕt+1(w * t+1 ) ≤ 1 + √ µα 4 wt+1 -w * t wt -w * t ⊤ P wt+1 -w * t wt -w * t + 1 + 4 √ µα w * t -w * t+1 w * t -w * t+1 ⊤ P w * t -w * t+1 w * t -w * t+1 + ϕt+1(wt+1) -ϕt+1(w * t+1 ) = 1 + √ µα 4 wt+1 -w * t wt -w * t ⊤ P wt+1 -w * t wt -w * t + ϕt(wt+1) -ϕt(w * t ) (A) + ϕt+1(wt+1) -ϕt+1(w * t+1 ) -1 + √ µα 4 (ϕt(wt+1) -ϕt(w * t )) (B) + 1 + 4 √ µα w * t -w * t+1 w * t -w * t+1 ⊤ P w * t -w * t+1 w * t -w * t+1 (C).Bounding (A). Let us define vector ω t ∈ R 2d and matrices A, B ∈ R 2d×2d asω t = ∇ϕ t (z t ) ε t , A = 1 + γ -γ 1 0 ⊗ I d , B = -α -α 0 0 ⊗ I d .By (9) we have θ t+1 = Aθ t + Bω t . Since ϕ t is µ-strongly convex, thenEQUATION)By L-smoothness of ϕ t and the fact that w t+1 = z t -αg t , we haveEQUATIONwhere we use g t = ∇ϕ t (z t ) + ε t in the second equality. Noting that by (9) we havew t -z t = -γ(w t -w * t ) + γ(w t-1 -w * t ), and combining ( 14) and ( 15) we obtainEQUATIONwhere matrix X 1 ∈ R 4d×4d is defined asX 1 = 1 2    µγ 2 -µγ 2 -γ 0 -µγ 2 µγ 2 γ 0 -γ γ α(2 -Lα) 0 0 0 0 0    ⊗ I d .Then applying the strong convexity of ϕ t again givesEQUATION)Noting thatw * t -z t = (w * t -w t ) -γ(w t -w * t ) + γ(w t-1 -w * t ) = -(1 + γ)(w t -w * t ) + γ(w t-1 -w * t ), and combining ( 15) and ( 17) we obtainEQUATIONwhere matrix X 2 ∈ R 4d×4d is defined asX 2 = 1 2    µ(1 + γ) 2 -µγ(1 + γ) -(1 + γ) 0 -µγ(1 + γ) µγ 2 γ 0 -(1 + γ) γ α(2 -Lα) 0 0 0 0 0    ⊗ I d .Next, we multiply ( 16) by ρ 2 and ( 18) by 1 -ρ 2 , then sum them up to getEQUATIONWith definition of γ, ρ 2 and P given in the statement of the theorem, by direct calculation we haveA ⊤ PA -ρ 2 P A ⊤ PB B ⊤ PA B ⊤ PB -(ρ 2 X 1 + (1 -ρ 2 )X 2 ) =       - µ(1- √ µα) 3 2 √ µα(1+ √ µα) µ(1- √ µα) 3 2 √ µα(1+ √ µα) 0 -1+µα 2(1+ √ µα) µ(1- √ µα) 3 2 √ µα(1+ √ µα) - µ(1- √ µα) 3 2 √ µα(1+ √ µα) 0 1- √ µα 2(1+ √ µα) 0 0 -α(1-Lα) 2 α 2 -1+µα 2(1+ √ µα) 1- √ µα 2(1+ √ µα) α 2 α 2       ⊗ I d := C,where matrix C is negative semidefinite [12] , i.e., C ⪯ 0. Then we use this fact to leveragew t+1 -w * t w t -w * t ⊤ P w t+1 -w * t w t -w * t -ρ 2 θ ⊤ t Pθ t = θ t ω t ⊤ (ρ 2 X 1 + (1 -ρ 2 )X 2 + C) θ t ω t ≤ θ t ω t ⊤ (ρ 2 X 1 + (1 -ρ 2 )X 2 ) θ t ω t ≤ -(ϕ t (w t+1 ) -ϕ t (w * t )) + ρ 2 (ϕ t (w t ) -ϕ t (w * t )) -α(1 -Lα)⟨∇ϕ t (z t ), ε t ⟩ + Lα 2 2 ∥ε t ∥ 2 ,where the last inequality follows by (19) . Rearrange the above inequality and by definition of the potential function V t , we obtainEQUATION)Now recall that in (8) our objective function has the form of ϕ t (w) = µ 2 ∥w -w * t ∥ 2 , hence ∇ϕ t (z t ) = µ(z t -w * t ). Plugging this into the above inequality givesEQUATION)where u t,1 and u t,2 are defined asu t,1 = w t -w * t ∥w t -w * t ∥ , u t,2 = z t -w t ∥z t -w t ∥ .Therefore, we conclude thatEQUATION)where the last inequality follows from the definition of ρ and simple calculation1 + √ µα 4 ρ 2 = 1 + √ µα 4 (1 - √ µα) ≤ 1 - 3 √ µα 4 .Bounding (B). Again, recall that under distributional drift, our objective function in (8) has the form of ϕ t (w) = µ 2 ∥w -w * t ∥ 2 , then we have(B) = ϕ t+1 (w t+1 ) -ϕ t+1 (w * t+1 ) -1 + √ µα 4 (ϕ t (w t+1 ) -ϕ t (w * t )) ≤ 1 + √ µα 4 (ϕ t+1 (w t+1 ) -ϕ t+1 (w * t+1 ) -ϕ t (w t+1 ) + ϕ t (w * t )) = µ 2 1 + √ µα 4 (∥w t+1 -w * t+1 ∥ 2 -∥w t+1 -w * t ∥ 2 ) = µ 2 1 + √ µα 4 ∥w * t -w * t+1 ∥∥w t+1 -w * t + w t+1 -w * t+1 ∥ ≤ µ 2 1 + √ µα 4 ∆ t (2∥w t+1 -w * t+1 ∥ + ∥w * t+1 -w * t ∥)Since ϕ t+1 is µ-strongly convex and matrix P is PSD, thenV t+1 = θ ⊤ t+1 Pθ t+1 + ϕ t+1 (w t+1 ) -ϕ t (w * t+1 ) ≥ ϕ t+1 (w t+1 ) -ϕ t (w * t+1 ) ≥ µ 2 ∥w t+1 -w * t+1 ∥ 2 =⇒ ∥w t+1 -w * t+1 ∥ ≤ 2 µ V t+1 .Plugging the above fact back into the upper bound for (B) givesEQUATIONBounding (C). For this part, we mainly handle the distributional drift. Recall that under this condition we consider one-dimensional case, that is, we let d = 1 in (8) such that w ∈ R. Then we use P ∈ R 2×2 to obtainEQUATIONwhere in the last equality we use the definition of P and the fact thatEQUATION)Final Bound for V t+1 . Now we are ready to derive the upper bound for V t+1 . Combining ( 22), ( 23) and ( 24) together yieldsEQUATIONFor simplicity, we define D as the followingD = 1 - 3 √ µα 4 V t + -2µα(1 -Lα)⟨u t,1 , ε t ⟩ µ 2 ∥w t -w * t ∥ - 2µ √ 2α 1 + √ µα α(1 -Lα)⟨u t,2 , ε t ⟩ 1 + √ µα 2 √ 2α ∥z t -w t ∥ + Lα 2 2 ∥ε t ∥ 2 + µ 1 + √ µα 8 + 2 √ µα ∆ 2 t .Hence ( 26) turns intoV t+1 -2µ∆ t 1 + √ µα 4 V t+1 -D ≤ 0.Solving the above inequality we getV t+1 ≤ 1 2   2µ∆ t 1 + √ µα 4 + 2µ∆ t 1 + √ µα 4 2 + 4D   ≤ 2µ∆ t 1 + √ µα 4 + √ DThen an application of Young's inequality revealsV t+1 ≤ 1 + √ µα 4 D + 1 + 4 √ µα 2µ∆ t 1 + √ µα 4 2 = 1 + √ µα 4 1 - 3 √ µα 4 V t + -2µα(1 -Lα)⟨u t,1 , ε t ⟩ µ 2 ∥w t -w * t ∥ - 2µ √ 2α 1 + √ µα α(1 -Lα)⟨u t,2 , ε t ⟩ 1 + √ µα 2 √ 2α ∥z t -w t ∥ + Lα 2 2 ∥ε t ∥ 2 + µ 1 + √ µα 8 + 2 √ µα ∆ 2 t + 1 + 4 √ µα 2µ∆ 2 t 1 + √ µα 4 2 ≤ 1 - √ µα 2 V t + 1 + √ µα 4 -2µα(1 -Lα)⟨u t,1 , ε t ⟩ µ 2 ∥w t -w * t ∥ - 2µ √ 2α 1 + √ µα α(1 -Lα)⟨u t,2 , ε t ⟩ 1 + √ µα 2 √ 2α ∥z t -w t ∥ + Lα 2 2 ∥ε t ∥ 2 + µ 1 + √ µα 4 1 + √ µα 8 + 2 √ µα ∆ 2 t + 1 + 4 √ µα 2µ∆ 2 t 1 + √ µα 4 2 ,where we plug in the definition of D for the first equality. Since the learning rate α ≤ 1/L and thus µα ≤ 1, then we have1 + √ µα 4 1 + √ µα 8 + 2 √ µα ≤ 125 32 √ µα , 2 1 + 4 √ µα 1 + √ µα 4 2 ≤ 125 8 √ µα .Fact (II) verification. By definition of matrix P and simple calculation we haveV t ≥ θ ⊤ t Pθ t = 1 2α ∥(w t -w * t ) + ( √ µα -1)(w t-1 -w * t )∥ = 1 2α ∥(1 - √ µα)(w t -w t-1 ) + √ µα(w t -w * t )∥ ≥ 1 2α (1 - √ µα)∥w t -w t-1 ∥ - µ 2 ∥w t -w * t ∥ ≥ 1 2α (1 - √ µα)∥w t -w t-1 ∥ -V t .Rearrange the above inequality, and recall the update rule of stochastic Nesterov accelerated gradient method, we have∥z t -w t ∥ = γ∥w t -w t-1 ∥ ≤ 2γ √ 2α 1 - √ µα V t = 2 √ 2α 1 + √ µα V t ,where for the last equality we use the definition of γ as in (10) . Rearrange it gives (30) .By Lemma C.3 we haveEQUATIONNote that under Assumption C.2, there exists an absolute constant c ≥ 1 such that for all t ≥ 0, ∥ε t ∥ 2 is sub-exponential conditioned on H t with parameter cσ 2 , and ε t is mean-zero sub-Gaussian conditioned on H t with parameter cσ [18, Theorem 30] . For convenience we simply let c = 1 here. Thus ⟨u t,1 , -ε t ⟩ is mean-zero sub-Gaussian conditioned on H t with parameter σ, and ∆ 2 t is sub-exponential conditioned on H t with parameter ∆ 2 by assumption. Hence, in light of (31), we apply Lemma A.1 withH t = H t , V t = V t , V ′ t,1 = µ 2 ∥w t -w * t ∥ 2 , V ′ t,2 = (1 + √ µα) 2 8α ∥z t -w t ∥ 2 , D t,1 = 1 + √ µα 4 2µα(1-Lα)⟨u t,1 , -ε t ⟩, D t,2 = 1 + √ µα 4 2µ √ 2α 1 + √ µα α(1-Lα)⟨u t,2 , -ε t ⟩, X t = 1 + √ µα 4 Lα 2 2 ∥ε t ∥ 2 + 20µ∆ 2 t √ µα , α t = 1 - √ µα 2 , κ t = 0, σ 1 = 1 + √ µα 4 2µα(1 -Lα)σ, σ 2 = 1 + √ µα 4 2µ √ 2α 1 + √ µα α(1 -Lα)σ, ν t = 1 + √ µα 4 Lα 2 2 σ 2 + 20µ∆ 2 √ µα ,yielding the following recursionE[exp(λV t+1 )] ≤ exp 1 + √ µα 4 Lα 2 2 σ 2 + 20µ∆ 2 √ µα E exp λ 1 - √ µα 4 V t (32) for all λ satisfying 0 ≤ λ ≤ min 2 125α √ µασ 2 , 1 5ασ 2 /4 + 40µ∆ 2 / √ µα .We then apply (32) recursively to deduceE[exp(λV t+1 )] ≤ exp λ 1 - √ µα 4 t V 0 + λ 1 + √ µα 4 Lα 2 2 σ 2 + 20µ∆ 2 √ µα t-1 i=0 1 - √ µα 4 i ≤ exp λ 1 - √ µα 4 t V 0 + 2 1 + √ µα 4 Lα 2 √ µα σ 2 + 80∆ 2 α for all λ satisfying 0 ≤ λ ≤ min 2 125α √ µασ 2 , 1 5ασ 2 /4 + 40µ∆ 2 / √ µα .Moreover, setting ν := 2ασ 2 + 80∆ 2 α and taking into account α ≤ min{1/2L, µ/2L 2 , 1/10 3 µ}, then we have2 1 + √ µα 4 Lα 2 √ µα σ 2 + 80∆ 2 α ≤ ν and 1 ν = 1 2ασ 2 + 80∆ 2 /α ≤ min 2 125α √ µασ 2 , 1 5ασ 2 /4 + 40µ∆ 2 /√ µα .Thus we obtainE exp λ V t -1 - √ µα 4 t V 0 ≤ exp(λν) for all 0 ≤ λ ≤ 1/ν.Taking λ = 1/ν and applying Markov's inequality and union bound completes the proof.The following result shows the second part of Lemma 4.3. Lemma C.6 (High-probability distance tracking, without drift). Under the same setting as in Lemma C.4 with α ≤ min{1/2L, µ/10 3 L 2 }, for any given δ ∈ (0, 1) and all t ∈ [T ], the following holds with probability at least 1 -δ over the randomness in H t :EQUATIONProof of Lemma C.6. First it is easy to verify that1 L µ 2 ∥∇ϕ t (w t ) -∇ϕ t (w * )∥ ≤ V t and 1 + √ µα 2L √ 2α ∥∇ϕ t (z t ) -∇ϕ t (w t )∥ ≤ V t .Then we apply Lemma A.1 to obtain the final result. We omit the detailed proof here since it follows the same procedure as in proof of Lemma C.5.We first present the following algebraic fact under suitable choice of parameters. Lemma D.1 (Parameter choice, informal). For any given δ ∈ (0, 1) and any small ϵ satisfyingEQUATIONif we set parameters α, β, η, T asEQUATIONwhere σ is defined in Lemma B.4, and d 0 and P are defined asEQUATION)Then the following holds for all t ∈ [T ]:4ασ 2 g,1 µ + 160η 2 l 2 g,1 µ 3 α ln eT δ ≤ ϵ 2 64L 2 0 . Proof of Lemma D.1. By Lemma B.1, we have ∆ t = ∥y * t -y * t+1 ∥ ≤ lg,1µ ∥x t -x t+1 ∥ = ηl g,1 /µ. Thus in our bilevel setting, we choose ∆ = ηl g,1 /µ, where ∆ is defined in Section 4.3.1. By choice of α, η, T as in (35) , we have4ασ 2 g,1 µ + 160η 2 l 2 g,1 µ 3 α ln eT δ = 4(1 -β)σ 2 g,1 µ 2 + 160η 2 l 2 g,1 µ 2 (1 -β) ln 4ed 0 δηϵ ≤ 164(1 -β)σ 2 g,1 µ 2 ln 4ed 0 δϵ(1 -β) max l g,1 σ g,1 , σ 4d 0 .Now we choose β to be1 -β = µ 2 ϵ 2 164 • 64L 2 0 σ 2 g,1 ln(P ), where P = 164• 64ed 0 L 2 0 σ 2 g,1 δµ 2 ϵ 3 max l g,1 σ g,1 , σ d 0 2 .Then we have164(1 -β)σ 2 g,1 µ 2 ln 4ed 0 l g,1 δϵσ g,1 (1 -β) = ϵ 2 64L 2 0 ln(P ) ln √ P ln(P ) ≤ ϵ 2 64L 2 0 ,where we use the fact that ln( √ P ln(P )) ≤ ln(P ) ≤ ln 2 (P ) for any P ≥ 4 by choice of ϵ as in (34) .In the rest of this section, we assume Assumptions 3.1 to 3.4 hold. In addition, the failure probability δ ∈ (0, 1) and ϵ > 0 are chosen in the same way as in Theorem 4.1.Lemma D.2 (Warm-start, Restatement of Lemma 4.4). Let {y init t } be the iterates produced by line 2 of Algorithm 2. Set α init = µα 2 = Θ(ϵ 4 ) with α defined in (35) and ϕ t (y) ≡ g(x 0 , y). Then ∥y init T0 -y * 0 ∥ ≤ µα 32 ϵ L0 holds with probability at least 1 -δ over the randomness in F init (we denote this event as E init ) in T 0 = O(ϵ -2 ) iterations, whereEQUATIONProof of Lemma D.2. By Lemmas C.6 and D.1 and µ-strong convexity of g in y, we have with probability at least 1 -δ over the randomness in F init that∥y init T0 -y * 0 ∥ 2 ≤ 2 µ 1 - µ 2 α 2 4 T0 U init 0 + 4µα 2 σ 2 g,1 µ ln eT 0 δ ≤ 2 µ 1 - µα 4 T0 U init 0 + µαϵ 2 64L 2 0 .By l g,1 -smoothness of g we haveU init 0 ≤ 2 -µ 2 α 2 + µ 2 α 2 2µα 2 ∥y init 0 -y * 0 ∥ 2 + g(x 0 , y init 0 ) -g(x 0 , y * 0 ) ≤ 3 2µα 2 ∥y init 0 -y * 0 ∥ 2 + l g,1 2 ∥y init 0 -y * 0 ∥ 2 ≤ 2 µα 2 ∥y init 0 -y * 0 ∥ 2 ,where the last inequality uses α ≤ 1/l g,1 . Now we set2 µ 1 - µα 4 T0 2 µα 2 ∥y init 0 -y * 0 ∥ 2 + µαϵ 2 64L 2 0 ≤ µαϵ 2 32L 2 0 , which gives T 0 ≥ ln µ 3 α 3 ϵ 2 256L 2 0 ∥y init 0 -y * 0 ∥ 2 ln 1 - µα 4 .By choice of α as in (35) and simple calculation we obtain T 0 = O(ϵ -2 ) when ϵ is small. Under event E init , let {y t } be the iterates produced by Option I. Set α = Θ(ϵ 2 ) as in (35) and ϕ t (y) = g(x t , y) = µ 2 ∥y -y * t ∥ 2 with y ∈ R. Then for any t ∈ [T ], Algorithm 2 guarantees with probability at least 1 -δ over the randomness in F 1T (we denote this event asE 1 y ) that ∥y t -y * t ∥ ≤ ϵ/2L 0 .Proof of Lemma D.3. By Lemmas C.5 and D.2 we have∥y t -y * t ∥ 2 ≤ 2 µ 1 - √ µα 4 t U 0 + 4ασ 2 g,1 µ + 160η 2 l 2 g,1 µ 3 α ln eT δ ≤ 2 µ 1 - √ µα 4 t 2 α ∥y init T0 -y * 0 ∥ 2 + 4ασ 2 g,1 µ + 160η 2 l 2 g,1 µ 3 α ln eT δ ≤ 4 µα ∥y init T0 -y * 0 ∥ 2 + ϵ 2 64L 2 0 ≤ ϵ 2 4L 2 0 .Thus we conclude that for all t ∈ [T ], we have ∥y t -y * t ∥ ≤ ϵ/2L 0 .Lemma D.4 (Option II, Restatement of Lemma 4.6). Under event E init , let {y t } be the iterates produced by Option II. Set α = Θ(ϵ 2 ) as in (35) and run SNAG in each update round forN = ln µα 128 ln 1 - √ µα 4 = O(ϵ -1 ) steps in every I = µϵ 2(1-β)L0σg,1 = O(ϵ -1) iterations, set ϕ t (y) = g(x t , y) when t is a multiple of I (i.e., x t is fixed for each update round of Option II so g can be general functions). Then for any t ∈ [T ], Algorithm 2 guarantees with probability at least 1 -δ over the randomness in σ(∪ t≤T F 2 t ) (we denote this event asE 2 y ) that ∥y t -y * t ∥ ≤ ϵ/L 0 .Proof of Lemma D.4. At the beginning of the first round, by Lemmas D.2 and D.3 we have ∥y 0y * 0 ∥ ≤ ϵ/2L 0 , then we do not update the lower-level variable until t = I-th iteration, then for t = I, we have∥y I -y * I ∥ = ∥y 0 -y * I ∥ ≤ ∥y 0 -y * 0 ∥ + I i=1 ∥y * i -y * i-1 ∥ ≤ ϵ 2L 0 + ηl g,1 µ I = ϵ L 0 ,where in the last equality we plug in the definition of η and I. By l g,1 -smoothness of g we haveU I = 2 -2 √ µα + µα 2α ∥y I -y * 0 ∥ 2 + g(x I , y I ) -g(x I , y * I ) ≤ 3 2α ∥y I -y * I ∥ 2 + l g,1 2 ∥y I -y * I ∥ 2 ≤ 2ϵ 2 αL 2 0Then for N steps update in the inner loops of t = I-th iteration, we setEQUATIONBy choice of α as in (35) and simple calculation we obtain N = O(ϵ -1 ) when ϵ is small. Now we have∥y I+1 -y * I ∥ 2 = ∥y N I -y * I ∥ 2 ≤ 2 µ 1 - √ µα 4 N 2ϵ 2 αL 2 0 + ϵ 2 64L 2 0 ≤ ϵ 2 16L 2 0 , which yields ∥y I+1 -y * I+1 ∥ ≤ ∥y I+1 -y * I ∥ + ∥y * I -y * I+1 ∥ ≤ ϵ 4L 0 + ηl g,1 µ ≤ ϵ 2L 0 ,where we choose 1 -β to be small (see (54) for details) such that η is small enough to make above inequality holds. Repeating the same process yields the result.D. =: ϑ.Proof of Lemma D.5. We will first show the following result by induction, i.e., for any t ≥ 0, the averaged sequence {ŷ t } satisfiesEQUATION)For t = 0, by Lemma D.2 we have∥ŷ 0 -y * 0 ∥ = ∥y init T0 -y * 0 ∥ ≤ µα 32 ϵ L 0 = 1 -β 32 ϵ L 0 ≤ ϵ L 0 ,thus the base case holds. Now suppose (38) holds for some t ≥ 0, then for time step t + 1 we have∥ŷ t+1 -y * t+1 ∥ = ∥(1 -τ )(ŷ t -y * t+1 ) + τ (y t+1 -y * t+1 )∥ = ∥(1 -τ )(ŷ t -y * t ) + (1 -τ )(y * t -y * t+1 ) + τ (y t+1 -y * t+1 )∥ ≤ (1 -τ )∥ŷ t -y * t ∥ + (1 -τ )∥y * t -y * t+1 ∥ + τ ∥y t+1 -y * t+1 ∥ ≤ (1 -τ ) (1 -τ )ηl g,1 τ µ + ϵ L 0 + (1 -τ )ηl g,1 µ + τ ϵ L 0 ≤ (1 -τ )ηl g,1 τ µ + ϵ L 0 ,where we use induction hypothesis in the second inequality. Therefore, we have that (38) holds for any t ≥ 0. Also, as a consequence, for any t ≥ 0 we have∥ŷ t+1 -ŷt ∥ = ∥τ (y t+1 -ŷt )∥ ≤ τ ∥y t+1 -y * t+1 ∥ + τ ∥y * t+1 -y * t ∥ + τ ∥y * t -ŷt ∥ ≤ τ ϵ L 0 + ηl g,1 µ + (1 -τ )ηl g,1 τ µ + ϵ L 0 = τ ηl g,1 τ µ + 2ϵ L 0 .Now we plug in the definition of α, β, τ, η as in (35) to obtainEQUATIOND. 1 T T -1 t=0 E∥ϵ t ∥ ≤ σ T (1 -β) + 1 -β σ + L0 √ 1 -β 2(η 2 + ϑ 2 ) S + L1 2(η 2 + ϑ 2 ) S(1 -β) 1 T T -1 t=0 E∥∇Φ(x t )∥.Proof of Lemma D. By definition of ϵ t , εt and εt , we have the following recursion for any t ≥ 0:EQUATIONThen we apply (39) recursively to obtainϵ t = β t ϵ 0 + β t i=1 β t-i εi + (1 -β) t i=1 β t-i εi ,which by triangle inequality and total expectation givesEQUATIONBounding Err 1 . By definition of ϵ 0 and Lemma B.4, along with Jensen's inequality, we haveEQUATIONBounding Err 2 . We apply Lemma B.4 and follow the similar procedure as in [35, Lemma D.9] to obtainEQUATIONwhere we use Jensen's inequality for the first step.Bounding Err 3 . We will first use induction to show that for 0 ≤ i ≤ t + 1, the following inequality holds:EQUATIONThen it's easy to check that by setting i = t + 1 we can obtain the bound. When i = 0, (43) holds obviously sinceE t j=1 β t-j εj ≤ E   t j=1 β t-j εj 2   = E t j=1 β t-j εj .Hence the base case stands. Now suppose (43) holds for some i ≥ 0, and we aim to show that (43) holds for i + 1. In fact, we haveEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONEQUATIONwhere (45) follows by law of total expectation, (46) follows by Jensen's inequality, (47) uses the fact that εj for j < t -i are F t-i-1 -measurable, and are uncorrelated with εt-i ; for (48) we use Lemmas B.6 and D.5 to deriveE t-i-1 [∥ε t-i ∥ 2 ] = E t-i-1 ∥ ∇f (x t-i , ŷt-i ; ξt-i ) -∇f (x t-i-1 , ŷt-i-1 ; ξt-i ) -E t-i [ ∇f (x t-i , ŷt-i ; ξt-i )] + E t-i [ ∇f (x t-i-1 , ŷt-i-1 ; ξt-i )]∥ 2 ≤ E t-i-1 ∥ ∇f (x t-i , ŷt-i ; ξt-i ) -∇f (x t-i-1 , ŷt-i-1 ; ξt-i )∥ 2 ≤ 2E t-i-1 ∥ ∇f (x t-i , ŷt-i ; ξt-i ) -∇f (x t-i , ŷt-i-1 ; ξt-i )∥ 2 + 2E t-i-1 ∥ ∇f (x t-i , ŷt-i-1 ; ξt-i ) -∇f (x t-i-1 , ŷt-i-1 ; ξt-i )∥ 2 ≤ 2 S ( L2 0 + L2 1 ∥∇Φ(x t-i )∥ 2 )(∥ŷ t-i -ŷt-i-1 ∥ 2 + ∥x t-i -x t-i-1 ∥ 2 ) = 2 S ( L2 0 + L2 1 ∥∇Φ(x t-i )∥ 2 )(η 2 + ϑ 2 ).And (49) follows from the fact that x t-i is F t-i-1 -measurable, (50) uses√ a + b ≤ √ a + √ b for a, b ≥ 0.Hence the induction proof is completed. We set i = t + 1 to obtainEQUATIONFinal Bound. Combining ( 41), ( 42) and ( 51) yieldsE∥ϵ t ∥ ≤ β t σ + 1 -β σ + L0 √ 1 -β 2(η 2 + ϑ 2 ) S + 2(η 2 + ϑ 2 ) S L1 t i=0 β t-i E∥∇Φ(x i )∥.Taking summation and dividing 1/T on both sides gives the final result1 T T -1 t=0 E∥ϵ t ∥ ≤ σ T (1 -β) + 1 -β σ + L0 √ 1 -β 2(η 2 + ϑ 2 ) S + L1 2(η 2 + ϑ 2 ) S(1 -β) 1 T T -1 t=0 E∥∇Φ(x t )∥.E Proof of Theorem 4.1Before starting the proof of main results, i.e., Theorem 4.1, we first need the following lemma.Lemma E.1. Suppose that Assumptions 3.1 to 3.4 hold. For any η satisfyingη ≤ 1 2(1 + l 2 g,1 /µ 2 )(L 2 x,1 + L 2 y,1 ),it holds that 1 - 1 2 ηL 1 -2L 1 ∥ŷ t -y * t ∥ 1 T T -1 t=0 E∥∇Φ(x t )∥ ≤ Φ(x 0 ) -Φ(x T ) T η + 2 T T -1 t=0 E∥ϵ t ∥ + 2l g,1 l f,0 µ 1 - µ l g,1 Q + 2L 0 T T -1 t=0 ∥ŷ t -y * t ∥ + 1 2 ηL 0 .Proof of Lemma E.1. Define h t = m t -∇Φ(x t ). Then we apply Lemma B.3 to obtainEQUATIONwhere for the last two lines we use Cauchy-Schwarz inequality and ∥h t ∥ = ∥∇Φ(x t ) + h t ∥ ≥ ∥∇Φ(x t )∥ -∥h t ∥. Now expanding h t by triangle inequality, we have∥h t ∥ = ∥m t -∇Φ(x t )∥ ≤ ∥m t -E t [ ∇f (x t , ŷt ; ξt )]∥ + ∥E t [ ∇f (x t , ŷt ; ξt )] -∇f (x t , ŷt )∥ + ∥ ∇f (x t , ŷt ) -∇Φ(x t )∥ ≤ ∥ϵ t ∥ + l g,1 l f,0 µ 1 - µ l g,1 Q + (L 0 + L 1 ∥∇Φ(x t )∥)∥ŷ t -y * t ∥,Proof of Theorem E.2. By Lemmas D.6 and E.1, we have1 - 1 2 ηL 1 -2L 1 ∥ŷ t -y * t ∥ 1 T T -1 t=0 E∥∇Φ(x t )∥ ≤ Φ(x 0 ) -Φ(x T ) T η + 2 T T -1 t=0 E∥ϵ t ∥ + 2l g,1 l f,0 µ 1 - µ l g,1 Q + 2L 0 T T -1 t=0 ∥ŷ t -y * t ∥ + 1 2 ηL 0 ≤ d 0 T η + 2l g,1 l f,0 µ 1 - µ l g,1 Q + 2L 0 T T -1 t=0 ∥ŷ t -y * t ∥ + 1 2 ηL 0 + 2σ T (1 -β) + 2 1 -β σ + 2 L0 √ 1 -β 2(η 2 + ϑ 2 ) S + 2 L1 2(η 2 + ϑ 2 ) S(1 -β) 1 T T -1 t=0 E∥∇Φ(x t )∥.Rearranging the above inequality givesEQUATION.Bounding (LHS). By Lemma D.5, we haveEQUATIONBounding (RHS). By choice of parameters, we haveEQUATIONCombining ( 59) and ( 60) finally yields 1 T T -1 t=0 E∥∇Φ(x t )∥ ≤ 20ϵ.∥ ∇f (x, y) -∇Φ(x)∥ ≤ ( L + L x,1 ∥∇Φ(x)∥)∥y -y * (x)∥,where constant L is defined asL := L x,0 + L x,1 l g,1 l f,0 µ + l g,1 µ (L y,0 + L y,1 l f,0 ) + l f,0 µl g,2 + l g,1 l g,2 µ 2 ≤ L 0 .Proof of Lemma B.5. Recall that the exact expressions of ∇f (x, y) and ∇Φ(x) are∇f (x, y) = ∇ x f (x, y) -∇ 2 xy g(x, y)[∇ 2 yy g(x, y)] -1 ∇ y f (x, y) and ∇Φ(x) = ∇ x f (x, y * (x)) -∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y * (x))] -1 ∇ y f (x, y * (x)). Then by Assumption 3.2 we have ∥ ∇f (x, y) -∇Φ(x)∥ ≤ ∥∇ x f (x, y) -∇ x f (x, y * (x))∥ + ∥∇ 2 xy g(x, y)[∇ 2 yy g(x, y)] -1 ∇ y f (x, y) -∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y * (x))] -1 ∇ y f (x, y * (x))∥ ≤ (L x,0 + L x,1 ∥∇ x f (x, y * (x))∥)∥y -y * (x)∥ + ∥∇ 2 xy g(x, y)[∇ 2 yy g(x, y)] -1 ∇ y f (x, y) -∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y)] -1 ∇ y f (x, y)∥ + ∥∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y)] -1 ∇ y f (x, y) -∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y * (x))] -1 ∇ y f (x, y)∥ + ∥∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y * (x))] -1 ∇ y f (x, y) -∇ 2 xy g(x, y * (x))[∇ 2 yy g(x, y * (x))] -1 ∇ y f (x, y * (x))∥ ≤ L x,0 + L x,1 l g,1 l f,0 µ + ∥∇Φ(x)∥ ∥y -y * (x)∥ + l f,0 µ l g,2 ∥y -y * (x)∥ + l f,0 l g,1 µ 2 l g,2 ∥y -y * (x)∥ + l g,1 µ (L y,0 + L y,1 ∥∇ y f (x, y * (x))∥)∥y -y * (x)∥ = L x,0 + L x,1 l g,1 l f,0 µ + l g,1 µ (L y,0 + L y,1 l f,0 ) + l f,0 µl g,2 + l g,1 l g,2 µ 2 + L x,1 ∥∇Φ(x)∥ ∥y -y * (x)∥.By definition of L we conclude the proof.F. ∈ R dx , Eξ∥ ∇f (x 1 , y; ξ) -∇f (x 2 , y; ξ)∥ 2 ≤ ( L2 0 + L2 1 ∥∇Φ(x 1 )∥ 2 )∥x 1 -x 2 ∥ 2 .(ii) For any fixed x ∈ R dx and any y 1 , y 2 ∈ R dy , Eξ∥ ∇f (x, y 1 ; ξ) -∇f (x, y2 ; ξ)∥ 2 ≤ ( L2 0 + L2 1 ∥∇Φ(x 1 )∥ 2 )∥x 1 -x 2 ∥ 2 .In the above expressions, we define L0 and L1 asL0 = 4 L x,0 + L x,1 l g,1 l f,0 µ + L x,0 + L x,1 l g,1 l f,0 µ ∥y 1 -y * 1 ∥ 2 + 6Q 2µl g,1 -µ 2 l 2 g,1 (L y,0 + L y,1 l f,0 ) 2 + l 2 f,0 l 2 g,2 + l 2 f,0 l 2 g,1 l 2 g,2 K 2 (l g,1 -µ) 2 1/2 , L1 = 2L x,1 (1 + L x,1 ∥y 1 -y * 1 ∥).Proof of Lemma B.6. We show statement (i) of the lemma, and (ii) follows by similar arguments. For any fixed y ∈ R dy and any x 1 , x 2 ∈ R dx , by definition of ∇f (x, y; ξ) we have∥ ∇f (x1, y; ξ) -∇f (x2, y; ξ)∥ 2 ≤ 2∥∇xF (x1, y; ξ) -∇xF (x2, y; ξ)∥ 2 + 2 ∇ 2 xy G(x1, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) ∇yF (x1, y; ξ) -∇ 2 xy G(x2, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) ∇yF (x2, y; ξ) 2 ≤ 2(Lx,0 + Lx,1∥∇xf (x1, y)∥) 2 ∥x1 -x2∥ 2 + 2 ∇ 2 xy G(x1, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) ∇yF (x1, y; ξ) -∇ 2 xy G(x2, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) ∇yF (x2, y; ξ) 2 .For the second term above, we have∇ 2 xy G(x1, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) ∇yF (x1, y; ξ) -∇ 2 xy G(x2, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) ∇yF (x2, y; ξ) 2 ≤ 3l 2 g,1 Q 2 l 2 g,1 1 - µ lg,1 2q ∥∇yF (x1, y; ξ) -∇yF (x2, y; ξ)∥ 2 + 3l 2 f,0 Q 2 l 2 g,1 1 - µ lg,1 2q ∥∇ 2 xy G(x1, y; ζ (0) ) -∇ 2 xy G(x2, y; ζ (0) )∥ 2 + 3l 2 g,1 l 2 f,0 Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) - Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) 2 ≤ 3Q 2 1 - µ lg,1 2q (Ly,0 + Ly,1∥∇yf (x1, y)∥) 2 ∥x1 -x2∥ 2 + 3 l 2 f,0 Q 2 l 2 g,1 1 - µ lg,1 2q l 2 g,2 ∥x1 -x2∥ 2 + 3l 2 f,0 Q 2 q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) - q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) 2 .Then we take expectation with respect to q and obtainEq ∇ 2 xy G(x1, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) ∇yF (x1, y; ξ) -∇ 2 xy G(x2, y; ζ (0) ) Q lg,1 q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) ∇yF (x2, y; ξ) 2 ≤ 3Q 2 (Ly,0 + Ly,1∥∇yf (x1, y)∥) 2 ∥x1 -x2∥ 2 + 3 l 2 f,0 Q 2 l 2 g,1 l 2 g,2 ∥x1 -x2∥ 2 Eq 1 - µ lg,1 2q + 3l 2 f,0 Q 2 Eq q i=1 I - 1 lg,1 ∇ 2 yy G(x1, y; ζ (i) ) - q i=1 I - 1 lg,1 ∇ 2 yy G(x2, y; ζ (i) ) 2 ≤ 3Q 2 (Ly,0 + Ly,1∥∇yf (x1, y)∥) 2 ∥x1 -x2∥ 2 + 3 l 2 f,0 Q 2 l 2 g,1 l 2 g,2 ∥x1 -x2∥ 2 • l 2 g,1 Q(2µlg,1 -µ 2 ) + 3l 2 f,0 Q 2 • l 2 g,1 l 2 g,2 Q (lg,1 -µ) 2 (2µlg,1 -µ 2 ) ∥x1 -x2∥ 2 ≤ 3Q 2µlg,1 -µ 2 l 2 g,1 (Ly,0 + Ly,1l f,0 ) 2 ∥x1 -x2∥ 2 + l 2 f,0 l 2 g,2 ∥x1 -x2∥ 2 + 3l 2 f,0 l 2 g,1 l 2 g,2 Q 3 (lg,1 -µ) 2 (2µlg,1 -µ 2 ) ∥x1 -x2∥ 2 = 3Q 2µlg,1 -µ 2 l 2 g,1 (Ly,0 + Ly,1l f,0 ) 2 + l 2 f,0 l 2 g,2 + l 2 f,0 l 2 g,1 l 2 g,2 Q 2 (lg,1 -µ) 2 ∥x1 -x2∥ 2Finally, taking expectation on both sides yields Eξ∥ ∇f (x 1 , y; ξ) -∇f (x 2 , y; ξ)∥ 2 ≤ 2(L x,0 + L x,1 ∥∇ x f (x 1 , y)∥) 2 ∥x 1 -x 2 ∥ 2 + 6Q 2µl g,1 -µ 2 l 2 g,1 (L y,0 + L y,1 l f,0 ) 2 + l 2 f,0 l 2 g,2 + l 2 f,0 l 2 g,1 l 2 g,2 Q 2 (l g,1 -µ) 2 ∥x 1 -x 2 ∥ 2 .Since for any y ∈ R dy , we have ∥∇ x f (x 1 , y) -∇ x f (x 1 , y * 1 )∥ ≤ (L x,0 + L x,1 ∥∇ x f (x 1 , y * 1 )∥)∥y 1 -y * 1 ∥ ≤ L x,0 + L x,1 l g,1 l f,0 µ + ∥∇Φ(x 1 )∥ ∥y 1 -y * 1 ∥ = L x,0 + L x,1 l g,1 l f,0 µ + L x,1 ∥∇Φ(x 1 )∥ ∥y 1 -y * 1 ∥,• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: This paper presents work whose goal is to advance the field of Machine Learning from algorithmic and theoretical aspects. We do not see any direct paths to negative societal impacts. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).The notation in Section 4.3.1 is independent of that in other sections, although there may be incidental overlaps in terminology.µ 2 (σ 2 f,1 + l 2 f,0 )(σ 2 g,2 + 2l 2 g,1 ) + σ 2 f,1 l 2 g,1 := σ2 .3 Here we define ∥v∥ Q := v ⊤ Qv for any vector v ∈ R d and PSD matrix Q ∈ R d×d .Note that Assumptions C.1 and C.2 are more concrete than that in Section 4.3.1. 


Decision: There was some disagreement among the reviewers about this submission. Based on re-reading the whole discussion (taking into account reviewer's levels of expertise) and my own evaluation of the paper, I recommend to accept.

Reviews: review: paper_summary: paper_summary: This paper considers a class of stochastic bilevel optimization problems where the upper-level function is nonconvex with potentially unbounded smoothness, and the lower-level problem is strongly convex. Their novel algorithms achieve an oracle complexity of $O(1/\epsilon^3)$ to find an $\epsilon$-stationary point.
strength_weakness: strength_weakness: 1. In the algorithm, for the update of the normalized stochastic gradient, it seems that it can be reformulated to tune the learning rate by dividing by the norm of the stochastic gradient. Therefore, the novelty and reasonableness of this step may be unclear.
2. In the experiments, it is common to use the same learning rate for baselines and the proposed algorithms for a fair comparison. However, in the experimental details, the authors provide the best learning rate pairs for different baselines, which seems unusual.
questions: questions: 1. Regarding the main challenge, the authors state why previous algorithms and analyses do not work. However, this explanation is not clear to me. Can you elaborate on how your work addresses the problem of hypergradient bias? Also, how does your potential function argument differ from previous work?
limitations: limitations: N/A
paper_summary: paper_summary: This paper investigates a class of stochastic bilevel optimization problems where the upper-level function is nonconvex with potentially unbounded smoothness, and the lower-level problem is strongly convex. To improve the convergence rate, the authors propose a new Accelerated Bilevel Optimization algorithm named AccBO. This algorithm updates the upper-level variable using normalized stochastic gradient descent and the lower-level variable using the stochastic Nesterov accelerated gradient descent algorithm with averaging. The proof shows that AccBO has an oracle complexity of $ \tilde{O}(1/\epsilon^3) $, relying on a novel lemma describing the dynamics of the stochastic Nesterov accelerated gradient descent algorithm under distribution drift. Experimental results demonstrate that AccBO significantly outperforms baseline algorithms in various tasks, achieving the predicted theoretical acceleration.
strength_weakness: strength_weakness: The convergence analysis for Option I of the algorithm is restricted to the lower-level problem being a one-dimensional quadratic function. This limitation reduces the generality and applicability of the theoretical results for more complex lower-level problems.
questions: questions: In reference [33], Assumption 2 states $|| \nabla_y f(x, y^*(x)) || \leq M$, while in this paper, Assumption 3.2 states$ || \nabla_y f(x, y) || \leq l_{f,0} $ . Can this assumption be relaxed to a less stringent form?
limitations: limitations: The convergence analysis for Option I of the AccBO algorithm assumes that the lower-level problem is a one-dimensional quadratic function.  This limits the algorithm's applicability to scenarios where the lower-level function deviates significantly from this simple form, especially in higher dimensions or more complex optimization landscapes.
paper_summary: paper_summary: The authors consider bilevel optimization problems under the unbounded smoothness assumption of the outer function. They propose AccBO, an AID-based method that uses Stochastic Nesterov Accelerated Gradient to approximate the lower-level solution and Neumann approximations to estimate the inverse Hessian-Vector product involved in the hypergradient expression. Then, a normalized STORM-like update is performed for the outer variable. They show a $\tilde{\mathcal{O}}(\epsilon^{-3})$ complexity of their algorithm. Finally, numerical experiments are provided on the deep AUC maximization and data hypercleaning problems.
strength_weakness: strength_weakness: ### Major

- **W1**: Some errors in the code may invalidate the experimental results. See the comment **C2** in the Code section for details.

- **W2**: AccBO comes with many hyperparameters, which can be a concern for practical applications.

- **W3**: l.219: "is nearly optimal in terms of the dependency of $\epsilon$". The authors should be careful when talking about optimality of an algorithm. In particular, the algorithm AccBO does not belong to the algorithm class considered in [1] which does not take into account the biasedness of the hypergradient estimate due to the approximation of the inverse Hessian-vector product and the solution of the inner optimization problem. To my knowledge, the only works considering lower bounds for bilevel problems are [2, 5].

- **W4**: The works [2, 3, 4] also use Nesterov acceleration for the inner optimization problem, but in the deterministic setting. They should be mentioned in the paper.

### Code
Even though it is a theoretical paper, I have several comments regarding the experiments code:

- **C1**: In the checklist, the authors indicate "The code and data are attached as a supplement with instructions for reproducibility." However, the code documentation is insufficient to run the experiments. The data are not present in the folders and there is no instruction indicating where to download them. As a consequence, when I launch the commmand `python main.py --methods accbo ` as indicated in the `README.md`, I get an error.

- **C2**: It seems that the implementation of the different methods do not correspond to the one described in their respective original paper, which might invalidate the empirical results:
    - For `TTSA`, `SUSTAIN` and `AccBO`, the estimation of the inverse Hessian-vector product is computed by the formula 
    $$
        v = \eta\sum_{q=-1}^Q\prod_{j=Q-q}^Q(I-\eta \nabla^2_{yy}G(x,y;\zeta^{(j)}))\nabla_y F(x,y;\xi)
    $$
    while in the original paper (and even the present paper for AccBO), the authors use the formula
    $$
        v = \eta\prod_{j=1}^{q(Q)}(I-\eta \nabla^2_{yy}G(x,y;\zeta^{(j)}))\nabla_y F(x,y;\xi)
    $$
    with $q(Q)\sim\mathcal{U}{\{1,\dots,Q-1\}}$.

    - For `SABA`, the SAGA-like variance reduction is absent in the code while being the core of the algorithm. Moreover, `SABA` uses constant step sizes.

    - For `VRBO`, the oracles in the case `step % self.args.update_interval == 0` should be computed with a larger batch-size (or even with full batches) than in the inner for-loop.

- **C3**: The code of `VRBO` for the datacleaning task is missing. The method to compute the hypergradient in `ma_soba.py` is missing as well for the datacleaning task (but not AUC experiment, thus I assume it is almost the same).

### Minor
- l.41: "transforms" -> "transformers"
- In `README.md` (both): `python main.py --medthods [algorithm] ` -> `python main.py --methods [algorithm] `

I am inclined to raise my score if my concerns, particularly regarding the experiments, are addressed.
questions: questions: N/A
limitations: limitations: N/A
paper_summary: paper_summary: This work develops the first algorithm that achieves the near-optimal oracle complexity $\widetilde{O}(\epsilon^{-3})$ (the number of evaluations to gradient or Hessian/Jacobian vector product) to achieve an $\epsilon$-stationary point of a stochastic bilevel optimization problem with unbounded smooth nonconvex upper-level function and strongly-convex lower-level function. The algorithm design is novel which updates the upper-level variable by normalized stochastic gradient descent with recursive momentum and updates the lower-level variable by the stochastic Nesterov accelerated gradient (SNAG) method.
strength_weakness: strength_weakness: The major weakness is Assumption 4.2 as mentioned in the question (1) below. Some points could be clarified as mentioned in other questions below.
questions: questions: (1) **My major concern** is that Assumption 4.2 is made on algorithm-generated variables not on problem setting, which cannot be verified. Is it possible to bound $\Delta_t$ in terms of $\|x_{t+1}-x_t\|$ to get rid of Assumption 4.2? If not, at least Assumption 4.2 should also be mentioned in the main result (Theorem 4.1). **I would like to raise my rating if all the assumptions become problem-related.**

(2) In line 48 "each realization of the stochastic oracle calls is Lipschitz with respect to its argument", do you mean $F(x,y,\xi)$ and $G(x,y,\xi)$ are Lipschitz continuous functions of $(x,y)$ for each $\xi$? You might express in formula to make it more clear. 

(3) You could also cite [1] as a paper on relaxed smoothness:  
[1] Convex and Non-convex Optimization Under Generalized Smoothness

(4) In line 127, do you mean "... is a good approximation of $\nabla \Phi(x)$ **if $y\approx y^*(x)$**"? 

(5) Is there an intuition why line 21 (averaging of y) is needed in Algorithm 2? Is it possible to implement $N=\mathcal{O}(1)$ steps of SNAG such that $\mathbb{E}\|y_{t+1}-y^*(x_t)\|\le \frac{1}{2}\mathbb{E}\|y_t-y^*(x_t)\|$ and remove averaging of y for every outer iteration t, and obtain convergence result purely in expectation instead of high probability? This yields the same oracle complexity per round since $I,N=\widetilde{\mathcal{O}}(\epsilon^{-1})$ are of the same order, and may get rid of Assumption 4.2 by bounding $\Delta_t$ in terms of $\|x_{t+1}-x_t\|$.

(6) You said your upper-level algorithm design is inspired from [17,49]. However, in line 23 (STORM variance reduction) of Algorithm 2, your second and third additive terms use the same samples $\overline{\xi}_t$ while [49] uses different samples. The algorithm design differs from [17] more. Such difference still exists even if Algorithm 2 is reduced to single-level. Why do you use such different design? 

(7) Also in line 23, does $\overline{\nabla}f(x_t,\hat{y}_t;\overline{\xi}_t)$ denote  the average of the expression right below Eq. (2) over samples $s=1,\ldots,S$? You could define more clearly. 

(8) At the beginning of Section 5 experiments, $ph(w; x)I_{[y=-1]}$ appears in your AIC maximization. Do you mean $rh(w; x)I_{[c=-1]}$? If not, what do $p$ and $y$ mean?
limitations: limitations: The conclusion section mentions the limitation that the convergence analysis for the Option I of Algorithm 2 relies on the lower-level problem being a one-dimensional quadratic function. Actually I think this is fine since Option II works for more general strongly convex lower-level function. The major limitation from my perspective in mentioned in my above question (1) about the undesired Assumption 4.2. 
The authors said in the checklist that they do not see any negative societal impacts in their theoretical work, which I agree.

Comments: title: Thank you!
Thanks for your reply!
Looking forward to post-rebuttal feedback!
Thank you for your feedback!
If you haven't received the anonymized code link, please reach out to AC for details.
Thank you!
Just raised rating
Thank you for your feedback!
Reviewer 2WtS agrees with the authors and will keep rating 5.
comment: Dear Reviewer wLtJ,

We are glad that our answer addressed your concerns. Thank you for reviewing our paper.

Best,

Authors
Thank you for your answer. I will raise my score accordingly. 

Best.
Thank you very much for your reply, below we provide further explanations to your additional concerns.

- First, we would like to clarify that the learning rate should be $\eta$ instead of $\frac{\eta}{\\|m_t\\|}$ when considering normalized stochastic gradient descent. For example, see Algorithm 1 of [(Cutkosky and Mehta, 2020)](https://arxiv.org/pdf/2002.03305) (Input: $\dots$, learning rate $\eta$, $\dots$), Algorithm 1 of [(Jin et al., 2021)](https://arxiv.org/pdf/2110.12459) (Input: $\dots$, learning rate $\gamma$, $\dots$), and Algorithms 1 and 2 of [(Liu et al., 2023b)](https://arxiv.org/pdf/2302.06763) (middle of page 6, "the step size $\eta$ is chosen by balancing every other term to get the right convergence rate") for details. We follow the same practice and refer the learning rate of the update $x_{t+1} = x_t - \eta\frac{m_t}{\\|m_t\\|}$ (i.e., line 24 of Algorithm 2) as $\eta$, rather than $\frac{\eta}{\\|m_t\\|}$. Additionally, our choice of learning rate $\eta$ is independent of the magnitude of $\\|m_t\\|$, please see Theorem 4.1 and Eq. (55) in Appendix E on page 31 for more details.

- Second, please note that we study the same problem setting as [(Hao et al., 2024)](https://arxiv.org/pdf/2401.09587.pdf), where the upper-level problem is unbounded smooth (i.e., $(L_{x,0}, L_{x,1}, L_{y,0}, L_{y,1})$-smooth as illustrated in Assumption 3.1). Both our proposed Algorithm 2 (i.e., AccBO) and Algorithm 1 in [(Hao et al., 2024)](https://arxiv.org/pdf/2401.09587.pdf) use normalization and choose a fixed learning rate $\eta$. We mention allowing a larger learning rate because, for any given small $\epsilon>0$ (i.e., the target gradient norm $\epsilon$), AccBO's learning rate is $\eta=\widetilde{\Theta}(\epsilon^2)$, while that of [(Hao et al., 2024)](https://arxiv.org/pdf/2401.09587.pdf) is $\eta=\widetilde{\Theta}(\epsilon^3)$. In other words, we say that we choose a larger learning rate $\eta$ specifically for comparison with [(Hao et al., 2024)](https://arxiv.org/pdf/2401.09587.pdf). The large learning rate is crucial for achieving acceleration: it improves the complexity from $\widetilde{O}(\epsilon^{-4})$ as in [(Hao et al., 2024)](https://arxiv.org/pdf/2401.09587.pdf) to $\widetilde{O}(\epsilon^{-3})$ as in this paper.

Please let us know if you have further questions. Thanks!

Best,

Authors
I appreciate the authors' reply. I understand the rationale for a larger learning rate from the convergence rate perspective. However, I still find the gradient normalization on line 24 in Algorithm 2 somewhat problematic in claiming this larger learning rate. From my perspective, when the norm of $m_t$ is larger, it would allow for a higher learning rate. Please correct me if I'm mistaken.
Dear Reviewer wLtJ,

Thank you for reviewing our paper. We have carefully addressed your concerns regarding the novelty and reasonableness of the learning rate choice, the fairness of comparisons between our proposed algorithm and other baselines, and how our work tackles the challenge of hypergradient bias.

Please let us know if our responses address your concerns accurately. If so, we kindly ask you to consider raising the rating of our work. We appreciate your time and efforts and are open to discussing any further questions you may have.

Best,

Authors
Dear Reviewer o6Pr,

We are glad that our responses addressed your concerns. Thank you for reviewing our paper.

Best,

Authors
I thank the authors for fixing their code and modified my score accordingly.
Dear Reviewer o6Pr,

Thank you for reviewing our paper. We have carefully addressed your concerns regarding the optimality of our algorithm, the implementation of the inverse Hessian-vector product and SABA, as well as the learning rate and batch size choices for the baselines SABA and VRBO.

We have updated our code and rerun the experiments. The new results, along with the hyperparameter settings, are included in the global rebuttal PDF file. **Please reach out to the AC for the anonymized code link if needed**, as we were notified that "If you were asked by the reviewers to provide code, please send an anonymized link to the AC in a separate comment."

Please let us know if our responses address your concerns accurately. We appreciate your time and efforts and are open to discussing any further questions you may have.

Best regards,

Authors
Dear Reviewer 2WtS,

Thank you for raising the rating. 

Best regards,

Authors
Dear Authors,

Reviewer 2WtS just raised rating from 5 to 6, and soundness from 2 to 3. 
Best regards,

Reviewer 2WtS
Dear Reviewer 2WtS,

We are glad that our responses addressed your major concerns. Thank you for reviewing our paper. You mentioned earlier that "I would like to raise my rating if all the assumptions become problem-related." We have now verified that Assumption 4.2 holds in the bilevel optimization setting (under Assumptions 3.1-3.4) and it is indeed problem-related. Could you please consider raising your score?

Please let us know if you have any further questions or concerns.

Best regards,

Authors
Reviewer 2WtS agrees with the authors and will keep rating 5.