Faster Gradient Methods for Highly-smooth Stochastic Bilevel Optimization
Authors: Lesi Chen
Junru Li
El Mahdi Chayti
Jingzhao Zhang
Venue: ICLR 2026 Poster (2026)
Keywords: bilevel optimization
stochastic acceleration
Primary area: optimization

Abstract: This paper studies the complexity of finding an $\epsilon$-stationary point for stochastic bilevel optimization when the upper-level problem is nonconvex and the lower-level problem is strongly convex.
Recent work proposed the first-order method, F${}^2$SA, achieving the $\tilde{\mathcal{O}}(\epsilon^{-6})$ upper complexity bound for first-order smooth problems. This is slower than the optimal $\Omega(\epsilon^{-4})$ complexity lower bound in its single-level counterpart. 
In this work, we show that faster rates are achievable for higher-order smooth problems. We first reformulate F$^2$SA as approximating the hyper-gradient with a forward difference. Based on this observation, we propose a class of methods F${}^2$SA-$p$ that uses $p$th-order finite difference for hyper-gradient approximation and improves the upper bound to $\tilde{\mathcal{O}}(p \epsilon^{-4-2/p})$ for $p$th-order smooth problems. Finally, we demonstrate that the
$\Omega(\epsilon^{-4})$ lower bound
also holds for stochastic bilevel problems when the high-order smoothness holds for the lower-level variable, indicating that the upper bound of F${}^2$SA-$p$ is nearly optimal in the highly smooth region $p = \Omega( \log \epsilon^{-1} / \log \log \epsilon^{-1})$.

Published as a conference paper at ICLR 2026
FASTER GRADIENT METHODS FOR HIGHLY-SMOOTH STOCHASTIC BILEVEL OPTIMIZATION
Lesi Chen and Junru Li Tsinghua University & Shanghai Qizhi Institude, China
(chenl jr-1124 @mai 1s. ts taiinnguaaaa.. edu.ccredccc
El Mahdi Chayti
Bioole Polytechnique Fédérale de Lausanne (EPFL., Switzerland el el-mahdi. chayti@epfl...
Jingzhao Zhang f
Tsinghua University Shanghai Qizhi Institude, China singzhaczzmail tsinghua. edu.cn cn
ABSTRACT
This paper studies the complexity of finding e-stationary point for stochastic bilevel optimization when the upper-level problem is nonconvex and the lower- level problem is strongly convex. Recent work proposed the first-order method, F*SA, achieving the ( upper complexity bound for first-order smooth prob- lems. This is slower than the optimal (2(=*)) complexity lower bound in its single- level counterpart. this work, we show that faster rates are achievable for higher- order smooth problems. We first reformulate F²SA approximating the hyper- gradient with forward difference. Based on this observation, we propose a class of methods SAp that uses pth-order finite difference for hyper-gradient ap- proximation and improves the upper bound to (pp--12/pp) for pth-order smooth problems. Finally, we demonstrate that the (?(=1) lower bound also holds for stochastic bilevel problems when the high-order smoothness holds for the lower- level variable, indicating that the upper bound of FASA-p is nearly optimal the region p Q(log e 1 log log e 1
INTRODUCTION
Many machine learning problems, such as meta-learning (Rajeswaran al. 2019). hyper-parameter
tuning (Bao et al. 2021 Franceschi el al., 2018 Mackay et al., 2019), adversarial training (Good-
fellow al.) 2020),
reinforcement learning Yang al., 9 Hong et al. 2023 Zeng al.,
2024 Zeng & Doan 2024) can be abstracted as solving the following bilevel optimization problem:
min (xx) zcR
f(x, y (æ)), y (x) arg min g 3), yER",
(1)
We call f and g the upper-level and lower-level functions, respectively, and call the hyper- objective. In this paper, we consider the most common sencorrvesttroonbb-oovvey setting where d + R is smooth and possibly nonconvex, and Rd - RR is smooth jointly in (æ, y) and strongly convex y. Under lower-level strong convexity assumption, the implicit function theorem (Dontchev & Rockafellar 2009) indicates the following closed form of the hyper-gradien
(0x Vp(x) (m, y (x)) V2 (a y (x)) ( y (z))) 'V (xx y (x)).
Following the works in nonconvex optimization (Carmon et al. 2020 2021 Arjevani et al. 2023 we consider the task of finding an e-stationary point of s. i.e., point E R such that
The corresponding author.

---

Published as a conference paper at ICLR 2026
[[V((z)] Σ Motivated by many real machine learning tasks, study the stochastic setting, where the algorithms only have access to the stochastic derivative estimators of both f and g.
The first efficient algorithm BSA Ghadimi & Wang 2018) for solving the stochastic bilevel problem leverages both stochastic gradient and Hessian- -ector--ro-ott (HVP) oracles to find an e-stationary point of (.)). Subsequently. Ji et al. (2021) proposed stocBiO by incorporating multiple enhanced designs in BSA to improve the complexity. Yang et al. (2023)) proposes FdeHBO that uses finite- differences to estimate HVP vectors. However, all these methods require the stochastic Hessian assumption 55 on the lower-level function, which is stronger than the standard SGD assumption.
Kwon al. 2023 proposed the first fully first-order method F2S, SA that works under standard SGI) assumptions on both f and (Assumption [.I]. The main idea is to solve the following penalty
problem (Liu et al. 2022 2023 Shen & Chen 2023] Shen et al. 2025b Lu & Mci 2024 2026):
min xr #E RERR
f(x, y) + x g(x, y) min g( x)
(3)
ZER"
where À is taken to be sufficiently large such that À = 1)( If we interpret À as the Lagrangian multiplier, then Problem (3) can be viewed as the Lagrangian function of the constrained opti- mization gen ye f(z, s),..t. g(x, y) S g(z,, v2 (z)). Thanks to Danskin's theorem, the gradient of the Problem (3) only involves gradient information. Therefore, F²SA does not require the stochastic Hessian assumptions (5). More importantly, by directly leveraging gradient oracles instead of more expensive HVP oracles, the F²SA is more efficient in practice (Shen et al. 2025a Xiao & Chen 2025 Jiang et al., 2025) and is also the only method that can be scaled to 32B sized
large language model (LLM) training (Pan et al. 2024).
Kwon al. (2023) proved that the F2SA method finds e-stationary point of (a())) with Ö(c-3) first-order oracle calls in the deterministic case and Ō(c=7) stochastic first-order oracle (SFO) calls in the stochastic case. Recently, [Chen et al.] (2025b) showed the two-time-scal stepsize strategy improves the upper complexity bound of F*SA method to Ō(c-2) in the deterministic case, which optimal up to logarithmic factors. However, the direct extension of their method in the stochastic case leads to the O(c-9) SFO complexity (Chen et al. 2025b Kwon al., 2024a) which still has a significant gap between the ((c-*) lower bound for SGD (Arjevani et al. 2023 It remains open whether optimal rates for stochastic bilevel problems can be achieved for fully first-order methods.
In this work, we revisit FSSA and interpret it as using forward difference to approximate the hyper- gradient. Our novel interpretation turn leads to straightforward algorithm extensions for the FASA method. Observing that the forward difference used by FSSA only has first-order error guarantee, a natural idea to improve the error guarantee is to use higher-order finite difference methods. For instance, we know that the central difference has an improved second-order error guarantee. Based on this fact, we can derive the F2SA-2 method that solves the following symmetric penalty problem:
min zg yER" f(x. v) + )g(xx, y) min (-f(x, 2) + Ag(xx
ZER
Compared with Eq. (3). this new penalty problem perturbs the lower-level variables y and z in the opposite direction to better cancel out the approximation errors Problem ( The connec- tion between bilevel optimization and finite difference approximations was recently established by Chayti & Jaggi (2024) in the context of meta-leaming, but their results were limited to symmetric approximations. We extend their findings beyond meta-learning and to general finite difference ap- proximations, addressing their conjecture about broader applicability. F2SA-2 provably improves the SFO complexity of FASA from Ō(c-0) to Ō(c-5) for second-order smooth problems. Moreover, our idea is generalizable for any pth-order smooth problems. It is known in numerical analysis there exists the pth-order central difference that uses p points to con- struct an estimator to the derivative of a unitary function with pth-order error guarantee, as recalled in Lemma 3.1, Motivated by this fact, we propose the FSSA-p algorithm and show that it enjoys the improved O 4-2/p) SFO complexity, as formally stated in Theorem 3.1
this work, we show that
Moreover, as formally stated in Theorem 4.1 it is easy to extend the Q(e-4) lower bound for SGD (Arjevani et al. 2023) from single-level optimization to bilevel optimization using a fully separable construction that automatically satisfies all our additional smoothness conditions in Definition 2.2 shows that F2SA-p is optimal up to logarithmic factors when p = M(oo //lo log e (see
Remark 3.4p. We summarize our main results in Table I and discuss open problems in the following.
2

---

Published as a conference paper at ICLR 2026
Open problems. Our upper bounds improve known results for high-order smooth problems, but our result still has a gap between the l2(e lower bound for = Ollo e log log c -1) Recently, Kwon al. (024a) obtained some preliminary results towards closing this gap for p 1, where they showed an (2(c lower bound holds under a more adversarial oracle. But it is still open whether their lower bounds can be extended to standard stochastic oracles. Another open problem is the tightness of the condition number dependency, for which the current upper and lower bounds have gap of (xxx)) demonstrated in Table [ Two recent concurrent works (Ji 2025 Chen & Zhang 2025) proposed tighter lower bounds for = 1: Ji (2025) showed a lower bound of ( )) and Chen & Zhang (2025) showed that of {2 ). However, it is still open to fully
close the gap in condition number dependency for both p 1 and p 2 2.
Method
Smoothness
Reference
Complexity
FSA
1st-order
Kwon
2023
O(poly(n))
FSA
1st-order
Kwon al.
2024a
((poly(k)e 6) Ö(st22
FSSA
1st-order
Chen et al.
2025b)
PSSA-p
Ist-order
Theorem 3.1
Ö(pn (4-2/1))))
+ pth-order in y
Lower Bound
Arjevani et al.
2023
((-
Table I: The SFO complexity of different methods to find e-stationary point for pth-order smooth first-order bilevel problems with condition number under standard SGD assumptions.
Notations. We use ||- || denote the Euclidean norm for vectors and the spectral norm for matrices and tensors. For any set S and functions g. h S + [0, 2o) we write = O(h) or h (gg) equivalently if there exists >> 0 such that y(s) S ch(s) for every E S. We use O and Ω -
hide logarithmic factors (.) and 22 ). We alsouse hi S hi to mean h (h2)), (.12)..1 hn Z h2 to mean h1 12hh2). and hı f2 to mean that both h1 $ h2 and h1 2 h2 hold. Additional notations
for tensors are introduced Appendix
2 PRELIMINARIES
The goal of bilevel optimization is minimize the hyper-objective p(x), which is in general non- convex. Since finding a global minimizer of general nonconvex function requires exponential complexity in the worst case (Nemirovskij Yudin 1983 1.6), we follow the literature (Carmon
et al. 2020 2021) to consider the task of finding an approximate stationary point.
Definition 2.1. Let i Rd, IR be the hyper-objective defined in Eq. variable & E de is an e-hyper-stationary point if V S
1
We say random
Next, we introduce the assumptions used in this paper, which ensure the tractability the above hyper-stationarity Compared to (Kwon et al. 2023 Chen al. 2025b), additionally assume the high-order smoothness in lower-level variable y to achieve further acceleration.
2.1 PROBLEM SETUP
First of all, we follow the standard assumptions on SGD (Arjevani et al. 2023) to assume that the stochastic gradient estimators satisfy the following assumption.
Assumption 2.1. There exists stochastic gradient estimators F(x. y) and G(æ, y) such that
EF(x, E) f y). E| y) f( r)ll < <²² EG(x, (. y: = Vg(n,y) El y) Vg(æ,, y) 2 < 0²
where o >> is the variance of the stochastic gradient estimators. We also partition F (Fa, and G (Ga, Gy) such that Fa, Fy, Gz, Gy are estimators VzJ, vyf, Vag, V Vyg. respectively.
3

---

Published as a conference paper at ICLR 2026
Second, we assume that the hyper-objective p(x) = f(z. y (z)) is lower bounded. Otherwise, any algorithm requires infinite time to find a stationary point. Note that the implicit condition infzeRd: ((x) > -oo below can also be easily implied by more explicit condition on the lower boundedness of upper-level function inf weR'", (ax, y) > -00.
Assumption 2.2. The hyper-objective defined in Eq. 7 is lower bounded, and we have
<0(x)) - inf WER
(xx)) S
where A 0 is the initial suboptimality gap and we assume To 0 without loss of generality.
Third, we assume the lower-level function g( y) is strongly convex in y. It guarantees the unique- ness of y (x: and the tractability of the bilevel problem. Although not all the problems in applica- tions satisfy the lower-level strong convexity assumption, it is impossible to derive dimension freeee upper bounds when the lower-level problem is only convex (Chen et al. 2024 Theorem 3.2). Hence, we follow most existing works to consider strongly convex lower-level problems.
Assumption 2.3. (z, y) is p-strongly convex in y. i.e., for any yl: y/2 E R we have
g a. ya) 2 g(x, y1) + (Vyg(z, yn). U2 y1) + 1 ll y2ll?
2
where μ > 0 is the strongly convex parameter.
Fourth, require the following smoothness assumptions following (Ghadimi & Wang 2018). Ac- cording to Eq. (2). these conditions are necessary and sufficient to guarantee the Lipschitz continuity of Vy(zz) which further ensure the tractability of an approximate stationary point of the nonconvex hyper-objective sp(x)) (Zhang et al., 2020 Kornowski & Shamir 2022).
Assumption 2.4. For the funer-looersunntionn I and lower-level funcriieveliunntioon g, we assume that
[(x, y) is Lo Lipichicii in y-
2. fx, 3/) and Vg(æ, y) are L ---iipschttzjoinnyy in (a, y).
3. Vayg((: y) and (2, y) are Lzz-Lipschitz in (æ: y)-
We refer to the problem class that jointly satisfies all the above Assumption 2.1 2.2 2.3 and 2.4 as first-order smooth bilevel problems, for which (K et al.) [202aa] Chen et al. 2025b) showed the F'SA method achieves the Ō(c-²) upper complexity bound. In this work, we show an improved bound under the following additional higher-order smoothness assumption on lower-level variable y.
Assumption 2.5 (High order smoothness y). Given p N4, we assume that
1.
Sy V (a( is -ips--it for all q =. .P-
g-+1 2. T-.-T Vg(æ, y) is Lips-hii in y for all = 1,.., ..- 1
We refer to problems jointly satisfying all the above assumptions as pth-order smooth bilevel prob- lems, and also formally define their condition numbers as follows.
Definition 2.2 (pth-order smooth bilevel problems). Given E N A 0, Lo. L1. Lp.1 >> O, and p S L1. we use Fre-sc (Lo. Lpt1, μ, A) denote the set of all bilevel instances satisfying Assumption 2.2 2.3 2.4 and 2.5 For this problem class, we define the largest smoothness constant
max L. and condition number K L/p.
All our above assumptions align with (Chen et al. 2025b except for the additional Assumption 2.5 A classic example a highly smooth function is the softmax function (Garg et al.] [021 Lemma 2(3)). Therefore, many hyper-parameter tuning problems for logistic regression provably highly smooth, such that our theory can be applied. We give two examples from (Pedregosa 2016): the first one aims learn the optimal weights of each sample in a corrupted training set, and the second one aims to learn the optimal regularizer of each parameter.
Example 2.1 (Data syper-cleaning Let x E R" parameterize the per--errappppe per-sample weight of a training set with samples via (3r.) exp(x:)) E exp(x;) and y E Rdd be the parameters of a linear
4

---

Published as a conference paper at ICLR 2026
model. Let Eval be the logistic loss of the linear model on the validation set and Eir be the logistic loss on the training sample i. The problem aims to solve
min Eaal øER
s.t. y E arg min
Σ (a)
yER
Example 2.2 -eeaa-to-regullriie) Let x E R.d parameterize the regularization matrix via W.. diag(exp(z)). and E Rd be the parameters of linear model. Let leaa and Eor be the logistic loss of the linear model on the validation set and training set, respectively. The problem aims solve
min (val WERd
s.t.
arg min lor + +y WWy. yERd
2.2 COMPARISON TO PREVIOUS WORKS
Before we show our improved upper bound, we first give a detailed discussion on other assumptions made in related works.
Stochastic Hessian assumption. Ghadimi & Wang 2018): Ji et al. 2021 assumes the access to a stochastic Hessian estimator (a. y) such that
(H y) V2g(a, y), EE E| y) V y)l || S 0
2020 Observation and 2), such an assumption is stronger than standard SGD assumptions and equivalent to the meea-ssuuared-ssoootnnss assumption (6) on the lower-level gradient estimator G under the mild condition of VG(z, y) (xx, y). Under this
According (Arjevani
assumption, in conjunction with Assumption 2.2 2.3 and 2.4 Ghadimi & Wang (2018) proposed
the BSA method that can find an E stationary point of (a()) with O SFO complexity and
Ö(c-4) stochastic HVP complexity. Later, Ji et al. (021) ) further improved the SFO complexity
term Ö(E Compared to them, we consider the setting where the algorithms only have access to stochastic gradient estimators, and make no assumptions on the stochastic Hessians.
Mean-squared smoothness assumption. Besides Assumption 2. 2.2 2.3 2.4 and the stochastic
Hessian assumption (5) Khanduri et al. (2021); Yang et al. (2021; ; 20236 further assumes that the
stochastic estimators gradients and Hessians are mean -quurre smooth:
F(x. y) - F(x', y'')l < [, y) - (æ', 3/ 3/)11
E G(x,, y) G(m', v')11 < (1) y) - (m', v/)
(6)
E H (z: y) HH', (('','')\l" < Lal y) - -x'''')11
Under this additional assumption, they proposed faster stochastic methods with upper complexity bound of Ō(E-3) via variance reduction (Fang et al. 2018:; Cutkosky Orabona 2019). In this paper. we only consider the setting without mean-squared smoothness assumptions and study different acceleration mechanism from variance reduction.
Jointly high-order smoothness assumption. Huang al. (2025) introduced second-order
smoothness assumption similar to but stronger than Assumption 2.5 when = 2. Specifically, they assumed the second-order smoothness jointly in (xx, y) instead of only:
2 (xx y) is L2-Lipschitz jointly in (x,
ga y) is Ly-Lipschitz jointly in y).
The jointly second-order smoothness 7 ensures that the hyper-objective (p(x) has Lipschitz con- tinuous Hessians, which further allows the application of known techniques in minimizing second- order smooth objectives. [Huang al. (2025) applied the technique from (Jin et al. 2017 202T Xu et al. 2018 Allen-Zhu & Li 2018) to show that an HVP-based method can find a second-order stationary point Ō(e-2) complexity under the deterministic setting, and in ( under
stochastic Hessian assumption ( Yang al. (2033a) applied the technique from (Li & Lin 2023)
to accelerate the complexity HVP-based method to O(751.75) in the deterministic setting. Chen
et al. 2025b) also proposed fully first-order method achieve the same Ō complexity.
Compared to these works, our work demonstrates a unique acceleration mechanism in stochastic bilevel optimization that only comes from the high-order smoothness y.
5

---

Published as a conference paper at ICLR 2026
3 THE FASA-p METHOD
To introduce our method, we first recall the prior F²SA method (Kwon et al. 2023 and establish their relationship between finite difference schemes, which further motivates us to design better algorithms by using better finite difference formulas.
3.1 HYPER-GRADIES GRADIENT APPROXIMATION VIA FINITE DIFFERENCE
The core idea of F2SA (won et al. 2023) the gradient of the penalty function to approximate the true hyper-gradient. To make connections of SASA to the finite difference method, let us introduce the extra notation Iv as the perturbed lower- lever problem with y. (x: and (xx)) being its optimal solution and optimal value, respectively:
solve the reformulated penalty problem e and use
(x, y) == fa y) + g( y).
V (z) = arg min 9u (x, y) VER
Eo (z)) == min go (z. y)
yER
Then we have de (z))) = lim, (v (o))) = lim,- fax yd (x))) [ ((m) g(m.ay (a))
In our notation, we can rewrite (Chen et al. we can also rewrite (Kwon et al.,
2025b Lemma B.3) as By Lu (x(x)) = (xx)) Similarly,
8
2023 Lemma 3.1) (l)))
8²
ly (=))===
dzdv (=))== Vvp(x)
(8)
Ōvdx
Let 1/A Eq. 3
Then the fully first-order aeper-gradient estimator (won et al. 2023 Chen et al. 2025b) is exactly using forward difference to approximate Vip(z), that is,
e. x
(x De
0²
ly (x))== = Vip(x))
(9)
Ovãx
However, the forward difference not the only way to approximate derivative. Essentially, falls into general class of pth-order finite difference (Atkinson & Han 2005) that can guarantee an (uv)) approximation error. We restate this known result (with generalization to vector-valued func- tions) in the following lemma and provide self-contained proof in Appendix for completeness.
Lemma 3.1. Assume the function sb : R + Rd has C-Lipschitz continuous pth-order derivative. There exist coefficients {ai} such that
V Σ
agy(iu) - (0() O(CUP).
If p even, the indices run - ... p/2. If odd, they run - 1)/2.
(p-11//2. ....
1)+1)/2.
Furthermore, all the coefficients satisfy رjayl $ for all j + and |ool S T/p is odd).
The explicit formulas for aj can be found in Appendix B
When p= we have 0g = 11 1,,00 01 = 1, and obtain the forward difference estimator (D)) (0)/v; When p
have a
-1 -1/2, a1 = 1/2 and we obtain the central difference estimator (w(v) + (-v))/ Lemma 3.1
tells us that in general we can always construct a finite difference estimator ( error with p points
for even or + 1 points for odd p under the given smoothness conditions. Inspired by Lemma 3.1 and Eq. combination of Pro (x) to achieve O(z²) approximation error to V((()) given that Bur0g (x) Lipschitz continuous v. further leads to Algorithm T that will be formally introduced the next subsection.
8 that 8 (x)
Vip(x),
propose fully first-order estimator via linear ar+1
3.2 THE PROPOSED ALGORITHM
Due to space limitations, we only present Algorithm [ designed for even p in the main text. The algorithm for odd p be designed similarly, and we defer the concrete algorithm to Appendix D

---

Published as a conference paper at ICLR 2026
Algorithm F2S (Zo: Vo), even p
1: y0 = -oo Vj E
2: for 0, 1.... .T. 1
3: parallel for j = -D/2 1, ,2/2
4:
Vt 3.0 y
5:
for k = 0, 1, ..... --1
6:
Sample random i..d indexes {(E). 23)3.
1. y By jiuP "/; Gy vb vi (c)))
7:
8:
end for
7/+1 = -
9:
10: end parallel for
Sample random iid indexes {(E F-
12/2 Let (a,))=-p/2 be the pth-order finite difference coefficients defined in Lemma 3.1
12:
4 1/2 Σ Σ31² 2z=-p/2 aj jFalar, v(+ + Gz (Ze: v(+
13:
14:
Ec++ The 1111111111
15:
end for
Algorithm follows the double-loop structure of FASA (Chen et al. 2025b Kwon et al. changes the hyper-gradient estimator to one introduced the previous section. Now, we give more detailed introduction to the procedures of the two loops of F-SA-p
2024a and
1. In the outer loop, the algorithm first samples a mini-batch with size S and uses Lemmal 3.1 construct t via the linear combination of Za (ax.) for j -p/2. ... p/ every iteration. After obtaining 4 as an approximation to V Vy((;;,, the algorithm then performs normalized gradient descent step # 2, 1h. / / dell with total T iterations.
2. In the inner loop, the algorithm returns an approximation to & the (xxe) for all - .2/. Note that Danskin's theorem indicates (xx) du Iu Y (x:)). It suffices to approximate (xx) to sufficient accuracy, which is achieved by taking K K-step single-batch SGD subroutine with stepsize ny on each function Ijv (x.
Remark 3.1 (Effect of normalized gradient step). Compared to (Chen et al. 2025b Kwon et al. 2023), the only modification make to the outer loop change the gradient step normalized gradient step. The normalization can control the change of (xx)) and make the analysis of inner
loops easier. We believe that all theoretical guarantees also hold for the standard gradient step via a more involved analysis.
3.3
COMPLEXITY ANALYSIS
This section contains the complexity analysis of Algorithm 1 We first derive the following lemma from the digh-dimensional Faà di Bruno formula (Licht) 2024).
Lemma 3.2. Let E (0, 12x)) For any instance in the pth-order smooth bilevel problem class Fre (Lo Lp+1, ,4, A) as Definition 2.2
Buroz (x is O(k2p+1 L)-Lisschitzc continuous in v.
gell l
Our result generalizes the prior result for the prior bounds for p = (Chen et al. 2025b) as we remark in the following.
2023 to any p N4 and also tightens (Kwon al.
Remark 3.2 (Tighter bounds for p 2). Note that the variables x and play equal roles in our analysis. Therefore, our result in = essentially implies that lv))) is 15 L>-Liscchitz continuous in D around zero, which tightens the ( L) bound of Hessian convergence in Cheen al. 20255b Lemma 5.1a) and is of independent interest. The main insight is to avoid the direct cal- culation of (2)((x) -o (x())))))) which involves third-order derivatives and makes the anal- ysis complex, but instead always analyze it through the limiting point lim Boogz (x)
(x))
83

---

Published as a conference paper at ICLR 2026
Recall Eq. ( that 8² (x(x))
âvōæ audzfu(z))
Vp(z)) Then Lemma 3.2 in conjunction with Lemma 3.1
indicates that the pth-order finite difference used in FASA-p guarantees Happroxmimiioitoii
error to (x), which always improves the (vr)-rrror guarantee of F2SA (Kwon et al. 2023 Chen
et al. 2025b) for any p 2 2. This improved error guarantee means that we can set v = ((et/2)
obtain an (ac)aucee hyper-gradient estimator Vip(ar), which further leads to the following improved complexity of our algorithm.
Theorem 3.1 (Main theorem). For any instance in the pth-order smooth bilevel problem class set the haperrpprameters.t
Fuess (Lo, Lpt11 n, A) per Definition 2.2
R
1/D
v²c²
min
Tx22+1
The
L
Fy LIRO²
(10)
x²
R²O²
RRL1
A
S x
K
log
T
V2e2'
v2e2
ve
ThE
where R = 11y0 y" (zo))l. Run Algorithm iif p is even Algorithm (in Appendix D if is odd. Then we can provably find an e-stationary point of ((x) with the total SFO calls upper bounded by
( 4+2/p
p
log (RILLK) (Hl)
=)= =
The above theorem shows that the F²SA-p method can achieve the 0pp log log((/))) SFO complexity for pth-order smooth bilevel problems. In the following, we give several remarks on the complexity in different regions of p.
Remark 3.3 (First-order smooth region). For p 1, our upper bound becomes ( which improves the ( bound in (Chen et al. 2025b) by a factor of K. The improvement comes from a tighter analysis in the lower-level SGD update and careful parameter setting.
Remark 3.4 (Highly-smooth region). For p ((log(+/e)/ log log(*/+) in Definition 2.2 we can run F SA-q with q x log(*/e)/ log (*))) and the O (x//2)22/ (og(4/))))g(((())) complexity in Theo- rem 3. simplifies to O "oog" (*/e)/ log log((/+)) = Ó(ne which matches the best-known complexity for HVP-based methods (Ji et al. [222] under stochastic Hessian assumption (5)
In the upcoming section, will derive an 2(c=*) lower bound prove that the FSSA-p is near- optimal in the above highly-smooth region the condition number λ constant. We leave the study of optimal complexity for non-constant future work.
Comparison of results for odd p and even Note that by Lemma 3.1 when is odd, we need to use p+1 points construct the estimator, which means the algorithm needs to solve p+ lower-level problems in each iteration to achieve an ((PP) еттог guarantee. In contrast, when is even, p points are enough since the pth-order central difference estimator satisfies that ag
It suggests that even when p is odd, the algorithm designed for odd may still be better. For instance, the F2SA-2 may always be a better choice than F2SA since its beneñit almost come for free: (1) still only needs to solve 2 lower-level problems as the FASA method, which means the per-iteration complexity remains the same. (2) Although the improved complexity of F"SA-2 relies on the second-order smooth condition, without such a condition, its error guarantee in hyper-gradien, estimation will only degenerate to a first-order one, which means it at least as good as F²SA.
4 AN AN 4) LOWER BOUND
In this section, we prove an 2(E-4) lower bound for stochastic bilevel optimization via reduction to single-level optimization. Our lower bound holds for any randomized algorithms A, which can be defined as a sequence of measurable mappings {Ar)T that is defined recursively by
(w 1(+1) = Av (r, F(2o, yo), G(oo. yo)), ... F( y:), G( y:))), t N4.
(1I)
where r is random seed drawn the beginning produce the queries, and F, G are the stochas- tie gradient estimators that satisfy Assumption 2.1 Without loss of generality, we assume that (æp, yo) (0, 0). Otherwise, can prove the same lower bound by shifting the functions.
8

---

Published as a conference paper at ICLR 2026
The construction. We construct separable bilevel instance such that the upper-level function f(x y) = (z) and its stochastic gradient align with the hard instance in (Arjevani et al. 2023), while the lower-level function is the simple quadratic ( y) g(y) 2yy with deterministic gradients. We defer the concrete construction to Appendix F For this separable bilevel instance,
show that for any randomized algorithm defined in Eq. (11) that uses oracles (Fu, G), the
progress in æ can be simulated by another randomized algorithm that only uses FU, meaning that the single-level lower bound (Arjevani et al. 2023) also holds.
Theorem 4.1 (Lower bound). There exist numerical constants > O such that for all A O, Li,.... Lp+1 > O and é < LvL, there exists a distribution over the function class Foo (Lo (Lo,.... 1,, A) and the stochastic gradient estimators satisfying Assumption 2.7 such that any randomized algorithm A defined as Eq. 477 can not find an e-stationary point of (p(x)= f(x.y' (z)) in less than ((AL SFO calls.
Below, we give detailed discussion on the constructions in related works.
Comparison to other bilevel lower bounds. Dagréou et al. (2024) proved lower bounds for finite- sum bilevel optimization via a similar reduction to single-level optimization. However, the direct extension of their construction in the fully stochastic setting gives f(x y) fo (3) and g(x, y) = (x - ") where the high-order derivatives of f(x y) not (1)-Lipsshitz in y and thus violates our assumptions. Kwon et al. (2024a) also proved an S2(c-1) lower bound for stochastic bilevel optimization. However, their construction f(x,,) v) = and g(æ, y) = (fr((r) - y)² violate the first-order smoothness of g(zz, y) in when y is far way from fo(x). In this work, we use fully separable construction to avoid all the aforementioned issues in other works.
5 EXPERIMENTS
In this section, we conduct numerical experiments to verify our theory. Following (Grazzi et al.) 2020 Ji et al. 202I). we consider the "learn-to- rergulaarre problem of logistic regression (Exam- ple 2.2) on the "20 Newsgroup" dataset, which provably satisfies the highly smooth assumption of any order. The dataset contains 18,000 samples, each sample consists of a feature vector in dimen- sional 130, 107 vector and label that takes value in {1. 20). We compare our proposed method F SA-p with E [2, 3. 5. 8, 10] with both the previous best fully first-order method F2SA (Kwon al. 2023 Chen et al. 2025b) and other Hessian- ceesiinnveetto-preee cector-product-based methods stOcBiO jji et al. 1021). MRBO and VRBO (Yang et al. 202T). We also include a baseline "w/o Reg" that means the training result of SGD without tuning any regularization. For all the algorithms, we search the other hyperparameters (including na Nly 2) in logarithmic scale with base 10. We run the al- gorithms with K = 10 iterations the inner loop, and T = 1000 iterations in the outer loop, and report the test loss/accuracy v.s. the number of outer-loop iterations / in Figure [ To demonstrate the potential of our methods on nonsmooth nonconvex problems, we also provide additional experi- ments 5-layer multilayer perceptron (MLP) network with ReLU activation in Appendix F The codes to reproduce our experiments are available online
6 CONCLUSIONS AND FUTURE WORKS
This paper proposes a class of fully first-order method FPSA-p that achieves the Oppe-4-2/p) SFO complexity for pth-order smooth bilevel problems. Our result generalized the best-known Ö(c-6) result (Kwon et al., 2024a] Chen et al. 20256) from p
to any p E N+. We also com- plement our result with $2(e lower bound to show that our method is near-optimal when p = M(ogg E 1 log log e ))). Nevertheless, a gap still exists when p is small, and how to fill it even is open problem. Another open problem is whether our theory can be extended our theory to structured nonconvex nonconvex-nonccover bilevel problems studied by many recent works (Kwon et al., 2024b; Chen et al.) 2024| 2025a; Diang et al., 2025 Xiao et al. 2023 Xiao & Chen 2025). In addition, it will also be interesting to further improve the convergence rate of methods by combining them with variance-reduction (Fang et al. 2018 Cutkosky & Orabona
for the basic setting
2019 or momentum techniques (Fang et al. 2019 Cutkosky & Mehta
2020).
httg:://gthhh.ccmm//aaaalllllttn

---

Published as a conference paper at ICLR 2026
0.90
0.82
0.85
Test Test Accuracy
0.80
Test Loss
0.80
0.78
0.75
0.76
0.70
0
250
500
750
1000
250
500
750
1000
#Iterations
#lterations
5A
FSA-3
FSSA-8
VRBO MRBO
stocBiO
F2SA-2
F2SA-5
F²SA-10
w/o Reg
Figure 1: Performances of different algorithms Example 2.2
ACKNOWLEDGMENTS
Lesi Chen thanks Jeongyeol Kwon for helpful discussions and Sanyou Mei for pointing out related references. Jingzhao Zhang supported by Shanghai Qi Zhi Institute Innovation Program, Tsinghua Dushi Funds, and Xiongan AI Institute.
REFERENCES
Zeyuan Allen-Zhu and Yuanzhi Li. Neon2: Finding local minima via first-order oracles. In NeurIPS, 2018.
Yossi Arjevani, Yair Carmon, John C. Duchi, Dylan J. Foster, Ayush Sekhari, and Karthik Sridha- ran. Second-order information in non-convex stochastic optimization: Power and limitations. In COLT, 2020.
Yossi Arjevani, Yair Carmon, John C. Duchi, Dylan Foster, Nathan Srebro, and Blake Woodworth. Lower bounds for non-convex stochastic optimization. Mathematical Programming, 199(1):165- 214, 2023.
Kendall Atkinson and Weimin Han. Finite difference method. Theoretical Numerical Analysis: A Functional Analysis Framework, pp. 249-271, 2005.
Fan Bao, Guoqiang Wu, Chongxuan Li. Jun Zhu, and Bo Zhang. Stability and generalization of bilevel programming in hyperparameter optimization. In NeurlPS, 2021.
Yair Carmon, John C. Duchi, Oliver Hinder, and Aaron Sidford. Lower bounds for finding stationary points i. Mathematical Programming, 184(1):71-120, 2020.
Yair Carmon, John C. Duchi, Oliver Hinder, and Aaron Sidford. Lower bounds for finding stationary points ii: first-order methods. Mathematical Programming, 185(1):315-355 2021.
El Mahdi Chayti and Martin Jaggi. A new first-order meta-learning algorithm with convergence guarantees. arXiv preprint 200.:22409.33682. 2024.
He Chen, Jiajin Li, and Anthony Man-cho So. Set smoothness unlocks clarke hyper-stationarity in bilevel optimization. NeurIPS, 2025a.
Lesi Chen and Jingzhao Zhang. On the condition number dependency in bilevel optimization. arXiv preprint vrXiv:2511.2233/1 2025.
Lesi Chen, Jing Xu, and Jingzhao Zhang. On finding small Hardness results and improved analysis. In COLT, 2024.
bilevel optimization:
10

---

Published as a conference paper at ICLR 2026
Lesi Chen, Yaohua Ma, and Jingzhao Zhang. Near-optimal nonconvex-sttrnnglyyoomee bilevel op- timization with fully first-order oracles. JMLR, 2025b.
Ashok Cutkosky and Harsh Mehta. Momentum improves normalized SGD. In ICML, 2020.
Ashok Cutkosky and Francesco Orabona. Momentum-based variance reduction in non-convex SGD. In NeurlPS, 2019.
Mathieu Dagréou, Thomas Moreau, Samuel Vaiter, and Pierre Ablin. A lower bound and a near- optimal algorithm for bilevel empirical risk minimization. AISTATS, 2024.
Asen L. Dontchev and R. Tyrrell Rockafellar. Implicit functions and solution mappings, volume 543. Springer, 2009.
Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang. Spider: Near-optimal non-convex optimization via stochastic auth-integrated differential estimator. In NeurlPS, 2018.
Cong Fang, Zhouchen Lin, and Tong Zhang. Sharp analysis for nonconvex SGD escaping from saddle points. In COLT, 2019.
Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In ICML, 2018.
Ankit Garg, Robin Kothari, Praneeth Netrapalli, and Suhail Sherif. Near-optimal lower bounds for convex optimization for all orders of smoothness. In NeurlPS, 2021.
Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. arXiv preprint arXiv: /802.0224, 2018.
Ian Goodfellow, Jean Pouget-Abadie, Mchdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial networks. Communications of the ACM, 63(11):139914) 2020.
Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration com- plexity of hypergradient computation. In ICML, 2020.
Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A two-timescale stochastic algorithm framework for bilevel optimization: Complexity analysis and application to actor-critic. SIAM Journal on Optimization, 33(1):147-180, 2023.
Minhui Huang, Xuxing Chen, Kaiyi Ji, Shiqian Ma, and Lifeng Lai. Efficiently escaping saddle points bilevel optimization. JMLR, 26(1):1-61, 2025.
Kaiyi Ji. Lower complexity bounds for nonconnee--ttronllyoooon bilevel optimization with first- order oracles. arXiv preprint arXiv:25/1.19656. 2025.
Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In ICML, 2021.
Liuyuan Jiang, Quan Xiao, Lisha Chen, and Tianyi Chen. Beyond value functions: Single-loop bilevel optimization under flatness conditions. In NeurlPS, 2025.
Chi Jin, Rong Ge, Praneeth Netrapalli, Sham M. Kakade, and Michael Jordan. How escape saddle points efficiently. In /CML, 2017.
Chi Jin, Praneeth Netrapalli, Rong Ge, Sham M. Kakade, and Michael I. Jordan. nonconvex optimization for machine learning: Gradients, stochasticity, and saddle points. Journal of the ACM, 68(2):1-29, 2021.
Ishtiaq Rasool Khan, Ryoji Ohba, and Noriyuki Hozumi. Mathematical proof of closed form pressions for finite difference approximations based on taylor series. Journal of Computational and Applied Mathematics, 150(2)3003-309 2003.

---

Published as a conference paper at ICLR 2026
Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A near-optimal algorithm for stochastic bilevel optimization via double-momentum. In NeurlPS, 2021.
Tamara G. Kolda and Brett Bader. Tensor decompositions and applications. SIAM review, 51(3): 455-500, 2009.
Guy Kornowski and Ohad Shamir. Oracle complexity nonsmooth nonconvex optimization. JMLR, pp. 1-44, 2022.
Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert D. Nowak. fully first-order method for stochastic bilevel optimization. In ICML, 2023.
Jeongyeol Kwon, Dohyun Kwon, and Hanbaek Lyu. On the complexity of first-order methods in stochastic bilevel optimization. In ICML, 2024a.
Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert D. Nowak. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. In ICLR, 2024b.
Huan Li and Zhouchen Lin. Restarted nonconvex accelerated gradient descent: No more polyloga- rithmic factor in the in the O(c-7/4 complexity. JMLR, 2023.
Martin Licht. Higher-order chain rules for tensor fields, generalized bell polynomials, and esti- mates in orticz-sooboeer-slooddkkk and total variation spaces. Journal of Mathematical Analysis and Applications, 534(1):128005, 2024.
Bo Liu, Mao Ye. Stephen Wright, Peter Stone, and Qiang Liu. Bome! bilevel optimization made easy: simple first-order approach. NeurIPS, 2022.
Risheng Liu, Xuan Liu, Shangzhi Zeng, Jin Zhang, and Yixuan Zhang. fulteefunctioo-baseed se- quential minimization for bi-level optimization. TPAMI, 15112))15530 2023.
Zhaosong Lu and Sanyou Mei. First-order penalty methods for bilevel optimization. SIAM Journal on Optimization, 34(2):1937-1969 2024.
Zhaosong Lu and Sanyou Mei. Solving bilevel optimization via sequential minimax optimization. Mathematics of Operations Research, 2026.
Luo Luo, Yujun Li, and Cheng Chen. Finding second-order stationary points in nonconvex-strongly. concave minimax optimization. In NeurIPS, 2022.
Matthew Mackay, Paul Vicol, Jonathan Lorraine, David Duvenaud, and Roger Grosse, Self-tuning networks: Bilevel optimization of hyperparameters using structured best-response functions. In ICLR, 2019.
Arkadij Semenovič Nemirovskij and David Borisovich Yudin. Problem complexity and method efficiency in optimization. 1983.
Rui Pan, Jipeng Zhang, Xingyuan Pan, Renjie Pi, Xiaoyu Wang, and Tong Zhang. ScaleBiO: Scal- able bilevel optimization for LLM data reweighting. arXiv preprint arXh::2406.19976 2024.
Fabian Pedregosa. Hyperparameter optimization with approximate gradient. In 2016.
Aravind Rajeswaran, Chelsea Finn, Sham M. Kakade, and Sergey Levine. Meta-learning with im- plicit gradients. In NeurlPS, volume 32, 2019.
Han Shen and Tianyi Chen. On penalty- based bilevel gradient descent method. In ICML, 2023.
Han Shen, Pin-Yo Yu Chen, Payel Das, and Tianyi Chen. Seal: Safety-enhanced aligned Ilm fine-tuning via bilevel data selection. In ICIR, 2025a.
Han Shen, Quan Xiao, and Tianyi Chen. On penalty-based bilevel gradient descent method. Math- ematical Programming, pp. 1-51, 2025b.
12

---

Published as a conference paper at ICLR 2026
Quan Xiao and Tianyi Chen. Unlocking global optimality in bilevel optimization: A pilot study. In ICLR, 2025.
Quan Xiao, Songtao Lu, and Tianyi Chen. An generalized alternating optimization method for bilevel problems under the polyak-tojasewwczz condition. In NeurIPS, 2023.
Yi Xu, Rong Jin, and Tianbao Yang- First-order stochastic algorithms for escaping from saddle points in almost linear time. In NeurlPS, 2018.
Haikuo Yang. Luo Luo, Chris Junchi Li, and Michael I Jordan. Accelerating inexact hypergradient descent for bilevel optimization. arXiv preprint arXiv:330.00112, 2023a.
Junjie Yang, Kaiyi Ji, and Yingbin Liang. Provably faster algorithms for bilevel optimization. In NeurlPS, 2021.
Yifan Yang, Peiyao Xiao, and Kaiyi Ji. Achieving ( complexity in hessian/jacobian-Fie free stochastic bilevel optimization. In NeurlPS, 2023b.
Zhuoran Yang, Yongxin Chen, Mingyi Hong, and Zhaoran Wang. Provably global convergence of actor-critic: A case for linear quadratic regulator with ergodic cost, 2019.
Sihan Zeng and Thinh Doan. Fast two-time-scale stochastic gradient method with applications in reinforcement learning. In COLT, 2024.
Sihan Zeng, Thinh Doan, and Justin Romberg, A two-time-scale stochastic optimization framework with applications in control and reinforcement learning. SIAM Journal on Optimization, 34(1): 946-976, 2024.
Jingzhao Zhang, Hongzhou Lin, Stefanie Jegelka, Suvrit Sra, and Ali Jadbabaie. Complexity of finding stationary points of nonconvex nonsmooth functions. In ICML, 2020.
13


Reviews:
[{"review_id": "BlARFgVWf2", "reviewer": "Reviewer_Rwzc", "summary": "This paper proposes a variant of the first-order method, F2SA, which uses higher-order central/forward finite difference to approximate the hypergradient. Therefore, the proposed method is called ‌F2SA-p, which uses p points to construct an estimator to the derivative of a unitary function (e.g., p = 1 or 2).", "questions": "1. What's the difference between the proposed method with F2SA in (Kwon et al., 2024a) and (Chen et al., 2025b)? \n\n2. The ‌problem setup‌ in Section 2 aligns with prior works like Chen et al. (2025b) and Kwon et al. (2024a), with the addition of ‌Assumption 2.5‌ (high-order smoothness).\n\n3. Fully first-order methods are ‌more scalable‌ than Hessian-vector-product-based approaches, particularly for large-scale problems like LLM training. What's the advantage of the proposed method with p= 2 or 3?\n\n4. The experimental results are not convincing. For instance, the authors should report the performance of the proposed method when p= 3 or 4.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "Strength:\n\n1. The work is ‌well-motivated‌, and provides both upper and lower bound analysis. \n\n2. For ‌pth-order smooth bilevel problems‌, the authors sated that the proposed method achieves an improved lower bound.", "weaknesses": "Weakness:\n\n1. What's the difference between the proposed method with F2SA in (Kwon et al., 2024a) and (Chen et al., 2025b)? \n\n2. The ‌problem setup‌ in Section 2 aligns with prior works like Chen et al. (2025b) and Kwon et al. (2024a), with the addition of ‌Assumption 2.5‌ (high-order smoothness).\n\n3. Fully first-order methods are ‌more scalable‌ than Hessian-vector-product-based approaches, particularly for large-scale problems like LLM training. What's the advantage of the proposed method with p= 2 or 3?\n\n4. The experimental results are not convincing. For instance, the authors should report the performance of the proposed method when p= 3 or 4.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "QQStodqeBl", "reviewer": "Reviewer_MZzi", "summary": "In this paper, the authors study a stochastic bilevel optimization problem where the upper-level objective is non-convex and the lower-level objective is strongly convex. Instead of relying on Hessian-vector products, they focus on first-order methods that use only gradient information. The key insight is that the state-of-the-art first-order method, $\\textsf{F}^2\\textsf{SA}$, can be viewed as approximating the hypergradient via a forward difference of a univariate function. Building on this idea, the authors employ higher-order finite difference schemes and, under proper smoothness assumptions, establish a complexity bound of $\\tilde{O}(p \\kappa^{9+2/p} \\epsilon^{-4-2/p})$. When $p=1$, their result improves the previous analysis by a factor of $\\kappa$, and when $p = \\Omega(\\log \\epsilon^{-1} / \\log\\log \\epsilon^{-1})$, the complexity improves to $\\tilde{O}(\\kappa^9 \\epsilon^{-4})$, matching the optimal rate for single level problems.", "questions": "Please see the \"Weaknesses\" Section above.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The $\\textsf{F}^2\\textsf{SA}$ method was originally proposed as a penalty-based approach. In this paper, the authors offer a novel interpretation by viewing it as a forward difference approximation of the hypergradient. This perspective naturally extends to higher-order finite-difference schemes and provides a clear pathway to improve first-order methods.  \n- The authors also perform a tight analysis to bound the Lipschitz constants of the higher-order derivatives of the perturbed function in terms of those of the original objectives. Under additional smoothness assumptions, they derive a sharper complexity bound than existing first-order methods.", "weaknesses": "- **Dimension dependence in Remark 3.1**. The remark states that Lemma 3.1 can be applied to each dimension and thus the overall error can be controlled in Euclidean norm. However, I think a direct coordinate-wise argument would incur an extra $\\sqrt{d}$ factor in the approximation error, which could be problematic. Could you please clarify whether your analysis avoids this factor, e.g., by analyzing the error in different coordinates jointly?\n- **Lower bound in Section 4**. The construction effectively considers a bilevel problem where the upper-level and lower-level objectives are independent, thus reducing to a single-level non-convex problem. While I am not aware of a prior formal statement in exactly this form, the resulting conclusion feels fairly direct.\n- The presentation can be improved, and several points would benefit from additional detail: \n    - In the context of the penalty-based method, $\\nu$ plays the role of the inverse of the penalty parameter $\\lambda$, so I find it somewhat counterintuitive that the algorithm considers negative $\\nu$ when $p$ is even. This also raises the question whether $\\min_{y} g_{\\nu}(x,y)$ can be solved efficiently. My understanding is that taking $\\nu$ sufficiently small preserves smoothness and strong convexity in $y$, so the subproblem can still be solved in a standard way.\n   - On a related note, the relationship between the symmetric penalty problem in (4), the penalty problem in (3), and the original bilevel objective is not immediately clear to me. A brief discussion clarifying their connections will be useful. \n    - Typos:\n         - Page 5, Line 262. The identity $\\frac{\\partial}{\\partial \\nu} \\ell_{\\nu}(x) \\vert_{\\nu=0} = \\phi(x)$ is cited to Chen et al. (2025b), Lemma B.3. However, that lemma appears to bound the gap between the penalized function and the hyperobjective.  \n         - In Lemma 3.1, the coefficients for even $p$ should satisify $\\alpha_j = -\\alpha_{-j}$. Also, for the $p$th-order forward difference estimator, the indices should run $j=0, \\dots, p$. \n         - In Lemma D.2, the SGD update is written as the normalized SGD by mistake.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "SOKnUIehPH", "reviewer": "Reviewer_pQay", "summary": "This paper studies the stochastic bilevel optimization problem in which the upper-level objective is nonconvex and the lower-level objective is strongly convex. Building on the first-order F2SA algorithm, which achieves $\\tilde{O}(\\epsilon^{-6})$ stochastic oracle complexity, the authors propose a family of methods F2SA-p that leverage p-th order finite-difference hyper-gradient approximation. They prove that under p-th order smoothness in the lower-level variable, the algorithm attains a tighter complexity bound $\\tilde{O}(p \\kappa^{9+2/p} \\epsilon^{-4 - 2/p})$. Moreover, they establish a matching lower bound $\\Omega(\\epsilon^{-4})$ under similar smoothness assumptions, showing near-optimality when $p = \\Omega(\\log(1/\\epsilon)/\\log\\log(1/\\epsilon))$.", "questions": "1. Regarding the convergence guarantee in Theorem 3.1, is the result stated in expectation or with high probability?\n2. I wonder whether the proposed algorithm can be combined with variance-reduction or momentum techniques to further improve the convergence rate.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The finite-difference modification of F2SA is both novel and conceptually elegant, offering a unified view that also tightens and refines previous analyses.\n2. The paper provides an explicit comparison of smoothness assumptions with related works and includes a thorough and contextualized literature review.\n3. A new lower-bound result demonstrates that the upper bound is nearly optimal in the regime $p = \\Omega(\\log(1/\\epsilon)/\\log\\log(1/\\epsilon))$.", "weaknesses": "1. The experiments are limited to a single small-scale dataset, with no evaluation on modern large-scale or deep bilevel tasks (e.g., hyperparameter tuning or meta-learning).\n2. The high-order smoothness assumption on the lower-level variable may be restrictive or difficult to verify in practical ML settings. It would strengthen the paper to discuss or cite works examining when such smoothness assumptions hold in real-world problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "jS5it1n9wI", "reviewer": "Reviewer_mbBx", "summary": "This paper provides a conceptually elegant and mathematically rigorous advancement in stochastic bilevel optimization. By linking finite-difference theory with $p$th-order first-order hyper-gradient estimation, it achieves improved convergence without Hessian information. The analysis is solid and the experiments validate both theory and practicality, though assumptions on high-order smoothness and complexity constants limit immediate applicability to all real-world problems. Nonetheless, it represents a notable theoretical milestone toward closing the gap between first-order bilevel optimization and its single-level counterpart.", "questions": "Could the authors clarify how much Assumption 2.5 (high-order smoothness) contributes to the improvement in the convergence rate? Additionally, could the authors provide justification or empirical evidence regarding the strength or practicality of this assumption?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "This paper provides a novel and elegant finite-difference interpretation of first-order bilevel optimization, enabling the F2SA-p method to achieve near-optimal convergence without any Hessian computation. It successfully bridges theory and practice by offering a fully first-order, scalable, and provably faster framework for highly smooth stochastic bilevel problems.", "weaknesses": "This paper relies on a new assumption of high-order smoothness in the lower-level variable $y$. The improved convergence rate largely stems from this assumption; however, since the method requires a relatively large $p$, the assumption may be quite strong in practice.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)