Achieving $\mathcal{O}(\epsilon^{-1.5})$ Complexity in Hessian/Jacobian-free Stochastic Bilevel Optimization
Authors: Yifan Yang
Peiyao Xiao
Kaiyi Ji
Venue: NeurIPS 2023 poster (2023)
Keywords: Stochastic bilevel optimization
Hessian-free algorithms
near-optimal complexity

Abstract: In this paper, we revisit the bilevel optimization problem, in which the upper-level objective function is generally nonconvex and the lower-level objective function is strongly convex. Although this type of problem has been studied extensively, it still remains an open question how to achieve an $\mathcal{O}(\epsilon^{-1.5})$ sample complexity in Hessian/Jacobian-free stochastic bilevel optimization without any second-order derivative computation. To fill this gap, we propose a novel Hessian/Jacobian-free bilevel optimizer named FdeHBO, which features a simple fully single-loop structure, a projection-aided finite-difference Hessian/Jacobian-vector approximation, and momentum-based updates. Theoretically, we show that FdeHBO requires $\mathcal{O}(\epsilon^{-1.5})$ iterations (each using $\mathcal{O}(1)$ samples and only first-order gradient information) to find an $\epsilon$-accurate stationary point. As far as we know, this is the first Hessian/Jacobian-free method with an $\mathcal{O}(\epsilon^{-1.5})$ sample complexity for nonconvex-strongly-convex stochastic bilevel optimization.

Achieving O Complexity in Hesian/JJcobian-frre Stochastic Bilevel Optimization
Yifan Yang, Peiyao Xiao and Kaiyi Ji
Department of Computer Science and Engineering University at Buffalo
Buffalo, NY 14260
(yyang99, peiyaoxi, (alyiij))buufllll edu
Abstract
In this paper, revisit the bilevel optimization problem, in which the upper-level objective function is generally nonconvex and the lower-level objective function is strongly convex. Although this type of problem has been studied extensively, it still remains an open question how to achieve an O(e sample complexity in Hessian/Jacobian-Fiee stochastic bilevel optimization without any second-order derivative computation. To fill this gap, we propose novel Hessian/Jacobian free bilevel optimizer named FdeHBO, which features a simple fully single-loop structure, projection-aide -idde finite-differee mation, and momentum- based updates. Theoretically, we show that FdeHBO requires O(e -15) iterations (each using (11) samples and only first-order gradient information) to find an e-accurate stationary point. far as we know, this is the first dessian/Jaccoian-fiie method with an ((c-1.) sample complexity for nonconvex-strongly-oooe convexx stochastic bilevel optimization.
Hessian/Jacobian veetoor approxi-
Introduction
Bilevel optimization has drawn intensive attention due its wide applications in meta-learning
18
50). hyperparameter optimization [8] 52 14] reinforcement learning [3]] (7), signal process 36 16
and communication [31] and federated learning (59) In this paper, we study the following stochastic bilevel optimization problem.
min ((x) = f(x (x)) [f (x))
#ERP
s.t. y" v²(r
arg min g(x, v) = E.g (x);,
MERY
where the upper- and lower-level objective functions f(x y) and g(x, y) take the expectation form the random variables & and Ç, and jointly continuously differentiable. In this paper, we focus on the sonconve--trronlyy-cowee convex and the upper-level function ((x) is nonconvex. This class of bilevel problems has been studied extensively from the theoretical perspective in recent years. Among them, [19 30 3 62]
bilevel setting, where the lower-level function g :) is strongly
proposed bilevel approaches with double-loop structure, which update and in nested manner. Single-loop bilevel algorithms have also attracted significant attention recently [27 62 34 25 g
40 due the simple updates on all variables simultaneously. Among them, the approaches in 34 25] have been shown to achieve an O(E 1.55 sample complexity, but with expensive evaluations of Hessian/Jacobian matrices Hessian/Jacobian-veee vector products.
Hessian/lacohian-free bilevel optimization has received increasing attention due to its high efficiency and feasibility in practical large-scale settings. In particular, [
48 61 directly ignored the
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Algorithm
Samples
Batch size
of iterations
Loops per iteration
PZOBO-S 58
3
O(c=1)
F2SA 37
O(e-3.))
((1)
((c-3,5)
SASA 37
O
O(1)
O(e-2.5)
FdeHBO (this paper)
O(E -15)
O(1)
Table 1: Comparison of stochastic Hessian//acooian-fee free bilevel optimization algorithms.
computation of all second-order derivatives. However, such climinations may lead to performance degeneration 2 [ and can vanish the hypergradient for bilevel problems with single-variable upper- level function, i.c., (x) = (x)) [56 23 proposed zeroth-order approaches that approximate the hypergradient using only function values. These methods do not have convergence rate guarantee.
Recently, several Hessian//acobian-fie free bilevel algorithms were proposed by [42] 57
53 g by reformulating the lower-level problem into the optimality-baseer constraints such as g y) &
miny g(x, 1). However, these approaches all focus on the deterministic setting, and their extensions to the stochastic setting remain unclear. In the stochastic case, [58] proposed evolution strategies based bilevel method, which achieves high sample complexity of O where is the problem dimension. Most recently, [37] proposed two fully first-order (i.e., Hessian/Jacobian -free)) value- function-based stochastic bilevel optimizer named F*SA and its momentum- -base version FSSA with a single-loop structure, which achieves sample complexities of O(e ((5-3.5 3.5 and O 2255 respectively. However, there is still large gap of E compared to the optimal complexity of O(e (.55., Then, an important open question, as recently proposed by [37. is:
Can we achieve (ce-15) sample/gradient complexity for noncomves-strroglyycomees bilevel optimization using only first-order gradient information?
1.1 Our Contributions
In this paper, we provide an affirmative answer to the above question by proposing a new Hessian/lacobian-feee stochastic bilevel optimizer named FdeHBO with three main features. First, FdeHBO takes the fully single-loop structure with momentum-based updates on three variables y, and x for optimizing the lower-level objective, the linear system (LS) of the Hessian- iessianninverse-vee -ecctr approximation, and the upper-level objective, respectively. Second, FdeHBO contains only single matrix-vector product at each iteration, which admits simple first-order finite-difference estimation. Third, FdeHBO involves an auxiliary projection Jessian-vectorr approximation error, the variance on momentum-based iterates, and the smoothness of the LS loss function. Our detailed contributions are summarized below.
updates to ensure the boundedness of the
Theoretically, we show that FdeHBO achieves a sample/gradient complexity of O(c-1.5) and iteration complexity of O to achieve an c-accurate stationary point, both of which outperforms existing results by a large margin. As far as we know, this is the first-known method with an O 1.5) sample complexity for vonconvex-trronnlyyconexxx stochastic bilevel optimization using only first-order gradient information.
Technically, we show that the auxiliary projection can provide more accurate iterates on v in solving the LS problem without affecting the overall convergence behavior, and in addition, provide a novel characterization of the gradient estimation error and the iterative progress during the updates, well as the impact of the y and e updates the momentum-basse hypergradient estimation, all which do not exist previous studies. addition, the finite-different approxi- mations make the unbiased assumptions in the comnentuu-based gradients no longer hold, and hence a more careful analysis is required.
As a byproduct, we further propose a fully single-loop momentum-bass based method named FMBO in the small-dimensional case with matrix-vector-beee hypergradient computations. Differently from existing momentum- based bilevel methods with (log :) Hessian-vest -esttor- vector evaluations per iteration, FMBO contains only single Hessian-vector computation per iteration with the same (c5-1-5) sample complexity.
We also want emphasize our technical differences from previous works below.

---

Comparison to existing momentum-bans bassd methods. Previous momentum- based methods [62 34 solve the linear system (LS) a high accuracy of (ee), whereas our algorithm includes new estimation error by the single-step momentum update on LS, and this error is also correlated with the lower-level updating error and the hypergradient estimation error. In addition, due to the finite- difference approximation, the stochastic gradients in all three updates on y, e, 2 are no longer unbiased. Non-trivial efforts need to be taken to deal with such challenges and derive the optimal complexity.
Comparison to existing fully single-loop methods. The analysis of the single-step momentum update in solving the LS requires the smoothness of the LS loss function and the boundedness of LS gradient variance, both of which may not be satisfied. To this end, we include an auxiliary projection and show not only guarantees these crucial properties, but also, in theory, provides an improved per-iteration progress. As comparison, existing works on fully single-loop stochastic bilevel optimization such as SOBA/SABA [1] and FL.SA 40) with a new time scale to update the LS problem often assume that the iterates on are bounded during the process. We do not require such assumptions. In addition, an (c--1i 1.5) complexity has not been established for fully single-loop bilevel algorithms yet.
1.2 Related Work
Bilevel optimization methods. Bilevel optimization, which was first introduced by [6]. has been
studied for decades. By replacing the lower-level problem with its optimality conditions, [26 20 [5]
55
reformulated the bilevel problem to the single-level problem. Gradient-baat bilevel methods have shown great promise recently, which can be divided into approximate implicit differentiation
(AID)
12
49
3 and iterative differentiation (ITD) 47
17
[52, 21] based approaches. 15]
Recently, bunch of stochastic bilevel algorithms has been proposed via Neumann series ( (0),
recursive momentum [62 28 25 and variance reduction 62
I1. Theoretically, the convergence
of bilevel optimization has been analyzed by 18. 52 45
30
27
3 [I]. Among them, [29]
provides the lower complexity bounds for deterministic bilevel optimization with (strongly-)convex
upper-level functions. [25] 62 34] achieved the near-optimal sample complexity with second-order derivative computations. Some works studied deterministic bilevel optimization with convex or
Dolyak-Lojasiewicz (PL) lower-level problems via mixed gradient aggregation [5T] 46 39], log-
barrier regularization 45 primal-dual method [57] and dynamic barrier 66 More results and
details can be found in the survey by 44.
Dessian/Jacobian-rres bilevel optimization. Some Hessian/Jacobian-ff free bilevel optimization methods have been proposed recently by [58, 43
15
23 56 48). Among them, FOMAML [5 48
and MUMOMAML [61] directly ignore the computation of all second-order derivatives. Several Hessian/Jacobian-fe
by replacing the lower- level problem with the optimality conditions as the constraints. However, these approaches focus only on the deterministic setting. Recently, zeroth-order stochastic approaches have been proposed for the hypergradient estimation [56] 23 58 Theoretically, (58) analyzed the convergence rate for their method. (37 proposed fully first-order stochastic bilevel optimization algorithms based on the value- based lower-level problem reformulation. This paper proposes new Hessian/Jacobian free stochastic bilevel algorithm that for the first time achieves an O(c sample complexity.
bilevel algorithms were proposed by 42 57
53
S
Momentum-bass based bilevel approaches. The recursive momentum technique was first introduced by 100 60 for minimization problems to improve the SGD-based updates in theory and in practice.
This technique has been incorporated in stochastic bilevel optimization [34] 9 24 25 62] These
approaches involve either gressian-inverss matrix computations or a subloop of number of iterations in the Hessian- -eetor approximation. As comparison, our proposed method takes the simpler fully single-loop structure, and only uses the first-order gradient information.
Finite-differenn matrix-vector approximation. The lnite-diiffrence matrix vector estimation has been studied extensively in the problems of escaping from saddle points ① 7 (some other related works can be found therein), neural architecture search (NAS) 43 and meta-learning [3]. However, such finite-different estimation can be sensitive to the selection of the smoothing constant, and may suffer from some numerical issues in practice 32 33 such as rounding errors. It is interesting but still open to developing fully first-order stochastic bilevel optimizer without the finite-different matrix-vector estimation. We would like to leave it for future study.
3

---

Algorithm 1 Hessian/lacobian-free Bilevel Optimizer via Projection-aided Finite-difference Estima-
tion
1: Input: {at, 8r, de A= and Tv-
2: Initialize:
3: fort== 1, .... --1 do
4: Compute the gradient estimator 63 by eq.
and update you1 yt Br11..
5: Compute the gradient estimator A. by eq. and update w+++1 delie.
We+1,
|| || S ros
6: Set U++1
[iwe+1ll [| > Tu-
7: Compute the gradient estimator bof by eq. (10) and update F++1 = -ochi
8:
end for
Algorithms
In this section, we first describe the hypergradient computation in bilevel optimization, and then present the proposed Hessian/Jacobian- bilevel method.
2.1 Hypergradient Computation
One major challenge in bilevel optimization lies in computing the hypergradient ($(x) due to the implicit and complex dependence of the lower-level minimizer y" on æ. To see this, if is twice differentiable, v, is continuously differentiable and the Hessian V Vy (2r, y (x)) is invertible, using the implicit function theorem (IFT) 22 5
the hypergradient ( takes the form of
( = V v²) - 22 y²) (yyg((
(z,
(2)
Note that the hypergradient in eq. requires computing the exact solution y" and the expensive Hessian inverse y y²) To approximate this hypergradient efficiently, we define the following (stochastic) hypergradient surrogates
(f =V (x, V2 22,
ff v.f y: v2 (1)
where E R8 is an auxiliary vector to approximate the Hessian-inverse-vectt product in eq. (2) and fx y/, u; {) can be regarded stochastic version of fz, y, v). Based on eq. (3), one needs to find an efficient estimate of y" e.g., via an iterative optimization procedure, as well as feasible estimate of the solution e* [V g y)] ( y) of linear system (LS) (equivalently quadratic programming) whose generic loss function is given by
1 V
(Linear system loss:) R(I, 2
V y)u V.f
(4)
where the gradient of R(x, y, v) w.r.t. is given by
V R(r, V Vigg(r 1) V
Similarly to eq. 3 we also define R,,,
= 1 v) Vyf(2, v: over any sample a as a stochastic version of V. V.,(Ir, y. v) in eq. 5 It can be seen from eq. (3). eq. (4) and
eq. (5) that the updates on the LS system involve the Hessian- and Jacobian-vectte vector products, which can be computationally intractable in the high-dimensional case. In the next section, we propose novel stochastic Hessian/lacobian-free bilevel algorithm.
2.2 Hessian/Jacobian-free Bilevel Optimizer via Projection-aided Finite-difference Estimation
As shown in Algorithm T we propose a fully single-loop stochastic Hessian/Jacobian free bilevel optimizer named FdeHBO via projection-aided liiit--differenc estimation. It can be seen that FdeHBO first minimizes the lower-level objective function g(r, y) w.r.t. by running a single-step momentum-basee update as Vt+1 = yt - B,hg where B, is the stepsize and hg is the momentum-based gradient estimator that takes the form of
h = V Vygg(E, ()) + +11- n7") (M
Veg(zi, Mes V.) - V V00(21-1.91-1)
(6)

---

where of E [0, is a tuning parameter. The next key step is to deal with the LS problem via solving the quadratic problem eq. w++1 Us Xehe with the momentum-based gradient h given by
h - A Ri,, ye vn, de v.) + - Th TT ( R. yes 2t, 8,s v/)
M0 u-,, 4/)),
where V.R is a Hessian-ffre version of the LS gradient V R in eq. ( given by
(First-order LS gradient:) V. R(r,, Ut. uz dz; 1/t) H (2:, Ut, vt, deg w/aa) VyJ ye; thx).
(8)
Note that in the above eq. (8). H (ar, yr, or, de: vtr) is the finite differdiffe estimation of the Hessian- vector product V ( yo! Ur )u, which takes the form of
(I Mc, U2, 8.1 V Vyg(rs, y dev v)) Vug(E Ve 5,44; v/2a)
26.
where 5. > 0 is a small constant. Note that in eq. ( if the iterative 2, unbounded, the approx- imation error between H and Vy yr: thy Ju can be uncontrollable as well. We further prove lemma
5 in appendix B that the bound of this gap relies on ||w:]| and but it is independent of the dimension of yt- To this end, after obtaining W our key step in line introduces an auxiliary projection on a ball (which can be generalized to any convex and bounded domain) with radius of ry
Wr+1s Pu [
| ||| S Tei |we++1| || > Tg-
(Auxiliary projection)
V++1
This auxiliary projection guarantees the boundedness ue, O,... T 1, which serves three important purposes. First, it ensures the smoothness of the LS loss function R(x, v) in eq. (5) w.r.t. all x. y and u, which crucial in the convergence analysis of the momentum based updates. Second, the boundedness of 0, also ensures that the estimation variance of the stochastic LS gradient V,dl(a,, Ne, 0ri 1/2) does not explode. Third, it guarantees the error of the finite-difference Hessian- vector approximation to be sufficiently small with proper ò,. We will show later that under proper choice of the radius Tu. this auxiliary projection provides better per-step progress, and the proposed algorithm achieves stronger convergence performance. Finally, for the upper-level problem, the momentum-based hypergradient estimate hig is designed as
Tf =m V (It, Ye, U2 8,i E.) + (1 T/ (h/_1 Vf Ur, Vf U1- v1-1: E)),
(10)
where Vf(z, E.) is the fully first-order hypergradient estimate evaluated at two consecutive iterates (rr, ye. 02) and (x:-1, Ut-1, 11-1) is given by
Vf y. &, =Va((m v; - Jxx y, v. de)
and J(z y, Et) the finite-differen Jacobian- approximation given by
i(z v. v. Vag(x,y 8.v V Vag(z, -0:0.v-660::5) 6,v;
(11)
28.
Note that R.R and f f are biased estimators of the gradients V. R and f. which further compli- cates the convergence analysis on the momentum-basee updates because the conventional analysis the recursive momentum requires the unbiased gradient estimation to ensure the variance reduction ef- fect. By controlling the perturbation 8, properly, will show that FdeHBO can achieve O(e 1.5) convergence and complexity performance without any second-order derivative computation.
2.3 Extension Small-Dimensional Case
As byproduct of our proposed FdeHBO, we further propose fully single-loop momentum-basse bilevel optimizer (FMBO), which more suitable in the small-dimensional case without finite- FMBO first takes the same lower-level updates
difference approximation. As shown in Algorithm 2
5

---

Algorithm 2 Fully Single-loop Momentum-based Bilevel Optimizer (FMBO)
1: Input: {au, Sr, Ae )= and Tu.
2: Initialize:
3: for = 0, 1, ..., do
Compute the gradient estimator hy by eq.
(6) and update y/+1 yt+ -.hh
5: Compute the gradient estimator by by eq. (2) and update Wc+1 = Ur Xehi
W1+1,
[| II S ruj
6: Set U2
1e [ > ru-
Compute gradient estimator by by oq. (33 and update 2x++1 = Jz achi
8: end for
on yr as in eq. 6 Then, it solves the LS problem as Wx+1 = Ut dehe where the momentum-based gradient estimator given by
R R V R(2,, yo, Uti 1h) + (1 mn (h. + V 26 h T
(12)
where differently from FdeHBO, we here use the precise gradient V, R without finite-difference approximation. Similarly to FdeHBBO, we add an auxiliary projection on the is updates to ensure the LS smoothness and bounded variance. Finally, for the upper-level problem, we optimize I based a momentum-based update as 24+1 = Jt oghy with the hypergradient estimator
bf =nd V ( Me, uzs Ex) (1 (6) Vf Ue, vei E. V Mt- 16-1 (13)
where nee E [0. is tuning parameter. Similarly, we directly use the hypergradient estimate in eq. (3) without the finiteddifference estimation. We note that compared to existing momentum- algorithms [62 34 that contains Ollog 4) steps in solving the LS problem, FMBO takes fully single-loop structure with a single-step nomentuu-baseed acceleration on the LS updates.
Main Results
3.1
Assumptions and Definitions
We make the following standard assumptions for the upper- and lower-level objective functions, as also adopted by (30 9 34]. The following assumption imposes the Lipschitz condition on the
upper-level function (x.y).
Assumption 1. For any x E Rd and Rd,, there exist positive constants L L Ofz and Cj such that V Vx y) and V v,,(I, y) are Ly-- and Lj, -Lipschit continuous w.r.t. (x,, y), and (v, < Ct., [V,[(a, < CC,,
The following assumption imposes the Lipschitz condition on the lower-level function g( y). Assumption 2. For any 8d and y Rd,, there exist positive constants Hg. Lg Lgny' L Lays
Cgep Capy such that
Function g( y) is twice continuously differentiable;
Function g(z,,) -) is coos:::
The derivatives V Vyg(I,8) ÿ). V2 V2,g(e,)) and V V2 are Lyz, Lazv" and Lgv. -Lipschitz continuous w.r.l. (x, y);
1V 1) < Cge,, and V 10z y)1 < Cam.*
The following assumption is adopted for the stochastic functions f(x, ys E) and g( C).

---

Assumption 3. Assumptions and hold for fz, Vi and g(z, y: f) V and Ç. Moreover we
assume that there exist positive constants o fz- fy og gay and 0 guv such that
E a v) VnI( y: 11 - i.
E yI V u (r, y: 0) S oj
E [IV (x - V Veg v: 11 _ og
E [IIV (0)( 22 V: 9 S 0z
E V gy y) Vz, y: O)I - "uu""
Definition 1. We say & is an e-accurate stationary point of a function ((x) if ( S where I is the output of an optimization algorithm.
3.2 Convergence and Complexity Analysis of FdeHBO
We further provide the convergence analysis for the proposed Hessian/Jacoiaa-ffree FdeHBO algo- rithm. We first characterize several estimation properties of FdeHBO. Let ef := Tif V (1 Mts Vt) A (xx, Nt. vt) denote the hypergradient estimation error.
Proposition 1. Under Assumption
3 the iterates of the outer problem by Algorithm 7 satiisy
E| 211e/111 2 S [11 12(11)² + 4 4L :. |lle/ + 441l411 + (4Lg., r. + 161 sax r2)))
+ 6(1 - )
3 / 1| + 217 [B|| + V v (
+ +0 yov x [ R2 || + r. - ni w.]13]]
for all t E tO....,T 1} with L3 = 2L3 + 1 r
The hypergradient estimator error O(EI|<(+111 ) contains three main components. The first term [11 - + + 11² indicates the per-iteration improvement induced by the momentum- based update, the error term E |12 is caused by the 2, updates, the error term (B + [V,,( {)) is caused by solving the lower-level problem, and the new error term ( L. - vi [1=) is induced by the one-step momentum update on the LS problem, which does not exist in previous momentum-based bilevel methods 62
34 25 that solve the LS problem to a high accuracy. Also note that the errors A Blle 1|² and 4.g.vre0 + 16L2 16
r252 are caused by the finiteddiferrenc approximation error. Fortunately, by choosing the perturbation level 6, in these two terms to be properly small, it can guarantee the descent factor (1 n )2 4.g.,, 7856 to be at
an order of (1 - O(T(+1)2 and hence the momentum-based variance reduction effect is still applied.
Proposition 2. For v/, define ef == Tif - V R.. Vt, vt). Under Assumptions 7 21 < <|( E
2
we have
m (2(11 + 966 + 4L y.u 25 EBle + (4L yuv -.. + SL $L9.0 (6)
+ 8 2 10= 2/ oJ,) + 96(1 / )a )
+ (
of 11²)
Dvv
+ 96
R
)2 ( Suy
2 + 12;,) a / 11 + 23 (E e E V Vyg(a, Vx( }11}
for all L E 10.,1,
As shown Proposition the LS gradient estimation error c contains an iteratively improved error component (1-m)t+1>² ( + 2.5. factor 5, sufficiently small, a finiteddifference approximation error ((6)) as well as an approximation error Ox};Bl1e 122) for solving the LS problem. The next step is to upper-bound Ellu; of 112
for the stepsize X and the approximation
Proposition 3. Under the Assumption
7 the iterates of the LS problem by Algorithm satisfy
El|ez+1
11
(Ln + L3)0g
$ (1 (1+ 700 (1++:))( [:))((1 21,
x2T.2) - v2[13]
Mg + Lg
A3Bllef112 + + + 7 ( 1+ 8,
+1+ )
1
22
207,
L2
[aRE 112 83 (2Elle/11 + 2E Vyg(ze, 31) 13)
71
11
7

---

for all t E T 1) with some % > 0 and 8, > 0.
provide the general convergence theorem for Theorem 1. Suppose Assumptions 1 and Lemma are satisfied. Choose Iv 2 Ctz and set
Based on the above important properties, FdeHBO.
a: (w + tw
P. cg
X exas,
Veé Cay ai The = Cyn
T/ = Cna a2
min min(eng Cog)
and 6, S B(Looy in the appendix. Then, the iterates generated by Algorithm 2 satissf
(-taa-yjary where the constants w, cB. Cx. Cnyy Can and Cyp are defined in eq. (67
(xro)
"
113o y" (zx))
1100 2" (xo, ( vo
E V
T2/3
T2/3
T2/3
of
2
o²
+
9
+
T2/3
T2/3
T2/3
T2/3
Corollary Under the same setting of Theorem FdeHBO requires (c--1.b) samples and gradient
evaluations, respectively, to achieve an c-accurate stationary point.
can be seen from Corollary that the proposed FdeHBO achieves an ( 1.5) sample complexity without any second-order derivative computation. As far we know, this is the first Hessian/Jacobian- free stochastic bilevel optimizer with an O ((c-1.5) sample complexity.
3.3 Convergence and Complexity Analysis of FMBO
In this section, we analyze the convergence and complexity of the simplified FMBO method. Theorem 2. Suppose Assumptions
Ciy and set paramelers
/ and 3 are satisfied. Choose Pa
0t
Bt = cgat, (w + 1)2/3* v/( Cay oi
d cx00t
v/ cna The = Cn, of
where w, CR. C)2 Cafs CaR and On, are defined in eq. by Algorithm 2 satisfy
(33 in the appendix. The iterates generated
(#xo) +
Ilyo
y" (xo))1
E ( <ō
T2/3
T2/3
v² (zo. vo)
o7
R
T23
72/3
T2/3
T2/3
Theorem 2 shows that the proposed fully single-loop FMBO achieves convergence rate of 22/3 which further yields the following complexity result.
Corollary 2. Under the same setting of l'heorem 2 FMBO requires totally (c-1.5) data samples, gradient and matrix-vector evaluations, respectively, achieve e-accurate stationary point.
Corollary shows that FMBO requires total number (c-11 of data samples, which matches the best sample complexity 34 62 28). More importantly, each iteration of FMBO contains only one Hessian-vector computation due to the simple fully single-loop implementation, whereas other momentum-based approaches require Oflog 4) Hessian-vector computations nested manner per iteration. Also, note that FMBO is the first fully single-loop bilevel optimizer that achieves the (c-11.5) sample complexity.
4 Experiments
In this section, we test the performance of the proposed FdeHBO and FMBO on two applications: hyper-representation and data hyper-cleaning, respectively.

---

PSA
PSA
Outer
PSA
PSA
P2080-5
Figure 1: Comparison syyer-representation with the LeNet neural network. Left plot: outer loss V.s. running time; right plot: accuracy v.s. running time.
4.1 Hyper-represennttion MNIST Dataset
We now compare the performance of our Hessian/Jacobian-fre FdeHBO with the relevant Hessian/Jacobian-fres methods PZOBO-S (58. FSSA (37 and FASA [37]. We perform the hyper-
representation with the 7-layer LeNet network problem.
38 which aims to solve the following bilevel
min (A) |
Σ
(E
s.t. w² (A) arg min Lin (A, w), Lin u)) w)
E
S
(T.Y))
where Low denotes the cross-entropy loss, S, and S, denote the training data and validation data, and 24) denotes the features extracted from the data r,. More details of the experimental setups are specified in Appendix A
As shown in Figure our FdeHBO converges much faster and more stably than PZOBO-S, FSA and SAS while achieving a higher training accuracy. This is consistent with our theoretical results, and validates the momentum-based approaches in reducing the variance during the entire training.
4.2 Hyper-cleaning on MNIST Dataset
We compare the performance of our FMBO various bilevel algorithms including AID-FP 21]. reverse 17 SUSTAIN 34), MRBO and VRBO (22, BSA [9), stocBiO [30], FSLA 40 and SOBA M. on low-dimensional data hyper-cleanine problem with a linear classifier on MNIST dataset, which takes the following formulation.
min Tap (A, wx)
Σ (Z(.):)S
Lor((w") (cce(w")Txxaaul 2x,,
|Sall
s.t. arg min L(A, w) = ST (H..V:)ES,
Σ a(A))C r..3) vi) + C/|
(14)
where LoE denotes the cross-entropy loss, S, and S, denote the training data and validation data, whose sizes are set to 20000 and 5000, respectively, À {Aihiss and C are the regularization
parameters, and (-) the sigmoid function. AmIGO [3] is not included the figures because
performs similarly to stocBiO. The experimental details can be found in Appendix A.2
As shown in Figure ((a) FMBO, stocBiO and AID-FP converge much faster and more stable than other algorithms. Compared to stocBiO and AID-FP, FMBO achieves a lower training loss. This demonstrates the effectiveness of momentum- based variance reduction in finding more accurate iterates. It can be seen from Figure 2(b) that FMBO converges faster than existing fully single-loop FSLA and SOBA algorithms with lower training loss.
5 Conclusion
In this paper, we propose a novel Hessian/Jacohian-feee bilevel optimizer named FdeHBO. We show that FdeHBO achieves an (c-51.) sample complexity, which outperforms existing algorithms of the

---

85A
PILA
VIDOO MEDO
VIDO
HBSC
PABO
PHBO
Rasning Time
Running Times
Running Time
noise p 0.1
(h) Single-loop methods
Figure (a) Comparison different algorithms on data hyper-cleaning with noise p = 0.1. Left plot: test loss v.s. running time; right plot: train loss v.s. running time. (b) Comparison among different single-loop algorithms: training loss v.s. running time.
same type by a large margin. Our experiments validate the theoretical results and the effectiveness of the proposed algorithms. We anticipate that the developed analysis will shed light on developing provable Hessian/Jacobian-free bilevel optimization algorithms and the proposed algorithms may be applied to other applications such as fair machine learning.
Acknowledgement
The work supported in part by NSF under grants 2326592 and 2311274.
References
[1] Z. Allen-Zhu and Y. Li. Ncon2: Finding local minima via first-order oracles. Advances in Neural Information Processing Systems, 31, 2018.
[2] A. Antoniou, H. Edwards, and A. Storkey. How to train your MAML. In International Conference on Learning Representations (ICLR), 2019.
[3] M. Arbel and J. Mairal. Amortized implicit differentiation for stochastic bilevel optimization. International Conference on Learning Representutions (ICLR), 2022.
[4] L. Bertinetto, J. F. Henriques, P. Torr, and A. Vedaldi. Meta-learning with differentiable closed-form solvers. In International Conference on Learning Representations (ICLR), 2018.
[5] M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert. Efficient and modular implicit differentiation. arXiv preprint arXiv:2105.1518. 2021.
[6] J. Bracken and J. T. McGill. Mathematical programs with optimization problems in the constraints. Operations Research, 21(1):37 44, 1973.
[7] Y. Carmon, J. C. Duchi, O. Hinder, and A. Sidford. Accelerated methods for nonconvex optimization. SIAM Journal on Optimization, 28(2):151-11772 2018.
[8] L. Chen, J. Xu, and J. Zhang. On bilevel optimization without lower-level strong convexity, arXiv preprint arXiv:2301.00712. 2023.
[9] Chen, Sun, and W. Yin. single-timescal stochastic bilevel optimization method. arXiv preprint 1121022102.0467,, 2021.
[10] A. Cutkosky and F. Orabona. Momentum-basee variance reduction in non-convex SGD. In Advances in Neural Information Processing Systems (NeurIPS). 2019.
[11] M. Dagréou, P. Ablin, S. Vaiter, and T. Moreau. A framework for bilevel optimization that enables stochastic and global variance reduction algorithms. arXiv preprint arXiv:2201. 13409, 2022.
10

---

[12] J. Domke. Generic methods for optimization based modeling, In Artificial Intelligence and Statistics (AISTATS), pages 318-326, 2012.
[13] A. Fallah, Mokhtari, and Ozdaglar. On the convergence theory of gradient gradient-based model- agnostic meta-learning algorithms. In International Conference Artificial Intelligence and Statistics (AISTATS), pages 1082-1092. PMLR, 2020.
[14] M. Feurer and Hutter. Hyperparameter optimization. In Automated Machine Learning, pages 3-33. Springer, Cham, 2019.
[15] C. Finn, P. Abbeel, and S. Levine. Modcl-agnostic meta-learning for fast adaptation of deep networks. Proc. International Conference on Machine Learning (ICML), pages 1126-1135, 2017.
[16] R. Flamary, A. Rakotomamonjy, and G. Gasso. Learning constrained task similarities in graphregularized multi-task learning. Regularization, Optimization, Kernels, and Support Vector Machines, page 103, 2014.
[17] L. Franceschi, M. Donini, P. Frasconi, and M. Pontil. Forward and reverse gradient-based hyperparameter optimization. In International Conference on Machine Learning (ICML), pages 1165-1173, 2017.
[18] L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In International Conference on Machine Learning (ICML), pages 1568-1577, 2018.
[19] S. Ghadimi and M. Wang. Approximation methods for bilevel programming. arXiv preprint 1.XXv:1802.02255 2018.
[20] S. Gould, B. Fernando, A. Cherian, P. Anderson, R. S. Cruz, and E. Guo. On differentiating parameterized argmin and argmax problems with application bi-level optimization. arXiv preprint 1607.055477 2016.
[21] R. Grazzi, L. Franceschi, M. Pontil, and S. Salzo. On the iteration complexity of hypergradient computation. In Proc. International Conference on Machine Learning (ICML), 2020.
[22] A. Griewank and A. Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. SIAM, 2008.
[23] B. Gu, G. Liu, Y. Zhang, X. Geng, and H. Huang. Optimizing large-scale hyperparameters via automated learning algorithm. arXiv preprint arXiv:2102.090066, 2021.
[24] Z. Guo. Y. Xu, W. Yin, R. Jin, and T. Yang. On stochastic moving- average estimators for non-convex optimization. arXiv preprint arXiv:2104. 14840, 2021.
[25] Z. Guo and T. Yang. Randomized stochastic variance-reduced methods for stochastic bilevel optimization. arXiv preprint arXiv:2 (05.02266, 2021.
[26] P. Hansen, B. Jaumard, and G. Savard. New branch-and-bound rules for linear bilevel program- ming. SIAM Journal on Scientific and Statistical Computing, 13(5):1194-12) 1992.
[27] M. Hong, H.-T. Wai, Z. Wang, and Yang. A two-tiesscaee framework for bilevel optimization: Complexity analysis and application to actor-critic. arXiv preprint arXiv:2007.05170, 2020.
[28] F. Huang and H. Huang. Biadam: Fast adaptive bilevel optimization methods. arXiv preprint arXiv:2106. 1139. 2021.
[29] K. and Y. Liang. Lower bounds and accelerated algorithms for bilevel optimization. arXiv preprint arXiv: 2021.
[30] K. Ji, J. Yang, and Y. Liang. Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning (ICML), pages
PMLR, 2021.
[31] K. and I.. Ying. Network utility maximization with general and unknown utility functions: distributed, data-driven bilevel optimization approach. Submitted, 2022.
11

---

[32] C. Jin, R. Ge. P. Netrapalli, S. M. Kakade, and M. I. Jordan. How to escape saddle points efficiently. In International conference on machine learning, pages 1724-1732. PMLR, 2017.
[33] N. Jorge and J. W. Stephen. Numerical oplimization. Spinger, 2006.
[34] P. Khanduri, S. Zeng, M. Hong, H.-T. Wai, Z. Wang. and Z. Yang. A near-optimal algorithm for stochastic bilevel optimization via Processing Systems (NeurIPS), 34:30271 001:300771-02885 300883, 2021.
Advances in Neural Information
[35] V. R. Konda and J. N. Tsitsiklis. Actor-critic algorithms. Advances in neural information processing systems (NeurIPS), pages 1008-1014, 2000.
[36] G. Kunapuli, K. Bennett, Hu, and J.-S. Pang. Classification model selection via bilevel programming. Optimization Methods & Software, 23(4):475 489, 2008
[37] Kwon, D. Kwon, S. Wright, and R. Nowak. A fully first-order method for stochastic bilevel optimization. arXiv preprint arXiv:2301. 10045,, 2023.
[38] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner. Gradient- based learning applied document recognition. Proceedings of the IEEE, 86(11) (1:(11)))7822) 2224,, 1998.
[39] J. Li, B. Gu, and H. Huang. Improved bilevel model: Fast and optimal algorithm with theoretical guarantee. arXiv preprint arXiv:2009. 0.00::0000900690, 2020.
[40] J. Li, B. Gu, and H. Huang. A fully single loop algorithm for bilevel optimization without hessian inverse. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 7426-7434, 2022.
[41] R. Liao, Y. Xiong. Fetaya, L. Zhang, K. Yoon, X. Pitkow, R. Urtasun, and R. Zemel. Reviving and improving recurrent back- poo-proopaattio. In Proc. International Conference on Machine Learning (ICML), 2018.
[42] B. Liu, M. Ye, S. Wright, P. Stone, and Q. Liu. Bome! bilevel optimization made easy: A simple first-order approach. Advances in Neural Information Processing Systems, 35:17248-17262. 2022.
[43] H. Liu, K. Simonyan, and Y. Yang. Darts: Differentiable architecture search. In International Conference on Learning Representations (ICLR), 2018
[44] R. Liu, J. Gao. J. Zhang, D. Meng, and Z. Lin. Investigating bi-level optimization for learning and vision from a unified perspective: A survey and beyond. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.
[45] R. Liu, X. Liu, X. Yuan, Zeng, and Zhang. A cellue-functionbarn based interior-point method for non-convex bi-level optimization. In International Conference on Machine Learning (ICML), 2021.
[46] R. Liu, P. Mu, X. Yuan, S. Zeng, and J. Zhang. A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton. In International Conference on Machine Learning (ICML), 2020.
[47] D. Maclaurin, D. Duvenaud, and R. Adams. Gradient-based hyperparameter optimization through reversible learning. In International Conference on Machine Learning (ICML), pages 21 2122, 2015.
[48] A. Nichol, J. Achiam, and J. Schulman. On first-order meta-learning algorithms. arXiv preprint arXiv:1803.02999 2018.
[49] Pedregosa. Hyperparameter optimization with approximate gradient. In International Conference Machine Learning (ICML), pages 737-746, 2016.
[50] A. Rajeswaran, C. Finn, S. M. Kakade, and S. Levine. Meta-learning with implicit gradients. In Advances in Neural Information Processing Systems (NeurIPS), pages 113-124, 2019.
12

---

[51] S. Sabach and S. Shtern. A first order method for solving convex bilevel optimization problems. SIAM Journal on Optimization, 27(2):640-660) 2017.
[52] A. Shaban, C.-A. Cheng, N. Hatch, and B. Boots. Truncated back-propagation for bilevel optimization. In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 1723-1732, 2019.
[53] H. Shen and T. Chen. On penalty-based bilevel gradient descent method. arXiv preprint a3Xiv:2302.05185, 2023.
[54] C. Shi, J. Lu, and Zhang. An extended kuhn-tucker approach for linear bilevel programming. Applied Mathematics and Computation, 162(1):51-63, 2005.
[55] A. Sinha. P. Malo, and K. Deb. A review on bilevel optimization: from classical to evolutionary approaches and applications. IEEE Transactions on Evolutionary Computation, 22(2):276-295, 2017.
[56] X. Song, W. Gao, Y. Yang, K. Choromanski, A. Pacchiano, and Y. Tang. ES-MAML: Simple Hessian- meta learning. In International Conference on Learning Representations (ICLR), 2019.
[57] D. Sow, K. Ji, Z. Guan, and Y. Liang. A constrained optimization approach bilevel optimiza- tion with multiple inner minima. arXiv preprint arXiv:2203.01 123, 2022.
[58] D. Sow, K. Ji, and Y. Liang. On the convergence theory for Hessian-free bilevel algorithms. Advances Neural Information Processing Systems (NeurIPS). 2022.
[59] D. A. Tarzanagh, M. Li, C. Thrampoulidis, and Oymak. Fednest: Federated bilevel, minimax, and compositional optimization. arXiv preprint arXiv:2205.02215. 2022.
[60] Q Tran-Dinh, N. H. Pham, D. T. Phan, and L. M. Nguyen. Hybrid stochastic gradient descent algorithms for stochastic nonconvex optimization. arXiv preprint arXiv:I90S. 190.009900. 2019.
[61] R. Vuorio, S.-H. Sun, H. Hu, and J. J. Lim. Multimodal model-agnostic meta-learning via task-aware modulation. Advances in Neural Information Processing Systems (NeurIPS), 32, 2019.
[62] J. Yang, Ji, and Y. Liang. Provably faster algorithms for bilevel optimization. Advances in Neural Information Processing Systems (NeurIPS), 34:13670-1368,, 2021.
[63] M. Ye, B. Liu, S. Wright, P. Stone, and Q Liu. Bome! bilevel optimization made easy: A simple first-order approach. In Conference on Neural Information Processing Systems (NeurIPS), page 1, 2022.
13


Reviews:
[{"review_id": "g9obPgrEQs", "reviewer": "Reviewer_Xzpd", "summary": "The paper studies stochastic bilevel optimization, where the first level objective function depends on the optimizer of the second level objective function. The paper provides an optimization procedure involving gradient computations and projections. The paper provides rates of convergence for the method. Since convexity is not assumed, convergence is expressed in terms of the gradient becoming small (not in terms of converging to a global optimizer). The benefit of the algorithm is demonstrated on various data sets. ", "questions": "Algorithm 1: Hessain -> Hessian.\n\nConclusion, why is it written that the sample complexity is O(epsilon^-2), while Corollary 1 writes O(epsilon^-1.5)?\n\nAppendix A.1. What are lambda and w? What are the dimensions? \n\n", "limitations": "The authors adequately discuss limitations.\nI think there is not concern about societal impact. ", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The technical contribution of the paper appears to be impressive. The algorithm suggested is involved. Its novelties and benefits compared to the state of the art are discussed extensively. The convergence proofs appear to be involved and innovative. The numerical performances are very good, especially when compared to the state of the art. ", "weaknesses": "The content is very technical, and it is not very easy to follow, especially given the space constraint in the main text. \nIn the appendix, the description of the experiments is also difficult to follow. There there are no space constraints, so I think the description could have been more pedagogical there.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank reviewer Xzpd for the time and valuable feedback! \n\n**Q1: Suggestion on presentations.**\n\n**A:** We will definitely follow your suggestions to improve the presentation.\n\n**Q2: Conclusion, why is it written that the sample complexity is** $O(\\epsilon^{-2})$**, while Corollary 1 writes** $O(\\epsilon^{-1.5})$**?**\n\n**A:** Sorry about this typo. The sample complexity should be $O(\\epsilon^{-1.5})$ in the conclusion. We will revise it.  \n\n**Q3: In Appendix A.1, what are $\\lambda$ and $w$? What are the dimensions?**\n\n**A:** We perform the hyper-representation with the 7-layer LeNet network to solve the bilevel problem. We take the last two layers as the lower-level parameters $w$, and all remaining layers as the upper-level parameters $\\lambda$. In our experiments, the dimension of $w$ is 850, and the dimension of  $\\lambda$ is 60856. \n"}, {"review_id": "tAi1A40TWW", "reviewer": "Reviewer_HtEu", "summary": "This work focuses on addressing a stochastic bilevel optimization problem characterized by a strongly convex lower-level objective function and a nonconvex upper-level objective. The authors specifically investigate a Hessian-free (fully first-order) approach to bilevel optimization, which avoids the need for Hessian or Jacobian-vector computations. They introduce a new method called FdeHBO, which adopts a single-loop structure and incorporates momentum-based updates for the upper-level problem, linear system, and lower-level problem.\nFdeHBO achieves a sample and gradient complexity of $O(\\epsilon^{-1.5})$, which is the first result for nonconvex-strongly-convex Hessian-free bilevel optimization. Additionally, the authors propose a simplified and nearly optimal variant called FMBO, which is specifically designed for scenarios with small-dimensional problems. To validate the effectiveness of their proposed methods, the authors conduct several experiments.", "questions": "See Weakness. \n\n$\\nabla_y^2 g$ in Line 126 should be corrected to $\\nabla_{yy}^2 g$?", "limitations": "yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is overall easy to follow. The topic of Hessian-free bilevel optimization, particularly in the nonconvex-strongly-convex setting, has not been extensively explored. Therefore, the investigation into the development of simple and nearly optimal stochastic Hessian-free bilevel algorithms is both interesting and significant.\n2. To the best of my knowledge, the proposed FdeHBO algorithm achieves near-optimal sample/gradient complexity for the first time. The utilization of the projection-aided finite-difference estimator in bilevel optimization is interesting. The experimental results depicted in Figure 1 provide substantial evidence of FdeHBO's superior performance compared to other approaches of the same kind.\n3. From a technical standpoint, addressing issues such as momentum-based errors, biased estimation caused by the finite-difference estimator, and auxiliary projection errors requires nontrivial treatments, further highlighting the sophistication of this work.", "weaknesses": "1. The paper's projection technique is limited to projecting onto a ball. Is it possible that  this approach can be generalized to encompass projections onto more general convex and closed sets?\n\n2. While this work primarily focuses on the momentum-based estimator, it appears that the approach can be simplified to a more general stochastic gradient descent (SGD)-type method by setting $\\eta_t^f = \\eta_t^g = \\eta_t^R = 0$. Does the analysis presented in the paper adequately cover this specific case? It is worth considering the importance of simpler SGD updates with reduced tuning requirements, as they often hold practical significance.\n\n3. The authors should provide further elaboration on why the auxiliary projection does not impact the final convergence guarantee. Although I have not thoroughly examined all the proofs, a clearer explanation on this aspect would be beneficial. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer HtEu for the time and valuable feedback! \n\n**Q1: The paper's projection technique is limited to projecting onto a ball. Is it possible that this approach can be generalized to encompass projections onto more general convex and closed sets?**\n\n**A:** Great question! The answer is yes. Our proof only requires this projection set to be convex and contain the optimal solution $v_t^*$. Since $v_t^*$ is bounded by $r_v$ (see Lemma 1), the proof still holds if we choose the convex set large enough. We will comment on this in the revision. \n\n**Q2: While this work primarily focuses on the momentum-based estimator, it appears that the approach can be simplified to a more general stochastic gradient descent (SGD)-type method by setting $\\eta_t^f = \\eta_t^g = \\eta_t^R = 0$. Does the analysis presented in the paper adequately cover this specific case? It is worth considering the importance of simpler SGD updates with reduced tuning requirements, as they often hold practical significance.**\n\n**A:** Good question! As the reviewer noted, our framework reduces to SGD-type updates if we set $\\eta_t^f = \\eta_t^g = \\eta_t^R = 0$, but the analysis may require proper adjustments in learning rate selection and variance characterization of gradient estimations. We would like to leave this for future study.\n\n**Q3: The authors should provide further elaboration on why the auxiliary projection does not impact the final convergence guarantee.**\n\n**A:** The projection appeals in two key steps in bounding the iterative gaps $E\\\\|v_{t+1} - v_t\\\\|^2$ in (25) and $E\\\\|v_{t+1} - v_t^*\\\\|^2$ in (28). In specific, by using the non-expansive property of projection that $\\\\|P_B(x)-P_B(y)\\\\|\\leq \\\\|x-y\\\\|$ for a projection $P$ on a set $B$, we bound these two gaps by $\\mathbb{E}\\\\|v_{t+1} - v_t\\\\|^2 \\leq \\mathbb{E}\\\\|w_{t+1} - v_t\\\\|^2 = \\mathbb{E}\\\\|\\lambda_t h_t^R\\\\|^2$ and  \n$\\mathbb{E}\\\\|v_{t+1} - v_t^*\\\\|^2 \\leq \\mathbb{E}\\\\|w_{t+1} - v_t^*\\\\|^2 = \\mathbb{E}\\\\|v_t - \\lambda_t h_t^R - v_t^*\\\\|^2$. The remaining steps are free from projection and the convergence will not be affected. We will clarify this in the revision. \n\n**Q4: $\\nabla_y^2 g$ in line 126 should be corrected to $\\nabla_{yy}^2 g$?**\n\n**A:** Yes, we will correct it in the revision. "}, {"review_id": "KKlFAe1rU2", "reviewer": "Reviewer_ySP3", "summary": "This paper provides a novel Hessian-free bilevel optimization algorithm called FdeHBO, which improves upon previous state of the art F3SA with a better iteration complexity. The author provide convergence analysis and empirical results to show the strength of the proposed method. However, I have some concerns regarding the proof, which is detailed in the weaknesses part. ", "questions": "See weakness part above.", "limitations": "I do not see the discussion of limitations, but I think this is not a major problem. I do not see potential negative societal impact as this paper mainly focuses on the theoretical side.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The algorithm is novel, as far as I am concerned, and the improved convergence rate seems promising. The paper is well organized and easy to follow. The experiment seems convincing. ", "weaknesses": "The proof of Lemma 5 is incorrect. In the proof, \nthe authors attempt to use the mean value theorem for a function $\\mathbb{R}\\to \\mathbb{R}^n$, but this is in general wrong. The estimation error of the finite difference of in the lemma is a key step, without which the proof of the main theorem is incorrect. \n\nIn general, a zeroth order optimization algorithm should have a dimension dependency, but this paper magically removes this. I believe the one of the reason is the mistake in Lemma 5. I will be glad to change my evaluation if I am wrong. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer ySP3 for the time and valuable feedback! \n\n**Q1: The proof of Lemma 5 is incorrect. In the proof, the authors attempt to use the mean value theorem for a function $R \\rightarrow R^n$, but this is in general wrong. The estimation error of the finite difference in the lemma is a key step, without which the proof of the main theorem is incorrect.**\n\n**A:** We thank the reviewer for pointing this out for us! For rigorous derivations, we provide a very detailed proof below for Lemma 5, which includes a generalized mean value theorem for vector-valued function, as well as a refined upper bound on the finite-difference estimation error. It can be seen that our characterizations on the estimation error $\\\\|e_t^J\\\\|^2$ are still bounded by a constant $C_J$ multiplied with a tunable sufficiently small parameter $\\delta_\\epsilon$ (similarly for the error  $\\\\|e_t^H\\\\|^2$). Thus, the remaining steps will not be affected. We will add such details in the revision. \n\n**Q2: In general, a zeroth order optimization algorithm should have a dimension dependency, but this paper magically removes this. I believe one of the reasons is the mistake in Lemma 5.**\n\n**A:** We would like to explain this dimension-free result from three perspectives. \n\nFirst, intuitively speaking, a typical zeroth-order method has such dimension dependence because it uses a one-dimension function value to approximate all $n$ coordinates of a $n$-dimension gradient vector. However, in our case, the finite-difference approach uses the $n$-dimension gradient information to approximate the $n$-dimension Hessian-vector product rather than the $n^2$-dimension Hessian matrix itself. More technically, one key property here is that the smoothing direction vector in our estimator is exactly the same as the vector of the Hessian-vector product. This is why the dimension dependence disappears in our results. \n\nSecond, rigorously speaking, we provide detailed proof below. It can be seen that the upper bound does not contain any dimension dependence. \n\nThird, we also note that similar results have been derived in previous studies, e.g., see the inequality right after eq.(140) in [1], where it can be seen that the finite-difference estimation error bound does not contain any dimension dependence either. \n\n[1] Fallah, Alireza, Aryan Mokhtari, and Asuman Ozdaglar. \"On the convergence theory of gradient-based model-agnostic meta-learning algorithms.\"\n\n---\n\n[**The mean-value theorem for vector-valued function**]\nLet $S \\subseteq \\mathbb{R}^n$ be open and let $f: S \\rightarrow \\mathbb{R}^m$ be differentiable on all of $S$. Let $x,y \\in S$ be such that the line segment connecting these two points is contained in $S$, i.e., $L(x,y) = \\\\{ (1-t)x+ty: t\\in [0,1]\\\\} \\subseteq S$. Then for every $\\mathbf{a} \\in \\mathbb{R}^m$, there exists a point $z \\in L(x,y)$ such that $\\mathbf{a} \\cdot [f(y)-f(x)] = \\ \\mathbf{a} \\cdot [f'(z)(y-x)]$. \n\nProof: Let $\\mathbf{a} \\in \\mathbb{R}^m$ and define a new function $F:[0,1] \\rightarrow \\mathbb{R}$ for all $t \\in [0,1]$ by: $F(t) = \\mathbf{a} \\cdot f(x+t(y-x))$. Since $f$ is differentiable on $S$, we have \nthat $F$ must be continuous on $[0,1]$. Furthermore, $F$ is differentiable on $(0,1)$ by the chain rule: $F'(t) = \\mathbf{a} \\cdot f'(x+t(y-x))(y-x)$. So by the Mean Value Theorem for single-variable real-valued functions, there exists a number $h \\in (0,1)$ for which: \n$$F(1) - F(0) = F'(h)(1-0),$$ \nwhere the left-hand side is \n$$F(1) - F(0) = \\mathbf{a} \\cdot f(y) - \\mathbf{a} \\cdot f(x)$$\nand the right-hand side is\n$$F'(h)(1-0) = \\mathbf{a} \\cdot f'(x+h(y-x))(y-x).$$\nSet $z = x+h(y-x)$, then $z \\in L(x,y)$ and we have from the equality that \n$$\\mathbf{a} \\cdot [f(y) - f(x)] = \\mathbf{a} \\cdot f'(z)(y-x).$$\n\n[**Rewrite Lemma 5**] \nUnder Assumption 2, for any $\\xi_t$ and $\\psi_t$, define $$e_t^J:= \\widetilde{J}(x_t, y_t, v_t, \\delta_{\\epsilon};\\xi_t) - \\nabla^2_{xy}g(x_t,y_t;\\xi_t)v_t, \\quad e_t^H:= \\widetilde{H}(x_t, y_t, v_t, \\delta_{\\epsilon};\\psi_t) - \\nabla^2_{yy}g(x_t,y_t;\\psi_t)v_t,$$\nthen we have the bounds of $e_t^J$ and $e_t^H$ as\n$$\\\\|e_t^J\\\\|^2 \\leq C_J\\delta_{\\epsilon}, \\quad \\\\|e_t^H\\\\|^2 \\leq C_H\\delta_{\\epsilon}, $$\nwhere $C_J := 4L_gL_{g_{xy}}r_v^3$, $C_H :=4L_gL_{g_{yy}}r_v^3$. \n\nProof: Here, in our case, we have \n$$\\mathbf{a} \\cdot \\bigg[\\frac{\\nabla_y g(x_t,y_t+\\delta_{\\epsilon}v_t;\\zeta_t) - \\nabla_y g(x_t,y_t-\\delta_{\\epsilon}v_t;\\zeta_t)}{2\\delta_{\\epsilon}} - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)v\\bigg] \\ = \\mathbf{a} \\cdot \\Big[\\nabla_{yy}^2 g(x_t,y_t+\\delta_{\\epsilon}'v_t;\\zeta_t) - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)\\Big] $$\nLet $\\mathbf{a} = \\bigg[\\Big(\\nabla_y g(x_t,y_t+\\delta_{\\epsilon}v_t;\\zeta_t) - \\nabla_y g(x_t,y_t-\\delta_{\\epsilon}v_t;\\zeta_t)\\Big)/(2\\delta_{\\epsilon}) - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)v_t\\bigg]^T$, then we have \n\n\\begin{align}\n\\\\|e_t^H\\\\|^2 &= \\bigg[\\frac{\\nabla_y g(x_t,y_t+\\delta_{\\epsilon}v_t;\\zeta_t) - \\nabla_y g(x_t,y_t-\\delta_{\\epsilon}v_t;\\zeta_t)}{2\\delta_{\\epsilon}} - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)v_t\\bigg]^T  \\cdot \\Big[\\nabla_{yy}^2 g(x_t,y_t+\\delta_{\\epsilon}'v_t;\\zeta_t) - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)\\Big]\n\\end{align}\n\\begin{align}\n\\leq \\bigg\\\\|\\frac{\\nabla_y g(x_t,y_t+\\delta_{\\epsilon}v_t;\\zeta_t) - \\nabla_y g(x_t,y_t-\\delta_{\\epsilon}v_t;\\zeta_t)}{2\\delta_{\\epsilon}} \\bigg\\\\| \\cdot \\Big\\\\|\\nabla_{yy}^2 g(x_t,y_t+\\delta_{\\epsilon}'v_t;\\zeta_t) - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)\\Big\\\\| \n\\end{align}\n\\begin{align}\n \\+ \\big\\\\|\\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)v_t\\big\\\\| \\cdot \\Big\\\\|\\nabla_{yy}^2 g(x_t,y_t+\\delta_{\\epsilon}'v_t;\\zeta_t) - \\nabla_{yy}^2 g(x_t,y_t;\\zeta_t)\\Big\\\\|\n\\end{align}\n\\begin{align}\n\\overset{(a)}{\\leq} 4L_gL_{g_{yy}}\\\\|v_t\\\\|^3\\delta_{\\epsilon}' \\leq 4L_gL_{g_{yy}}r_v^3\\delta_{\\epsilon} ,\\nonumber\n\\end{align}\nwhere (a) uses the Lipschitz continuity of $\\nabla_y g(x,\\cdot)$ and $\\nabla^2_{yy} g(x,\\cdot)$ in Assumption 2.\nSimilarly, we have $\\\\|e_t^H\\\\|^2 \\leq 4L_gL_{g_{xy}}r_v^3\\delta_{\\epsilon}$."}, {"review_id": "CEXtLgc0ke", "reviewer": "Reviewer_GD5q", "summary": "This paper proposed faster first-order stochastic algorithms for bilevel optimization, based on the idea that a Hessian-vector product can be implemented by finite difference of the gradient along the given direction.\n", "questions": "I have a question about the claim of \"near-optimal''.\nAssumption 3 supposes that \"Assumptions 1 and 2 hold for $f(x,y;\\xi)$ and $g(x, y;\\zeta)$ for $\\forall \\xi$ and $\\forall \\zeta$.\"\nBut according to \"Lower Bounds for Non-Convex Stochastic  Optimization\"(https://arxiv.org/pdf/1912.02365.pdf), the lower bound function is constructed under the mean-square-smooth Assumption. Therefore it seems that the mean-square-smoothness should be imposed for $f$, which is slightly stronger than the smoothness of $f(x,y;\\xi)$ for $\\forall \\xi$?\n\nI think it may be a fixable bug. But the paper may need some modifications.\n\nPlease let me know if I am wrong. I will adjust my ratings based on the rebuttal.", "limitations": "N/A", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper studies a very important problem and the provided complexity upper bound strictly improves previous works.", "weaknesses": "1. It is well-known that the Hessian-vector product can be implemented by the finite difference of the gradient along the given direction.\nAnd the computation of a Hessian-vector product is known to be as easy as the computation of gradients. Therefore the results and algorithms seem not to be novel compared to the existing methods (AID/ITD) that use Hessian-vector products as oracles.\n\n2. The paper uses three different names for the same baseline, including F${}^2$SA, F2SA, $F{}^2$SA  (Line 442, Table 1, Figure 1). Should use a unified name.\n\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer GD5q for the time and valuable feedback! \n\n**Q1: The paper uses three different names for the same baseline, including FSA, F2SA, FSA (Line 442, Table 1, Figure 1). Should use a unified name.**\n\n**A:** Many thanks! We will use a unified name in the revision.\n\n**Q2: I have a question about the claim of \"near-optimal''. Assumption 3 supposes that \"Assumptions 1 and 2 hold for $f(x,y; \\xi)$ and $g(x,y; \\zeta)$ for $\\forall \\xi$ and $\\forall \\zeta$. But the lower bound function is constructed under the mean-square-smooth assumption. Therefore it seems that the mean-square-smoothness should be imposed for $f$, which is slightly stronger than the smoothness of $f(x,y;\\xi)$ for $\\forall \\xi$?**\n\n**A:** Good question! Note that if our assumption on the smoothness of $f(x,y;\\xi)$ for $\\forall \\xi$ holds, then the mean-square-smooth assumption (eq. (4) in [2] pointed out by the reviewer) must hold. This means that the problem class considered in our paper is a subset of that in [2]. Thus, if there exists a case in our problem class that achieves a lower bound smaller than $\\Omega(\\epsilon^{-1.5})$, then since this case also belongs to that of [2], the lower bound in [2] should also be smaller than $\\Omega(\\epsilon^{-1.5})$. This contradicts the result in [2]. Thus, by contradiction, we can show that the lower bound under our assumptions is larger than or equal to $\\Omega(\\epsilon^{-1.5})$. Since our proved upper bound of $O(\\epsilon^{-1.5})$ matches this $\\Omega(\\epsilon^{-1.5})$ lower bound, we can claim that our algorithm is “near-optimal”. \n\n[2] Arjevani, Yossi, et al. \"Lower bounds for non-convex stochastic optimization.\" Mathematical Programming 199.1-2 (2023): 165-214.\n\n"}]

Author rebuttal:

Meta-review: This paper proposes a bilevel optimization algorithm that do not use the Hessian of the outer objective and that it has the best complexity known complexity. It was overall well received by the reviewers, and the rebuttal allowed to fix some important concerns. 

It is **absolutely crucial** that the camera-ready acknowledges the fact that FdeHBO is not proved to be near-optimal by changing the title and the core of the text. It is also important that the authors discusses the fact if being Hessian-free is a *good* property. Indeed the computation of an HVP is marginally higher than a JVP. The paper does not really discuss this aspect, either from a theoretical aspect or from an experimental point of view.
Decision: Accept (poster)