Efficient Curvature-Aware Hypergradient Approximation for Bilevel Optimization
Authors: Youran Dong
Junfeng Yang
Wei Yao
Jin Zhang
Venue: ICML 2025 poster (2025)
Keywords: bilevel optimization
hypergradient
inexact Newton method
computational complexity
Primary area: optimization
TL;DR: We propose a new algorithmic framework for bilevel optimization that incorporates curvature information into hypergradient approximations, achieving improved computational complexity and demonstrating significant practical performance.

Abstract: Bilevel optimization is a powerful tool for many machine learning problems, such as hyperparameter optimization and meta-learning. Estimating hypergradients (also known as implicit gradients) is crucial for developing gradient-based methods for bilevel optimization. In this work, we propose a computationally efficient technique for incorporating curvature information into the approximation of hypergradients and present a novel algorithmic framework based on the resulting enhanced hypergradient computation. We provide convergence rate guarantees for the proposed framework in both deterministic and stochastic scenarios, particularly showing improved computational complexity over popular gradient-based methods in the deterministic setting. This improvement in complexity arises from a careful exploitation of the hypergradient structure and the inexact Newton method. In addition to the theoretical speedup, numerical experiments demonstrate the significant practical performance benefits of incorporating curvature information.

Efficient Curvature-Aware Hypergradient Approximation for Bilevel Optimization
Youran Dong'
Junfeng Yang'
Wei Yao 23 Jin Zhang 324
Abstract
architecture search (Liu et al., 2018; Liang et al., 2019).
Bilevel optimization is a powerful tool for many machine learning problems, such hyperparam- eter optimization and meta-learning, Estimating hypergradients (also known as implicit gradients) is crucial for developing gradient- gradient-based methods for bilevel optimization. this work, we propose computationally efficient technique for incor- porating curvature information into the approx- imation of hypergradients and present novel algorithmic framework based on the resulting enhanced hypergradient computation. We pro- vide convergence rate guarantees for the proposed framework both deterministic and stochastic scenarios, particularly showing improved compu- tational complexity over popular gradient-based methods in the deterministic setting. This im- provement in complexity arises from careful exploitation of the hypergradient structure and the inexact Newton method. In addition to the theoret- ical speedup, numerical experiments demonstrate the significant practical performance benefits of incorporating curvature information.
Bilevel optimization amounts to solving an optimization problem with a constraint defined by another optimization problem. this work, we focus on the following bilevel optimization problem:
min (x))= == f( (x))
JE Rm
(1)
8.1. y (x) arg min ( y);
yER"
where the upper- and lower-level objective functions f and are real-valued functions defined on Ih x RR. We assume that g is strongly convex w.r. the lower-level variable y. which guarantees the uniqueness of the lower-level solution. However, solving this bilevel optimization problem remains challenging, as yy (x:) can typically only be approximated through iterative schemes (Pedregosa, 2016; Grazzi et al., 2020; Dagréou et al., 2022: Chu et al., 2024).
In large-scale scenarios, gradient based bilevel optimization methods have gained popularity due their effectiveness (Ghadimi & Wang, 2018: Ji et al., 2021). When solving the bilevel problem (1), it is essential to estimate the hypergra- dient (also known as implicit gradient), which represents the gradient of (xr. By the implicit function theorem, the hypergradient can be expressed as:
1. Introduction
Bilevel optimization has been widely applied to solve enor- mous machine learning problems, such as hyperparameter optimization (Pedregosa, 2016; Franceschi et al.. 2018), meta-learning (Franceschi et al., 2018; Rajeswaran et al., 2019; Ji al., 2020), adversarial training (Bishop et al., 2020; Wang et al., 2021; Zhang al., 2022), reinforcement learning (Yang al., 2019; Liu et al., 2021), and neural
(z), (2) Based the structure in (2), estimating hypergradients requires solving lower-level problems and computing the Hessian inverse-vector products. Several studies have emerged to effectively address this challenge, such as those by Ghadimi & Wang (2018); Ji et al. (2021; 2022); Arbel & Mairal (2022a); Dagréou et al. (2022). However, most existing studies primarily focus on estimating the Hessian inverse- vector products, including Neumann series approximations (Chen et al., 2021; Hong et al., 2023), conjugate gradient descent (Pedregosa, 2016; Ji et al., 2021; Arbel & Mairal, 2022a), gradient descent (Arbel & Mairal, 2022a; Dagréou et al., 2022), and subspace techniques (Yang et al., 2025). In approximating y (a="), Pedregosa (2016) employs the L-BEGS method to solve the lower-level problem up to
= Vi
(x)) -V))))))VV((aa)
where u² ux(x) (=agl
'Shool of Mathematies, Nanjing University, Nanjing, China NNational Center for Applied Mathematics Shenzhen, South- ern University of Science and Technology, Shenzhen, China "Department of Mathematics, Southern University of Science and Technology, Shenzhen, China *Detection Institute for Advanced Technology Longhua-Shenan Shenzhen (DIATLHSZ), Shenzhen, China. Correspondence to: Jin Zhang <zhangj9 @ sustech. anstech.edu.cn>..
Proceedings of the ndnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Efficient Curvature-Aware Hypergradient Approximation
specified tolerance. Recently, two commonly used ap- proaches have emerged. The first approach involves per- forming single (stochastic) gradient descent step (Ji et al., 2022; Dagréou et al., 2022; Hong et al.. 2023). The second approach entails executing multiple (stochastic) gradient descent steps (Ghadimi & Wang. 2018; Chen et al., 2021; Ji et al., 2021; 2022; Arbel & Mairal, 2022a).
deterministic setting.
We conduct numerical experiments to compare the pro- posed algorithms with popular gradient- based methods, demonstrating the significant practical performance benefits of incorporating curvature information.
1.2. Notation
From the above, is clear that solving lower-level problems and computing Hessian inverse- products typically treated as separate tasks. However, it is important to note that the Hessian V229 in Hessian inverse-vector products originates from the lower-level objective function. Hence, hypergradient approximation has an intrinsic structure in which solving lower-level problems and computing Hes- sian inverse-vector products share the same Hessian. This intuitive structure just deseribed has been demonstrated in Ramzi et al. (2022), which introduces SHINE, novel method that solves the lower-level problem using quasi- Newton (qN) methods and employs the associated qÑ ma- trices, along with refinement strategies, to compute Hessian inverse-vector products.
We refer the optimal value of problem (1) as क*. The gradient of g w.r.f the variables and y are denoted by V "1 y) and V2g(a, y), respectively. The Jacobian ma- trix of Vig and the Hessian matrix of g w.r.f y denoted by V22 y) and V229(1, y). respectively. Unless other- wise specified, the notation || - denotes the Ez norm for vectors and the Frobenius norm for matrices. Furthermore, the operator norm of a matrix Z is denoted by /
2. Curvature-Aware Bilevel Optimization Framework
2.1. Hypergradient Approximation
However, the theoretical analysis of SHINE is hindered by the mixing of quasi-Newton recursion schemes and com- plex refinement strategies. As a result, Ramzi et al. (2022) focuses on the asymptotic convergence analysis of hyper- gradient approximation in deterministic setting. The con- vergence rate and computational complexity are lacking. Therefore, the intuitive benefit of the intrinsic hypergradient structure has not been fully realized in Ramzi al. (2022) or in the existing literature.
The hypergradient given by (2) is intractable in practice because it requires knowing y (z). To address this, the seminal work (Ghadimi & Wang, 2018) approximates the hypergradient by
(o(x,)) = Vis(a,v) - Vigg(,, v)u (I,
where y is an approximation of (²)) and w wx == V220(z, (z. v) -V)- (2,( y). Under appropriate hypotheses,
Lemma 2.2 of Ghadimi & Wang (2018) provides an error
bound:
1.1. Contributions
( v) v = Cllv"()) - yll.
(3)
This paper aims to explore and exploit the benefits of lever- aging the hypergradient structure. Our contributions are summarized below.
where C a constant. Observe that V22g in u² u'(x, y) is the Hessian of g. Hypergradient approximation has a struc- ture in which solving lower-level problems and computing Hessian inverse-vector products share the same Hessian.
We propose simple Newton-based framework, NBO, for bilevel optimization that integrates solving the lower-level problems with computing Hessian inverse- vector products. This framework is built on a new curvature-aware hypergradient approximation, which utilizes the hypergradient structure and inexact Newton methods. When each subproblem in NBO is approx- imately solved using single gradient descent step initialized at zero, NBO simplifies to the well-known single-loop algorithm framework presented in Dagréou et al. (2022).
To exploit the hypergradient structure, it is natural to utilize the Hessian V22g. which provides curvature information, in solving the lower-level problem. The canonical second- order optimization scheme for this purpose is Newton's 's method. However, Newton's 's method is impractical for large- scale problems; thus, we instead consider inexact Newton methods and propose approximation, consisting of two steps:
technique for hypergradient
(i) Given æ and y, compute an inexact solution (v, u) of the linear system:
We establish the convergence rate and computational complexity for two specific instances of NBO in both deterministic and stochastic scenarios. In particular, we demonstrate improved computational complexity compared to popular gradient-based methods in the
V22g(2, ()) u) (V2g(1, v). Vef(a, (4)
ii) Compute the approximated hypergradient:
d, := = VIf(z,y v) Vi2g(r, y v)u. (5)
2

---

Efficient Curvature-Aware Hypergradient Approximation
One can easily observe that v is inexact ap- proximation of the Newton direction o (x, [222((,,)) y)] y), just as is an inexact approx- imation of u² (z,
Algorithm 1 Newton-based framework for Bilevel Opti- mization (NBO)
1: Input: Initialize y², u. xx;; step size ar..
2: for k = 0,1. K - 1o
Compute inexact Newton directions v², w by mini-
The motivation is twofold: (1) Denote yt == y [V22(1,)))
mizing the quadratic functions in (8) and (10).
1) as a single Newton step. Since (2. y (a)) = O, and assuming that V22g(z,) -) is gg Lipschitz continuous, by Lemma 1.2.4 in Nesterov (2018),
4: Update yh and w
g y² - vh: w+t u² w
(22((, v) VVg(2, y² (x)
Compute the approximated hypergradient:
V2g -²) V2g(z, y) (220(1, ")l("" - b)
Vif w8" Vi2g(=", v/
2.2 2/
Ily (x) y11?
(6)
Update *+1;
where we denote y" := yd (a) to save space when there is no ambiguity. Hence, we obtain:
x k+1 = x² and
7: end for
Cl ( (x)
y/1..
24
slightly different from (5), we use y and u² instead of t and ++1 in de to facilitate potential parallel computation. We can also update these variables in an alternating manner.
This inequality shows that a single classical Newton step can accelerate the hypergradient estimation, leading to quadratic decay. (2) Observe that u (x) and u² vI)
share the same Hessian inverse.
Since the subproblems are strongly convex quadratic pro- gramming problems, the proposed NBO has the potential to adapt to stochastic scenarios and be implemented using various techniques from both deterministic and stochastic optimization, similar to the frameworks outlined in Arbel & Mairal (2022a); Dagréou et al. (2022). In the following, we study two specific examples.
2.2. Description of the Algorithmic Framework
Building the new hypergradient approximation, we in- troduce our algorithmic framework for solving the bilevel optimization problem in (1).
First, motivated by Arbel & Mairal (2022a); Dagréou et al. (2022), solving the linear system (4) at (a*, y² is reformu- lated as minimizing the following two quadratic problems, which share the same Hessian Hş V22g(" ,y*):
2.3. Examples
the NBO-GD algorithm. In the first example, we compute inexact Newton directions in NBO using gradient descent (GD) steps with pre-defined number of iterations 1'+ The resulting algorithm, referred to as NBO-GD, is detailed in Algorithm 2, with its subrou- tine GD(a*, ,v,, u²: outlined in Algorithm 3, where
Deterministic Setting:
min
V229(> y²) ju, e) (V2g(0 ,,,,), 2
(8)
min- a 2 (222(<" yM Ju, u) (V2 w/t) ),u)) u).
dt = V20(2 )²) dt = (x )/) V2f(=*.y*)
Second, using warm-start procedure to initialize the solver for u, write w² w, where w² is the current iterate of u. Then w minimizes the following quadratic function:
Remark 2.1. (i) When T = O, that is, when each subprob- lem in NBO is approximately solved using a single GD step initialized at zero, NBO-GD reduces to the single-loop algorithm framework in Dagréou al. (2022).
|Iwlllr.. (V22000 )a (a (."),,"") w), (10)
where ||wl| := (Hyw, w). Since computing the exact minimizers may be computationally demanding, instead seek inexact solutions. We denote these inexact solutions by w and & for the quadratic problems (8) and (10).
(ii) The subproblems in NBO can also be approximately solved using the conjugate gradient method (Nocedal & Wright, 2006), gradient descent with a Chebyshev step size (Young, 1953), or subspace techniques (Yang et al., 2025).
Third, we compute the approximated hypergradient and update r using an inexact hypergradient descent step. The overall procedure is summarized in Algorithm 1. Here,
(iii) It is worth noting that for a fixed k, the subroutine GD( (r.k yt w²;T) does not involve gradient computations
3

---

Efficient Curvature-Aware Hypergradient Approximation
Algorithm 2 NBO-GD
Algorithm 4 NSBO-SGD
1: Input: Initialize y²: u², _ number of iterations K, T;
1: Input: Initialize y", u", x:: number of iterations K, T;
step size OFk-
step size Ork-
2: for k k= 0, K 1 do
2: for k = 0, 0........ K do
3: Compute inexact Newton directions:
Compute the inexact subsampled Newton directions:
vh th = GD(x "0" T),
w w = SGD(2 y ;T)
4: Update y + and uht
4: Update y and wk+1;
g y² ut; w++1 wh w
g = y - ut: whtl w² -w
5: Compute the approximated hypergradient:
Update *+1,
ViS "("(a",,"))) Vi2g(>" y² Jua.
k+1 x =
6:
Update xh+1
where D is unbiased estimator of de in (15).
6: end for
x k+1 = x² - axd
Algorithm SGDD(* v², , T)
7: end for
1: Initialize --1 = 0) and w 1,, = O: step size Tk-
2: Sample batches B By. B and compute D D using (13) and (14), respectively.
Algorithm 3 GD(x", y², TT
1: Initialize e -11 = and w 1,b = 0; step size "/k-
3: for 1, -,, T do
2: for -,, T do
4: Sample batch BiAh and update
3: Update
[utti, t
[ur+ u 1,kk
(11)
=11 - H "IO wh + 7a DD' D.l, (16)
=[1 - 7/k V 122 y [utt, wot + nild (e....). d.gl.
where H := (12 n²: Bl.r.). end for
4: end for
5:
Return u = o.. w = w.rkk
6: Return enth Tik wk = T.k
but relies solely on Ressian-vectoon product computations, sharing the same Hessian (22 .
Remark 2.2. Inspired by recent works (Dagréou et al., 2022; Chu al., 2024), broad class of stochastic gradient esti- mation techniques, such as SAGA (Defazio et al., 2014), STORM (Cutkosky & Orabona, 2019), and PAGE (Li et al., 2021), can be incorporated into the NSBO framework.
Stochastic Setting: the NSBO-SGD algorithm. We con- sider our framework in the stochastic setting, where
f(z v) Ee [z E),,, g(z, y) Ez[G(z, yi ().. (12)
3. Convergence Analysis
For simplicity, we focus on basic stochastic gradient descent (SGD) algorithms. In Algorithm 4, NSBO-SGD, an adapta- tion of SGD to the newton-based framework for Stochastic Bilevel Optimization (NSBO). is presented, with its sub- routine SGD(2* "², wk T) outlined in Algorithm 5. The descent directions D D D are unbiased estimators of d di, dy, given by
3.1. Assumptions
Before presenting the theoretical results, we introduce the assumptions that will be used throughout this paper.
Assumption 3.1. (a) For any x, g(x, -) is strongly convex
with parameter μ 0.
(b) Vg is Lipschitz continuous with a Lipschitz constant Lg.1, and V22g and Vi2g are Lipschitz continuous with Lipschitz constant T.2.2-
D* V2G(a*, ,*; BB),
(13)
D = ViO s²: BB)'" - 12F(a" . BB),
(14)
D = ViE(at, ::BB) VigG(at, s/ Bl)a". (15)
(c) V Lipschitz continuous with a Lipschitz con- stant ITg,1 and there exists a constant I, such that

---

Efficient Curvature-Aware Hypergradient Approximation
(x, (x))1I s o, for all æ.
in BOX 1. This one-lime cost requires (x: log k) gradient computations and Hessian-vector products and is included in the total computational complexity.
(d) There exists a constant C5.0 such that [ (f y)1l S C,o0 for all x and
Remark 3.4. If we replace GID with the conjugate gradient method (CG) in line of NBO-GD for solving the sub- problems (8) and (10), we obtain NBO-CG. Our analysis in
Assumptions 3. are standard in the bilevel optimiza- tion literature (Ghadimi & Wang, 2018; Chen et al., 2021; Khanduri et al., 2021; Arbel & Mairal, 2022a; Dagréou et al., 2022; Ji et al., 2022; Hong et al., 2023; Chu et al., 2024). Under Assumptions 3. (C), the hypergradient v is Lipschitz continuous with a constant given by:
Appendix E. shows that NBO-CG requires fewer Hessian- vector product computations than NBO-GD, specifically O1e55 log k/e) Hessian-verctor products.
Theorem 3.2 provides theoretical complexity guarantee for NBO-GD. As shown in Table 1, the computational com- plexity of NBO-GI improves upon that of AID-BiO (Ji et al., 2021) and AmIGO (Arbel & Mairal, 2022a). Specif- ically, the gradient complexity of NBO-GD surpasses the ofatc-of-the-art result by an order of x log..
La= == L5,1 + 2L + Lg,2L1,00
(17)
1
2.g L + L2 1 +
L..2L.g11 5.0 n
+
112
as established in Lemma 2.2 of Ghadimi & Wang (2018). As
in Ji et al. (2021; 2022), we define == max (T.g,1, }
3.3. Convergence Analysis for NSBO-SGD
and k = L/p. Then L4 O(s)).
The following assumptions are made regarding the stochas- tic oracles, which will be used to analyze the convergence rate and sample complexity of NSBO-SGD Algorithm 4.
Note that Assumption (d) is not employed in the afore- mentioned literature but has been adopted in Ji et al. (2021); ); Liu et al. (2022); Kwon et al. (2023); Yang et al. (2025) for their respective purposes. In this work, this assumption is utilized to ensure that all iteration points y remain within a predefined neighborhood of y" (x*), provided that the initial point y lies within a predefined neighborhood of v² (x²)). This condition aligns with the local quadratic convergence rate of Newton's 's method.
Assumption 3.5. There exist positive constants f,1, Gg.1+ g. such that
E VF () Vf(x,)) - o),,,
E (IV - Vg(a, v)1² S og
V'G(I, - V29(1, w/]""] < o
E
3.2. Convergence Analysis for NBO-GD
Assumption 3.6. There is a constant r 2 1 such that for any iterate (x*, y²² generated by Algorithm 4, we have
Under the above assumptions, we establish the conver- gence properties of NBO-GD. The detailed proof is pro- First, we define BOX
E [lly* - y' (P S [1o*
y (=+)))) (18)
vided in Appendix E.I.
to represent the set of initial points satisfying 1ly0 y" (20)] S min { " 1 2/1,) } and l|u0 -x(x")l S
This assumption, which imposes bounded moments iter- ates, is commonly used in the stochastic Newton literature to establish convergence results in expectation form (Bollapra- gada et al., 2019; Meng et al., 2020; Berahas et al., 2020). We now establish the convergence properties of NSBO-SGD as follows. The proof details can be found in Appendix E.2. Below, we define BOX to represent the set of initial points
5L, min| "bg.2' ALz,, where L1 == L5.1 + Lg.2 Lea #
Theorem 3.2. Under Assumption 3.1, choose initial iterate (v", u,, xx in BOX 1. Then, for any constant step size Yik 7 S L/L,, there exists proper constant step size ag = a = On=33 and T 2 (K) such that NBO-GD
hus the following properties:
satisfying EE [[li] y² (a (1a) < min { 20 2,2 42, } and
a²²
E [Ilu* u (x.))1 S min (5
417 4.
(a) For all integers K3 2 1, mino 1/V )|2 < 20(a")224' (*²) That is, NBO-GD can find ak an e-optimal solution &. (i.e., [[] < é) in K=
Theorem 3.7. Under Assumptions 3.1, 3.5 and 3.6, choose an initial iterate (y u,,",, x²) in BOX 2. Then, for any con- stant step size 7Kk = Y < 1/L there exists a proper constant step size a = a = O(K and T 2 O(k) such that NSBO-SGD has the following properties:
1 steps.
(b) The computational complexity of NBO-GD is: O(K3/a) gradient computations and Jacobian-vector products, and O(k²) /c) destian-vector products.
(a) Fix K 2 1. For samples with batch sizes B 2
Remark 3.3. One can apply the gradient descent (GD) method, as described in Appendix D, to select initial points
((x²). B > 0(nK + x²), |B 2 O(k3 K + kt))
B 2
|B|| Z B(k-1K) K), it holds
5

---

Efficient Curvature-Aware Hypergradient Approximation
Table Comparison of the computational complexities of two NBO implementations with state-of-the-ar methods in the deterministic setting, including AID-BiO (Ji al., 2021), AmIGO (Arbel & Mairal, 2022a), No-loop AID (Ji et al., 2022), and F'SA (Kwon et al.. 2023; Chen al., 2023a). Note that the dependence on log R not explicitly stated AID-BiO and AmIGO, it can derived from equation (30) in Ji al. (2021) and Proposition 10 in Arbel & Mairal (2022a).
Algorithms
Convergence Rate
Gradient Computations
Hessian-Vector Vector Products
AmIGO-GD (Arbel Mairal, 2022a)
O(k²e"))
O((x* x)g
O(Kn log x)E =1)
AID-BiO (Ji al., 2021)
O(x²c-1)
O(K** log
O((13 O((c3-00ogg ()
AmlGO-CG (Arbel Mairal, 2022a)
No-loop ATD (Ji et al., 2022)
O(w"E-t)
O(w%
FASA (Kwon et al., 2023; Chen et al., 2023a)
O(k2c-))
0 (W"C
NBO-GD (this paper)
O(w(e-t)
O(x46-t)
NBO-CG (this paper)
O(x"e-))
O((x2"* k)e
that minoskçk E [[]V(2>))]] NSBO-SGD can find an e-optimal solution in K- (1)²c=1) steps.
O That is,
respectively. Unlike existing bilevel optimization literature, by leveraging a one-step inexact Newton method, we prove that after taking T 2 ((x) iterations the subroutine, llye* - n" [x*)]] satisfies the descent condition:
(b) The computational complexity of NSBO-SGD is: ( gradient complexity for F, 2())e=) gradi- ent complexity for G, O (xxec2) Jacobian-vecte prod- uet complexity, and O(K² 22) Hessian-vector product complexity.
T.g.2
1l of (x
11 (x*)x 1/112
24r
+ Tg..1 | x11| +
1
$ l1 3 (x(*))]
Next, by selecting an appropriate constant step size a and applying induction while leveraging the initialization strat- egy, we establish that lly" y²Il I 2 Hence, it follows that
Note that the number T of inner-loop steps remains at constant level. Moreover, as shown in Table 2, the computa- tional complexity of our stochastic algorithm, NSBO-SGD, improves upon the state-of-the-ar result in AmIGO by a factor of log K; for comparison, refer to Proposition 10 and Corollary 4 in Arbel Mairal (2022a). Indeed, in each outer iteration, AmIGO requires (xx (x x) gradient com- putations with batch size of O1) for each gradient, while NSBO-SGD performs single gradient computation with batch size of ( 1). As result, the complexity improvement Theorem 3.7 is not as significant as that in Theorem 3.2.
for all k.
1 A+1 y "2
2 1110" (x* s/1l
Lg | /
x²
Combining this with the descent properties of faa and ||u* == (***)] we derive the descent property of the Lyapunov function:
(ax,, yy" (a*))
VK+1 V S an V (2 |12 Ail k+1 20112
3.4. Proof Sketch
Aelly² - y [1²) Asllu - u -====""") (1²
In this section, we present a proof sketch for the determinis- tie setting as an illustrative example. The proof strategy for the stochastic case follows a similar approach. The detailed proofs of Theorem 3.2, Theorem 3.7, and Remark 3.4 are provided in Appendix E.
where A1, A2. Ag are positive constants.
Finally, by summing over iterations using a telescoping argument and estimating an upper bound for Vo. we derive the results presented in Theorem 3.2.
First, we define a Lyapunov function: V4 f(**, y where by and b, are constants that depend on и. Then, following a classical Lyapunov analysis, we establish the descent of fa y² (a*)), l V" (a*) || u" (x**))
4. Experiments
+b,llle"y"* ( 11?)
In this section, we present experiments to evaluate the prac- tical performance of the proposed NBO framework. Specifi- cally, we compare our NBO-GD and NSBO-SGD methods
6

---

Efficient Curvature-Aware Hypergradient Approximation
Table Comparison of the sample complexities of NSBO-SGD with state-of-the-art results in the stochastic setting. including those reported in ALSET (Chen et al., 2021). SOBA (Dagréou et al., 2022; Huang et al., 2025), AmIGO (Arbel & Mairal, 2022a). FSLA (Li al., 2022), FAA (Kwon al., 2023; Chen et al., 2024a). Note that p(x) indicates that the explicit dependence on r not provided in the corresponding references. For AmIGO, the dependence on log is derived from Proposition 10 in Arbel & Mairal (2022a).
Algorithm
Sample Complexity
Stochastic Estimators
ALSET (Chen et al., 2021)
SGD
SOBA (Dagrćou et al., 2022; Huang et al., 2025)
((p(k)c-2)
SGD
AmIGO (Arbel & Mairal, 2022a)
O(11)) log (s)
SGD
F'SA (Kwon et al., 2023; Chen et al., 2024a)
O
SGD
FSLA (Li et al., 2022)
Olp(s)E-3
SGD+r-momentum
NSBO-SGD (this paper)
O(1' 7)
SGD
with several widely used gradient-based algorithms, includ- ing SOBA, SABA (Dagréou et al., 2022), StoBiO (Ji al., 2021), AmIGO (Arbel & Mairal, 2022a), SHINE (Ramzi et al., 2022), F2SA (Kwon et al., 2023), and MA-SABA (Chu et al., 2024). Details of the experimental setup and additional results are provided Appendix C.
4.1. Synthetic Problem
First, we consider a synthetic problem to study NBO in a deterministic and controlled scenario. This problem focuses on hyperparameter optimization, where the upper-level and lower-level objective functions are defined as follows:
Figure I. Experimental results on synthetic data with r' 1.
to solve lower-level problems and quadratic subproblems related to Hessian inverse- vector products. Since we will later compare NBO-type algorithms with other gradient- based algorithms, we restrict the comparison in the synthetic problem NBO-GD and AmIGO for simplicity.
[(A,w)
1
Σ
D
g(A,w)
Σ
vllareee)
+
m w
D
As shown in Figure 1(a), our NBO-GD, which employs single step of the inexact Newton method, performs com- parably to AmIGO, which uses Q = 10 steps of gradient descent. Notably, the inner loop in NBO-GD is used to ap- proximate the inexact Newton direction. More importantly, Figure ((a) demonstrates that NBO-GD maintains strong performance even when the Newton direction is approxi- mated with T = 1. In contrast, AmIGO's performance degrades significantly when only one step of gradient de- scent is used. Finally, when running time is considered, Figure ((b) shows that NBO-GD (T = outperforms other methods. Based on these findings, we use 1 for NBO- type algorithms in subsequent experiments.
(EDW
where Rº represents the hyperparameter, w R x denotes the model parameter, and (+) log(1 + is the logistic loss function. Here, D' and D represent the validation and training datasets, respectively. The synthetic data is generated following a procedure similar to that de- scribed in Chen et al. (2023b; 2024b); Dong et al. (2023). Specifically, the distribution of x, follows a normal distribu- tion N (²)) and ye. = wx, + 0.1 where is sampled from (O, 1). Subsequently, ye transformed into binary labels: if ye exceeds the median of the dataset, it is set to 1; otherwise, it is set to 1.
The experimental results for r' = are shown in Figure 1, while the results for = 0.5 and r' = are provided in Appendix C. AmIGO is a representative method that employs multiple (Q) (stochastic) gradient descent steps
4.2. Hyperparameter Optimization
In this section, we evaluate the empirical performance of our NSBO-SGD algorithm on hyperparameter selection prob-

---

Efficient Curvature-Aware Hypergradient Approximation
lems, a typical bilevel optimization task (Franceschi al., 2018; Ji et al., 2021; Dagréou al., 2022). The upper- and lower-level objective functions are defined as follows:
solving lower-level problems and quadratic subproblems in AmIGO.
4.3. Data Hyper-CCeaning
[(A,w)
Σ
I (Ile: Veiw),
|D'|
We also conduct data hyper-cleaning experiments (Franceschi et al., 2017: Dagréou et al., 2022) two datasets: MNIST and FashionMNIST (Xiao al., 2017). Data hyper-c involves training a multiclass classifier while addressing training samples with noisy labels. It can be formulated as a bilevel optimization problem with the following objective functions:
g(A.w)= 1 (D]
Σ
l (Fe: W) + r(),w),
(ED.
where D' and D denote the validation and training datasets, respectively. Here, À is the hyperparameter, w E Rex denotes the model parameter, I is the loss function, and r(A,w) w) is a regularizer.
f(A, w) = [D'
Llwad' ve)s
We conduct experiments on two datasets: IJCNNI and Cov- type. The IJCNN1 dataset corresponds to binary classifica- tion problem using logistic regression, with p 22, 1. The loss function is defined as E( ye:w) := (wreBee). where & is the logistic function. The regularizer is given
((A,w) 1
Σ [Xee)[(wTer Ve) crllol? (ED:
D
where C denotes the cross-entropy loss, o is the sigmoid function, and G is regularization parameter. For this experiment, we set the corruption probability to p' = 0.5. as in (Dagréou et al., 2022).
by 2 :i= e>im"" w. The Covtype dataset corresponds to a multi-class classification problem using logistic regres- sion, with p- 54, 17. The loss function is the cross- entropy function, and the regularization term is given by 1 e² El=1 wf.-
505A
2254
S1INN
N1000-11
30.74
Asa00
SudtHO
50DA
F25A
SUINZ
NSBO-ED
33,0%
2
20.7%
31,07
24074
16/90
22.54
20/9
130
Tixai
Thxat
16*
2500%
130 2700g
100
Figure 3. Comparison between NSBO-SGI and other algorithms on data hyper-cleaning. Left: MNIST dataset; Right: FashionM- NIST dataset.
Figure 2. Comparison between NSBO-SGD and other algorithms on hyperparameter optimization. Left: IJCNN1 dataset; Right: Covtype dataset.
In Figure 3, we report the test errors for each method w.t running time on MNIST and FashionMNIST. We observe that both AmIGO and NSBO-SGD outperform all other methods by reaching the smallest error. Meanwhile, NSBO- SGD is the fastest, demonstrating the efficiency of incorpo- rating curvature information in bilevel optimization.
In Figure 2, we present the suboptimality gap on IJCNNI and the test error on Covtype as functions of time. Our first observation is that, among all methods, NSBO-SGD ex- hibits the fastest convergence both datasets, with its per- formance advantage being particularly evident on IJCNN1. Second, the gap between SOBA and AmIGO on IJCNN1 highlights the benefits of performing multiple SGD steps when solving lower-level problems, as well as handling quadratic subproblems involving Hessian rnverse-vector products. In comparison, NSBO-SGD highlights the practi- cal advantages of incorporating curvature information from the subroutine that computes an approximation of Newton directions. Notably, in Figure 2, NSBO-SGD uses T' = 1, whereas AmiGO employs Q = 10. This demonstrates that the inner loop for approximating Newton directions in NSBO-SGD is more effective than the inner loop used for
4.4. Other Implementations of NBO
Variance reduction and momentum. Similar to the gradient-based framework in (Dagréou et al., 2022), the simplicity of our NBO approach allows us to easily incorpo- rate reeducee gradient estimators and momentum techniques (Chen et al., 2024c; Chu et al., 2024). We con- duct numerical experiments to test the versatility of NBO by comparing it with the framework in (Dagréou et al., 2022), using the same auriance-reducee gradient estimator and momentum technique.
8

---

Efficient Curvature-Aware Hypergradient Approximation
SABA (Dagréou al., 2022) is an adaptation of the vari- ance reduction algorithm SAGA (Defazio et al., 2014) for bilevel optimization. MA-SABA (Chu et al., 2024) extends the SABA algorithm by incorporating an additional stan- dard momentum (also referred to as moving average) into the update of the upper-level variable. We take a similar approach and extend SAGA to NBO by replacing SGD in NSBO-SGD with SAGA, resulting in a new implementation of NBO, referred to ASBO-SAGA. If we further incorpo- rate an additional standard momentum into the update of the upper-level variable, we obtain MA-NSBO-SAGA, another new implementation of NBO.
is straightforward to derive that the optimal solution is (e, e, e). The results, shown in Figure 5, demonstrate that BAMM+NBO significantly outperforms BAMM, highlight- ing the effectiveness of the NBO framework.
BAND
BAMM
In Figure 4, we present the comparison results between SABA and NSBO-SAGA, as well as between MA-SABA and MA-NSBO-SAGA. We observe that both NSBO-SAGA and MA-NSBO-SAAGE achieve better performance.
Figure 5. Performance of the NBO framework on a toy example with a non-strongly convex lower-level problem.
5. Conclusion
SABA
NEBO-FAGA
MA-SABA
MANSBO-5AGA
10.7%
In this work, we design and analyze a simple and efficient framework (NBO) for bilevel optimization, leveraging the hypergradient structure and inexact Newton methods. The convergence analysis and experimental results for specific examples of NBO demonstrate the benefits of incorporat- ing curvature information into the optimization process for bilevel problems. However, many additional benefits and extensions could be explored, such as the exploitation of parallel and distributed computation, and the integration of noise reduction techniques.
15.0%
7
$4.99
1409:
18.
13.90
Figure 4. Comparison between the NBO framework and other al- gorithms incorporating variance reduction and the moving average technique. Left: Hyperparameter optimization on IJCNNI; Right: Data hyper-cleaning on MNIST.
A Acknowledgements
Extension to lower-level non-strongly convex structures. The NBO framework is primarily designed for bilevel op- timization problems in which the lower-level objective is strongly convex. For problems with a non-strongly convex lower-level structure, existing methods often reformulate the original problem to induce strong convexity. Once this condition is met, the NBO framework can be applied. For in- stance, in the BAMM method (Liu et al., 2023b), when y is merely convex, an aggregation function e,, + -4 is introduced. This function becomes strongly convex if f is strongly convex. An approximate hypergradient de can then be computed by substituting g with dpp-
Authors listed in alphabetical order. This work is supported by the National Natural Science Foundation of China (12431011, 12371301, 12371305, 12222 10622106, 12326605), Natural Science Foundation for Distinguished Young Scholars of Gansu Province (22JR5RA223), Guang- dong Basic and Applied Basic Research Foundation (No. 2022B 515000082) and the Longhua District Science and Innovation Commission Project Grants of Shenzhen (Grant 2o.202501 13G43468522). We thank the anonymous re- viewers for their valuable comments and constructive sug- gestions on this work.
We compare the BAMM method and BAMM+NBO, where the NBO framework is used to compute di within BAMM, on the toy example proposed in Liu et al. (2023b). The example is defined as follows:
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
1 min || 22 + 2 zER" 2
1i ell
s.t. y (v1,22) E arg min (v..
1 |1 x² yi.
TTy,
References
2
Arbel, M. and Mairal, J. Amortized implicit differentia- tion for stochastic bilevel optimization. In International
where is a vector with all components equal to It

---

Efficient Curvature-Aware Hypergradient Approximation
Conference on Learning Representations, 2022a.
Cutkosky, A. and Orabona, F. Momentum-based variance
reduction in non-convex sgd. In Advances in Neural Information Processing Systems, volume 32, 2019.
Arbel, M. and Mairal, J. Non-convex bilevel games with critical point selection maps. In Advances in Neural Information Processing Systems, volume 35, pp. 8013- 8026, 2022b.
Dagréou, M., Ablin, P., Vaiter, S., and Moreau, T. frame- work for bilevel optimization that enables stochastic and global variance reduction algorithms. In Advances in Neural Information Processing Systems, volume 35, pp. 26698-2671 2022.
Berahas, A. S., Bollapragada, R., and Nocedal, J. An inves- tigation of newton-sketer and subsampled newton meth- ods. Optimization Methods and Software, 35(4):661-68) 2020.
Dagréou, M., Ablin, P., Vaiter, S., and Moreau, T. How to compute hessian-vector products? ICLR Blogposts, 2024.
Bishop. N., ing from verified training data. In Advances in Neural Information Processing Systems, volume 33, pp. 9520- 9529, 2020.
L., and Gerding. E. Optimal leamn
Dagréou, M., Moreau, T., Vaiter, S., and Ablin, P. A lower bound and a near-optimal algorithm for bilevel empiri- cal risk minimization. International Conference on Artificial Intelligence and Statistics, pp. 82-90, 2024.
Bollapragada, R., Byrd, H., and Nocedal, J. Exact and in- exact subsampled newton methods for optimization. IMA Journal of Numerical Analysis, 39(2):545-578 2019.
Defazio, A., Bach, F., and Lacoste-Julien, S. Saga: A fast in- cremental gradient method with support for non-strongly convex composite objectives. In Advances in Neural Information Processing Systems, volume 27, 2014.
Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D.. Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne S., and Zhang, Q. JAX: composable transformations of Python+NumPy programs, 2018.
Dong, Y., Ma, S., Yang. J., and Yin, C. A single- loop algorithm for decentralized bilevel optimization. arXiv:2311.0894 2023.
Chen, L., Ma, Y., and Zhang, J. Near-optimal nonconvex- strongly-convex bilevel optimization with fully first-order oracles. arXiv:2306. 14853, 2023a.
Franceschi, L., Donini, M., Frasconi, P., and Pontil, M. Forward and reverse gradient-based hyperparameter opti- mization. In International Conference Machine Learn- ing, pp. 1175, 2017.
Chen, L., Xu, J., and Zhang. On finding small hyper- gradients in bilevel optimization: Hardness results and improved analysis. In Annual Conference on Learning Theory, pp. 947-980. PMLR, 2024a.
Franceschi, L., Frasconi, P., Salzo, S., Grazzi, R., and Pontil, M. Bilevel programming for hyperparameter optimiza- tion and meta-learning. In International Conference on Machine Learning, volume 80, pp. 1568-1577, 2018.
Chen, T., Sun, Y., and Yin, W. Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems. In Advances in Neural Information Processing Systems, volume 34, pp. 25294-25307, 2021.
Ghadimi, S. and Wang, M. Approximation methods for bilevel programming. arXiv:1802.02246. 2018.
Chen, X., Huang, M., Ma, S., and Balasubramanian, K. De- centralized stochastic bilevel optimization with improved per-iteration complexity. In International Conference Machine Learning, volume 202, pp. 4641-4671, 2023b.
Grazzi, R., Franceschi, L., Pontil, M., and Salzo, S. On the iteration complexity of hypergradient computation. In International Conference on Machine Learning, volume 119, pp. 3748-3758, 2020.
Chen, X., Huang, M., and Ma, Decentralized bilevel optimization. Optimization Letters, pp. 1-65, 2024b.
Hong, M. Wai, H.-T., Wang. Z., and Yang. Z. A two- timescale stochastic algorithm framework for bilevel op- timization: Complexity analysis and application actor- critic. SIAM Journal on Optimization, 33(1):147-180, 2023.
Chen, X., Xiao, T., and Balasubramanian, K. Optimal algo- rithms for stochastic bilevel optimization under relaxed smoothness conditions. Journal of Machine Learning Research, -51 2024c.
Huang, Optimal Hessian/lacobian-free nonconvex-PL bilevel optimization. In International Conference Machine Learning, 2024.
Chu, T., Xu, D.. Yao, W., and Zhang, J. SPABA: A single- loop and probabilistic stochastic bilevel algorithm achiev- ing optimal sample complexity. In International Confer- ence on Machine Learning, 2024.
Huang, Y., Wu, Z., Ma, S., and Ling, Q. Single- multi-sequence stochastic approximation without fixed
10

---

Efficient Curvature-Aware Hypergradient Approximation
point smoothness: Theories and applications. IEEE Trans- actions Signal Processing, 2025.
Liu, R., Gao, J.ang. J., Meng, D., and Lin, Z. Investigat- ing bi-level optimization for learning and vision from unified perspective: A survey and beyond. IEEE Trans- actions on Pattern Analysis and Machine Intelligence, 2021.
Ji, K., Lee, J. D., Liang, Y., and Poor, H. V. Convergence of meta-learning with task-specific adaptation over partial parameters. In Advances Neural Information Process- ing Systems, volume 33. pp. 14490-11500, 2020.
Liu, R., Liu, X., Zeng, S., Zhang, J., and Zhang, Y. Value- function-based sequential minimization for bi-level op- timization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023a.
Ji, K., Yang, J., and Liang, Y. Bilevel optimization: Con- vergence analysis and enhanced design. In International Conference on Machine Learning, pp. 4882-4892, 2021.
Liu, R., Liu, Y., Yao, W., Zeng, S., and Zhang. J. Averaged method of multipliers for bi-level optimization without lower-level strong convexity. In International Conference on Machine Learning, pp. 11839-21866, 2023b.
Ji, K., Liu, M., Liang. Y., and Ying, L. Will bilevel optimiz- ers benefit from loops. In Advances in Neural Information Processing Systems, volume 35, pp. 3011-3023, 2022.
Khanduri, P., Zeng. S., Hong, M., Wai, H.-T., Wang. 7 and Yang, Z. A near-optimal algorithm for stochastic In Advances
Maclaurin, D., Duvenaud, D., and Adams, R. Gradient- based hyperparameter optimization through reversible learning. In International Conference on Machine Learn- ing, 2015.
bilevel optimization via in Neural Information Processing Systems, volume 34, pp. 30271-30283, 2021.
Meng, S. Y., Vaswani, S., Laradji, I. H., Schmidt, M., and Lacoste-Julien, S. Fast and furious convergence: Stochas- tic second order methods under interpolation. In Interna- tional Conference on Artificial Intelligence and Statistics, pp. 1375 1386, 2020.
Kwon, J., Kwon, D., Wright, S., and Nowak, R. D. A fully first-order method for stochastic bilevel optimiza- tion. In International Conference on Machine Learning, pp. 18083-18113, 2023.
Kwon, J., Kwon, D., Wright, S., and Nowak, R. D. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. In International Conference on Learning Representations, 2024.
Moreau, T., Massias, M., Gramfort, A., Ablin, P., Bannier, P.-A., Charlier, B., Dagréou, M., Dupre la Tour, T., Durif, G., Dantas, C. F,, et al. Benchopt: Reproducible, efficient and collaborative optimization benchmarks. In Advances in Neural Information Processing Systems, volume 35, 2022.
Lam, S., Pitrou, A., and Seibert, S. Numba: a llym-based python jit compiler. In Workshop on the LLVM Compiler Infrastructure HPC, pp. 1-6,1 2015.
pp.
Nesterov, Y. Lectures on convex optimization, volume 137. Springer, Cham, 2018
Li, J., Gu, B., and Huang, H. A fully single loop algo- rithm for bilevel optimization without hessian inverse. In AAAI Conference on Artificial Intelligence, volume 36, pp. 7426-7434, 2022.
Nocedal, J. and Wright, S. J. Conjugate gradient methods. Numerical optimization, pp. 101-134, 2006.
Li, Z.,, Bao, H., Zhang, X., and Richtárik, P. Page: A simple and optimal probabilistic gradient estimator for nonconvex optimization. In International Conference on Machine Learning, pp. 6286-6295, 2021.
Pearlmutter, B. A. Fast exact multiplication by the hessian. Neural computation, 6( 6(1):147-160, 1994.
Pedregosa, F. Hyperparameter optimization with approxi- mate gradient. In International Conference on Machine Learning, volume 48, pp. 737-746, 2016.
Liang, H., Zhang, S., Sun, J., He, X., Huang, W., Zhuang, K., and Li, Z. Darts+: Improved differentiable architecture search with early stopping. arXiv: 00030, 2019.
Rajeswaran, A., Finn, C., Kakade, S. M., and Levine, S. Meta-learning with implicit gradients. In Advances in Neural Information Processing Systems, volume 32, pp. 113-124, 2019.
Liu, B., Ye, M., Wright, S., Stone, P., and Liu, Q. Bome! bilevel optimization made easy: A simple first-order ap- proach. In Advances in Neural Information Processing Systems, volume 35, pp. 77724-17266 2022.
Ramzi, Z., Mannel, F., Bai, S., Starck, J.-L.,, Ciuciu, P., and Moreau, T. Shine: Sharing the inverse estimate from the forward pass for bi-level optimization and implicit models. In International Conference on Learning Repre- sentations, 2022.
Liu, H., Simonyan, K., and Yang, Y. Darts: Differentiable ar- chitecture search. In International Conference on Learn- ing Representations, 2018

---

Efficient Curvature-Aware Hypergradient Approximation
Shen, H. and Chen, T. On penalty-based bilevel gradient descent method. In International Conference Machine Learning, pp. 30992-31015, PMLR, 2023.
Wang, J., Chen, H., Jiang, R., Li, X., and Li, Fast algorithms for Stackelberg prediction game with least squares loss. In International Conference on Machine Learning, volume 139, pp. 10708-10716, 2021.
Xiao, H., Rasul, K., and Vollgraf, R. Fashion-mnist: novel image dataset for benchmarking machine learning algorithms. 10X.v:7708.0777, 2017.
Xiao, Q., Lu, S., and Chen, A generalized alternat- ing method for bilevel optimization under the polyak- tojasiewicz condition. Advances in Neural Information Processing Systems, pp. 63847-63873, 2023.
Yang, J., Ji, K., and Liang, Y. Provably faster algorithms for bilevel optimization. Advances in Neural Information Processing Systems, volume 34, pp. 13670-13682, 2021.
Yang, Y., Gao, B., and Yuan, Y.-x. Lancbio: Dynamic lanczos-aided bilevel optimization via krylov subspace. In International Conference on Learning Representations, 2025.
Yang, Z., Chen, Y., Hong, M., and Wang, Z. Provably global convergence of actor-critic: A case for linear quadratic regulator with ergodic cost. In Advances Neural Infor- mation Processing Systems, volume 32, 2019.
Ye, J. J. and Zhu, D. Optimality conditions for bilevel programming problems. Optimization, 33(1):9-27, 1995.
Young, D. On richardson's method for solving linear sys- tems with positive definite matrices. Journal of Mathe- matics and Physics,
1953.
Zhang, Y., Zhang, G., Khanduri. P., Hong. M.. Chang, S., and Liu, S. Revisiting and advancing fast adversarial training through the lens of bi-level optimization. In International Conference Machine Learning, volume 162, pp. 26693-26712, 2022.
Zhang, Y., Khanduri, P., Tsaknakis, I., Yao, Y., Hong, M., and Liu, S. An introduction to bilevel optimization: Foun- dations and applications in signal processing and machine learning. IEEE Signal Processing Magazine, 41(1))38- 59, 2024.
12

---

Efficient Curvature-Aware Hypergradient Approximation
A. Appendix
The appendix is organized as follows:
More related work is provided in Section B.
Experimental details and additional experiments are provided in Section C.
The initialization strategy of the proposed algorithms is provided in Section D.
Detailed proofs of the main theorems are provided in Section E.
B. More Related Works
Bilevel optimization. Bilevel optimization addresses the challenges associated with nested optimization structures commonly encountered in various machine learning applications, as discussed in recent survey papers (Liu et al., 2021; Zhang et al., 2024). Although numerous methods have emerged, the numerical computation of bilevel optimization remains a significant challenge, even when the lower-level problem has a unique solution (Pedregosa, 2016; Ghadimi & Wang, 2018; Kwon al., 2023). Assuming further that the Hessian of the lower-level objective function w.r. the lower-level variables invertible, the hypergradient well-defined and can be derived from the implicit function theorem. Recently, several strategies have been proposed solve bilevel optimization. For instance, the iterative differentiation (TTD)-bascd approach (Maclaurin et al., 2015; Franceschi et al., 2018; 8: Grazzi et al., 2020; Ji et al., 2021) estimates the Jacobian of the lower-level solution map by differentiating the steps used to compute an approximation of the lower-level solution. The approximate implicit differentiation (AID)-based approach (Ghadimi & Wang, 2018; Chen et al., 2021; Ji et al., 2021; 2022; Arbel & Mairal, 2022a; Dagréou et al.. 2022: Hong et al., 2023) is directly based on equation (2). It performs several (stochastic) gradient descent steps in the lower-level problem, followed by the estimation of the Hessian inverse-vecte vector product, which can be computed using Neumann approximations (Chen et al., 2021; Hong al.. 2023), solving a linear system (Pedregosa, 2016; Ji et al., 2021), or solving quadratic programming problem (Arbel & Mairal, 2022a; Dagréou al., 2022).
The ITD- and AlD-based methods involve numerous Hessian- and Jacobian-vector products, which can be efficiently computed and stored using existing automatic differentiation packages (Pearlmutter, 1994; Dagréou et al., 2024). Another class of methods, which rely solely on the first-order gradients of the upper- and lower-level objective functions, is based on the value function reformulation (Ye & Zhu, 1995) and referred to as the value function-based approach (Liu et al., 2022; Kwon et al., 2023: Chen et al.. 2023a; Liu et al., 2023a; Kwon et al., 2024). From the hypergradient perspective, all of the aforementioned gradient-based bilevel optimization algorithms are inexact hypergradient methods. The distinction between these methods lies in whether the hypergradient approximation is directly derived from equation (2). For instance, the AID-based approach is directly based on equation (2), whereas the ITDD-baeed and value function-based approaches are not. However, what these methods do is to approximate the hypergrudient. For example, Proposition 2 in Ji et al. (2021) provides the explicit form of the ITDD-bseed hypergradient estimate, while Lemma 3.1 in Kwon et al. (2023) offers a hypergradient estimate for the value function sunction-based approach. In stochastic bilevel optimization, various stochastic techniques (e.g., momentum and variance reduction) from single-level optimization have been employed to improve the convergence rate (Yang et al., 2021; Dagréou al., 2022; Khanduri al., 2021; Chen al., 2024c; Dagréou et al., 2024; Chu al., 2024; Huang, 2024).
Although our paper focuses on the lower-level strongly convex case, the lower-level non-strongly convex case is also of great interest. In some works, the hypergradient can still be approximated in non-strongly convex case, by combining the pseudo inverse of the Hessian (Arbel & Mairal, 2022b; Xiao et al., 2023) or aggregation functions (Liu et al., 2023b). Then this approximate hypergradient can be used to iteratively update the upper-level variable. There are also works that employ penalty function methods, including but not limited to Kwon al. (2024); Shen & Chen (2023); Chen et al. (2024a).
C. Experimental Details and Additional Experiments
Our experiments were conducted using the Bilevel Optimization Benchmark framework (Dagréou et al., 2022) and the Benchopt library (Moreau al., 2022). All experiments were performed on system equipped with an Intel(R) Xeon(R)
13

---

Effiiiett Curvature-Aware Hypergradient Approximation
Gold 5218R CPU running at 2.10 GHz and NVIDIA H100 GPU with 80 GB of memory. Each experiment repeated ten times.
C.1. Synthetic Data
For the synthetic dataset, we use 16,000 training samples and 4,000 validation samples. The dimension size is set p 50. For both AmIGO and NBO, the outer step size is set to I, and the inner step size is set to 0.03. Both algorithms deterministic, employing full-batch updates. The experimental results for r' = 0.5 and r' = 2 are presented in Figures 6 and 7, respectively.
(1+10)
(I-1)
AmtOO (Q-12)
AstGO (Q-1)
10"
10
00
Figure 6. Experimental results on synthetic data with 0.5.
NBO-CID (1w10)
NBO-GD (1+1)
Amilo (Q=10)
Amtio
In"!
6.
5.0
7.5
10.0
125
Derelion
Thooes
Figure 7. Experimental results on synthetic data with r' = 2.
The results align with our discussion in Section 4.1. Figures and 7 demonstrate that NBO-GD maintains strong performance even when the Newton direction is approximated with when only single step of gradient descent is used. When considering runtime, the results show that NBO-GD (T' outperforms other methods. Moreover, evaluate the scalability of NBO synthetie data while progressively increasing the problem dimension in Figure 8. We record the time it takes to achieve (/(A*, - (X*)) /')/ w (2")) - x*) 10 when the dimension increasing, The results after 5000 iterations is used as 4*.
In contrast, AmIGO's performance deteriorates significantly
C.2. Hyperparameter Optimization and Data
Datasets. We provide detailed information regarding the datasets used in our experiments. For IJCNNI', we employ 49.990 training samples and 91,701 validation samples. For Covtype', we utilize 371,847 training samples, 92,962 validation samples, and 16,203 testing samples. For 3 and FashionMNISTT*. we use 20,000 training samples, 5,000 validation
wwt.saaaa.nattttt.fetettse covtype.html
'htpps//github.co hhtpps//gttuh..co//aaaaaaat htt.g//wwh...mmmmmmmmdmmea///iii
14

---

Efficient Curvature-Aware Hypergradient Approximation
XB0-DD
35
10
100
200
303
400
340
Figure 8. Scalability NBO to large dimensional problems on synthetic data.
samples, and 10,000 testing samples, and we corrupt the labels of these two datasets with a probability of p' = 0.5. In terms of computational frameworks, we use JAX (Bradbury et al., 2018) for MNIST, FashionMNIST, and Covtype. For IJCNNI, we use Numba (Lam et al., 2015), as demonstrates faster performance compared JAX for this dataset.
Algorithm Settings. We provide a detailed description of the batch size, step size, and other settings for the algorithms compared in Section 4.2. The batch size for all algorithms is set to 64, except for NSBO SGD and SHINE. For NSBO-SOC since the size of B is significantly larger than that of other batches in theory, set B 256, while the other batches remain at 64. For SHINE, we set the batch size to 256, as it is a deterministic algorithm by design. In the benchmark, the double-loop algorithms, including AmIGO, StocBiO, and F2SA, are configured with 10 inner loop iterations, and we set the number of iterations for BFGS in SHINE to 5. The step sizes are tuned via grid search, following a method similar to (Dagréou et al., 2022). The grid search procedure is described in detail below. For algorithms that employ decreasing step size, the step sizes we search for correspond their initial values. For SHINE, the inner step size we search for is related to the initial step size of the ttrong-Wolfe line search. The outer step size is computed inner step size outer ratio
The inner step size is chosen from 6 values between 2-5 and 1, spaced on a logarithmic scale. The outer ratio is chosen in {10-² 10 1.5 10 10-0.5, 1].
IJCNN1:
0.5
Covtype: The inner step size is chosen from 7 values between 2-6 and spaced logarithmic scale. The outer ratio is chosen from 4 values between 10-2 and 10, spaced on logarithmic scale.
MNIST and FashionMNIST. The inner step size chosen from values between 10-3 33 and 1, spaced on logarithmic scale. The outer ratio is chosen from values between 10 10-6 and 10-3, spaced on logarithmic scale.
Moreover, for MA-SABA and MA-NSBO-SAGA, we set the moving average coefficient to 0.99 for MNIST and FashionM- NIST, 0.9 for IJCNNI, 0.8 for Covtype. Last but not least, we do not follow the initialization in Algorithm 2 and Algorithm 4. Instead, we set the initial points for our framework in the same way as other algorithms in the benchmark. Other settings not mentioned are kept consistent with the benchmark.
Additional results. We present the additional experimental results for Covtype and FashionMNIST. In Figure 9 (a). the additional results of hyperparameter optimization on Covtype are shown. In Figure (b), the additional results of data hyper-cleaning on FashionMNIST are displayed. These results demonstrate the effectiveness of our framework under different implementations. In Figure (b), the performance of MA-SABA and MA-NSBO-SA SAGA is similar, likely because MA-SABA is fast enough on MNShionMNIST, leaving little room for further acceleration.
15

---

Effiiient Curvature-Aware Hypergradient Approximation
SABA
NSBO-SAGA
MA-SABA
MA-NSBOSAGA
28.5%
28.45
23. 536
22.9% 11
28
21.5%
31.0%%
Troels
opticiiaaionn
Figure 9. Comparison between NBO framework and other algorithms equipped with variance reduction and moving average technique. Left: Hyperparameter optimization on Coviype: Right: Data hyper-cleaning on MNSST.
D. Initialization Strategy of the Proposed Algorithms
In this section, we introduce two initialization boxes for our proposed algorithms. These initializations ensure that the initial points y", 2 are close y" (xxx), u² (z²), respectively, thereby guaranteeing that y², u² remain within neighborhood of y" (z*)), (ax)) for all 0, 0,-.. K with controlled step size ag. important note that the primary purpose of these initialization requirements is to ensure theoretical validity. In practice, strict adherence to these constraints is not necessary, as the proposed initialization process does not introduce significant computational overhead. Some constants in this section, such as ru and L1, are defined in (29).
Choose α² 000 000 stepsize Bo S 1.S11
1
min { \ 2t
Take No 2 2 log
log(1 - BB04) = (1s log x) and compute y = No.0.
|| \o ny اللl
for = 0, 1,.1,,... ,.,...,No. No
y n+1.0 1m - 3o V2g(z 00
end for
Take Qo > log
min
52, 2k
2m
log(1 3014 = Oxx log k) and compute v = 40.00
010000 + 4 | 0= .0)
for
g = 0, ,0
" 90+00 u.00 - Bu V22 22 - V2f(o",
end for
Output x0 w0,y",000 u.
BOX 1: Initialization of Algorithm 2.
16

---

Efficient Curvature-Aware Hypergradient Approximation
Choose x0, 000 0,0 stepsize 3o 5 Lya" 1. batch size Bo | 2 2(x"), 2 (1)) according to (24).
min हoं tor
Take No 2 log yy
/log(1 Boll) = O log k) and compute y = No..
for n 0, ,,,1,,,,,,, .ooNoo1
y -+1,0 = 1/ Bo v2 Bo,a)
end for
213 11
min 6-12 2/12
Take Qo 2 log
log(1 - Bo//2) Ox re) and compute u 20.00
4 [1100.01+ VVV 2
for = 0,1...,,0 1
u5-1,0 = - 3o V220(xx,, ; Bu, )b..) (2f Bb,))
end for
Output ,²,y" u.
BOX 2: Initialization of Algorithm 4.
The following two lemmas describe the properties of the output sequences from BOX and BOX 2, respectively. Lemma D.1. Suppose Assumption 3.1 holds, then x.. v u in BOX I satisfy 14ll v=(x)] s min { 2 2v 1
||u0 (29)) S min {
5LL vLa }.
2L,
Proof. For y = ymh - Bu V2g yy we have
11 y" (2 = -1 8o 22 ym. - - (1
= -11y y (1 - 280 (y - y² (x"), V29(2 y(mlo,r,)))) )) 2 V2g(2 ( )1
Since V2g(2 y² "eg(a"."""""" )x² 0, Assumption 3.1 and Theorem 2.1.12 in (Nesterov, 2018) implies that
(y - n²² (a.), T2 (y - (x²),, V2g(2 v00) V2g(<" y² (xxx())) (x)1))
|L mm 3" 112 + # Ln
1
V2g(20) ,ym00)
# + Lg
Since 8o I 1 I #+h111' we have
L,
(1- 2,00 1bbgg11
- (2))) < (1 Bou)lly - (12
(19)
11 Tzg.1
Then we can deduce that
|| - (1 s (1 No/2 - (1')") (1 (8)) No | Vg(a" 1. < min (
No/2
1 P
2//1 2T.
where the second inequality follows from
0.0 - (2xx)) s - -VVooe 000) Vg(x", ,3"" (xx)) = = [VVg(z ا
(20)
17

---

Efficient Curvature-Aware Hypergradient Approximation
by /-Strron convexity of g(x, -) and the optimality of y (x²). Next, recall that u" (x(,,)) (22 - (2, v). the difference between u" (2)) and u" (x, is
a||u" (x, v) u² (x)]l 1V220a v)(" (z, -)
< [2]((, - V2f(0, (x))1l + (x))ll 22 y (x)) [229(1,))11
5 BL f, + galgg))lly y (=) Lilly y (x)ll,
(21)
where ||u* (x)ll S Lr Tu by definition, the second incquality is because V22g(2, y (x)))²(x) V (x)) and V22g(1, y)u" (v. = 22](,,v) v). For uo+1.0 = w - 30 V V20 ,0 nn.0 (25(x" 1) by p-strong convexity of
g(Ø,-) and 8o S Lg. we have
| | 1.0 - u (5 v0)
|| - 3V (000((""))001 )0))=000 - -V2s(a, x0)) - n "''x",1' v
$ T 3o V V220(0 0 y y>11 ||000 - u² (x.), vn) | S _o1))lue.0 u² (x) )l
(22)
So we can deduce that
|| - - (x(x))11 = -110 (x"", + 1| (x" - vi' (x9))
(20)(21) (1 Bot) 20 1000 2² u²(xx, x0)l + L
-11000 y" (x")|l
< (1 Bolt) (112 11 + (1 (2",0))1 )))
L
min {
p } min
5L1 vhn
}.
2
2L,
2L..2 p
where
11:0 ( S [2"" n²' (x9)) +l=""
21) L (=9))
S
11100 -
+ ru
L S Vg(2 ,00 +/0+++00 + tru.
(23)
122
Hence we complete the proof.
Lemma D.2. Suppose Assumptions 3.1, 3.5 hold, when we take the bachsize
8300 9.2
22
32 2 ( Vo(a 19)l + +ru)² O(x²),
BBll 2 max
-(xx),
j mmmnn(2 (120/400222 (20) 1,)
(5 L L
11
16:3007
|Bb 2 0(1)
2
(24)
mmnn Li/(n 22,), L1 (4//))
then x0, 20,8, ud in BOX 2 satisfy E [111] - v" S min
4² 1 and [||e0 (29))1] < min {
4L Lg
}.
20r aL,
5r
Proof. Similar to Lemma D.I, we can deduce that
mly Veg(a" .vm00 v²(x'"))))')))"") (2²))²² s (1- 1 y" (2
Then denote a conditional expectation Fy = [. |x.h, y0. for 10 ym - 30 22G .ym, Bo), we have
EB 1110 "y [2²) = || Bo 22 00 12" (x (x0) 112 + BEE || 12 ) - v2G(20 000 n. B
< (1 - (.30))l 3" 5² )12 + 88 2.11
(25)
Bol
lake expectation and by induction we obtain
22
112
E. No.D y² [1 s (1 No. 22²")) + 30
< min
4L, }
/4|
Bo
20r 1 9.2
18

---

Efficient Curvature-Aware Hypergradient Approximation
where 1,000 - (((2") A / |
from (20). Next, from (22) we have
|10:00 - Po (V (x0 5²) va) V2f(2" - u² (x" , s s
(1- - (Boo)))llu - -a"""" 30)),
denote Eg = B 1=0, v², , we have
Ey 11 u (z) . || = =1|ug00 - 30 (3(V222 3u (V20(22) ,2 no -
- V2f .v0)) - u (x", 2/
= ||| 1220 2/ 22000 - V2f(=> v0)) G2
y::
+ +2F(2=) :o)
,y Bu)l
<(1 Bo14) 11 - u¹ (x) )2 + 2,8
111
+ 4,6
9.2
11000
u u
","""''''' y0 + 4,88 9.2
|| (xx ,
B
Bpl
S(1 - (o-oop)) u² (x )2 2
f11
+ 43
9.1
||"
(2", I1?
B
[B
8 Then taking conditional expetation E' = =E|-|2"; y on both side of the
where the last inequality is because inequality and by induction,
|Bll 2
of 1 B
2.11
E' [1 u² (x", - S (1 - S01l2) 10 u² (2).
+ 43u-
+ 830
1 (x", y11
B
Combining with (21) and (23), we obtain
E' [Il u² (=²")))"" [[]] < 2E' [I - -u((x" ",''"")1" + + 211u0(("))
un (x",, ,,,²
2 - 8/ 2 (11u 0 || + (xx ,y0)1)) + 880 0
+ 16,30 0g11
21
(x y² +
E [100 - y (2
BBB
B
112
Ia Vg(2"; y0 + ra)²
$2 8ou//) 0 000 (12001
+ 2
+ 8/30
f.1
+ 16,80
Li Vg(a ا0) + +rm)²+ +
2+
2LF 112
E [al] - (x < min
4L3 L 5r : 112 J:
Bo
2
then we complete the proof.
Remark D.3. From Lemma D.1, we deduce that the initialization in BOX 1 requires ( log x) gradient computations and Hessian-vectus vector products. Similarly, Lemma D.2 indicates that the initialization in BOX 2 has a gradient complexity of O(* log x) for g. Ox x) for f. and Hessian-vector product complexity of O(55 log
E. Detailed Proofs
In this section of the Appendix, we provide detailed proof' of Theorem 3.2, Theorem 3.7, and Remark 3.4. Recall that
d = Vif(x , - Vizo(=* )²)²", D = 12F y BB) d = V2g(x* ,²1,, D3 = v2O(2" BB). d = V22g(a", "2")"² Vzl(a * ,²), D : Glz y B+iu* V2F(1*, (B))
; - GizG v Bl)u",
(26)
(27)
(28)
For convenience, we further define
Tu := I.S0//, M1 T. = (T.p1 + 27 (1 + 1,,1/1,, 1= = I. + I. Ln == Cy,0 + I.
(29)
wa* ",,u") - - -IV220(0** )0 VV ,²²),
di,kk = "200l=",'')0 V2g(2* )² D../A = G22 s
- - ,,;; BB))
(30)
dt.k Vieg(at "²th )ut V220(1" *, .// ))(u²² + V2s(=","' y²²)
Dit = G2aG( B+ )wn" -V22 ,3;; B+)u² + 22F(2*, ,y;;; Ba),
(31)
E* [ 1o , ",A,, [1]]]E[-aae y², ",ttt biik w
=
S = 2
|Bj 2 0022
07
Sg.2 =
9.2
29.2
9.2 B
S :=
(32)
B
B
B
19

---

Efficient Curvature-Aware Hypergradient Approximation
The following lemma is fundamental and will be repeatedly invoked throughout the proof. Lemma E.1. Suppose Assumption 3.1 holds, we have
11y" (a') - y- (-())] = Lg.1 || - x|l and I|u (xx) - u² (x)] & - - xx
(33)
Proof. From the definition of u² (ar), by Assumption 3.1, we have [x)]l - 11V (z, y (x))] (x))) (l2 (z)) =
41.0 Due to the optimality of y (a:), we have V2g(x, = O. Since g(z, -) is j-strrnglyy convex and g(: y) is -smooth,, we can deduce that
ully" (0) - y (x')ll - V2g(s, )(()) - V2g(a, (a'))l
VVgg(,, (z') V2g(2' "y (z')) L. 11 2'11,
which implies the first inequality in (33), Next, by the strong convexity of gzz, )) and smoothness of g, f. we have
(22g(,, "x()ii 2 Allw" (x) (z('')l
(34)
and
(vffftto(()))-Vff(f(f())(())))( (V Vi2((, y (((((,,((()))))))))))( (z))] $ (L111 Lg,27u) L -
(gy)))))) y (x))))
(35)
By the definition of w" (x)), the following equality holds
V22 y (x))(u((() - (²())) = Ves(a,y - Veel(E, (x))x²(z) + Venu(', ( (z')) (2') v=l(a',v'(z'))
Taking the norm on both sides of the above inequality and combining (34), (35), we can deduce the second inequlity in (33).
E.1. Detailed Proof of Theorem 3.2 and Remark 3.4
We first establish the descent of (<²"), Ily y" (a*))l and ||ue - u (1²))1 respectively.
Lemma E.2. Suppose Assumption 3.1 holds, the sequences generated by Algorithm 2 satisfy
(1*++1)
<2(x*)
az (2 2
at 2
Loak 2
||
+ Laxlle" y" (=x))) + L2g axlle* u 11
(36)
Proof. From Lemma 2.2 in (Ghadimi & Wang, 2018) we know that (xx) is Lg smooth, and L.p is defined (17), then we have
Lap
((xht1) S ( (V k+1 2²)
[| 111
2
(tx*tt a des) + Tu
|
= ( + as
( 12
ak
V
akk
||ab!1? +
Laa 2
||ot
(37)
2
2
2
2
where
|V(x**) dل'l < [15(t* " - Vi , + 11 w/a)
Vielatt w2 (x*1))u""(>
Vi2((> yu² (1*)]]
$ (T,y,1 + I lly 3/" (s*) + I.gallu* - u (a*))))
(38)
Then combining (38) with (37), we can obtain the desired result.
20

---

Efficient Curvature-Aware Hypergradient Approximation
Lemma E.3. Suppose Assumption 3.1 holds. When we take T 2 ((K), the sequences generated by Algorithm 2 satisfy
1.9.2 11 (x*1 - 1²11² + 1. || x011 + 7 +11/* 32 (x=x(((*)) (x*))ll 2/
11y - y² [1]+
(39)
||a/t1= - ui (2++ - 1 I|ut ui (x**)|
5/1 11 32 (2*))]
+ I'
Lim || - xz1||
(40)
Aye
Proof. From the iteration of Algorithm we know that
y² c 1+ - [I no V229(2 (0))))))""())))))))) bit - no Ver(("," )²)
Combining V2g(2*, 10 (x = 0, we have
y - 1l1 (00))) = -
n V20( }²]] (0+ - r -
- )²)). (V y (x*)) - V2g(e* (oole","")) Ve2g(t o (2)
Since V229(1, -) is - -iipschiz continuous, Lemma 1.2.4 in (Nesterov, 2018) implies that
Vag( ,y (x*)) V2g(a", Viag(a (y -\") s Lg 2
11 (x* - (x*))=)))) 112
then we can deduce that
18 - 1 - (=(x")ll S || 7/h V220(0 )/))) "lppll1" with - 2² (2**)] I11 +
1. (x* y²112
2
$ (1 - ml)))ll** wtk - (1*))) +1+ 12227 (x*t) y²1
(41)
2
where the second inequality is because -(x,,) is p-stroogly convex and 7k S L...
Then we can calculate that
11g+1 1² (1("')] 1105* T. y² (1**\|] S ()xa) )² "10* o.o yy (xx")] + Lg.2 I1y" (a*) 111
24
Laa
(1 - - 7 220(=* .,²) - (x*
lly (a*) - 111
2/c
L Ily" (2*ttt) 1²112
(1 ++ v (a.
2/1
S I1 - y² y² + Lg
11 (x*)) y 11?
(42)
2μ
where the second inequality is derived based on the same reason as in (19), and the last ineauqlity holds when T 1/2 (k)). Finally we have
>
1/1
[n(1-19))
11
++1 - 3²
51xx )| S 11
k+1 1/" (2+1)] - 11 (x (1 - 10" (a"ti) I|
(42)
1
11 (>(*>\]1- + Igg.2
11 (x=*) - 8r 112 + 1110(( (1=+1)= - 3¹ (-^*))
4
2p
(33) 1
L (a*k) - 11112 + L 1 - -2*11 2*11.
Σ 4
1
3² (a+t)] +
2p
Similarily, from the update of wi we know that
||u² w +1. - u (x*, y || =11 [T -
- V V220 ,00] (ur² with - ui (2 (
<(1 - (177))))l++ with - 2i ,
(43)
21

---

Efficient Curvature-Aware Hypergradient Approximation
This inequality implies
||u u- (x )² = |/"² Th u (x* "²) S (1- p/)) w u² (x 3))
= (1 - mmx) (1-1ys)² 11 Y V22 Dr + nhe V2f u² (zt y")
= (1 mmx))1 11 [T 7/ (22 ,r²)] (u* - - (a w"""
2 (1 wrys) ) 1|² u" (x*) ) (1 + (z (x* yt) ) u2" (2*)
1
$ 11u* u ('(x") || + 1u (x" 5²¹) 2i (2*))]
(44)
where the last inequality holds when T 2 In pap
In 1/4
1. Then we have
||u - u²' (1*++1||
|| <
k+1 u²¹ (x* ) + 11ue (x*) - u² (x* y² + [(r - u² (x*++1))
5
- - 110 - 2² (a*>)l || + Ill0(x* (at, )² u² (x*)| + |||0"( (2ltt)) - 22 (2*)))
4
(21)
S 7 |100 - u (x + 5L1
llu" 5" (xxh)l + Lu I -x r+11-
Ap
p
Hence we complete the proof.
Next, we prove that llv* - (²(x))l and I|u* - u- (x*")| have uniform upper bounds.
Lemma F.4. Suppose Assumption 3.1 holds. Take T 2 (1:) and
0k S min {
112
1²2
5pL,
1. O(1-8),
مو 2012 n.1 Li
BLabz 4gg, Lu
we have 111*~ 10' (x*))| S 2h and 1|ue - u² w²(x*))l )1 S zhg
3h1 for all S.k S K, then
1 "ll0* - v² (x*))1 + Lg || || 1/
11 - 3² (x s 2
(45)
Proof: From (26) and || (x*)) 11 S ru we know that
||aa|| = fa* )²t Viz , )n s Cy.0 + T. 11* - u -x()a\] + I.g. S To + I.gallu* u (x*)))] (46)
For ||a. - x$1 || = ax ||||dt ll. combine (46) with (39) and (40), we have
1 3¹ (1)+ -
L.. l1l" (x*) 112 L0,10k (L2 + Lg 11t - -u Illef= - w2 (2*))]
240
1
|lutt - u1" (2++1)]1 S - - u²' (x + 511
Luok (L2 + Lo.lllu* - u² (2)))).
-11 - n² (a*))] + 4u
11
a If Ily* - (x*))" S 2hyy , Ilut u² (x*)))) S 5h2 and ax I min 8hgghghog 201 g. Li 5L, Shgbo : TZg. Ln } we can 2h.y,
4²
4²
Sul,
deduce that ||g - y² (x S z ||u u² (1*+1)]] 1 shn By induction, since 11 v=("")] S =
2ha.²
||u0 n² (x*))l < 5Lr (Lemma D.I), we have 11y - 1 (r(*))]
5L1 for all 0 < k < K. S Zt. I| - u (x*")l s 21,
2h,
Moreover, we can get (45) by taking 1 - - (*+))] V 22
|_ in to (39).
Then, can obtain the results in Theorem 2.
Theorem E.5 (Restatement of Theorem 3.2). Under Assumption 3.1, choose an initial iterate (y u. xx)) in BOX 1 that 5h
satisfies 11y" y" (x*))l < min { ( 2 and ||uu u² (x")]l min 7h = 7 S 1/ /.g. there exists a proper constant step size
LLn
}. Then, for any constant step size
2bg
a S min
5pL1 Lo LL.L2Lgg 4Lg La
P
L
- 0(k
4L. 4L2,
12
9.1
and T 2 (k)) such that NBO-GD, as described in Algorithm 2, has the following properties:
22

---

Efficient Curvature-Aware Hypergradient Approximation
(a) For all integers K 2 1, minos [1 (V 2=(=")-24"' 240) +1 = =0 ). That is, NBO-GD can find an e-optimal solution ā (i.e., [+] in O steps.
ak
(b) The computational complexity of NBO-GD is: O /c) gradient computations and Jacobian-vector products, and ( /e) Hessian-vector products.
Proof. Define a Lyapunov function
Vi = f(=" (z*))) 4- + bbyll** - y (x*))) + bullu² u' (x^))
where by = b = SL, When ag satisfies the bound in Lemma E.4, because of Cauchy- inequality and xkt1 x² and inequalities (45) and (40) become
+ g o
1 82 (2 (2)
1|2 - 11oe - vx"
||
2
8L7
4L3a2
- -i | (x
< 8 110
2s )1 +
11 - y (12)
+
12
v²
Combining (36), we obtain
V4+1 Vx 5
03k a6 [V((=))
akk
Loak
4b,Liaaz
2by 122 o2
|1
2
2
112
422
2/4 1,
by
8b,L9
110/ - (x²)))1
7bb
1
1
u² (2*)))
2
122
8
Further controlling on S min|
1 L1
n²
}. we have
Tho :
TIZ 644, L3
O%
Loa
4b, 12
2by
g.g10 0k
2
by
Daax
1.
2
by
Thu
1 9,100
2
Thu 16
2 2
2²
122
8
2
122
8
8
Then, by telescoping sum, we have ming / (2*))) S 1 ER-
S
2V
0 where Vo
Kk
(x=)) 6- 4L111y" y [12 + ش²= |u0 - u² (2 S 2x)") d because from Lemma D.I we know that llu0 v (x")) - के and | u' (1 4.
Additionally, we calculate the complexities of Algorithm to achieve a stationary point. To achieve min (| S e, we need gradient computations: Ge(c) = O(K + No) = O ( Matrix-vector products: MV(() = OKT + Qo) = 1t).
Moreover, we give a brief proof of Remark 3.4.
Proof. (Proof of Remark 3.4) For y and nt if we use CG solve (8) with initial point w11k 0, we have
11 y² (1* S 1 h+1 3² (**\)) + 11y (a (a*+1) y (x**)]
(33)
11y* ("(x". - x² (a*t)) + 11o vs (a* 2/ + Lg || =x
L 11 y" (1* + 2//l(
V 1 T+1 1100 1,k 1¹ (a*ttx)) + Lg.1 ||a x*
V + 1
1
2/
< T.g.2 111" 12 + Vhllgg.1
VRi 1
T+
y² (a*))) + I.g.1 11.2111/1111 x
2/
/
V +
where the third inequality is due to (6) and the linear rate of CG (refer (17) in (Grazzi al., 2020)), and the last inequality
folloes from -1k and ||o" (a+ ) = 11(V229(0* (V y*) - V2g(x y² (²(a***))))))))"))" (a*))) || s La Ily* y (a: )|
#
23

---

Efficient Curvature-Aware Hypergradient Approximation
When we take T = (og/log((s))) the above inequality becomes the same (39). For we and wx, if we use CG solve (10) with initial point w w-1k = 0, similarily, we obain
1 | - u² "x"", ,1) ||
w T.k: - u² "x"","""" ,)""'l
-u"
|w w" (x* " " x()) S 22/1l(
v ) ) 110 - w w"(x"",,"," (a",,", w²)ll
Vk 1 T+
=2 Vn- - 1 T+1
11 )
u² (x* ,/
V +
22
Vk - 1 T+1
1| 1 u² (x*))\ 2
-k 1 )
Vk -
|| ||" -
n- (x*))-
Vk +
vk VK+11
When we take T = logkk))))) we can get (44).
The proof process following (39) and (44) is the same as that of NBO-GD. Therefore, if we replace GD with CG in Algorithm 2, the convergence rate is minosksk-1 [[ S k E (>* < 2// Additionally, in order to achieve mingose ( I| S e, we need gradient computations: Gc(e) - O(K + No) = O(x2 ()) Matrix-vector products: MV(e) OKT + Qo) (13 log Re t).
E.2. Detailed Proof of Theorem 3.7
Here, we provide a detailed proof of Theorem 3.7 for the stochastic setting. Firstly, similar to the deterministic setting, we establish the descent of E[ E [lly* y" (x*))1 and E [||u² u² )13)
Lemma E.6. Suppose Assumption 3.1 and 3.5 hold, then the sequences generated by Algorithm 4 satissf
La
B[0(x || S B[0(2*)) oh [1] @
E [1 + 22 E [llx* - v² (a*))13)
2
+ (212 + 45 )a |||u* n² (x (1²) + (25 + 45 cn.
(47)
Proof. Recall that E+ = El-- / )² from Assumption 3.5 we can deduce that
E* [1I - [3]
<2R^ [IV] F(=*, n² B3:) - VIs("", [()",,"))]]] + 2F G . BB) Vingls "....//)) )1 |/
5 + 25 |11 < 28, + 25 (2|1e* n² (1 + 2||u² (a*)
$25 + 48g 13 4 u² (I²)
(48)
where the last inequality follows from ||u S For x = w cna DJ. replace d with D in (37) and take expectation, we have
EE (2 k+]]] = E [(*x))
+ ak E |1/ -
DEP - o [V
(
ak
L
[ 11),
(49)
2
2
2
where
E [[[ - 1 < 2E [|] - 1 + 2E [|I - d
(48)
4L2E [11"*
(2 || + (4LY, + 858 E [||u* - u (x+) 45,,1 2
(50)
Then we can get the desired inequality by putting (50) into (49).
Lemma F.7. Suppose Assumption 3,3.3. and 3.6 hold, if Sq.2 s afe then the sequences generated by Algorithm 4 satisfy
5r
20.9.2
E [1 - y² (1
<
9.2 2
(E
[11>* (a** - y²th 1112) + r5(5(1
10tK 271
[11** 2² (2*
2
12
1r
9.
2
E
||P + 5r(1
176
S 9.1
20r
(51)
2
2
12
24

---

Efficient Curvature-Aware Hypergradient Approximation
Proof. From the iteration of ckk we know that
11 1 32 (x") || S 11 uik + txy - (1(x") + 7M >al1 Di - d::11
(41) S (1 pm))ll - olt - (x*x)) + L
Lg 11 (10x(x*) y0 + 7h Dtt - d7,71l
(52)
2
For D
d ||, we have
(IEL [IID - d/hh [1]1]) [1]2 S B [IDA - de:
|1 ,, Bl Ve2g(0*
- - . BB²) + V29(2 " . {1))
<E
G122 / Bl V22g(0* (xt /¹ )|] |1a + VV y BB) V29( )0))
<S%72110441) + Ga . Bb) V2g(0 y1
Then by trangle inequality, we obtain
Eitk |1D de*1Il S v S9 1l || y² B3) 229 w²)]| (x*)) y + V2G , ,²; Bb) Vag(e* ,0)"
<V 110 gtt y² (x*)] + / sp - yy
(53)
Take expectation on both side of the above inequality and combine (52), we have
E [llv* t 3" [1+)
(52)
3 $1 - (1+)))BBB [lle** - wik (x*))[]] + Lg
E [11x"(x*))) r*1 + TB [ID - -d,+11 [1]
2
(51)
S (1 - w7 + va 3 |||e* - giik - y² (a*))] + 1 E
y 11
2
+ 7/ V Sq. [llu** * of [a*))]] + + [I] G( .v:; Bh) V2g(a* y w0
For Sy, S 2² 4, we have 1 - wry + nk / s492 s - mT this implies
E [|1y h+1 - y² (2 [|]] = E [110* T. v (0 + R [1ly" (1lly"(( (1 - y [1*))]]
1922 E [|l"" (xx²) - 10 2 2/ s4
<(1 17 T |||s* o. yy (x*))] +
E [II - y² [x+)
1
2
+ [lly* (x (1) - y² (1="²²()))))) [1] + μ
V V2 v ; BB} - V29(2* |}]
and from (33) we know that IE [Ily* (a*+)) - yn [z+))]] - fui EE [|| - I|] Following (Ein a,)a,) Σm= af and (E[X))² v 2] we can deduce that
5L2 (E [llv* (I 11???
(E [|l y (2 (a*+))I)) )² < 5(1
17
FTT EE [lle* a - y 1²] +
9.2 12
2
5L2
E [|| - 1011 + 20.9.2 .
20
+
E |||a* v/" (0 || +
E || * .vAB
V2g(a* "or((=".""") 11
n12
n12
112
5L3 2 1²2
(L [1ly" (xx - 11 2 + 20.5922
<5(1
1ik 22T+11
E |||y* - y¹ (x**)]
122 I1l** - 11
2
+
5
g.1
E [1|a
x²1 111 + 5(1
1ilk
2T 25
+
+
20S
112
2
g.1
112
where the last inequality is because
E [|Va ,/:: BB) Vag(a* ..),t,,)]]]]] EE [1+
[IVZG(=* w/: Veg(a*, v < S.
25

---

Efficient Curvature-Aware Hypergradient Approximation
and
E I1ye -o --=xx))11 = 7/8 V2G( BB) y" (x 112)
similar (25) S
=E [E* [lla** - 7h V2G ,/: Bb) - v (s [|
(1 w7x) E [lle** - y" (1 (x=) + 5
Then we can complete the proof by combining xk+1 al = Da and (E [lly - 3 y0 (2a++))]) )²
E [[l y" (**+))) from Assumption 3.6.
Lemma E.8. Suppose Assumption 3.1 and 3.5 hold, if S9:2
< e.e g.
we have
IE ||1u u0" =x(x**+)11
485
4852 112
<66(1 (61
10th 3
+
+ 12(1
STSsgsgg 3 7k +
E [1lu** u= (x**))11])
2
112
2
485 122
3L
3 +6(1
E
[1ly** th - . (x ||a]
E [[la
- 11
n12
2
112
+ 6(1 yT
s
+
24
+ (12
M7lk yT S. 2k
4859
(54)
2
2
12
Proof. The iteration of witk implies that
1
11u² w u²( (x" )/ = (1 + mrx))ll² utt + nh u "a""
22 + (1
+
-222112
(43)
< (1 + +1)1)))1111 - (1)x) Pl1"" with - u (xt sn)) )11 + 27k Dit - dt.h1p?.
(55)
For D - d'all we have
Ek [1] Dr.k dt., = E* [ELL D1 - |:]
< Lt. [1 G .ed,t." ; B V22 3 |)1 || ||]0 + [1] DE + d
SS.."
4 []] + E² [1D d
SS9.1
1B [1 1/2] + Eett [1[ n² Bt) V220(0" (oo(==",~)))0 M1 ( J4 22f (2((="."'))11
/ | 2 1 112 + 2
<25
E*
|||u* - = u (z /² >| + 25 11u2 - (=*,,"²)"1 + 2/ + 2S
Take expectation and combine (55),
E [1|u** w0 1-1,k u² (x* 2/²² )|2]
5(1 - u7 + 15 7h )E I||u* - wi ui (x* (1))))1) + 4 72 E |1|ue - u² (x* (02)
+ 45% 7k
E [1] + S
Since Sg, S n² 8 we obtain
E |||uA11 n" (x* y )11c = E [llut T.t. "x )13]
TE [||u² wi - w² (x" >/? + 85
8.5
85
<(1
112 EE |||u** - 21 1x y |?] +
E
[1|u [|]]
+
2
n12
112
<(1 <(11
1 T+2 +
lrlk
8.091
E [[lut ui (a* )² + (2
17kk TT S 7k
8.59
E
[|| |12
2
112
2
112
17 yT TinyySS,,3 5 + 2
8811
+
(56)
122
26

---

Efficient Curvature-Aware Hypergradient Approximation
where the last inequlity is becuase
E || wh u² (x" ./)) E [E* f Dx u² (z " [[]]]
S [||u* - nnd u (x ) [|]] + + | [|| - d=]
SE ||[[ 7/ V220(0" (x - -=""""""""
+ 2 |||u~ + 25
<11 - -uxx) [||u* - u (x" ,/a",//"'' - + 2 [|] + 25
Then by (Ei= a,)² < n -m= a. we have
E [ [| - 2 (x + ||a]
S3E [Ilu* h- - u² (x" v/ ]1
+ 3E (1" -
- u5 (x*))) + 3E [[lu~ (a
u²
[1]
(6)
$
(3(1-
M74 yT42 +
245
E [1Iu* u (x" (x".""))) ) + (6(1 - m ST $9
7k +
245g.2
[|| [|]
2
42
2
112
312 +
E [lly** 1/ (1(x²") ]1(]]
+
3L2 112
E [|a - 111 + 6(1
17ik YT f +
24
112
2
11
<(0(1- 6(1
47kk 2
T+2 +
485
+ 12(1
1Ikk
T SS2227
485 2
E [||u* - w (x²))
112
2
112
3 + 6(1
1ik TT2
+ 485
E |11y* 3² (z M1 +
3
a
2
E |11
2
122
12
+ 6(1 Hiik y s +
24
+ (12(1- 10l yT
y sg: 7k
485
2
n.2
where the last inequality follows from
E [llu* u "=*,,M"""""
|
s 2E [||u* - - (1*x)\))
+ 21E [x* .x²) - u=" (=*"))1
2L9
s 2E |||e* - u )1*))
E ||le" v=(x*)1111
n
and E [[]]~[[]]
2I [||u** u² (x*) + 2r
Next, we prove that there exist constant bounds for both E [lle* y [1+ and E [llu* - u (==))
Lemma F.9. Suppose Assumptions 3. 1, 3.5 and 3.6 hold, choose T, stepsize on and batch sizes satisfy the conditions
T max(In(1/ (00(11/400)), (n(99)6) [K)),
40r)
((11-yrk/2)
on S ā = min {
1
112
1.2
6 L0 } O(x -3),
L1b
16V2LuLgg, 8/330 L3
80r L22
112
Sy,1 S min| {
12], Sq,1 S 3360r 1,2,2"
37
9.2
s < min!
112
n22
}.
2
< min 735r 2,2""
171²
1ra
}, sg.² S min (Lg L3
160r'
768
768
/ J.
then we have E. [1ly** y² [1+))|]]]
2²
E [||u**
n² (x+
11 for all 0 < k < K.
20r 20
Proof. Take sg² < min (L2,1: 1 ). Sy I L3 from (46) and (48) we know that
E [[ [[] < E [[ + E [IE* D ][] < (2L3, + +5g - u² 22+x+)112+ 213 + +25 + +48g.222 r
.1 1| - u² (x*))12 + ST3.
(57)
27

---

Efficient Curvature-Aware Hypergradient Approximation
Take T 2
max(In(1/vv0r In( (/966) ST.1 & min {
J.
y²²
combine 7/k S Z,
and r 2 1, then (51) and
1111-17//)
160r 768
768
(54) becomes
E [lly +1 - y¹ (z [1
5r Lr 9.2 TE [1ly" (xx*) (x=*) "²]))) 1 [1l y" 3
Lr
o
21r
9.
E
[1D
(58)
112
2²
1
and
E |||u w² (x/t111)))11)
4
f,
49.5%
$ 4 [llu* u- (x*) [13] +
E [lle* - (2*)) +
E
[1D
+
+
ru
(59)
4²
Put (57) into the above inequalities, we obtain
L2 2
E [1l k+1 3/ (2 (x )| v
EE (1²² - )
122
30r L
a
40r I
Lza
21r
E [lle** - 50 (=xx))))
+
9.1
E
[llu* - (x+))11
+
(60)
2
1
and
E [lluAt1
w¹ -²(x"t1)11 }13
4EG 112
E Ille* - v²(x))1 13
25S 4
49.5%
<(
E |||e* u² (x*))
r
(61)
4
42
122
42
For s u outputed by BOX 2 satisfy E [1lv0 y²xx") S
وي and E [[|ux
u² u²'x
S
4L3
then we complete
20r
5r
the proof by induction: If E [|ly* y" (a* [13]
20rL2,1 _
E I| ur [++
41.2
and
5r
1
n²2
4²
0k min
(62)
6622LLLg1 L,330L,,' L2
80r L L
6v110L Lu
9.1
17
1²²
Sy11
Sy, S 3360r L2
S
(63)
375
735r
1,2
r2:
8.2
9.2
9.2
we can deduce that E [[ly*+1 y (a*ht))11
S
2
and E [|lu ut (a*
S
1L 5y 6.2
from (60) and (61).
20
622
Finally, we can prove Theorem 3.7. Theorem E.10 (Restatement of Theorem 3.7). Under Assumptions 3.1.3.5 and 3.6, choose an initial iterate (g) u. x²) in x²2
BOX 2 that satisfies E | y² (1)) I min { 23, 41, and E [|lua u² [x]) I min { 5r n2a }- Then, for
41,3 L1
any constant step size 7 7 S Lg there exists proper constant step size an a O(x and T O(x) such that NSBO-SGD, described in Algorithm 4. has the following properties:
(a) Fix K 1. For samples with batch sizes |B** 0(K²), B 2 OKK x2), B 0 + k²), |B (K) [B1 2 Oxx=1K), ), it holds that mino E
x(x That is, NSBO-SGD can find an V = 0
e-optimal solution in K = O -1) steps.
(b) The computational complexity of NSBO-SSD is: = gradient complexity for F, O(Ka gradient complexity for G. O Jacobian-vectes product complexity, O 22) Hessian-vector product complexity.
Proof. Define a Lyapumov function
Vx = E [f(a 4² (=+1)] x- 4-+ b, [11x* v² (at))1²) + b,, [llu² u (1*))]
28

---

Efficient Curvature-Aware Hypergradient Approximation
where by = 41, by = 821" the inequality (58) becomes
If the step sizes and T satisfy the conditions in Lemma E.9, for E
[lle* y (x*)
S 20
5r Lr a |lly y² y" (x - 22 I1 - y [2
21r
+
11|
(64)
112
1r²
Combining (47) and (59), we have
ak Vx+1 - Vi 1 2
E
|||
|| - (ax
L
5r Lr
3b
L [|]
2
2
4²2
122
by 2
4by T
3b
21
E [1la*
'/" "²(xa")13]]
21
1
45% 0kk
E [lIu*
2n (x^*)]²]
4²
4
9.
21r
25 1
496, S2 2,2
+ (25 + 453- ok +
)
(65)
n
112
On the basis of conditions in Lemma E.9, further add
112
L
n
og min 4L4 160r L, 12 9.1
O(x
Sg.2 <
(66)
3L2
128
then we can deduce that coefficients in (65) are all positive:
a3.
Laak
5rby 122 a
3b, Lb,L2a
2 8 2
as by
2L
4bnl2 11
2
by
3b
21 21, az
45g ax
2 by
2
2
2²
n2
4
4
4
Take og as constant step sizes a, have
? e
[[]10(>>>))9 < V/-
Vie11 + (28 + +45g
21rbyS 112
25h
49b,, 22.2
41
1
By telescoping we obtain
K-1
min 0sksK-1
E
s Σ k=0
< aK 4S 888. va 42rbyS
2V
50,,
98b 112
+
12a
12a
where Vo < E [f(x, y" (x*))] - + 2 because of Lemma D.2. Take Sy.1 = 0( Sy.1 = 0 ju S? ) )1 0
( ) we have minos E [|]v((>*)1*) = 0(*)) Combining the conditions in Lemma E.9,
we can deduce the batch sizes from (32) as follows |Bfx1| 2 (kx)). |B|| > O(sK + x²), |B|| > 0(* + x²,, | 2 (kx-K) |B 2 (<=-K))
Moreover, to achieve mino K- E (V<))) S need to choose K 2 e(=²), ) and B ( |B > O
B 2 0 B > 0 B 2 @ Then, we need the following complexities: Gradient complexity of
F: K 1 - Qal Qal B = ((xa--) Gradient complexity of G: KBB + Tol Bol = O 9 -2): Jacobian-vector product
complexity: K K|B
0 -2;;;; Hessian-vector product complexity: KT BT + K|B + Qol ( O
29


Reviews:
[{"review_id": "A3oq5XD7l7", "reviewer": "Reviewer_e2rS", "summary": "This paper introduces a Newton-based approach to efficiently compute hypergradients in bilevel optimization. Instead of directly inverting the Hessian, which is costly, the method approximates Hessian-inverse-vector products (HVPs) to improve computational efficiency. The proposed Newton-based bilevel optimizer (NBO) works for both deterministic and stochastic settings and provides theoretical convergence guarantees.\nThe main idea is to use curvature information from the inner problem to improve hypergradient estimation while keeping computation manageable. Compared to existing first-order methods, the approach reduces complexity and converges faster. Experiments on meta-learning and hyperparameter optimization show that NBO requires fewer iterations to reach comparable or better solutions.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable comments and suggestions. We will address each point in detail below. The referenced figures are compiled in a single-page PDF (containing only figures) available at the anonymous link: https://drive.google.com/file/d/1ZEzn2mKcwrPzlBeziFpC1mKGWDeyq6-C/view?usp=sharing\n\n## Methods And Evaluation Criteria:\n\n* **Combine with Bregman-based methods:** \n  * The improvements of our NBO framework and Bregman-based methods stem from two different perspectives: the NBO framework improves hypergradient approximation (related to updates of the lower-level variable $y$), while Bregman-based methods use mirror descent to update the upper-level variable $x$.\n  * Thank you for your constructive suggestion. The idea of combining our NBO framework with Bregman-based methods is both interesting and promising. To explore this, we conducted a numerical experiment where we used the NBO framework to approximate the hypergradient, followed by updating $x$ using the mirror descent method. We then compared SBiO-BreD in [1] and SBiO-BreD+NSBO (i.e., using our proposed framework to compute $w\\_t$ in SBiO-BreD) in the context of data hyper-cleaning. The result, presented in Fig. 1 of the [anonymous link](https://drive.google.com/file/d/1ZEzn2mKcwrPzlBeziFpC1mKGWDeyq6-C/view?usp=sharing), demonstrates that SBiO-BreD+NSBO significantly outperforms SBiO-BreD, highlighting the effectiveness of the NBO framework.\n\n## Experimental Designs Or Analyses & Relation To Broader Scientific Literature:\n\n* **Comparison against other curvature-aware algorithms:** We would like to clarify that we have compared our method with SHINE in Experiments (Section 4.2 and 4.3). SHINE is a popular curvature-aware algorithm that employs the quasi-Newton method to solve the lower-level problem.\n\n## Essential References Not Discussed:\n\n* **Comparison against exact Hessian methods:** Directly computing the Hessian and its inverse is both time-consuming and memory-consuming. In our setting, using the Hessian-vector product is a more practical choice, as it can be efficiently \n  computed and stored with modern automatic differentiation frameworks. To demonstrate this, we conducted a numerical experiment comparing NBO with the exact Hessian inverse (implemented using jax.hessian and jnp.linalg.inv) and NBO-GD for hyperparameter optimization on synthetic data. The result in Fig. 2 of the [anonymous link](https://drive.google.com/file/d/1ZEzn2mKcwrPzlBeziFpC1mKGWDeyq6-C/view?usp=sharing) shows that NBO with the exact Hessian inverse is significantly slower.\n* About comparison against Bregman-based methods, please see the response to \"Methods And Evaluation Criteria\".\n\n## Other Strengths And Weaknesses:\n\n* **Scalability to large dimensional problems:** Thanks for the suggestion. As shown in Fig. 3 of the [anonymous link](https://drive.google.com/file/d/1ZEzn2mKcwrPzlBeziFpC1mKGWDeyq6-C/view?usp=sharing), we evaluated the scalability of NBO by testing hyperparameter optimization on synthetic data while progressively increasing the problem dimension.\n* **Extension to non-strongly convex problems:** The NBO framework is primarily designed for bilevel optimization problems with a strongly convex structure. For non-strongly convex problems, existing methods typically involve reformulating the original problem to incorporate a strongly convex structure. Once this structure is established, our NBO framework can be applied. For instance, we compared the BAMM method [2] and BAMM+NBO (i.e., using NBO to compute $d_x^k$ in BAMM). The result, presented in Fig. 4 of the [anonymous link](https://drive.google.com/file/d/1ZEzn2mKcwrPzlBeziFpC1mKGWDeyq6-C/view?usp=sharing), shows that BAMM+NBO significantly outperforms BAMM, highlighting the effectiveness of the NBO framework. Due to space constraints, we primarily focused on strongly convex lower-level problems in this work, which limited our ability to fully showcase the versatility of NBO. We appreciate your valuable suggestions and will include additional discussion in the revised manuscript to further demonstrate the versatility of NBO.\n\n[1] Huang et al., Enhanced bilevel optimization via bregman distance, NeurIPS 2022.\n\n[2] Liu et al., Averaged Method of Multipliers for BiLevel Optimization without Lower-Level Strong Convexity, ICML 2023.", "claims_and_evidence": "The paper’s claims are mostly well supported by empirical evidence and theory.", "methods_and_evaluation_criteria": "Consider Bregman-Based Curvature-Aware Methods:\n- Would Bregman-based methods achieve similar improvements?\n- Could a hybrid approach combining Newton and Bregman techniques be better?\n\nSee, e.g.,\n- Enhanced Bilevel Optimization via Bregman Distance by Feihu Huang, Junyi Li, Shangqian Gao, Heng Huang\n- Online Nonconvex Bilevel Optimization with Bregman Divergences by Jason Bohne, David Rosenberg, Gary Kazantsev, Pawel Polak", "theoretical_claims": "I did not review the proofs.", "experimental_designs_or_analyses": "I suggest to compare against other curvature-aware algorithms.", "relation_to_broader_scientific_literature": "Lack of comparison against second order methods.", "essential_references_not_discussed": "It does not benchmark against full second-order methods (e.g., ones that directly compute Hessians rather than approximations).\nIt does not compare against Bregman-based methods.", "other_strengths_and_weaknesses": "The paper is well-structured with strong theoretical backing, but it could be improve by discussing the scalability of the method to large dimensional problems and how well it performs empirically when the strong-convexity assumption is violated.", "other_comments_or_suggestions": "No", "questions_for_authors": "No", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xQe00Euao3", "reviewer": "Reviewer_Y7jC", "summary": "This paper proposes a novel method for bilevel optimization, focusing on improving hypergradient estimation by incorporating curvature information. The key contributions include: (1) New Algorithmic Framework: An enhanced algorithm using an inexact Newton method, with improved computational complexity and convergence guarantees in both deterministic and stochastic settings. (2) Empirical Results: Numerical experiments demonstrate significant performance benefits compared to existing gradient-based methods.\n\n\"## update after rebuttal\"\nThrough the rebuttal, I believe this work presents a general and efficient gradient-based framework for bilevel optimization, which can be extended to the non-strongly convex setting using techniques such as the BAMM method. The motivation and approach are highly interesting and are supported by strong theoretical guarantees. Therefore, I recommend acceptance.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable comments and suggestions. We will address each point in detail below. The referenced figures are compiled in a single-page PDF (containing only figures) available at the anonymous link: \nhttps://drive.google.com/file/d/15xqtvUMRk7Ah7Gi5hnvBQWyuZ6W15zX8/view?usp=sharing  \n \n## Other Strengths And Weaknesses & Questions For Authors:\n* **Extension to lower-level non-strongly convex problems:** The NBO framework is primarily designed for bilevel optimization problems with a strongly convex structure. For non-strongly convex problems, existing methods typically involve reformulating the original problem to incorporate a strongly convex structure. Once this structure is established, our NBO framework can be applied. For instance, in BAMM method [1], when $g$ is merely convex, an aggregation function $\\phi_{\\mu} = \\mu f + (1 - \\mu) g$ is defined, which is strongly convex when $f$ is strongly convex, then an approximated hypergradient $d_x^k$ can be computed by replacing $g$ with $\\phi_{\\mu}$. We compared the BAMM method and BAMM+NBO (i.e., using NBO to compute $d_x^k$ in BAMM) on the toy example (13) in [1]. The result, presented in Fig. 1 of the [anonymous link](https://drive.google.com/file/d/15xqtvUMRk7Ah7Gi5hnvBQWyuZ6W15zX8/view?usp=sharing), shows that BAMM+NBO significantly outperforms BAMM, highlighting the effectiveness of the NBO framework. Due to space constraints, we primarily focused on strongly convex lower-level problems in this work, which limited our ability to fully showcase the versatility of NBO. We appreciate your valuable suggestions and will include additional discussion in the revised manuscript to further demonstrate the versatility of NBO.\n\n## Other Comments Or Suggestions:\nThank you for pointing this out. We will revise this citation typo.\n\n[1] Liu et al., Averaged Method of Multipliers for BiLevel Optimization without Lower-Level Strong Convexity, ICML 2023.", "claims_and_evidence": "Yes, the paper provides both theoretical and empirical evidence to support the claims made.", "methods_and_evaluation_criteria": "Yes. In the theoretical part, the authors use $\\| \\nabla \\Phi(x) \\|$ as the stationary measure, which is reasonable for bilevel problems with a strongly convex lower-level problem. In the experimental part, the benchmark used by the authors are also reasonable, various types of bilevel algorithms are involved in the benchmark they used.", "theoretical_claims": "I have checked the proofs of Theorems 3.2 and 3.7 and did not find any errors.", "experimental_designs_or_analyses": "I have checked the experimental designs and analyses, including the models, datasets, and step size selection. Specifically, the step sizes are tuned via grid search, which is fair.", "relation_to_broader_scientific_literature": "The hypergradient method is a commonly employed approach for solving bilevel optimization problems in the literature. However, most of previous works utilized first-order methods to estimate the hypergradient. The key contribution of this paper lies in its use of the curvature information, which is theoretically proven to be more efficient by the authors. This advancement represents an improvement over existing gradient-based methods.", "essential_references_not_discussed": "I do not find any missing references.  As far as I know, in terms of complexity, this paper compares its results with the state-of-the-art works.", "other_strengths_and_weaknesses": "Strengths: \n\n(1)The authors propose a novel method for estimating the hypergradient, which marks a significant departure from classical gradient-based methods. This method leverages the unique structure of the hypergradient, where a Hessian is shared when using inexact Newton method. This idea is interesting and I think it is a breakthrough. \n\n(2)Particularly, the authors provide theoretical results that the proposed algorithm achieves lower computational complexity compared to existing methods, thereby highlighting the advantages of using curvature information.\n\nWeaknesses:  \n\nThe authors only consider the lower-level strongly convex case. (Actually, I know that the non-strongly convex case is very challenging.)", "other_comments_or_suggestions": "There is a citation typo in Table 1. The reference for \"No-loop AID\" should be (Ji et al., 2022).", "questions_for_authors": "Please discuss the possibility of extending the algorithm proposed in this paper to bilevel optimization problems where the lower-level problem is non-strongly convex.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8OJdnFKDMZ", "reviewer": "Reviewer_92Dz", "summary": "This paper consider the bilevel problem $\\min_x \\Phi(x) = f(x, y^*(x))$ where $y^*(x) = \\arg\\min_y g(x, y)$ where the inner function $g$ is strongly convex w.r.t. the inner variable $y$. The paper proposes a new AID-based method where the inner variable $y$ and the linear system variable $u$ are updated by an approximate Newton step. These approximate Newton directions are computed by applying (S)GD to the quadratic functions associated with the two linear systems. Then, the paper shows that in the deterministic setting, the proposed method improves the gradient complexity by the factor $(\\kappa\\log(\\kappa))^{-1}$ in comparison with classical AID-based methods.  The computational complexity of the stochastic variant is also provided. Numerical expriments on a synthetic problem, hyperparameter optimization and data cleaning are provided.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thorough review and valuable feedback. We will address each point in detail below. The figures mentioned below are in this anonymous link: https://drive.google.com/file/d/1v6ftNYExUb_ClkgoS7b9wU2Q3nNY1IsP/view?usp=sharing\n\n## Other Strengths And Weaknesses:\n\n* **Concern about batch size:** The large batch size we choose helps achieve better sample complexity. If the batch size is set to $O(1)$, we can obtain a sample complexity of $O(\\kappa^{16} \\epsilon^{-2})$, which is worse than $O(\\kappa^9 \\epsilon^{-2})$ in Theorem 3.7. Since SOBA is a special case of our NBO framework (as noted in Remark 2.1 (i) of our work), their analysis also applies to our NSBO-SGD when $T=0$. Moreover, a large batch size related to $ \\epsilon $ or $ K $ is commonly used in the stochastic bilevel optimization literature (see, e.g., Ji et al., 2021; Arbel & Mairal, 2022).\n\n* **Our improvements in the stochastic setting and in practice:** In theory, as stated after Theorem 3.7, our NSBO-SGD improves upon the SOTA result of AmIGO by a factor of $ \\log \\kappa $ in the stochastic setting, where AmIGO also employs a large batch size. \n  \n  In practice, our algorithms show significant advantages in experiments, even when using small batch sizes (64 or 256). Additionally, we add experiments on meta-learning using the Omniglot and miniImageNet (about 3GB) datasets, with CNN4 networks. We compared our NBO with the SOTA algorithms PZOBO in [1] and qNBO in [2]. The results, presented in Fig. 1 and 2 of the [anonymous link](https://drive.google.com/file/d/1v6ftNYExUb_ClkgoS7b9wU2Q3nNY1IsP/view?usp=sharing), demonstrate that NBO consistently outperforms the other methods on both datasets.\n  \n  [1] Sow et al., On the convergence theory for hessian-free bilevel algorithms, NeurIPS 2022.\n  \n  [2] Fang et al., qNBO: quasi-Newton Meets Bilevel Optimization, ICLR 2025.\n\n## Other Comments Or Suggestions:\n\nThank you for your careful reading. We will address each point as follows:\n\n* **Line 39:** We will revise it.  \n* **Line 131:** This sentence represents the second aspect of our motivation. Although we have not fully utilized this point, we still believe it is important.  \n* **Line 150:** Here we use warm-start for $u$ and will make revisions to clarify this.\n* **Line 303:** We will revise it.  \n* **Algorithm 3 and 5:** The initialization point $w^{-1,k}=0$ is correct. With this initialization, we obtain $w^{0,k}=\\gamma_k d_u^k$ (where $d_u^k$ is defined in Line 150), ensuring that the second equality in (44) holds.  \n* **Theorem 3.2 and 3.7:** BOX 1 and 2 represent the initialization strategies of Algorithm 2 and 4, respectively. We will add an introduction before these theorems.  \n* **Summary:** Depending on the available space in the final manuscript, we will consider adding a summary table for the stochastic setting.  \n* **Box 2:** Thank you for pointing out the typo. $B_0$ and $B_0^{'}$ should be replaced by $B_{0,n}$ and $B^{'}_{0,q}$ respectively, while maintaining the same batch size.  \n* **Line 1049:** We will revise it.  \n* **Line 1142-1143:** The constant $\\frac{1}{2}$ arises from $2 \\times (\\frac{1}{2})^2$.  \n* **Section E.2 and the first line of section E.2:** We will revise them.  \n* **Equation (52):** We will revise it.  \n\n## Questions For Authors:\n\n* **Concern about warm-start:** The reason we use warm-start for $u$ but not for $v$ is that $v$ directly affects $||y^{k+1}-y^{k}||$, while $u$ does not. Indeed, in the $k$-th iteration of NBO, $v^k$ serves as an inexact approximation of the Newton direction $v^*(x^k,y^k):=[\\nabla_{22}^2g(x^k,y^k)]^{-1} \\nabla_2 g(x^k,y^k)$, whereas $y^{k}$ approximates $y^*(x^{k})$. We update $y^{k+1}= y^k - v^k$ using a single inexact Newton step. If we further compute $v^{k+1} = v^{k} - \\gamma_k d_v^{k} $ (**using one-step gradient descent with warm-start**), according to the Lyapunov function argument, we estimate $||v^{k+1} - v^*(x^{k+1},y^{k+1})|| \\leq (1-\\gamma_{k} \\mu)||v^{k} - v^*(x^{k},y^{k})|| + ||v^*(x^{k+1},y^{k+1})-v^*(x^{k},y^{k})||. (1) $ Similar to $y^*(x)$ and $u^*(x)$, $v^*(x,y)$ is Lipschitz continuous and we can get $||v^*(x^{k+1},y^{k+1})-v^*(x^{k},y^{k})||\\leq \\frac{2L_{g,1}}{\\mu}(||x^{k+1} - x^{k}||+||y^{k+1} - y^{k}||).          (2)$ \n  \n  Note that $||y^{k+1} - y^{k}||\\leq ||v^k-v^*(x^{k},y^{k})||+||v^*(x^{k},y^{k})||\\leq  ||v^k-v^*(x^{k},y^{k})||+\\frac{L_{g,1}}{\\mu}||y^{k}-y^*(x^k)||. (3)$ Substituting (2) and (3) into (1), we obtain $||v^{k+1} - v^*(x^{k+1},y^{k+1})|| \\leq (1-\\gamma_{k} \\mu+\\frac{2L_{g,1}}{\\mu})||v^{k} - v^*(x^{k},y^{k})|| + others. $\n  \n  Since $1-\\gamma_{k} \\mu+\\frac{2L_{g,1}}{\\mu}>1$ when $\\gamma_{k}\\leq 1/L_{g,1}$, the Lyapunov function argument fails to hold, indicating that the above warm-start strategy for $v$ is problematic. If multi-step gradient descent is performed for $v$, we choose to initialize at 0 because, in this setting, when $T=0$, NBO reduces to the single-loop algorithm framework proposed by (Dagréou et al., 2022) .", "claims_and_evidence": "* Convergence rates in the deterministic and stochastic settings are provided and sound.\n* NBO-GD achieves a convergence rate in $\\mathcal{O}\\left(\\frac{\\kappa^3}{K}\\right)$ leading to a gradient complexity of $\\mathcal{O}\\left(\\frac{\\kappa^3}{\\epsilon}\\right)$ and a HVP complexity in $\\mathcal{O}\\left(\\frac{\\kappa^4}{\\epsilon}\\right)$ (which can be reduced to $\\mathcal{O}\\left(\\frac{\\kappa^{3.5}}{\\epsilon}\\right)$ by using CG instead of GD).\n* NBO-SGD achieves a convergence rate $\\mathcal{O}\\left(\\frac{\\kappa^3}{\\epsilon}\\right)$ (by assuming batch sizes in $\\Theta(K)$) leading to a gradient complexity in $\\mathcal{O}(\\kappa^5/\\epsilon)$ for $F$, $\\mathcal{O}(\\kappa^9/\\epsilon)$ for $G$, and HVP complexity in $\\mathcal{O}(\\kappa^7/\\epsilon)$.", "methods_and_evaluation_criteria": "The method is numerically evaluated on a synthetic problem, hyperparameter optimization and data cleaning. This setting is classical for the numerical evaluation of bilevel optimization methods.", "theoretical_claims": "I checked the proof for the deterministic setting and the result are sound apart from some inconsequential typos (see **Other Comments Or Suggestions** section for detail).", "experimental_designs_or_analyses": "Experimental results show the practical interest of the proposed method.", "relation_to_broader_scientific_literature": "The paper proposes an adaptation of the single-loop algorithmic framework introduced in [1], where the directions for solving the inner linear system are approximated using Newton steps. It is worth noting that [2] and [3] explore the use of quasi-Newton steps to approximate the inner solution. Furthermore, [3] utilizes the quasi-Newton directions employed in solving the inner problem to enhance the resolution of the linear system. However, unlike this paper, [2] and [3] do not provide non-asymptotic results.\n\n[1] Dagréou, M., Ablin, P., Vaiter, S., and Moreau, T. *A framework for bilevel optimization that enables stochastic and global variance reduction algorithms*. NeurIPS 2022.\n\n[2] Pedregosa, F. *Hyperparameter optimization with approximate gradient*. ICML 2016.\n\n[3] Ramzi, Z., Mannel, F., Bai, S., Starck, J.-L., Ciuciu, P., and Moreau, T. *Shine: Sharing the inverse estimate from the forward pass for bi-level optimization and implicit models*. ICLR, 2022.", "essential_references_not_discussed": "To my knowledge, any essential reference is missing.", "other_strengths_and_weaknesses": "### Strengths\n* The paper is well-written\n* The method is novel to my knowledge.\n* Potential improvements with variance reduction and momentum are discussed.\n\n### Weaknesses\n\n* The stochastic result assumes large batch sizes which scale linearly with the number of iterations. This setting does not match the practice.\n* The gain in the stochastic setting is limited while this setting being important in practice.", "other_comments_or_suggestions": "* **Line 39**: when the citation is part of the sentence, it should not be in parentheses (i.e. it should be `\\citet` instead of `\\citep`).\n\n* **Line 131**: *\"(2) Observe that v∗(x, y) and u∗(x, y) share the same Hessian inverse.\"* It seems that this sentence is not supposed to be here. \n\n* **Line 150**: *\"we write $u = u^k - w$\"*, $u^k$ is used before being defined.\n\n* **Line 303**: The equation goes beyond the margin.\n\n* **Algorithms 3 and 5**: I guess by reading the proof that it is $w^{-1},k = u^k$ instead of $w^{-1},k = 0$ in the initialization of the algorithms.\n\n* **Theorem 3.2 and 3.7**: By reading the text, it is not clear what BOX 1 and BOX 2 refer to. These things should be introduced before the theorems.\n\n* It would be nice to have a summury of complexity results for the stochastic setting, as done in Table 1 for the deterministic setting.\n\n* **Box 2**: The same batches $B_0$ and $B_0'$ are used in all the iterations. Why?\n\n* **Line 1049**: *\"$\\nabla_2f$\"* -> *\"$\\nabla_2g$\"*.\n\n* **Line 1142-1143**: I think there is no $\\frac12$ in factor of $\\lVert y^k-y^*(x^k)\\rVert^2$ anymore if the inequality comes from $(a+b)^2\\leq 2a^2 + 2b^2$.\n\n* **Section E.2 and first line of section E.2**: It Theorem 3.7 and not Theorem 3.2. \n\n* **Equation (52)**: Isn't the equality an inequality actually?", "questions_for_authors": "* As it is common to warm start the subsolvers in bilevel optimization [1, 2], I wonder why Algorithm 3 and 5 not set $v^{-1,k} = y^k$ instead of $v^{-1,k} = 0$.\n\n[1] Ji, K., Yang, J., and Liang, Y. *Bilevel optimization: Convergence analysis and enhanced design*. ICML 2021\n[2] Arbel M. and Mairal J. *Amortized implicit differentiation for stochastic bilevel optimization*. ICLR 2022", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "c88PbULgYL", "reviewer": "Reviewer_ihhr", "summary": "The paper focus on bilevel optimization and incorporates curvature information into the approximation of hypergradients in bilevel optimization.  The authors propose a Newton-based framework (NBO) that solves lower-level problems with computing Hessian inverse-vector products.  They establish the convergence rate guarantees in both deterministic and stochastic scenarios, demonstrating improved computational complexity over popular gradient-based methods.  Numerical experiments validate the effectiveness of the proposed method.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable feedback. We will address each point in detail. The figures and the tables mentioned below are in this anonymous link: https://drive.google.com/file/d/1lCY1UF3isNnoujM8AGCPIdJlRqHctj-b/view?usp=sharing\n\n## Essential References Not Discussed:\n\n* **Comparison with fully first-order methods (numerical evaluation):** **We clarify that we have compared our method with the fully first-order method F2SA [1] in Experiments (Section 4.2 and 4.3).**  Our method consistently outperforms F2SA, even when using small batch sizes (64 or 256). It is worth noting that, based on the results from [2,3], F2SA is considered a SOTA fully first-order bilevel algorithm.\n  \n  Furthermore, we conducted **additional experiments on meta-learning** using the Omniglot and miniImageNet (approximately 3GB) datasets, with 4-layer convolutional neural networks (CNN4). We compared our NBO with the SOTA algorithms PZOBO [4] and qNBO [5]. The results, presented in Fig. 1 and 2 of the [anonymous link](https://drive.google.com/file/d/1lCY1UF3isNnoujM8AGCPIdJlRqHctj-b/view?usp=sharing), show that NBO consistently outperforms the other methods on both datasets, highlighting the effectiveness of our framework.\n  \n  Note that PZOBO is a widely used Hessian-free bilevel algorithm, while qNBO is a recently proposed curvature-aware algorithm that employs the quasi-Newton method to solve the lower-level problem.\n\n* **More references:** First, among the references mentioned by the reviewer, [R1, R2, R4] focus on nonconvex lower-level problems, whereas our work focuses on strongly convex lower-level problems. Second, the reference [R3] does not propose a new bilevel algorithm but rather provides an improved analysis for SOBA. Since SOBA is a special case of our NBO framework (as noted in Remark 2.1 (i) of our work), their analysis also applies to our NSBO-SGD when $T=0$.\n  \n  We appreciate the reviewer’s suggestions and will include these references in our work. Additionally, we have added a summary table for the stochastic setting that incorporates the references mentioned by the reviewer, as shown in Table 1 in the [anonymous link](https://drive.google.com/file/d/1lCY1UF3isNnoujM8AGCPIdJlRqHctj-b/view?usp=sharing) (the second page).\n\n## Other Strengths And Weaknesses:\n\n* **In the stochastic setting, we also establish a theoretical improvement:** As stated after Theorem 3.7, our NSBO-SGD improves upon the SOTA result of AmIGO by a factor of $\\log \\kappa$ in the stochastic setting, where AmIGO also employs a large batch size.\n* **Concern about batch size:** The large batch size we choose helps achieve better sample complexity. If the batch size is set to $O(1)$, we can obtain a sample complexity of $O(\\kappa^{16} \\epsilon^{-2})$, which is worse than $O(\\kappa^9 \\epsilon^{-2})$ in Theorem 3.7. It is also worth noting that SOBA is a special case of our NBO framework. \n* **Extension to non-strongly convex problems:** The NBO framework is primarily designed for bilevel optimization problems with a strongly convex structure. For non-strongly convex problems, existing methods often involve reformulating the original problem to incorporate a strongly convex structure. Once this structure is established, our NBO framework can be applied. For instance, we compared the BAMM method [7] and BAMM+NBO (i.e., using NBO to compute $d_x^k$ in BAMM). The result, presented in Fig. 3 of the [anonymous link](https://drive.google.com/file/d/1lCY1UF3isNnoujM8AGCPIdJlRqHctj-b/view?usp=sharing), shows that BAMM+NBO significantly outperforms BAMM, highlighting the effectiveness of the NBO framework. \n* **Comparison with fully first-order methods (convergence rate and complexity):** Please see Table 1 in the [anonymous link](https://drive.google.com/file/d/1lCY1UF3isNnoujM8AGCPIdJlRqHctj-b/view?usp=sharing) for the stochastic setting. For the deterministic setting, we will add a new row to Table 1 in our paper, presented in Table 2 of the [anonymous link](https://drive.google.com/file/d/1lCY1UF3isNnoujM8AGCPIdJlRqHctj-b/view?usp=sharing) (the second page).\n\n[1] Kwon et al., A fully first-order method for stochastic bilevel optimization, ICML 2023.\n\n[2] Chen et al., Near-optimal nonconvex-strongly-convex bilevel optimization with fully first-order oracles. arXiv preprint 2023.\n\n[3] Chen et al., On finding small hyper-gradients in bilevel optimization: Hardness results and improved analysis. COLT 2024.\n\n[4] Sow et al., On the convergence theory for hessian-free bilevel algorithms, NeurIPS 2022.\n\n[5] Fang et al., qNBO: quasi-Newton Meets Bilevel Optimization, ICLR 2025.\n\n[6] Arbel et al., Amortized implicit differentiation for stochastic bilevel optimization, ICLR 2022.\n\n[7] Liu et al., Averaged Method of Multipliers for BiLevel Optimization without Lower-Level Strong Convexity, ICML 2023.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "I did not fully check the correctness of the proofs, but the theoretical claims seems reasonable.", "experimental_designs_or_analyses": "Yes", "relation_to_broader_scientific_literature": "This paper achieves enhanced iteration complexity for bilevel algorithm with strongly convex lower-level problem in terms of the condition number in determinstic setting.", "essential_references_not_discussed": "The key contribution is in quadratic programming-based bilevel methods; however, the literature on this topic is incomplete, e.g. [R1]-[R3]. Most notably, [R3], which employs a constant batch size for the stochastic setting, in contrast to the increasing batch size required by the proposed method. \n\n[R1] Non-Convex Bilevel Games with Critical Point Selection Maps. Michael Arbel, et. al. NeurIPS 2022. \n\n[R2] A Generalized Alternating Method for Bilevel Optimization under the Polyak-Łojasiewicz Condition. Quan Xiao, et. al. NeurIPS 2023.  \n\n[R3] Single-Timescale Multi-Sequence Stochastic Approximation Without Fixed Point Smoothness: Theories and Applications. Yue Huang, et. al. \n\nFurthermore, since their method leverages second-order information, it is essential to compare it with fully first-order bilevel methods to demonstrate that the use of second-order terms does not impede numerical performance. However, the paper lacks references to the relevant literature in this area. \n\n[R4] On Penalty-based Bilevel Gradient Descent Method. Han Shen, et. al. ICML 2023.", "other_strengths_and_weaknesses": "Strengths:\n\n1. The paper presents a novel method for hypergradient approximation that efficiently incorporates curvature information with rigorous convergence rate guarantees and improves rate in deterministic setting.  \n2. Numerical experiments showcase the effectiveness of the proposed method, showing improved performance over existing gradient-based methods. \n\nWeaknesses:\n\n1. The theoretical improvements are established only in the deterministic setting. Moreover, the analysis for the stochastic setting requires an increasing batch size inversely proportional to $\\epsilon$, while other bilevel methods use a constant batch size (see [R3], [R5]–[R7]). This requirement is impractical for large-scale machine learning applications. Also the lower-level strongly convexity is somewhat restrictive. \n\n2. Although the paper primarily builds upon a fully single-loop second order method, it would be beneficial to include a more comprehensive convergence rate comparison with other bilevel approaches, such as fully first-order methods in discussion and Table 1. \n\n[R1] Non-Convex Bilevel Games with Critical Point Selection Maps. Michael Arbel, et. al. NeurIPS 2022. \n\n[R2] A Generalized Alternating Method for Bilevel Optimization under the Polyak-Łojasiewicz Condition. Quan Xiao, et. al. NeurIPS 2023.  \n\n[R3] Single-Timescale Multi-Sequence Stochastic Approximation Without Fixed Point Smoothness: Theories and Applications. Yue Huang, et. al. \n\n[R4] On Penalty-based Bilevel Gradient Descent Method. Han Shen, et. al. ICML 2023.\n\n[R5] A framework for bilevel optimization that enables stochastic and global variance reduction algorithms. Mathieu Dagréou, et. al. NeurIPS 2022.\n\n[R6] Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems. Tianyi Chen, et. al. NeurIPS 2021. \n\n[R7] A Fully Single Loop Algorithm for Bilevel Optimization without Hessian Inverse. Junyi Li, et. al. AAAI 2022.", "other_comments_or_suggestions": "No", "questions_for_authors": "See weakness", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper proposes a Newton-based algorithm (NBO) for bilevel optimization, with both deterministic and stochastic variants. The authors provide convergence rate guarantees for both settings. The experimental evaluations are sufficiently comprehensive to support the main claims.

All reviewers recommend acceptance after the rebuttal. I concur with their assessment and recommend acceptance. Please incorporate necessary changes requested by the reviewers in the final version, such as addressing the scalability of the method to high-dimensional problems, discussing its performance when the strong convexity assumption is violated, and including additional comparisons to other methods mentioned in the rebuttal.
Decision: Accept (poster)