LancBiO: Dynamic Lanczos-aided Bilevel Optimization via Krylov Subspace
Authors: Yan Yang
Bin Gao
Ya-xiang Yuan
Venue: ICLR 2025 Poster (2025)
Keywords: Bilevel Optimization
Lanczos Process
Krylov Subspace
Primary area: optimization

Abstract: Bilevel optimization, with broad applications in machine learning, has an intricate hierarchical structure. Gradient-based methods have emerged as a common approach to large-scale bilevel problems. However, the computation of the hyper-gradient, which involves a Hessian inverse vector product, confines the efficiency and is regarded as a bottleneck. To circumvent the inverse, we construct a sequence of low-dimensional approximate Krylov subspaces with the aid of the Lanczos process. As a result, the constructed subspace is able to dynamically and incrementally approximate the Hessian inverse vector product with less effort and thus leads to a favorable estimate of the hyper-gradient. Moreover, we propose a provable subspace-based framework for bilevel problems where one central step is to solve a small-size tridiagonal linear system. To the best of our knowledge, this is the first time that subspace techniques are incorporated into bilevel optimization. This successful trial not only enjoys $\mathcal{O}(\epsilon^{-1})$ convergence rate but also demonstrates efficiency in a synthetic problem and two deep learning tasks.

Published as a conference paper at ICLR 2025
LANCBIO: DYNAMIC LANCZOS-AIDED BILEVEL OPTI- MIZATION VIA KRYLOV SUBSPACE
Yan Yang
Bin Gao" & Ya-xiang Yuan
Academy of Mathematics and Systems Science Chinese Academy of Sciences University of Chinese Academy of Sciences yangyan&amss ac.
Academy of Mathematics and Systems Science
Chinese Academy of Sciences
(gaobin, yyx @lsec. cc ac..c.cc
ABSTRACT
Bilevel optimization, with broad applications in machine learning, has an intri- cate hierarchical structure. Gradient-based methods have emerged as common approach to large-scale bilevel problems. However, the computation of the hyper- gradient, which involves Hessian inverse vector product, confines the efficiency and regarded as a bottleneck. To circumvent the inverse, we construct a sequence of now-dimensional approximate Krylov subspaces with the aid of the Lanczos process. As a result, the constructed subspace able to dynamically and incremen- tally approximate the Hessian inverse vector product with less effort and thus leads to a favorable estimate of the hyper-gradient. Moreover, we propose a provable subspace-baset framework for bilevel problems where one central step is to solve a small-size tridiagonal linear system. To the best of our knowledge, this is the first time that subspace techniques are incorporated into bilevel optimization. This successful trial not only enjoys Occ=l) convergence rate but also demonstrates efficiency synthetic problem and two deep learning tasks.
INTRODUCTION
Bilevel optimization, in which upper-level and lower-level problems are nested with each other, mirrors multitude of applications, c.g., game theory (Stackelberg, 1952), hyper-parameter selection (Ye et al., 2023), data poisoning (Liu al., 2024), meta-learning (Bertinetto al., 2018), neural architecture search (Liu et al., 2018; Wang et al., 2022), adversarial training (Wang et al., 2021). reinforcement learning (Chakraborty et al., 2024; Thoma et al., 2024; Yang et al., 2025), computer vision (Liu al., 2021a). In this paper, consider the bilevel problem:
min xxx) (x,
t.t. t.. (xx) E arg min ( v). yER",
where the upper-level function and the lower-level function g are defined on R&d Rdv, op called the hyper-objective, and the gradient of sp(x) is referred to as the hyper-gradien (Pedregosa, 2016; Grazzi et al., 2020; Yang et al., 2023; Chen et al., 2024) it exists. In contrast to standard single-level optimization problems, bilevel optimization is inherently challenging due to its intertwined structure. Specifically, the formulation (1) underscores the crucial role of the lower-level solution y" (ar) in each update of or.
One of the focal points in recent bilevel methods has shifted towards nonconvex upper-level problems
coupled with strongly convex lower-level problems (Ghadimi and Wang, 2018; Ji et al., 2021: Dagréou et al., 2022; Li et al., 2022; Hong et al., 2023; Hu et al., 2024). This configuration ensures that y (x) single-valued function of x, i..., y (z) arg min men ( y). Subsequently, Vy((x) can be computed via the implicit function theorem following (Ghadimi and Wang, 2018).
V = Va. (I, (x)) V (x, y² (x)) [V (x, y (x)) vyd (x, (x))
(2)
"Corresponding author.

---

Published as a conference paper at ICLR 2025
10
1-10
1.20
6x 10
1650
1+152
Residual norm
10:2
10:1
SusBIO
AmioO-D
5086
1000
2000 3000
4000
5000
40
120
140
160
Time [sec]
Figure 1: Left: test loss for the method stocBiO (Ji al., 2021) with different inner iterations I to approximate the Hessian inverse vector product; Right: Estimation error of the Hessian inverse vector product in hyper-data cleaning task with corruption rate 0.5 for different methods: LancBiO and SubBiO (ours), AmIGO (Arbel and Mairal, 2022), and SOBA (Dagréou et al., 2022).
Gradient methods based on the hyper-gradient, #x+1 = 2x - -V are known as the approximate implicit differentiation (AID) based methods (Ji et al., 2021; Liu et al., 2023b; Huang, 2024). However, the computation of the hyper-gradien (2) suffers from two pains: 1) solving the lower-level problem
obtain y" (x); 2) assembling the Hessian inverse vector product
(x)) = Vygg (z, y² (x)) Vuf (2, (x)),
equivalently, solving large linear system in terms of
( y (z))u= y=(x
(4)
To this end, beneficial adopt few inner iterations to approximate yy (x) and (x)) within each outer iteration (i.e., the update of 2). Note that the approximation accuracy of u* (x() is crucial for AlD-based methods; see (Ji et al., 2022; Li et al., 2022). Specifically, the left of Figure confirms that the more inner iterations, the higher quality of the estimate of v", and the more enhanced descent of the objective function within the same number of outer iterations.
Approximation: Existing efforts dedicated approximating v" in different fashions by regulating the number of inner iterations, e.g.. the Neumann series approximation (Ghadimi and Wang, 2018; Ji al., 2021) for the inverse, gradient descent (Arbel and Mairal, 2022; Dagréou al., 2022) and conjugate gradient descent (Pedregosa, 2016; Yang et al., 2023) for the linear system.
Moreover, there studies aimed amortizing the of approximation through outer iterations. These methods include using the inner estimate from the previous outer iteration as warm start for the current outer iteration (Ji al., 2021; Arbel and Mairal, 2022; Dagréou et al., 2022; Ji et al., 2022; Li et al., 2022; Xiao et al., 2023), or employing a refined step size control (Hong al., 2023).
Amortization:
Subspace techniques, widely adopted in nonlinear optimization (Yuan, 2014), approximately solve large-scale problems in lower-dimensional subspaces, which not only reduce the computational cost significantly but also enjoy favorable theoretical properties in full space models. Taking into account the above two principles, it is reasonable to consider subspace techniques bilevel optimization. Specifically, and sequentially solve linear systems (4) in these subspaces to approximate e* accurately.
efficiently amortize the construction of low-dimensional subspaces
1.1 CONTRIBUTIONS
In this paper, taking advantage the Krylov subspace and the Lanczos process, we develop an inno- vative subspace-bai based Snnework-Lanchhe which features an efficient and accurate approximation of the Hessian inverse vector product e" in the eyper-gradientt for bilevel optimization. The main contributions summarized as follows.
Firstly, we build up a dynamic process for constructing now-dimensional subspaces that are tailored from the Krylov subspace for bilevel optimization. This process effectively reduces the large-scale
2

---

Published as a conference paper at ICLR 2025
subproblem (4) to the small-size tridiagonal linear system, which draws on the spirit of the Lanezos process. To the best of our knowledge, this is the first time that the subspace technique is leveraged in bilevel optimization.
Moreover, the constructed subspaces enable us to dynamically and incrementally approximate o across outer iterations, thereby achieving enhanced estimate of the hyper- sepee-gradiet; the right of Figure illustrates that the proposed LancBiO reaches the best estimation error for v" Hence, we provide new perspective for approximating the Hessian inverse vector product in bilevel optimization. Specifically, the number of Hessian-vector products averages at (1 + 1/m) per outer iteration with the subspace dimension m, which is favorably comparable with the existing methods.
Finally, we offer analysis to circumvent the instability in the process of approximating subspaces, with the result that LancBiO can profit from the benign properties of the Krylov subspace. We prove that the proposed method LancBiO is globally convergent with the convergence rate O(c-l), ). In addition, the efficiency of LancBiO is validated by a synthetic problem and two deep learning tasks.
1.2 RELATED WORK
A detailed introduction to bilevel optimization methods can be found in Appendix A.
Krylov subspace methods: Subspace techniques have gained significant recognition in the realm of numerical linear algebra (Parlett, 1998; Saad, 2011; Golub and Van Loan, 2013) and nonlinear optimization (Yuan, 2014; Liu et al., 2021c). Specifically, numerous optimization methods utilized subspace techniques to improve efficiency, including acceleration technique (Lii et al., 2020), diagonal preconditioning (Gao et al., 2023), and derivative-free optimization methods (Cartis and Roberts, 2023). Krylov subspace (Krylov, 1931), due its special structure,
XCx (A, := span {b. Ab, A AN-1b)
with the dimension N for a matrix A and a vector b, exhibits advantageous properties in convex quadratic optimization (Nesterov al., 2018), eigenvalue computation (Kuczyáski and Wozniakowski, 1992), and regularized nonconvex quadratic problems (Carmon and Duchi. 2018). Krylov subspace has been widely considered in large-scale optimization such trust region methods (Gould et al., 1999), trace maximization problems (Liu et al., 2013), and cubic Newton methods (Cartis et al., 2011; Jiang al., 2024). Lanezos process (Lanezos, 1950) is an orthogonal projection method onto the Krylov subspace, which reduces a dense symmetric matrix to tridiagonal form. Details of the Krylov subspace and the Lanczos process are summarized Appendix
Approximation of the Hessian inverse vector product: It is cumbersome to compute the Hessian inverse vector product in bilevel optimization. bypass it, several strategies implemented through inner iterations were proposed, e.g., the Neumann series approximation (Ghadimi and Wang, 2018; et al., 2021), gradient descent (Arbel and Mairal, 2022; Dagréou et al., 2022), and conjugate gradient descent (Pedregosa, 2016; Arbel and Mairal, 2022; Yang et al., 2023). Alternatively, the previous information was exploited in (Ji et al., 2021; Arbel and Mairal, 2022) as a warm start for outer iterations; Ramzi et al. (2022) suggested approximating the Hessian inverse in the manner of quasi-Newton; Dagréou et al. (2022) and Li et al. (2022) proposed the frameworks without inner iterations approximate the Hessian inverse vector product.
2 LUBSPACE-BASED ALGORITHMS
In this section, we dive into the development of bilevel optimization algorithms for solving (I), which dynamically construct subspaces to approximate the Hessian inverse vector product.
The (hyper-)gradien descent method carries out the k-th outer iteration as Fb+1 = 2s - AVp(())) where the dyper-gradions is exactly computed by Vi((zz) V v.f (æx, vi) V (Ix, vy) oi with V := y (Ix) and vt ;= v (Ix) defined in (3). In view of the computational intricacy of of and ot. is commonly concerned with the following estimator for the hyper-gradient
Vep (Nz, Uk: uk) Vz (ZR, 9K) V (IK; yk. Uk:
3

---

Published as a conference paper at ICLR 2025
where Ui is approximation of vi. Denote A V (Z Wt) and b (zk, yk). The vector 2g serves as the (approximate) solution of the quadratic optimization problem
1
min Ag of ba: "ES& 2
where Sk is the full space Rd and the exact solution is A Subsequently, in order to reduce the computational cost, is natural to ask:
Can we construct low-dimensional subspace S₁. such that the solution of the quadratic problem (6) satisfactorily approximates A
General subspace constructions introduced in the existing subspace methods (Yuan, 2014; Liu et al.. 2021c) are not straightforward and not exploited in the bilevel setting, rendering the exploration of appropriate subspaces challenging. In the following subsections, we construct approximate Krylov subspaces and propose an elaborate subspace-based framework for bilevel problems.
2.1 WHY KRYLOV SUBSPACE: THE SUBBIO ALGORITHM
= In light of the Neumann series for a suitable E IR, b- E - nA)'b, is observed from Appendix B that A -1 belongs to Krylov subspace for some N i...
A K.v (A, (N n),
Hence, it is reasonable to consider a Krylov subspace for the construction of Sp.
Given a constant & N, consider approximation of A-1b in lower-dimensional Krylov subspace K,, (A, b). i... Vn (A, b) K,,(1 .1,, and v, Σi= G (1 my) b A-'b. Note that the approximation en is composed of the set {(T 7)) by=nd in the sense of the Neumann series. Moreover, observe that (I n)A)),, E (x,1(A,1) and hence we recursively choose
Un+1 Snt1 := span (b, A)),, S KA,1(A,b)
since subspace span (b, (I 1))),,} includes the information of the increased set - {(--)4) o) In summary, we can construct sequence of two-dimensional subspaces (S,) that implicitly filters information from the Krylov subspaces. The rationale for this procedure can be illustrated in Figure 2.
In the context of bilevel optimization, seek the best solution Ux to the subproblem (6) in the subspace
A-16
A-1
SN
Sn
Sp41
A 1b
S = span (be, (T- (T-
(7)
Repeating the procedure is capable of dy- namically approximating the Hessian inverse vector product, i.c., 8/2 approximates A by- The Krylov Subspace-aided Bilevel Opti- mization algorithm (SubBiO) listed Al- gorithm 1.
(A, C. Kn-i(A,b) C C (A,
Figure 2: Illustration of approximating A- E KNN (A. b) by 2n in the two-dimensional sub- space S, C (A,b).
Algorithm SubBiO
Input: iteration threshold K, step sizes 0, X, 17, initialization F1, V/, Up
1: for k = 1.2,. 1.2...., K do
2: A V2 b (IK, m1k), by V v (rk,,
3: S₁ span (be, (I nAAs) Uk
4:
arg min 3 A&U bir
5: PE+1 =
(V (T yk) - Vay Mh )bk
6: yk+11 = Uk - oV g 1k+1. yk)
7: end for
Output: (UK+1,9K+1
4

---

Published as a conference paper at ICLR 2025
2.2
WHY DYNAMIC LANCZOS: THE LANCBIO FRAMEWORK
Notice that the subproblem in SubBiO (Algorithm 1) can be equivalently reduced to
1 TSp AgSx)= by SE2, zE 2
min
where Sx = [bx (1-))A)) vk-1] Rdyx2, The solution =" R²² results in es Sx.=". It is two- dimensional subproblem, whereas computing the projection S A₂ requires two Hessian-vector products, which dominate the cost of the subproblem. Therefore, is crucial to reduce or amortize the projection cost while preserving the advantages of the Krylov subspace. To this end, we find that the Lanezos process (Appendix B) provides an enlightening way (Lanczos, 1950; Saad, 2011; Golub and Van Loan, 2013) since it allows for the construction of a Krylov subspace and maintaining tridiagonal matrix as the projection matrix, which significantly reduces the computational cost.
In bilevel optimization, since the quadratic problem (6) evolves through outer iterations, it is difficult to leverage the Lanczos process to amortize the projection cost while updating variables as in SubBiO. Specifically, the Lanczos process is inherently unstable (Paige, 1980), and thus the accumulative difference among } and {ba} will make the Lanczos process invalid.
In order to address the above difficulties, we propose restart mechanism to guarantee the benign behavior of approximating the Krylov subspace and consider solving residual systems to employ the historical information. In summary, we propose dynamic Lanczos-aided Bilevel Optimization framework, LancBiO, which is listed in Algorithm 2. The only difference between LancBiO and SubBiO is solving the subproblem (line 3-4 in Algorithm 1).
we adapt the standard Lanczos process for tridiagonalizing Ak. by starting from QI / 1, Qo
Algorithm 2 LancBiO
Input: iteration threshold K, step sizes 0, X, initialization 21: V1, U1 initial correction AND = O, subspace dimen- sion m, initial epoch h -1
O, B1 = 0 and using the dynamic ma- trices A, for j = 1, ....,, k. the pro- cess as follows.
1: for A: = 1.2..... K do
2:
Ag V ( yk), by
V Vuf(aa Vk)
U Ag4; Bz9j-1;
3:
if (A: mod m) then
0 9j Uj.
4:
= h +1
W Wj a,ıs
5:
Th Uk
6:
= A&ū
9--1 = ||ررا|| ||,
7:
Q 1 (br. wx)/ |6 wall
4j+1 w;/ wu;//jt1
8:
T = Empty Ma Matrix
9:
3R = 0
and T3 is tridiagonal matrix recur- sively computed from
10:
end if
11:
(Tss, 1e.Bk-1) Diancooooiocc Qu-1, Ak, Bx)
12:
Tk bn w/
13:
Auk = Qx QT
0
T
14:
Us Ün + Avk
T =
15:
Jx+1 = Ik (V.I( (lyg(ok,,
8
16: Mc++1 = - OV v yg
0
17: end for
Output: (œK*I, BK41)
Restart mechanism:
contrast to SubBiO, we construct the subspace
S₁ span Q2+) span((g...........
Consequently, Q; approximates the basis of the true Krylov subspace KC (Ak. b&), and Tş approxi- mates the projection Aş onto it, i.c., Tr x Q Ax However, Qx will lose the orthogonality due to the evolution of Aj. and Ti will deviate from the true projection. Based on this observation, we restart the subspace spanned by the matrix Q for every m outer iterations (line to line 10 in Algorithm 2). The restart mechanism allows to mitigate the accumulation of the difference among (A)..., A& }. and hence, we can maintain a more reliable basis matrix to approximate Krylov subspaces. The above dynamic Lanezos subroutine, DLanczos, is summarized in Appendix C.
5

---

Published as a conference paper at ICLR 2025
Residual minimization: The preceding discussion reveals that the dimension of the subspace Sk should be moderate to retain the reliability of the dynamic process. However, the limited dimension stagnates the approximation of the subproblem (6) to the full space problem. Therefore, instead of directly solving the quadratic subproblem (6) in the subspace, we intend to find v E Sg span such that Agu by- To this end, after going through m outer iterations, denote the current approximation by U. In the subsequent outer iterations, we concentrate on a linear system with residual in the form
Agãe = by AR Au E Sk. Specifically, we inexactly solve the minimal residual subproblem
(8)
min AvES
[1(bx. A. AxAull
and use the solution Ar, correction to (line 12 to line 14 in Algorithm 2), c2 6 Ace.
Consequently, taking into account the two strategies above, we illuminate the framework LancBiO in Figure 3. It is structured into epochs, with each epoch built by outer iterations. Notably, each epoch restarts by incrementally constructing from one- mieedimesioonl space Q1 to an m-diimensional space Qm, aiming approximate the solution the residual system within these subspaces. The solution Auş to the subproblem serves as correction to enhance the hyper-gradient estimation, which facilitates the (xr,)) updating.
The combination of the two strategies, restart mechanism and residual minimization, not only controls the dimension of the subspace but also utilizes historical information to enhance the ap- proximation accuracy. By considering simpli- fied scenario, we reduce the two strategies into solving a standard linear system problem Ar = b with A and fixed. Note that, from the perspec- tive of Theorem in (Carmon and Duchi, 2018), the residual associated with solving a linear sys- tem in an mn-dimensional Krylov subspace decays faster than rate O (1/m²) after each restart. In other words, the estimation error of the Hessian inverse vector product experiences decay rate of O (1/m²) after every restart, i.e.,
Epochs
One Epoch
m Outer Iterations
Subspace Construction
|bb AUm(A+121
(b A AA m
Aumill 1|
1
Ava
Q1 Q2
m2
Residual Minimization
Remark 2.1. The classic Lanczos process is known for its capability to solve indefinite linear systems (Greenbaum et al., 1999). In similar fashion, the LancBiO framework can be adapted to the bilevel problems with nonconvex lower- level problem. Interested readers are referred to Appendix D for details.
by
A
Au An
Auy
A Amm
Hyper-gradient Estimation
V V2
V Pm
2.3 RELATION TO EXISTING ALGORITHMS
(x.y) Updating
The proposed SubBiO and LaneBiO have intrin- sic connections the existing algorithms. Gen- erally, in each outer iteration, methods such as BSA (Ghadimi and Wang, 2018) and TTSA (Hong et al., 2023) truncate the Neumann series at N, exploiting information from an N-dimensional Krylov subspace. In contrast, both SubBiO and LancBiO implicitly gather knowledge from a nigher-dimmesiooall Krylov subspace with less effort.
Figure 3: An overview of LancBiO.
SubBiO shares similarities with SOBA (Dagréou et al., 2022) and FSLA (Li et al., 2022). The update rule for the estimator of the Hessian inverse vector product in SOBA and FSLA is
Uk Uk-1 (AA
bb) 7 Uk-1 + 17b::

---

Published as a conference paper at ICLR 2025
while the proposed SubBiO constructs two-dimensional subspace, S₁ span {bx, (I ndR - defined in (7). It is worth noting that the updated 2k in SOBA and FSLA belongs to the subspace Sg. Furthermore, in the sense of solving the two-dimensional subproblem (6), SubBiO selects the optimal solution v in the subspace.
In addition, if the subspace dimension is set one, LancBiO is simplified scenario in which one conjugate gradient (CG) step with the warm start mechanism is performed each outer iteration, which exactly recovers the algorithm AmIGO-CG (Arbel and Mairal, 2022) with one inner iteration to the update of v. Alternatively. if the step size in Algorithm 2 is set to 0 within each m-steps (i.e., only inner iterations are invoked), LancBiO reduces the algorithm AmIGO-CG (Arbel and Mairal, 2022) with m inner iterations.
3 THEORETICAL ANALYSIS
In this section, provide non-asymptotic convergence analysis for LancBiO. Firstly, we introduce some appropriate assumptions. Subsequently, to address the principal theoretical challenges, we analyze the properties and dynamics of the subspaces constructed in Section Finally, we prove the global convergence of LancBiO and give the iteration complexity; the detailed proofs are provided in the appendices.
Assumption 3.1. The upper-level function f is twice continuously differentiable. The gradients Vef(z,)) y) and V, fr, y) are Lip-sihct and L -ipsciip and |V (z, y" vx)))] s Cfy. Assumption 3.2. The lower-level function is twice continuously differentiable. Vzg(z,)) y) and Vyg(z. y) are Lgæ-Lipschinn and Lip-Sipschitz. The derivative V V2yg(1, y) and the Hessian matrix V (2 3) are Lgay-Lipschite and Lips-iihch...
Assumption 3.3. For any Rd=, the lower-level function g(z,)) :) is -tr-nnoy convex.
The Lipschitz properties of f.g g and the strong convexity of the lower-level problem revealed by the above assumptions are standard in bilevel optimization (Ghadimi and Wang. 2018; Chen et al., 2021; Ji et al., 2021; Khanduri et al., 2021; Arbel and Mairal, 2022; Chen al., 2022; Dagréou et al., 2022; Li et al., 2022; Ji et al., 2022; Hong et al., 2023). These assumptions ensure the smoothness of yp and y"; see the following results (Ghadimi and Wang, 2018).
Lemma 3.4. Under the Assumptions 3.2 and 3.3, y" (x: is Lgz/ My -Lipschitz continuous, i.e., for any æ1, 22 E RR,, (x1) - (x2)) < Lug || Ø211-
Lemma 3.5. Under the Assumptions 3.1, 3.2 and 3.3, the hyper- gradient Vi((-) Lip-schh continuous, i.e. for any $1,42 E IRd=, [|V V S Lu || - x2I1, where L.p > is defined in Appendix E.
Assumption 3.6. There exists a constant Cfz so that fz, v)ll - Cja-
Assumption 3.6, commonly adopted in (Ghadimi and Wang, 2018; 8; Ji et al., 2021; Liu et al., 2022: Kwon et al., 2023), is helpful in ensuring the stable behavior of the dynamic Lanczos process; see Section 3.1.
3.1 SUBSPACE PROPERTIES DYNAMIC LANCZOS PROCESS
In view of the inherent instability of the Lanezos process (Paige, 1980; Meurant and Strakoš, 2006) and the evolution of the Hessian {AA} and the gradient {br} in LaneBiO, the analysis of the constructed subspaces is intricate. Based on the existing work (Paige, 1976; 1980; Greenbaum, 1997), this subsection sheds light on the analysis of the subspaces and the effectiveness the subproblem in approximating the full space problem in LancBiO.
An epoch is constituted of a complete m-step dynamic Lanczos process between two restarts, namely, after h epochs, the number of outer iterations is mh. Given the outer iterations A mh + j for j = 1, 2, ...,, m,, we denote
e : (1+
Lgz
|| - Fmattell + |mmk+ss - Vib4oll
Mg
7

---

Published as a conference paper at ICLR 2025
for s, s, 1, ..... and ei maxiga,te Eae, serving the accumulative difference. For brevity, we omit the superscript where there is no ambiguity, and we are slightly abusing of notation that at the current epoch, (Amb+j) and {6mb+t} } are simplified by {Aj) and (b,} for j = 1..... m. In addition, the approximations in the residual system (8) are simplified by & and b == b1 A1ū.
The following proposition demonstrates that the dynamic subspace constructed in Algorithm within an epoch is indeed approximate Krylov subspace.
Proposition 3.7. At the j-th step within an epoch (j = -1,2,... ..... 1), the subspace spanned by the matrix Qj+1 in Algorithm 2 satisfies
span(Q)+1) S span {A9) A2"
A7 b o, or 1.2.......
Specifically, when A = A2 = A₁ and Qj+1 of full rank, span
span(QQ+1) (Q++1) = Kjtt (A.
Denote A = V (2,((,)") and b' = vf fxr Notice that the dynamic Lanczos process in Algo- rithm 2 centers on A, instead of AT. The subsequent lemma interprets the perturbation analysis for the dynamic Lanczos process in terms of A,, which satisfies an approximate three-term recurrence with perturbation term óQ.
Lemma 3.8. Suppose Assumptions 3.1 to 3.3 hold. The dynamic Lanczos process in Algorithm with normalized qn satisfies
,2Q; = QiT; + 85+145+11 + +02, for 11,2...,,
where Qj = [0,, ............... V.). 8Q, löq1 892, ..,,,...,00 5gg] with [|6g;1| S Lyg and Tj is j x symmetric tridiagonal matrix with diagonal elements ( ; and subdiagonal elements (B2....,B)
3.2 CONVERGENCE ANALYSIS
To guarantee the stable behavior of the dynamic process, we need the subsequent assumption.
Assumption 3.9. The initialization of y1 in Algorithm satisfies|lyr vill
/3p
8(m 1) Lowv
Similar initialization refinement is used in (Hao al., 2024), which be achieved by implementing several gradient descent steps for the smooth and strongly convex lower-level problem. The following lemma reveals that the dynamic process yields improved solution for the subproblem (9).
Lemma 3.10. Suppose Assumptions 1, 3.2, 3.3, 3.6 and 3.9 hold. Within each epoch, if set the step size 0 ~ O(1 (1/mm) a constant for y, and the step size for x as zero in the first mg ~ (1(1) steps and the others as an appropriate constant X O(1/m²), we have the following inequality.
|rlll ||oll|
+ -
22VR())
Lggoe (J),
where & (j) = L 2 (j+1)
My 201 (1+1)%1.gu55
Theorem 3.11. Suppose Assumptions 3.1, 3.2, 3.3. 3.6 and 3.9 hold. Within each epoch, if we set the step size 9 ~ (l)mm) a constant for y, and the step size for as zero in the first my steps and the others as an appropriate constant (1/mm))) the iterates {ox} generated by Algorithm satisfy
md
[ =0 (m- mo)
K (m mg
k=O. mn)mom
where mo 2(log m) is constant and m is the subspace dimension.
In other words, we prove that the proposed LancBiO is globally convergent, and the average norm square of the hyper-gradient |/+a(())) achieves within ( outer iterations.
8

---

Published as a conference paper at ICLR 2025
4 NUMERICAL EXPERIMENTS
In this section, we conduct experiments in the deterministic setting to empirically validate the performance of the proposed algorithms. We test synthetic problem and two deep learning tasks. The selection of parameters and more details of the experiments are deferred to Appendix I. We have made the code available on https://gitubb.com/U.AA Yan Yang/LancBiO.
LancBio SubBiO
AmIGO-GD
SOBA
stocBiO
HJFBIO
AmIGO-CG
TTSA
F2SA
0.90
10
901
10-
20.86
7x
Resiual
0.8
60-10-
0.82
20-1
0.0
0.78
20
40
100
140
Time [swc]
[sec]
Time [sec]
Figure 4: Comparison of the bilevel algorithms data hyper-cleaning task when p 0.8. Left: test accuracy; Center: test loss; Right: residual norm of the linear system, || bull-
Table 1: Comparison of the bilevel algorithms on data hyper-cleaning task across two corruption rates p- = 0.5 and p = 0.8. The results are averaged over 10 runs and followed by the standard deviation. The results are conducted after 40 and 60 seconds running time.
Algorithm
p = 0.3
p 0.8
Test accuracy (%)
Test loss
Residual
Test accuracy (%)
Test loss
Residual
LaneBiO
90.35 = 0.1716
0.36 + 0.0028
1.20e-1 1 + 2.52e-5
89.45 0.24170
0.42 + 0, 0,0038
3.18e-3-
2.77e-5
SubBiO
90.21 + 0.2159
0.36 + 0.0035
2.22e-2
2.633-3
89.22
0.2387
0.42 0.0050
2.49e-2 229e-2210
1.93e-3
AmIGO-GD
90.16 1 0.2114
0.37 1. 0.0044
3.66c-2
4.00e-3
3
89.14 0.2722
0.43 0.0044
1.07e-2
8.25e-
AmGGO-CG
90.06 + 0.2305
0.38 4: 0.0053
5.89e-4 4 4 2.52e
88.57 = 0.5839
0.46 +: 0.0176
6.32e-4
4 3.74e-
SOBA
90.00 0.1811
0.37 + 0.0051
2.74e-2 2 8.52e-3
88.99 0.2661
0.42 0.0054
1.73e-2 -2 1.70e-3
TISA
89.35 - 0.2747
0.40 + 0.0103
82.91 = 0.4516
0.74 + 0,0072
StocBiO
89.20 0.1824
0.0033
86.44 0.2907
0.0064
F2SA
89.78 = 0.1969
0.40 + 0.0073
88.66 0.2828
0.51 0,0055
IUFBIO
90.21 + 0.2027
0.37 1: 0.0048
89.30 0.3394
0.43 1: 0.0040
The data hyper-cleaning task (Shaban et al., 2019) aims train a classifier in a corruption scenario, where the labels of the training data are randomly altered to incorrect classification numbers at probability results on the MNIST dataset are presented in Figure 4 and Table 1. Note that LaneBiO is crafted for approximating the Hessian inverse vector product v², while the two solid methods, TTSA and stocBiO are not. Consequently, with respect to the residual norm of the linear system, i.e., || by ll. we only compare the results with AmIGO-GD, AmIGO-CG, and SOBA. Observe that the proposed subspace- based LancBiO achieves the lowest residual norm and the best test accuracy, and subBiO is comparable to the other algorithms. Specifically, in Figure 4, the efficiency of LancBiO stems from its accurate approximation of the linear system. Additionally, while AmIGO-CG is also adept at approximating v*, the results in Table indicate that it tends to yield higher variance. Moreover, algorithms are also evaluated on the Rushion-MNIST and Kuzushiji-MNN datasets; see Figure 10 and Figure 11, respectively. The proposed LancBiO performs better than other algorithms and showcases robustness across various datasets.
Data hyper-cleaning on three datasets:
referred to as the corruption rate. The
Synthetic problem: We concentrate on synthetic bilevel optimization (1) with d, dy d and
[(x, y) = ei cos (x T Diy) + =/ D2x - y122,
d
g w) 02 1=1
sin(x, + ni) + log Σ minh
+
1 2 (Ds + G)y.

---

Published as a conference paper at ICLR 2025
LancBiO_ma m10
10P
107
LancBiO_m22
Hypergradiert narm
LancBio _ancBlo_m50 m50
2.0
LancBiO_men _m8
LancBiO_min m150
AmIGO-GD
10-3
AmIGO-GD |10
AmIGO-CG
AmIGO-CG _100
Figure 5: Influence of the subspace dimension m LancBiO. The post-fix of legend represents the subspace dimension m or the inner iteration Left: norm of the hyper-gradient; Right: residual norm of the linear system, ||Aguu bull.
It can be seen from Figure 12 that LaneBiO achieves the final accuracy the fastest, which benefits from the more accurate e* estimation. Figure 5 illustrates how variations m and influence the performance of LancBiO and AmIGO, tested across range from 10 to 150 for m, and from 2 to for 1. For clarity, we set the seed of the experiment at 4, and present typical results to encapsulate the observed trends. It is observed that the increase of accelerates the decrease in the residual norm, thus achieving better convergence of the of the classic Lanczos process. Under the same outer iterations, to attain comparable convergence property, for AmIGO-GG should be set to 10. Furthermore, given that the number of Hessian-vector products averages at + 1/m) per outer iteration for LancBiO, whereas AmIGO involves / calculations, it follows that LaneBiO is more efficient. Moreover, to illustrate how the methods scale with increasing dimensions, we present the convergence time and the final upper-level value under different problem dimensions d = 10', i = 1, 2. 3, 4 in Table 2. The results demonstrate the proposed methods maintain decent performance across different problem dimensions.
which aligns with the spirit
LancBiO
AmiGO-GD
SOBA
stocBiO
HJFBiO
SubBiO
AmIGO-CG
TTSA
F25A
0,475
10-1
6 20
10*2
0.50
Salddaiion
Sesidual
0880
10-1
0.755
0.750
Figure 6: Comparison of the bilevel algorithms the hyper-parameterrs selection task. Left: validation accuracy; Center: validation loss; Right: residual norm of the linear system, | byll
Logistic regression on 20Newsgroup: Consider the hyper-parameter selection task on the 20Newsgroups dataset (Grazzi et al., 2020). The goal is to train a linear classifier w and deter- mine the optimal regularization parameter Ç. As shown in Figure 6, AmIGO-CG exhibits slightly better performance in reducing the residual norm. Nevertheless, under the same time, LancBiO implements more outer iterations to update x, which optimizes the hyper-function more efficiently.
Generally, to solve standard linear systems, the Lanczos process is recognized for its efficiency and versatility over gradient descent methods. LiOc in a sense, reflects this principle in the context of bilevel optimization, underscoring the effectiveness of the dynamic Lanczos-aided approach.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
Bin Gao was supported by the Young Elite Scientist Sponsorship Program by CAST. Ya-xiang Yuan was supported by the National Natural Science Foundation of China (grant No. 12288201). The authors are grateful the Program and Area Chairs and Reviewers for their valuable comments and suggestions.
REFERENCES
Michael Arbel and Julien Mairal. Amortized implicit differentiation for stochastic bilevel optimization. In International Conference on Learning Representations, 2022.
Stefania Bellavia and Benedetta Morini. A globally convergent Newton-GMRES subspace method for systems of nonlinear equations. SIAM Journal on Scientific Compuring, 23(3):940-99) 2001.
Luca Bertinetto, Joao F Henriques, Philip HS Torr, and Andrea Vedaldi. Meta-learning with differen- tiable closed-form solvers. arXiv preprint 180.00810000.3666 2018.
Peter N Brown and Youcef Saad. Hybrid Krylov methods for nonlinear systems of equations. SIAM
Journal on Scientific and Statistical Computing, 11(3):450-481, 1990.
Yair Carmon and John C Duchi. Analysis of Krylov subspace solutions of regularized non-convex quadratic problems. Advances in Neural Information Processing Systems, 31, 2018.
Coralia Cartis and Lindon Roberts. Scalable subspace methods for derivative-free nonlinear least-
squares optimization. Mathematical Programming, 199(1-2):461-522 2023.
Coralia Cartis, Nicholas IM Gould, and Philippe L Toint. Adaptive cubic regularisation methods for
unconstrained optimization, part I: motivation, convergence and numerical results. Mathematical Programming, (27(2):245-295 2011.
Souradip Chakraborty, Amrit Bedi, Alec Koppel, Huazheng Wang, Dinesh Manocha, Mengdi Wang. and Furong Huang. PARL: unified framework for policy alignment in reinforcement learning. In International Conference Learning Representations, 2024.
Lesi Chen, Jing Xu, and Jingzhao Zhang. On finding small syper-gradients in bilevel optimization: hardness results and improved analysis. Annual Conference on Learning Theory, pages 947-980. PMLR, 2024.
Tianyi Chen, Yuejiao Sun, and Wotao Yin. Closing the gap: tighter analysis of alternating stochastic gradient methods for bilevel problems. In Advances in Neural Information Processing Systems, 2021.
Tianyi Chen, Yuejiao Sun, Quan Xiao, and Wotao Yin. A single-timescale method for stochastic bilevel optimization. In International Conference on Artificial Intelligence and Statistics, pages 1466-2488. PMLR, 2022.
Tarin Clanuwat, Mikel Bober- Asanobu Kitamoto, Alex Lamb, Kazuaki Yamamoto, and David Ha. Deep learning for classical Japanese literature. In Advances in Neural Information Processing Systems, 2018.
Mathieu Dagréou, Pierre Ablin, Samuel Vaiter, and Thomas Moreau. A framework for bilevel optimization that enables stochastic and global variance reduction algorithms. Advances in Neural Information Processing Systems, 35:26698 200888226600000 2022.
Stephan Dempe and Joydeep Dutta. Is bilevel programming special case of mathematical program with complemenaarity constraints? Mathematical programming, 131:37 -8, 2012.
Stephan Dempe and Alain B Zemkoho. The bilevel programming problem: reformulations, constraint qualifications and optimality conditions. Mathematical Programming, 138:447-473, 2013.
Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil. Forward and reverse gradient- based hyperparameter optimization. In International Conference on Machine Learning. pages 1165-1173 173. PMLR, 2017.
11

---

Published as a conference paper at ICLR 2025
Wenzhi Gao, Zhaonan Qu, Madeleine Udell, and Yinyu Ye. Scalable approximate optimal diagonal preconditioning. arXiv preprint arXiv:2312.15594, 2023.
Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. arXiv preprint /02.:0680..02246, 2018.
Gene H Golub and Charles F Van Loan. Malrix computations. JHU press, 2013.
Nicholas IM Gould, Stefano Lucidi, Massimo Roma, and Philippe L Toint. Solving the trust-region subproblem using the Lanczos method. SIAM Journal on Optimization, 9():5044-255 1999.
Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration com- plexity of hypergradient computation. In International Conference on Machine Learning, pages 3748 -358. PMLR, 2020.
Greenbaum, Druskin, and LA Knizhnerman. solving indefinite symmetric linear systems by means of the Lanczos method. Comput. Muth. Math. Phys, 39(3):350-356 1999.
Anne Greenbaum. Iterative methods for solving linear systems. SIAM, 1997.
Jie Hao, Xiaochuan Gong, and Mingrui Liu. Bilevel optimization under unbounded smoothness: new algorithm and convergence analysis. In International Conference on Learning Representations, 2024.
Yutong He, Jie Hu, Xinmeng Huang, Songtao Lu, Bin Wang, and Kun Yuan. Distributed bilevel optimization with communication compression. In International Conference Machine Learning, 2024.
Magnus Hestenes, Eduard Stiefel, al. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6):409-435 406449699, 1952.
Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. two-timescale stochastic algorithm framework for bilevel optimization: complexity analysis and application to actor-critic. SIAM Journal on Optimization, 23(1):147-188) 2023.
Xiaoyin Hu, Nachuan Xiao, Xin Liu, and Kim-Chua Toh. An improved unconstrained approach for bilevel optimization. SIAM Journal on Optimization, 33(4):2801-222)) 2023.
Yifan Hu, Jie Wang, Yao Xie, Andreas Krause, and Daniel Kuhn. Contextual stochastic bilevel optimization. In Advances in Neural Information Processing Systems, volume 36, 2024.
Feihu Huang. Optimal Hessian/lacobian-fres nonconvex-PL. bilevel optimization. In International Conference on Machine Learning. 2024.
Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: convergence analysis and enhanced
design. In International Conference on Machine Learning, pages 4882-4892, PMLR, 2021.
Kaiyi Ji, Mingrui Liu, Yingbin Liang, and Lei Ying. Will bilevel optimizers benefit from loops. In Advances in Neural Information Processing Systems, volume 35, pages 3011-3022 2022.
Ruichen Jiang, Parameswaran Raman, Shoham Sabach, Aryan Mokhtari, Mingyi Hong, and Volkan Cevher. Krylov cubic regularized Newton: subspace second-order method with fimension-fres convergence rate. In International Conference on Artificial Intelligence and Statisties, pages 4411-4419. PMLR, 2024.
Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. near-optimal algorithm for stochastic bilevel optimization via double-momentum. In Advances in Neural Information Processing Systems, volume 34, pages 30271-30283 2021.
Boao Kong, Shuchen Zhu, Songtao Lu, Xinmeng Huang, and Kun Yuan. Decentralized bilevel optimization over graphs: loopless algorithmic update and transient iteration complexity. arXiv preprint 20X.::24022.0366... 2024.
12

---

Published as a conference paper at ICLR 2025
Aleksey Nikolaevich Krylov. On the numerical solution of the equation by which in technical questions frequencies small oscillations of material systems are determined. Javestija AN SSSR (News of Academy of Sciences of the USSR), Otdel. mat. estest. nauk, 7(4):491-539, 1931.
Jacek Kuczyński and Henryk Wožniakowski. Estimating the largest eigenvalue by the power and Lanczos algorithms with random start. SIAM journal Matrix Analysis and Applications, 13 (4):1094-1122, 1992.
Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert D Nowak. A fully first-order method for stochastic bilevel optimization. In International Conference on Machine Learning, pages 18083-18113. PMLR, 2023.
Jeongyeol Kwon, Dohyun Kwon, Steve Wright, and Robert Nowak. On penalty methods for noncon- vex bilevel optimization and first-order stochastic approximation. In International Conference on Learning Representations, 2024.
Comelius Lanezos. An iteration method for the solution of the eigenvalue problem linear differ- ential and integral operators. Journal of Research of the National Bureau of Standards, 45(4), 1950.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(1)) (6(11))2278-324, 1998.
Junyi Li, Bin Gu, and Heng Huang. A fully single loop algorithm for bilevel optimization without Hessian inverse. In Annual AAAI Conference on Artificial Intelligence, volume 36, pages 7426- 7434, 2022.
Yongfeng Li, Haoyang Liu, Zaiwen Wen, and Ya-xiang Yuan. Low-rank matrix iteration using polynomial-filten -illtered subspace extraction. SIAM Journal on Scientific Computing, 42(3):A 1686 A1713, 2020. doi: 10. 00.1337/1M1259444
Yuwei Li, Gui-Hua Lin, Jin Zhang, and Xide Zhu. A novel approach for bilevel programs based on Wolfe duality. arXiv preprint arXiv:2302.066 06838, 2023.
Gui-Hua Lin, Mengwei Xu, and Jane J Ye. On solving simple bilevel programs with nonconvex lower level program. Mathematical Programming, 144(1-2):2773355 2014.
Bo Liu, Mao Ye, Stephen Wright, Peter Stone, and Qiang Liu. Bome! bilevel optimization made easy: simple first-order approach. In Advances in Neural Information Processing Systems, volume 35, pages 17248-17262, 2022.
Hanxiao Liu, Karen Simonyan, and Yiming Yang. Darts: differentiable architecture search. arXiv preprint arXiv: 180.00880.,,055, 2018.
Risheng Liu, Pan Mu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton. In International Conference on Machine Learning, pages 6315-6315. PMLR, 2020.
Risheng Liu, Jiaxin Gao, Jin Zhang, Deyu Meng, and Zhouchen Lin. Investigating bi-level optimiza- tion for learning and vision from unified perspective: a survey and beyond. JEEE Transactions 2021a.
Patlern Analysis and Machine Intelligence, (44(11):00045-1000)
Risheng Liu, Xuan Liu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. value-function--baed interior-point method for non-convex bi-level optimization. In International Conference Machine Learning, pages 6882-6892. PMLR, 2021b.
Risheng Liu, Xuan Liu, Shangzhi Zeng, Jin Zhang, and Yixuan Zhang. Value-function-basee sequential minimization for bi-level optimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023a.
Risheng Liu, Yaohua Liu, Wei Yao, Shangzhi Zeng, and Jin Zhang. Averaged method of multipliers for bi-level optimization without lower-level strong convexity. In International Conference on Machine Learning, volume 202, pages
PMLR, 2023b.
13

---

Published as a conference paper at ICLR 2025
Shuang Liu, Yihan Wang, and Xiao- Shan Gao. Game-theoretic unlearnable example generator. In Annual AAAI Conference on Artificial Intelligence, volume 38, pages 13349-21358, 2024.
Xin Liu, Zaiwen Wen, and Yin Zhang. Limited memory block Krylov subspace optimization for computing dominant singular value decompositions. SIAM Journal on Scientific Computing, 35(3): A A1641-A1668, 2013.
Xin Liu, Zaiwen Wen, and Ya-Xiang Yuan. Subspace methods for nonlinear optimization. CSIAM Trans. Appl. Math, 2(4):585-663 6515, 2021c.
Zhaosong Lu and Sanyou Mei. First-order penalty methods for bilevel optimization. SIAM Journal on Optimization, 34(2):1937-1999 2024.
Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient- bbased hyperparameter optimization through reversible learning. In International Conference on Machine Learning, pages 2112-2122. PMLR, 2015.
Gérard Meurant and Zdeněk Strakoš. The Lanczos and conjugate gradient algorithms in finite precision arithmetic. Acta Numerica, 15:471-542, 2006.
Yurii Nesterov et al. Lectures on convex optimization, volume 137. Springer, 2018.
Jifí V Outrata. On the numerical solution of a class of Stackelberg problems. Zeitschrift für Operations Research, 34:255-2777 1990.
Chris C Paige. Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem. Linear Algebra and Its Applications, 34:235-2553 1980.
Christopher C Paige. Error analysis of the Lanczos algorithm for tridiagonalizing symmetric matrix. IMA Journal of Applied Mathematcc, 18(3):341-34, 1976.
Christopher C Paige and Michael A Saunders. Solution of sparse indefinite systems of linear equations.
SIAM Journal on Numerical Analysis, 12(4):617-629) 1975.
Christopher Conway Paige. The computation of eigenvalues and eigenvectors of very large sparse matrices. PhD thesis, University of London, 1971.
Beresford N Parlett. The symmetric eigenvalue problem. SIAM, 1998.
Fabian Pedregosa.
Hyperparameter optimization with approximate gradient. In International Conference on Machine Learning, pages 737-746. PMLR, 2016
Zaccharie Ramzi, Florian Mannel, Shaojie Bai, Jean-Luc Starck, Philippe Ciuciu, and Thomas Moreau. SHINE: sharing the inverse estimate from the forward pass for bi-level optimization and implicit models. In International Conference on Learning Representations, 2022.
Yousef Saad. Numerical methods for large eigenvalue problems. SIAM, 2011.
Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots. Truncated back- propagation for bilevel optimization. In International Conference Artificial Intelligence and Statistics, pages 1723-1732. PMLR, 2019.
Heinrich von Stackelberg. The theory of the market economy. Oxford University Press, 1952.
Vinzenz Thoma, Barna Pásztor, Andreas Krause, Giorgia Ramponi, and Yifan Hu. Contextual bilevel reinforcement learning for incentive alignment. Advances Neural Information Processing Systems, 2024.
Jiali Wang, He Chen, Rujun Jiang, Xudong Li, and Zihao Li. Fast algorithms for Stackelberg prediction game with least squares loss. In International Conference on Muchine Learning, pages 10708-10716. PMLR. 2021.
Xiaoxing Wang, Wenxuan Guo, Jianlin Su, Xiaokang Yang, and Junchi Yan. Zarts: zero-order optimization for neural architecture search. In Advances in Neural Information Processing Systems, volume 35, pages 12868-1222 -2880, 2022.
14

---

Published as a conference paper at ICLR 2025
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mnist novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXin:1008.07747, 2017.
Quan Xiao, Songtao Lu, and Tianyi Chen. An alternating optimization method for bilevel problems under the Polyak-tojasiewicz condition. In Advances Neural Information Processing Systems, 2023.
Mengwei Xu and Jane J Ye. A smoothing augmented Lagrangian method for solving simple bilevel programs. Computational Optimization and Applications, 59:353-377, 2014.
Haikuo Yang, Luo Luo, Chris Junchi Li, Michael Jordan, and Maryam Fazel. Accelerating inexact hypergradient descent for bilevel optimization. In OPT 2023: Optimization for Machine Learning. 2023.
Yan Yang, Bin Gao, and Ya-xiang Yuan. Bilevel reinforcement learning via the development of hyper-gradient without lower-level convexity. In International Conference on Artificial Intelligence and Statistics, 2025.
Wei Yao, Chengming Yu, Shangzhi Zeng, and Jin Zhang. Constrained bi-level optimization: proximal Lagrangian value function approach and Hessian-free algorithm. International Conference on Learning Representations, 2024.
Jane Ye and Daoli Zhu. New necessary optimality conditions for bilevel programs by combining the
MPEC and value function approaches. SIAM Journal on Optimization, 20(4):1885-1905 2010.
Jane Ye and DL Zhu. Optimality conditions for bilevel programming problems. Optimization, 33 (1):9-27, 1995.
Jane J Ye, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. Difference of convex algorithms for bilevel programs with applications in hyperparameter selection. Mathematical Programming, 198 (2):1583-1616, 2023.
Ya-xiang Yuan. A review on subspace methods for nonlinear optimization. In Proceedings of the International Congress of Mathematics, pages 807-827, 2014.
Lei-Hong Zhang, Chungen Shen, and Ren-Cang Li. On the generalized Lanezos trust-region method. SIAM Journal on Optimization, 27(3):210-21122, 2017.
Yihua Zhang, Prashant Khanduri, loannis Tsaknakis, Yuguang Yao, Mingyi Hong, and Sijia Liu. introduction to bilevel optimization: foundations and applications signal processing and machine learning. IEEE Signal Processing Magazine, 41(():38-59) 2024.
Shuchen Zhu, Boao Kong, Songtao Lu, Xinmeng Huang, and Kun Yuan. SPARKLE: unified single-loop primal-dual framework for decentralized bilevel optimization. In Advances in Neural Information Processing Systemx, volume 37, pages 62912-62987. 2024.
15

---

Published as a conference paper at ICLR 2025
Appendix
A RELATED WORK IN BILEVEL OPTIMIZATION
A variety of bilevel optimization algorithms are based on reformulation (Liu et al., 2021b; Yao et al., 2024). These algorithms involve transforming the lower-level problem into a set of constraints, such as the optimal conditions of the lower-level problem (Dempe and Dutta, 2012; Li et al., 2023), or the optimal value condition (Outrata, 1990; Ye and Zhu, 1995; 2010; Dempe and Zemkoho, 2013; Lin et al., 2014; Xu and Ye, 2014). Furthermore, incorporating the constraints of the reformulated problem as the penalty function into the upper-level objective inspires a series of algorithms (Liu et al., 2022; Hu et al., 2023; Kwon al., 2023; Lu and Mei, 2024; Kwon et al., 2024). Another category of methods in bilevel optimization is the iterative differentiation (ITD) based method (Maclaurin et al., 2015; Franceschi et al., 2017; Shaban et al., 2019; Grazzi al., 2020; Liu al., 2020; Ji et al.. 2021), which takes advantage of the automatic differentiation technique. Central to this approach is the construction of a computational graph during each outer iteration, achieved by solving the lower-level problem. This setup facilitates the approximation of the hyper- -rped-gradient through backpropagation, and is noted that parts of these methods share unified structure, characterized by recursive equations (Ji et al., 2021; Li et al., 2022: Zhang et al., 2024). The approximate implicit differentiation (AID) treats the lower-level variable as a function of the upper-level variable. calculates the hyper-gradient to implement alternating gradient descent between the two levels (Ghadimi and Wang, 2018; Ji et al., 2021: Chen et al., 2022: Dagréou et al., 2022; Li et al., 2022; Hong et al., 2023). Moreover, extending the spirits of AID to distributed settings has garnered increasing interest in recent years(Kong et al., 2024: He et al., 2024; Zhu et al., 2024).
B KRYLOV SUBSPACE AND LANCZOS PROCESS
Krylov subspace (Krylov, 1931) is fundamental in numerical linear algebra (Parlett, 1998; Saad, 2011; Golub and Van Loan, 2013) and nonlinear optimization (Yuan, 2014; Liu et al., 2021c), specifically in the context of solving large linear systems and eigenvalue problems. In this section, we will briefly introduce the Krylov subspace and the Lanczos process, and recap some important properties; readers are referred to Saad (2011), Golub and Van Loan (2013) for more details.
An N dimennsiional Krylov subspace generated by a matrix A and vector b is defined as follows,
Kx Kv b) == span { {b, Ab, A 1..b),,
A
and the sequence of vectors { (, widely acknowledged for its favorable properties various aspects, including approximating eigenvalues (Kuczyński and Wožniakowski, 1992), solving the regularized nonconvex quadratic problems (Gould et al., 1999; Zhang et al., 2017: Carmon and Duchi, 2018), and reducing computation cost (Brown and Saad, 1990; Bellavia and Morini, 2001; Liu et al., 2013; Jiang et al., 2024).
N forms the basis for it. The Krylov subspace
The Lanezos process (Lanczos, 1950) is an algorithm that exploits the structure of the Krylov subspace when A is symmetric. Specifically, in the j-th step of the Lanczos process, we can efficiently maintain an orthogonal basis Q, of (A, 6). so that Tj = - AQ, is tridiagonal, which means tridiagonal matrix T3 approximates A in the Krylov subspace. Consequently, it allows to solve the minimal residual problem or the eigenvalue problem efficiently within the Krylov subspace. There several equivalent variants of the Lanezos process (Paige, 1971; 1976; Meurant and Strakok, 2006), and we follow the update rule as shown in Algorithm 3. We now present several key properties of the Krylov subspace and the Lanczos Process from Saad (2011).
Definition B.1. The minimal polynomial of a vector E RD with respect to a matrix A E R defined as the non-zero monic polynomial p of the lowest degree such that p( (A) = 0, where a monic polynomial is a non-zero univariate polynomial with the coefficient of highest degree equal to 1.
Remark B.2. The degree of the minimal polynomial p does not exceed n because the set of n + vectors (Anv, An-lu,... v. Av, v} is linearly dependent.
Remark B.3. Suppose the minimal polynomial of a vector e with respect to a matrix A is
p(x)= th Cm 13 m-1
c202 c1 Co:
16

---

Published as a conference paper at ICLR 2025
Algorithm 3 Lanczos process
Input: dimension m, matrix A E [Raxn initial vector b I
I: Initialization: Q1 [[6]| = 0, B1 = 0, Qo To = Empty Matrix
2: for j 1, 2, ..... do
3: = Agj 35
4: aj = q} U
5: Wj Uj aj
6: Bi+1 = |
7:
gj+1 Wj/ Bj+1
8:
Q = Qj 4j
T3
9:
T3
Bg Bj Oj
10: end for
Output: Tms Qm: ||6||
and has degree of m. If co and is invertible, by Definition B.1,
A nm + Cm-1 1m + ....... + 12 + 6j Av + oy 0,
multiply both sides of the equation by A and rearrange the equation,
A v (A + Cm 1 + A2 + tcv).
In other words, L belongs to the Krylov subspace Km (A)
Proposition B.4. Denote the x j matrix with column vectors 91..... 9j by Q3 and the j X tridiagonal matrix by Tj, all of which are generated by Algorithm Then the following three-term recurrence holds.
AQ; Q,T; + B5+19541e].
AQ, = Tj.
Based on Proposition B.4, the Lanczos process is illustrated Figure 7.
A
Qj
Qj
Bit1di+te]
Figure 7: Classic three-term recurrence
C DYNAMIC LANCZOS SUBROUTINE
This section lists the DLanczos subroutine (Algorithm 4) invoked in the LancBiO framework (Al- gorithm 2). One of the main differences between Algorithm 3 and Algorithm 4 is that Algorithm represents the entire m-step Lanczos process, while Algorithm 4 serves Specifically, LancBiO invokes DLanczos once each outer iteration (line in Algorithm 2), expanding both 7 and Q by one dimension. Consequently, the inputs of Algorithm 4 are not in- dexed to avoid confusion, with their corresponding variables (T Qx-1, Ag. 3x) in Algorithm evolving across outer iterations indexed by k. Another difference lies in the dynamic property of the DLanczos subroutine, i.c., the matrix A passed during each invocation varies in Algorithm while the classic Lanczos process (Algorithm 3) employs static A.
one-step subroutine.
17

---

Published as a conference paper at ICLR 2025
Algorithm 4 Dynamic Lanczos subroutine for LancBiO (DLanczos)
Input: tridiagonal matrix T, basis matrix Q with j columns, Hessian matrix A, and /
1: if j - then
2:
[91] Q, 9-1 0
3: else
4:
191 92. ................................ [/} Q
5: end if
6: új Agj Bqj-1
7: aj 9 új
8: Jj 03
9: 81+1 ||
10: 4/+1 = w/ B1+1
21+1 IQ gj+1
12: TJ+1
Output: Tj+1. 13+1
D EXTENDING LANCBIO TO NON-CONVEX LOWER-LEVER PROBLEM
The Lanczos process is known for its efficiency of constructing Krylov subspaces and is capable of solving indefinite linear systems (Greenbaum et al., 1999). In this section, we will briefly demonstrate that the dynamic Lanczos-aided bilevel optimization framework, LancBiO, can also handle lower- level problems with the indefinite Hessian.
Suppose A is invertible, and consider solving standard linear system
Ax
with initial ponit Zo, initial residual ro = Arg and initial error = A- Zg- If the matrix is positive-definite the classic Lanezos algorithm equivalent to the Conjugate Gradient (CG) algorithm (Hestenes et al., 1952), both of which minimize the A-norm of the error in affine space (Greenbaum, 1997; Meurant and Strakoš, 2006), i.e., the m-th step.
arg min
||A- - ==|a-
#E
If the matrix A is not dositive-definite MINRES (Paige and Saunders, 1975) is the algorithm recognized to minimize the 2-norm of the residual in an affine space (Greenbaum, 1997; Mcurant and Strakoš, 2006), i.e., the m-th step,
arg min ||b - AzIl.
(10)
EEzo (. (A.6)
Additionally, based on Q.. as the basis of the Krylov subspace Km (A, b). Tm as the projection of A onto (A, b), and the three-term recurrence
AQm Q..T + 9m 11+++1mm44 1m
we can rewrite (10) as
Fm = xo Qm
with
Gn = arg min
Iro AQ
arg min
Iro 1
e
= arg min 1Qm+1 (l|roll 61 - Tm+l1m))
arg min |||rol| e Tm+
18

---

Published as a conference paper at ICLR 2025
where
T'm B mm+t Cy
T 1m+1.m
T
In the spirit of MINRES, address the bilevel problem where the lower-level problem exhibits an indefinite Hessian, the framework LancBiO (Algorithm 2) requires only a minor modification. Specifically, line 13 in Algorithm 2, which solves small-size tridiagonal linear system, will be replaced by solving a low-dimensional least squares problem
C arg arggimn min ||]|=// 61 Tst1.el|
and computing the correction
Av = QkCk.
where
The Bx
T++1 ==
E PROOF OF SMOOTHNESS OF yy AND 4
To ensure completeness, we provide detailed proofs for the preliminary lemmas that characterize the smoothness of the lower level solution y" and the oyper-objective s.
Lemma E.1. Under the Assumplions 3.2 and 3.3, y (ar) 21,22 Rd,,
Luz -Lisschitz continuous, i.e., for any Ma
Lgg Ug
(x:) 3" (x2))l S
|
x2 ll
Proof. The assunption that v.g(a, y) is -ipscchit reveals || Viyg(a w)ll S Lig. Then
I (x)] = Vay (x, y) V (z, w)]
V2 2 (z,
v\
Vinn (r, v) w
Tigr Fg
since g(x, )) is My strongly convex.
Lemma E.2. Under the Assumptions 3.1, 3.2 and 3.3, the hyper-grudient Vip(a)) is Li -ipscchc continuous, iee, for any Iı, x22 E Rdd,,
[V((z1) vip(x2)l1 I Lip || x2ll,
where Lip (1 Lige) Ma
(Lfz
Loge
+ Ly Cre 14
14
Proof. By combining
(A1) b7 - by
= (AA) -, (Ai) -1 bj + (Ai) b9 (A2) 65
=(A}) (bi (0)-))) + (Ai) (A} A[) (A5)- bū
with the properties revealed by Assumptions 3.1 3.2 and 3.3, we can derive
(()))* -A5) bz
S
Liy 1
Lagn || - -
Cty
(1 (1+
Legar
|| xzll
14g
In a similar way, the subsequent decomposition holds,
V - V (xz) = (Val (a1,8" (x1) Val (22, v" (xz)))
V2yy (I1, (x1 )) (47) bi + V2 (αi)) (A3) bg 22gg (12, nn (x2)) (At) -1 69 (agg(1 (a1)) (A:) -1 65.
19

---

Published as a conference paper at ICLR 2025
It follows that
Lag
(xi)) - (=())))V((z>)))
s Lig
1
|21 -22l
1lg
Lgz
( (Lyy
CiyLauy
+ Lyz
1
+
+
|| x2
14
Mg
Cfu
Ligz
L Lgzy
1
||:0
x2
14g
Mg
Lip
201
x22
where Lup = (1+ Lggz) Ms
(Lia
Lign higy
+ Toga Cy
Hg
F PROPERTIES OF DYNAMIC SUBSPACE SECTION 3.1
In this section, focus the properties of the basis matrix Q and the tridiagonal matrix T constructed within each epoch of the dynamic Lanczos process. Denote
A
V Bn (Ix, y2) and b% v,f (æt, vk)
An epoch is constituted of complete m-step dynamic Lanczos process between two restarts, namely, after epochs, the number outer iterations is mh. Given the outer iterations = mh + j for j = 1,2 2, MT, we denote
(1+ Lyz Mg
|| Imh+ell [mmm+t +mub+sll
Est
=
for 5,t = 1,2, m and
max 1Sa,5
serving as the accumulative difference. For brevity, we omit the superscript where there is no ambiguity, and we are slightly abusing of notation that at the current epoch, {Amh+j] and (bmt+t) are simplified by {A,] and {b,} for j 1....,,, m. In addition, the approximations in the residual system (8) are simplified by 0 and 6 = 61 A
We rewrite the dynamic update rule from Section 2.2
U A,g 839j-11
(1I)
aj qj Ujj
T
(12)
Wj Uj ajli+
(13)
8- = ||رب|||
(14)
f- wj/ wj/Bj+1 /Bj+1,
(15)
for 1,2.... with qo 0 8 0 and Q 91 b/ [6]. The following proposition characterizes that the dynamic subspace constructed in Algorithm 2 within an epoch indeed an approximate Krylov subspace.
Proposition F.1. At the j-th step within an epoch (j = =,,, 1,2,.... m 1), the subspace spanned by the matrix Qj+1 in Algorithm 2 satisfies
or A A ... A Vs 1, ...
}
span(Qj+1) S span
(16)
Specifically, when A A2 -.. A; A and Qi+1 is of full rank.
span 2;+1) = Kj+1 (A, b)
b satisfies (16). We will give a proof by induction. Proof. Note that Q1 sp) with q |ba
Suppose for = 1,2, ......,... j, it holds that
a, = 0 or 1
span(Q:+1) S span
A A ... A
vs 1, 1,2,....
20

---

Published as a conference paper at ICLR 2025
By the dynamic Lanezos process, it yields
45+2 =
(Ajt10jt1- - -B++10 a3+105+1)
(17)
Bj+2
Since
or 0 = =1
}
4j+1 span A/ 6
V/s 1, 1.2..... }
then we have
Aj+195+1 span
Aq A91 A2 4
02 = or O2
= = 1,2, ..... 1
It follows from (17) that
9j+2 E span A A7 ... ..A.1'
0, or as
}
= 1,2, j+1 }
By induction, complete the proof.
Although estimate the difference between the basis of the above two subspaces, it is noted that the Krylov subspaces can be very sensitive to small perturbation (Meurant and Strakoš, 2006; Paige, 1976; 1980; Greenbaum, 1997). The next lemma interprets the perturbation analysis for the dynamic Lanczos process in terms of Aj, which satisfies an approximate three-term recurrence with a perturbation term 8Q.
Lemma F.2. Suppose Assumptions 3.1 to 3.3 hold. The dynamic Lanczos process in Algorithum with normalized q and aj, B33 9j satisfies
A7Q; Qj Bs+1 6Q}
(18)
for j 1.2,...,
, m, where Q3 = [q1, 92. ......... Gil. 8Q, = 15q1,
01
B2
B2
022
Bs
B3
B3
B3
Okj
The columns of the perturbation 8Q, satisfy
||óqull S LoweEj for 1,2,.....
Additionally, we decompose Qj as
0, RR + Rí
(19)
with Ry as a strictly upper triangular matrix, then
R, B 17 90+1e} + &R,, where 8Ry is stricily upper triangular with elements Karl S E. s
(20)
Proof. From
a, = u, 93 T..gAAAp A,g; B,9 9--1
and
1
1
1/+19/
3541 w/ 91
(u, - -3 9i By+1
(A,4 - a,9j B, Aj.
3,41
we can derive
4/+19 =
(21)
21

---

Published as a conference paper at ICLR 2025
by induction. Then, we combine equations (11), (12), (13) and (15), and rewrite them in the perturbed form:
84+19+11 A.q. Bıqi-1 a,q, Ajq. B.qi-1 0gg; + ōq., for i = 1.2,..... ..... where ||óqil| I E due to Assumpstions 3.2 and 3.3. Specifically, (22) can be rewritten in compact form:
(22)
= + Then, we consider the orthogonality of matrix Qj. which is reflected by R2 in (19). Multiply on both sides of (18) by 05,
&02--
Q, BJ+1 gt11] + +Q,80,
=
Combining its symmetry with the decomposition (19), we obtain
T, (RJ + R₂) - [R] + R, - B,-1 (Q) ++1e] - ( + 03 -8Q Qy. (23)
Denote M, T,R, R,T, which is upper triangular. Since the consecutive q, is orthogonal as revealed by (21), we conclude that the diagonal elements of Mj are O. Furthermore, by extracting the upper triangular part of the right hand side of (23), can get
M₁ R2 R2T; BJ+ Q3 1itte + +RRj.
where R2 is strictly upper triangular with elements Çat satisfying: for = =.......... .......
{
(1-1,1 -- pq: 59.119 q
Sat u. 6g: 64 yt:
1,2, .......:4. 2
From the boundedness of [[ōg; |l. it follows that for 1 < 15al S 2..
L
Lemma F22 illustrates the influence of the dynamics in Algorithm 2 imposed on the three-term Lanezos recurrence, and (19) reveals, R reflects the loss of orthogonality of the basis Q. However, the following lemmas demonstrate that the range of eigenvalues of the approximate projection matrix T's indeed controllable.
To proceed, we establish the Ritz pairs of Ti (u/".u)) (1) for ..., such that (1) (5) T2Y y(u) diag ( 1....... ;...,,33)) -1) )
where the normalized (v. 1= form the orthogonal matrix y with the elements . s, 8 S j, and we arrange the Ritz values a specific order,
Set for 1 <
We define the j-th approximate eigenvector matrix
z = 2 .2
S
-/ Q2r().
and the corresponding Rayleigh quotients of A
v00) = (2())
2/2
(1) for = 1, =1,2,.... ....
(
The subsequent lemma describes the difference of eigenvalues between T, and some T,, constructed in preceding steps.
Lemma F.3. Suppose Assumptions 3.1 to 3.3 hold. For any eigenpair (v.)).1))) of there exists an integer pair (s, n)) where 1 < s S n < j. such that
(4(2) | s
2 LygvEz
(24)
v3
(v!5)
R R,y
22

---

Published as a conference paper at ICLR 2025
(1)
Proof. Multiply the extended eigenvector yr of T; by T3, where i <j
y/ (1)
(i) (
T
0
-
B B 10)
Mr y
(25)
Then multiply (u!") on the both sides of (25),
(n!))) 149)) (v!2)
3r
Be (0 ())
(26)
0
S141,87
and multiply the eigenvectors (v! T
and u on the left and right of equation (20), respectively,
(112) ( (u!)) R5yl (6) Sjx -/ B5+1 (z11) 9i+1 + + ())
(27)
where we define
2t = (1!0)) T
&Rzyt
(28)
( Sst < Lj because of Lemma F.2. Specifically, taking s = t in (27).
Note that
9j+1
3 B5+15
5)
(29)
can rewrite (29) in the matrix form
Q 27 9r+1 -r(rer
(30)
where for s,
c Gr ==
Ca8
By observing that O, gr+1 = (()) , is the (r + 1)-th column Ry, we can derive
(r)
(1) Ryyl)=
Σ
(
Ctt
yr
(31)
B
(r)
0
t=1
(r)
(32)
14
where (31) and (32) follow from (30) and (26) respectively. Consequently, the definition (28) reveals
i
2 = |8 1F S 2 Lowe
Σ
(33)
8,1=1
Based on (32) and the orthogonality of y(i),
61
<
u
(
1d 11)
min sddcl<j
i-1
Σ
(
(r)
min 1sdsk<
follows from the Cauchy- Schwarz inequality and (33) that
vJv E Di=EE
Ei c(tt
2
m0/ 5
min 15dsk<j
(v.0))
v3
(v)"))
23

---

Published as a conference paper at ICLR 2025
Lemma F.4. Suppose Assumptions 3.1 to 3.3 hold. Q2 and T the basis matrix and the approxi- mate tridiagonal matrix in the j-th step and R, is the strictly upper triangular matrix defined in (19) characterizing the orthogonality of Qg. Given ui the i-th eigenvalue of TT. then fhe i 1.
2/3 G + 1)² =g S 144) S L +
2/3 (j 1)3 L gg
Mg
(34)
Proof. By conducting the left multiplication with (1 2,) and the right multiplication with y/" ( on both sides of equation (18), we obtain the following equation,
A
Mi 2
(3)
(1)
21
-c + (2(11) 8Q2y1 (
Dividing by
(2())
T
2;
())
we have
AJ2
(5)
14
(1)
(
2
(
5 18
+
(>()) T
8Q ()
(
(
+ 2 (v!0)
Ryy.
6
+ L gon 1 (11+2
(35)
1 2 ) Ryy,
where the inequalities come from
(
1|
2
1 2 (1(
Rgv ())
(36)
Case I: If
(1150) T
Rgy
(5)
3 Ej
(37)
8 2
holds, then from (36),
Furthermore, (35) reveals that there exists a Rayleigh quotient v6 of A3 that satisfies
L00) - 14 < 2
+ V 22llgme S (4j 22 T -.
Case II: If the condition (37) does not hold, by applying Lemma F.3, we find an integer pair (s1,01) with S $1 S ny < such that
(( (11)))
2j² Lg
S L
v3
(2(1))
Rgy,") ()
By observing that ( (n1)
can also be categorized into Case I or Case II, we can repeat this process and construct sequence ((), n with S 7441 < m < B1 < no j until
y
(m1)
and
yo
(m1)
(m1) Re, yoy
(m141)
Em4s
<
8
2
24

---

Published as a conference paper at ICLR 2025
Inequality (37) holds when the superscript j = 1 since Ti
[an]; V1 (1) 1 and 21 q1- Therefore,
we can obtain {(st, nz) ) in finite steps, resulting in the following estimate.
b - -4 $ p(u)
( Ps41
+
Ma441
(me 1) - 14
(m4)
$ 4jj + 22j Lgoves
2v 2/ 3ng E
1=0
2 $
(j 1)3 LggnEs
for some T0).
Any Rayleigh quotient of A bounded by its eigenvalues (Parlett, 1998), i.e., for any vl1,
(D) S v()))
A mozt
where X and X the minimum and maximal eigenvalue of A3, respectively. Based Assumption 3.2 and Assumption 3.3, we complete the proof.
G PROOF OF LEMMA 3.10
We detail the proof for Lemma 3.10 in this section, starting with proof sketch.
G.1 PROOF SKETCH
The proof of Lemma 3.10 structured by four steps.
Stepl: extending E ē in the lemmas given in Appendix F. In Appendix F, we adopt A and bj as reference values with each epoch for the analysis, i.e., for
j = 1,2, ........... m,
A3 = V (Iz. v5) b7 = Vyl (1.. b = b, A1 To = b, rg = b - A3 Ang.
and
Lgz)
Ej = max 1
1 +
||Emitts Imh+tell + || Fmn++oll
Mg
In parallel, we can also view the A1, d2 as reference values. In this way, denote the similar quantities
Ay = 2 g (z), 85), by Vuf (z), B' := A1 Fo = b', ry = A. (38)
and
:
|amitt Rmh++ell + 11 ymhh+ell 1
max
(39)
Consequently, we extend the lemmas in Appendix F. The results listed in Appendix G.2 the recipe of Step2.
upper-bounding the residual In Appendix G.3, we then demonstrate that if the value of ě, is not too large, | can be bounded, which is an important lemma for Step3.
Step2:
||
Step3: controlling ěÿ and Ej by induction.
Since the expression of ē, does not involve y/", its magnitude can be controlled by adjusting the step size, implied by (39). With the help of the stability of Ej, can prove
2/3 (j+ 1)3 Lyuv > 0.
Mg
(40)
3
25

---

Published as a conference paper at ICLR 2025
Step4: proof of Lemma 3.10 Based on the conclusion (40) revealing the benign property of the dynamic process, we achieve the proof of Lemma 3.10.
G.2 EXTENDED LEMMAS FROM APPENDIX F
In this part, we view the A1, b reference values. In this way, we can institute A₃ with A₁ and Ej with ε in the lemmas from Appendix F, resulting in the extended version lemmas.
Lemma G.1. (Extended version of Lemma F.2) Suppose Assumptions 3. 3.3 hold. The dynamic Lanczos process in Algorithm 2 with normalized q1 satisfies
AiQ; - QjT 15+10++1e +0'
(41)
for j = 1,2, ........ m, where Q2 = 991, 92.
9zl. 60' = [ō0h o8h.........8
α1
B2
B2
0r2
T3 =
B3
Bg
The columns of the perturbation 8Q/ satisfy
| S L for i = 1,2,.....
(42)
Lemma G.2. (Extended version of Lemma F:4) Suppose Assumptions 3. to 3.3 hold. Q2 and T3 are the basis matrix and the approximate tridiagonal matrix in the j-th step. Take pl (3) eigenvalue of T then for i = I.....J.
as the i-th
2/3 (j (1+1)² Loug S (1) S Loy + Mg
2/3 (i+ 1) Lowe
(43)
3
G.3 PROOF OF STEP2
The following lemma demonstrates that the value of &, is not too large, Lemma G.3. Suppose Assumpsions 3.1 to 3.3 and
||| can bounded.
2/3 Ui + 1)3 Layy > 0 Bg
are satisfied within epoch. Then, it holds that
R )) +
S 2
ViL (i),
where K (j) == Lay 2 (j+1) L
By 243 11+ L.yay
Proof. Denote the solution in the dynamic subspace in the j-th step by
A = (T,))' || (Tz) e1.
(44)
where Tj is nonsingular because of Lemma G.2. By (38), (41), and (44), we have
7/ = 13 DE, AQ,
It follows that
F
< |80
- (t,)-tei t
(45)
26

---

Published as a conference paper at ICLR 2025
The first term the right side of (45) can be bounded (42) and (43):
80 (T3)
Vilguz
T.gy
S ViL (j).
(46)
Mg - (j + 1)² T. Tygyy
Recall that
T TJ+1J 1341e]
(47)
and that for any symmetric tridiagonal matrix T, where the upper left (j +1 x block is Tz+1.5- the
application of the classic Lanczos algorithm to T, starting with the initial vector e1 will result in the
matrix Tz+1, at the j-th step. To construct a suitable (j + 1) (j +11) symmetric tridiagonal matrix
T, consider virtual step with X 0 which leads to (4]+1) 15+1) (Ij, Uj). 8J+1 = Ej. and
Tz+1J
=
33+1
1āj+1
By Lemma G.2, given any eigenvalue μ of
2v3
(j+ + 1) LyuvEy Sus Lov + 2/3 (i+ (i+1)² Lowe
Mg 3
In this way, [ (T,)-² ej can be seen as the residual in the j-th step of the classic Lanezos process with the positive-defintte matrix T and the iniitial vector e1. Since the cigenvalues of T satisfy (47), follows from the standard convergence property of the Lanczos process (Greenbaum, 1997) that
(
[3412e (T,) e sais2ve $2V((()
which completes the proof.
G.4 PROOF OF My 2/3 ()+1)² L > 0
In this part, we will give the detailed proof of Pg p-- demonstrate that, with appropriate step sizes, the auxiliary variable 242 bounded. thereby showing that the hyper-gradient estimator (5) remains bounded. This ensures the stable behavior of the dynamic Lanezos process,
2/3 (j + 1)3 LyovE > 0. At the same time, we
Lemma G.4. Suppose Assumptions 3.1, 3.2, 3.3, 3.6 and 3.9 hold. If within each epoch, we set the step size 0 ~ (11/ /m) a constant for y. and the step size for x. as zero in the first ma ~ ((1) steps and the others as an appropriate constant X ((1/m²) then for any epoch,
2/3 (++ 1)3 Lyg > for j 1,2, ..... 1,
3
and there exists a constant C > so that ||we|| S C. for un generated by Algorithm
Proof. Consider the iterates within one epoch and the constants
Lgv + Mg E
0 << غ < M4g and =
It follows from Lemma G.3 that if Eg S
Jae
then
1)m+11um Laxu"
1b1 Aruyl
Vk
V3E
< 2
+
x < 3VR.
(48)
bb
VR +
2(m + 1)2
27

---

Published as a conference paper at ICLR 2025
Then, we give a proof by induction. At the beginning of the algorithm, Assumption 3.9 reveals
3Bpg
1 = 11 will S
(m + 1) Lovy
Combining it with ||6 - Holl < C, constructs the start of induction within an epoch and induction between epochs. Within an epoch, suppose the following statements hold for :.2....:J. J.
||b Alusll 3V/CC.
1 (3VkC,+ C
| S M9
|pp S C5z + 1 (3ViC, + Cia),
V3E
6, 2(m + 1)3T
/31lg
8, S 4 4(m + 1) Lowa
Then by setting the stepsizes
S
Mg my
(49)
3
Lggy
)((rr=
1 (3VKC, +c1)))
A A5 4(m + 1)
1 0
Crez
Ug
uz
and noticing that
V (0 Ui)1l
|| 1/ (I,+1, W1) - V Vuo(2442,(++1)
S Lou [/s yi+1
Lgg (llyr - yē l| + +1107 wi+all)
Lg 11 wēl 11 + LogLgg ||:s se41 11,
can get
EJt1 AA
Vyg (x)41.31) p +0
1=1
1
LogLoz
Vips
+ gg
11:: vfll
Mig
Lag Lez
1 +
+ 0L
E1
Pg
1=1
< 1 + 0 Lgy gz
V3
+ +0Lgu
Pg
4 (m + 1)² Loun
V3E
< 2(m + Lyvn
(50)
It follows from (48) that
|b, A
|1 S <3V/Cr,
(51)
1 (WKC,+C())
||ح+رر|| ||A (b1 - A11541- -bill S Mgg
(52)
p
SCzz 1
(3VRC, + +Crm).
(53)
28

---

Published as a conference paper at ICLR 2025
Additionally, the descent property of lly= v/ll and the Lipschitz continuity of y reveal that l v.ll _ (1 8)gg) & 1l5-1 - vill
(54)
$ (1 81gg) ) ||ue11 - y || + (1 - (1-0p
Ligg
||0
Æz-1l
||
Mg
$ (1 814)) + l1 vill + (1 - 01ty)
Laz
t=
x
Mg
1=1
|1 -vill + (1 024)8
Lgz
P
(55)
Mi
t=1
Setting
V3gg 8(m + 1 199
Csz + 1 (3VKC,+
(1
Lyr bg
Lgz
+
1
+ (1 - 0s1g)*
guu
Mg
yields from Assumption 3.9, (52), (55) that
(1 Lg
Ej+1 S
max
| vell +
max
11
Vzl
Pg
150,55+1
15e55+1
i=j
Legn
Lugz
S 11 vill +
+ (1 - -(1gg)
Hg
Mgg
Lig
+ (1 0gg)) )
(Lyz))
< 1y vil +
1
Hig
pgg
/3pg
2 4 (m 1)a Lovy
(56)
1 (3VRC, + CFz). have
As for the next epoch, denoting C. No
[|mm+t Am+ || - (Tefa + Co) Em- + 11b1 - -Aluml
V R2 -
(Lyo + Co) 22
+ VmLggyymtak Cr
VTc +1
(
(is + Logen
m
Cr
Cr Cr ) 5m41 + 2
emmggyymm 18:
VK+1
Cr:
(57)
by choosing m, E such that
Lis Loonn (Li++
Vk
m
) Em+i
2V/
e m+ i 1.
C Cr
VR+11
Moreover, since the step size for is set as zero at the first my steps in the next epoch, we obtain for i 1,2,..., mo.
v3μg 4 (m 1)3² Lggy
Emti = Ilymte Vmntall (1 0p,)) llym
Em vin
(58)
Specifically,
v3pg
1mmmmo 4m+mol || < (1 Optd) mg 2 11 wiall S 8 (m + 1)* Logyy
(59)
we choose mi so that (1 - 0pg) < J. Therefore, by induction within an epoch (50), (51), (52), (53), (56) and induction between epochs (57), (58), (59), we conclude the lemma.
29

---

Published as a conference paper at ICLR 2025
G.5 PROOF OF LEMMA 3.10
Lemma G.5. Suppose Assumptions 3.1, 3.2, 3.3, 3.6 and 3.9 hold. If within each epoch. we set the step size 0 ~ O(1/m) a constant for y and the step size for as zero the first me ((1) steps, and the others as an appropriate constant X ~ (11 fm²²), we have the following inequality,
Fلlll Foll
VR()= -10--1
22VR())
+ Wilo99933 (i)
(R()++1
Lout 2 L
where & (f)
My 2 ( Ey
Proof.
Lemma G.4 guarantees the condition
2/3 satisfied. The remaining proof can be directly adapted from Lemma G.3.
My (j 1) L Lg > 0
H PROOF OF THE MAIN THEOREM
In this section, provide proof of the main theorem presented in Section 3.2. Let A V (Tk: yis) and by = Vyf yk). and let the reference values be A = V ( 67 of (Ix. vt) and 2 = (At))16G. To begin with, short proof sketch is provided for guidance, which structured in four main steps.
Stepl: upper-bounding the residual I|0% - r
Appendix F and G lay foundation to bound the residual term | ufll (as corollary of Lemma G.5).
studying the descent property of || - v Lemma H.l reveals the descent property of the estimation error for y" follows,
Step2:
2 + (1+ 2) (1 Opo) (1-0pt)
[Lgur)
2 || xxll?
1 vi+1l (11
llue (1-)))))
Mg
Step3: Controlling the hyper-gradient estimation error [[vx(x))) vip(xx)l Defining
6x = fall + L2 C2
11
vEll
1|0 - o
L 9
and incorporating the results from the last two steps, then in Lemma H.2 we can establish the upper bound for [[()))))) - V<((x>)] and dy recursively.
Step4: Assembling the estimations above and achieving the conclusion Consider the descent property of La == S + 5g. Substituting the inequalities developed in stepl to step3, and telescoping the index from to K gives the convergence results.
The following lemma displays the descent property of the iterates {yk).
Lemma H.1. Suppose Assumptions 3.2 and 3.3 hold. Setting 0 < 0 S Hot
2
we have
|| vitill $ 0pg) 1lue- 1/ (1 [ 1+ ) (1- 0pp)
Legs fgg
| =xal??
for any 0.
30

---

Published as a conference paper at ICLR 2025
Proof. Algorithm 2 executes single-step gradient descent the strongly convex function (( during the outer iteration. Leveraging the established convergence properties of strongly convex functions (Nesterov et al., 2018), are thus able derive the following.
||++1 - y 11 - (1 0pr) 11 - 11+
By Young's s inequality that |a b/2 S (1 + o |aI + (1 + = 11612 with any > 0,
|uk+1 vetill = - 1lve= vEll++ [+ + )) (1 8(1g)
2 Log
|
xxll?
Mg
In the context of bilevel optimization, define the initial residual in the (h + 1)-th epoch as and the residual in k-th step
Th+1 == bmi+1 Ama+1 Ub.
Tk = (bmi+1 - Ama+10n) - A Auk.
Based on the boundness of U2 in Lemma G.4 we can estimate
Mg 110k - al - (bg A57a) 1gAnm
(bb Aāoa) AãAnz Tk rkll
< |6 bmn+ || + ||A - Amh+11 | + [|re/l
< fu (3 vā) - (Imh+1, 1, Mmmh+1))l
+ Tgy [1 uk) - (Imh +1 Vmnh [| |[rell
= (Ieru + T.gvv |1 11) s0 + ||rall
I (Lyv + Co) e/ + ||r=|l,
(60)
which comes from (A7)) < 02 Un Aok. and [u S C..
Lemma H.2. Suppose Assumptions 3.1, 3.2, 3.3, 3.6 and 3.9 hold. Within each epoch, we set the step size 0 ~ O(1/m) a constant for y and the step size for as zero in the first mo steps, and the O(1/m²) //mm)) then the iterates
others as an appropriate constant
{x+} for k = mh + h = 0,,, and j = mo + 1, mo + .....,
generated by Algorithum satisfy
(61)
and
8% < mmmm-m)) (0m(h-1)) + 2m 8mih-1) + 200 2008 8m(h-2)
+ 12m2 X2imp (/ (/V:p((mmmt)l" + V 1| + 1 (
Σ
( (62)
t=mp
t=mo
t=mo
where 0 <<< mo 2ilog m) constants, is the subspace dimension.
T2 T22 2 1.g Cig
6z =
11yk vill + || 2i
12201
31

---

Published as a conference paper at ICLR 2025
Proof. By (60), conclusion from Lemma G.3 and the Young's inequality, it holds that
Liy + Lggyy Co
2
110e
11 2
(e))
||xl1
Mg
2
L L + Lggy
(et)²
Mg
1
Vi( -
Vi guue (i) ())
+ 2
2VR()
||
(f( +
Liy + LOye Cv.
2 (e7)²
2
1gg
2/
1
Vi(i) -1 (R(+++
+ (V) )2 (E))) ||
+ 4
46 (j)
+1
An estimate can be made for ||r4+I||.
||a+1
[1mmt+1 Amh+18pll
[lbmb+ Uynh + A A mn Amh+1
S Ly [Imma Amn++)) (Embs Fmnt
Ciu Lgvy [2 Vrnl+1) (Eemttr Vmaa + Lgg 11itan Umuh
My
vmall Mg MoLsy by + Csyy Lggu (1
Lye Mg
Tgg |m
| Ima
Mg
and thus
24
22g
Vi( - [R(5++1 +1
11 uill < We (
+ (()
umall
R (j) +
( + CyyLooy
(1+
Lix
))
V K (ī) -1
||Imh+1 Imh 12
mg
Mi
(J( +1
(63)
where
Liw Lon 2
We
Mg
(Vmi (m)) Ligue??
02
(Hglisy + Crylggn
Lox (1+
Lexe
)) 02 O 0mm).
+
112
1g
By definition of e/, the update rule of x, and Young's inequality, it holds that
2 2j-1
(et)² 2 < 2m22
[+
Lig Mg
Σ
E (am)++))
|| + 2llumbhtte) - vin++a())
32

---

Published as a conference paper at ICLR 2025
for some my +1 S i(f) S j. recursively to derive
Then apply the descent property of 11 yill (Lemma H.1)
2 S 2m12 (ch)²
(1+
Lg
2 25-1
Σ
V (#ma+s)
||
+ + +0)()) (() (1 - -1p,) ((() 11 - Vmnnl 11
Pg
f())
7.gx
+ 2
1
(1 (pg)
Σ
(1 + +o (1 (1-)) 1| 1k-
Hg
S 22
1+
L
2
m +
++ σ
1
(1 - 0pg)
1-1
- (zma+t)))
+ l1ymm - Vin 11
Rig
$ 2X
1
Legn
m
_ (1+= (1 - 0pg)
1
4 (Imatt))
1gg
+ 2 (1 + a) (1 011g)
mh
llymnh - Yonk 11
(64)
Since Lemma G.4 reveals that under the appropriate step-size setting O (11 and 0 (()))
R(j) for R ~ Q Eex) ), also choose > 0 such that
Aay
1 V (1)
::= max
<1 (41)
(65)
VR +1
Additionally, set the warm-up steps mg to satisfy
2mo
max
{
fully
gay
12 2 + 1zz 19 +
L
L 5y R
Ug Lg + CiyLgve
1 Lgg
))
2
2wdpp
Pg
Mg
(9/122 u2 + 11.g
which means mo Ω (log m). In this manner, adding of (63) and incorporating the estimation for d (64) yield that
11 vE on both sides
d6 < (3()mmo) (6mh) + 11yma vmnill + 6m2 X20
T (ommate)) + +|
1=mB
< m20) mo) (6mi) ,. (1mm(b-))) Vin(h-11ll
+ 6m X
( + ( (
m- ( + (Vp(Em((-)))l?)
(66)
1=mo
where we apply Lemma H.l recursively to obtain the second inequality, and adopt the notation
We
Lgz
RL
[g MgLLy + +CyyLyu
(1+
(1+
Lgz
))
Wp == max
1 +
mm2ui
Mg
Mig
1 = (1 0p4))
T.7
1: C2
(Lgz 2
+ +m
m
1.2 12
Mg
33

---

Published as a conference paper at ICLR 2025
Consequently, expanding 6mb in (66) in the same way derives (62). Regrading the upper- bound of the hyper- hyper-gradient gradient estimation error, we have
( (xx)
3 zf (k)) Va (x Bill +3 V (I mwll
1|ux
+ + m / (It, Uk) Vang (at, v/lll ||
Ciu My
<3 Ja lluk vtll +T 110 rI 31
llue
vkll
= 33 2 5u. gz
Theorem H.3. Suppose Assumptions 3.1, 3.2, 3.3, 3.6 and 3.9 hold. Within each epoch, if we set the step size 0 ~ 0l1/m) a constant for y, and the step size for x as zero in the first mg steps and the others appropriate constant À ((1 (1/m²)) the ilerates {ax) generated by Algorithm satisfy
Σ 10
K
md-1 (m -mm)) ma))
K (m mo
mo)
(A mor vs mm
where mg (log m) is constant and m is the subspace dimension.
Proof. Consider the Lyapunov function C3 == ((xx) + 6R. According to Lemma E.2, a gradient
descent step leads to the decrease in the hyper-function:
4p (0x++1) (xxx) & (*+) Ix+1 - xx) + Le
| - xx.1?.
(67)
2
Then, telescoping Cx+1 - Ls over the index set == f 0 S k s K, (lmodmm) > mo}.
Σ keZ
Cx+1 L4
K+1
Σ p((x+1) - p(xx) +
Σ
8x
8x
KET
((N 1) mod m)) Somo
(Vp(zz), xk) Lip
h-1
+66 (om)
2x +
2mi
2( - mo 8me
Σ 6F
KET
c=0
KET.
(10)
Σ
(Tx + 72m)ww)) |[
2
KEI
(
+ >22 r2m2ww)) [Vp((x)) Fp((z))) + Σ
2me 2(m mo) 8mn Σ 8k
+
KEI
(112
Σ
(TTp + T2m))u) [[(p((>)))
KEI
(2) (2 + wp)
+ 3T2 gz
Σ KEI
5x
+
2m 2(m mo)
8mo
Σ 8k KEI
(68)
e=0
where (i) follows from inequalities (62) and (67); (ti) results from the update rule dh++1 ( and Young's inequality; (tii comes from (61). Taking the coefficients of 8, into account, the dimension parameters (m, mp) satisfying 2m-(maa)) < 1/4 and the step size À such that
1
1
A < min
(en)))))))) L2 6L. (12 72m
1/2 4 (L + 72m )
(69)
34

---

Published as a conference paper at ICLR 2025
In this way, obtain the following result from (68),
(2) -2-x0 X2 (Iip 72m2wp) V1 (xx))1
LK+1 Lo
Cx+1
Lx
mod m))mm
(Ae nod mi) 2-mg
Rearrange the above inequality and denote
4² minzeRd :0(x),
[ & 4 (10 (2o) - 40")
400
x +
(Rooad m)) 2-mo
which completes the proof by dividing both sides by m- K.
m
I DETAILS ON EXPERIMENTS
1.1 GENERAL SETTINGS
We conduct experiments to empirically validate the performance of the proposed algorithms. We test synthetic problem, a hyper-parameters selection task, and data hyper-cleaning task. We compare the proposed SubBiO and LaneBiO with the existing algorithms stocBiO (Ji et al., 2021), AmIGO-GD and AmIGO-CG (Arbel and Mairal, 2022), SOBA (Dagréou et al., 2022) and TTSA (Hong et al., 2023), F2SA (K (wonn& et al., 2023) and HJFBIO (Huang, 2024). produced
bilevel optimization:
The experiments
workstation that consists of two Intel@ Xeon@ Gold 6330 CPUs (total 2228 cores), 512GB RAM, and one NVIDIA A800 (80GB memory) GPU. The synthetic problem and the deep learning experiments carried out on the CPUs and the GPU, respectively. For wider accessibility and application, we have made the code available on http://gthhu..oo/UCCC.a Yan Yang/LancBiO.
For the proposed LancBiO, we initiate the subspace dimension at 1, and gradually increase it to m 10 for the deep learning experiments and to m 80 for the synthetic problem. For all the compared algorithms, we employ grid search strategy to optimize the parameters. The optimal parameters yield the lowest loss. The experiment results are averaged over 10 runs. Note that Assumption 3.9 for initialization is not used in practice, for which treat it as theoretical assumption rather than incorporating into Algorithm 2 in this paper.
In this paper, we consider the algorithms SubBiO and LancBiO in the deterministic scenario, initially compare them against the baseline algorithms with a full batch (i.e., deterministic gradient). In this setting, LancBiO yields favorable numerical results. Moreover, in the data hyper-cleaning task, to facilitate a more effective comparison with algorithms designed for stochastic applications, we implement all compared methods with small batch size, finding that the proposed methods show competitive performance.
1.2 DATA HYPER-CLEANING CLEANNNG
The data hyper-cleaning task (Shaban et al., 2019), conducted on the MNIST dataset (LeCun et al., 1998), aims to train classifier in corruption scenario, where the labels of the training data randomly altered to incorrect classification numbers at a certain probability referred to as the corruption rate. The task is formulated follows,
min A
Leat w²) = Prull EiennEED I(w v(, Vi)
8. w" = arg min Cep (w,
Σ
o()))(ww,, vi) - Cr|lw||3
Dul
(EL
where L(.) is the cross-entropy loss, a(-) is the sigmoid function, and C, is regularization parameter. In addition, serves linear classifier and (A:) be viewed as the confidence of each data.
35

---

Published as a conference paper at ICLR 2025
In the deterministic setting, where implement all compared methods with full-batch, the training set, the validation set and the test set contain 5000, 5000 and 10000 samples, respec- tively. For algorithms that incorporate inner iterations to approximate y² or v², select the inner iteration number from the set {5i = 1, 2,2,3, 3,4.. The step size of inner iteration is selected from the set {0.01,0.1, 1, and the step size of outer iteration is chosen from {5 10' i -3 2, - 1,0,,, Regarding the Hessian//acobian-Fe 10 to implement finite difference methods. The results presented in Figure 4. Note that LancBiO is crafted for approximating the Hessian inverse vector product v", while the solid methods stocBiO, TTSA, F2SA, and HJFBiO not. Consequently, concerning the residual norm of the linear system, i... || v2 bull, we only compare the results with AmIGO-GD, AmIGO-CG and SOBA. Observe that the proposed subspace-based LancBiO achieves the lowest residual norm and the best test accuracy, and SubBiO is comparable to the other algorithms. Specif- ically, Figure 4, the efficiency of LaneBiO stems from its accurate approximation of the linear system. Furthermore, we implement the solvers designed for the stochastic setting using mini-batch enable broader comparison Figure 8. is shown that the stochastic algorithm SOBA tends to converge faster initially, but algorithms employing full-batch approach achieve higher accuracy.
algorithm HJFBiO,
the step size
18169
0.30
LancBiO full SubBiO full AmIGO-GD 256 AmIGO-CG full SOBA 128
0.88
0.06
Tustioss
Nett
0.00
TTSA 256
0.82
4x10
stoc 256
[sec]
Time [soc]
Figure Comparison of the bilevel algorithms on data hyper-cleaning task with mini-batch when p 0.5. The training set, the validation set and the test set contain 5000, 5000 and 10000 samples, respectively. The post-fix of legend represents the batch size. Left: test accuracy; Right:test loss.
To explore the potential for extending our proposed methods to a stochastic setting, we also con- duct experiment with stochastic gradients. In this setting, where implement all compared methods with mini-batch, the training set, the validation set and the test set contain 20000, 5000 and 10000 samples, respectively. For algorithms that incorporate inner iterations approximate y" or . we select the inner iteration number from the set (3i iteration selected from the set {0.01, 0.1, 10}, the step size of outer iteration is chosen from -3, 2, 0,, 1, and the batch size is picked from {32 2
1, 2, 4}. The step size of inner
0,3, AmIGO-CG is not presented since it fails in this experiment in our setting. The results in Figure demonstrate that LancBiO maintains reasonable performance with stochastic gradients, exhibiting fast convergence rate, although the final convergence accuracy is slightly lower.
{1 10'
3x10
0.775
LancBiO 256
30.800
2x10
SubBiO 256
0.255
Text toxe
AmIGO-GD 256
0.800
SOBA 128
0.75
TTSA 256
0.50
0.725
Sx1111
stocBiO_
30
20
Figure 9: Comparison of the bilevel algorithms data hyper-cleaning task with mini-batch when p 0.5. The training set, the validation set and the test set contain 20000, 5000 and 10000 samples, respectively. The fost-fix of legend represents the batch size. Left: test accuracy: Right:test loss.
36

---

Published as a conference paper at ICLR 2025
Additionally, also evaluate the performance of bilevel algorithms Fashion-MNIST (Xiao et al., 2017) and Kuzushiji-MNIST (Clanuwat et al., 2018) datasets, both of which present more complexity compared MNIST. Specifically, Fashion- -NNIST serves as modern replacement for MNIST, featuring grayscale images of clothing items across 10 categories, and Kuzushiji-MNIS culturally rich dataset of handwritten Japanese characters. The results, reported in Figures 10 and 11, reveal that LancBiO performs better than other algorithms and showcases robustness across various datasets, and SubBiO delivers comparable convergence property.
LancBiO
AmIGO-GD AmlGO-CG
SOBA
stocBiO
HJFBIO
SubBiO
TTSA
F2SA
0815
2.70
10"
20805
8000
0.795
2.50
0.700
13-4
0.780
30
Time [sec]]
[sec]]
[sec]]
Figure 10: Data hyper-cleaning task tested on the Fashion- MNIST dataset when p = 0.5. Left: test accuracy; Center: test loss; Right: residual norm of the linear system, || ball-
LancBiO SubBiO
AmIGO-GD AmlGO-CG
SOBA
stocBiO
HJFBIO
TTSA
F2SA
2.50
201
0.66
3,45
2.40
10-
3.35
3.30
0.63
3.25
0.62
0.63
2:15
0.60
10
30
30
45
50
1.0
20
30
40
50
15
35
20
30 1ime [Se]]
25
35
45
45
50
Time [snc]]
Tirm [sec]]
Figure 11: Data byper-cleaning task tested on the Kuzushiji- MNNIST dataset when p = 0.6. Left: test accuracy; Center: test loss; Right: residual norm of the linear system, || bull.
1.3 SYNTHETIC PROBLEMS
We concentrate synthetic scenario in bilevel optimization:
min f(z,y") == cos (.VDiv") + 1/4 1Dar 11
t. =²= arg min g( v) y
(70)
sin(x; + y1) + log
Σ 2ap
+ アリ (Do+G)
where we incorporate the trigonometric and log-sum-exp functions to enhance the complexity of the objective functions. In addition, we utilize the positive-definit matrix to ensure a strongly convex lower-level problem, and diagonal matrices D; (i 1, 3) to control the condition numbers of both levels.
In this experiment, set the problem dimension d 10 and the constants C1 0.1, 0 0.5. G is constructed randomly with eigenvalues from to 10², We generate entries of D,, D2 and D from uniform distributions over the intervals [-5, 5. [0.1, 1.1] and [0, 0.5). respectively. Taking into
37

---

Published as a conference paper at ICLR 2025
10
101
10
pypergraeeent nom
101
LancBio
Resioul nom
101
SubBiO
AmIGO-GD
20-
10h
AmlGO-CG
10-4
SOBA
TTSA
10
stocBiO
40
50
50
100
120
140
10-3
20
43
80
100
120
140
160
Time
Figure 12: Comparison of the bilevel algorithms on the synthetic problem. Left: norm of the hyper-gradient; Right: residual norm of the linear system, ||Agu byll-
account the condition numbers dominated by D; (i 1,2 and G, we choose À 1 and 0 10 for all algorithms compared after manual search.
It can be seen from Figure 12 that LancBiO achieves the final accuracy the fastest, which benefits from the more accurate v" estimation. Figure illustrates how variations in m and I influence the performance of LancBiO and AmIGO, tested across range from 10 150 for and from 2 to 10 for I. For clarity, we set the seed of the experiment at 4, and present typical results to encapsulate the observed trends. It is observed that the increase of accelerates the decrease in the residual norm, thus achieving better convergence of the hyper-gradient, which aligns with the spirit of the classic Lanczos process.
When m convergence, which is demonstrated in Figure that for 50, further increases in m merely enhance the convergence of the residual norm. Under the same outer iterations, to attain comparable convergence property, for AmIGO-CC should be to 10. Furthermore, given that the number of Hessian- vector products averages at (1 + 1/m) outer iteration for LancBiO, whereas AmIGO requires 2, follows that LancBiO is more efficient.
50, the estimation of v" is sufficiently accurate to facilitate effective hyper-gradient
Table 2: Comparison the synthetic problem (70). The dimension of problems is denoted by d. Results are averaged over 10 runs.
d 10
d 100
d 1000
d = 10000
Algorithm
UI. Val.
Time (S)
UL Val.
Time (S)
UL Val.
Time (S)
UL Val.
Time (S)
LancBiO
4.52e-2
0.32
6.37e-2
0.53
5.29e
1.30
1.21e-2
16.36
SubBiO
3.73e-2
1.00
7.19e-2
1.17
4.63e-2
1.72
-2.91e22
21.26
AmIGO-GD
1.67e-1
2.44
1.05e-1
3.48
1.05e-1
1.46
4.96e-2
46.53
AmIGO-CG
5.56e-2
0.40
7.65e-2
1.90
5.73e-2
2.44
3.68e-2
25.00
SOBA
1.70e-1
0.60
1.28e-1
1.64
1.03e-1
2.52
3.49e-2
33.89
TTSA
5.66e-2
0.47
5.52e-2
0.89
6.52
2.81
1.87e-1
121.07
stocBiO
6.24e-2 2
0.29
6.02e
0.45
5.21e-2
1.34
1.30e-2
20.66
Moreover, to illustrate how the proposed methods scale with increasing dimensions, we present the convergence time and the final upper-level (UL.) value under different problem dimensions d 10', i 1, 2, 4 in Table 2. The results show the robustness of the proposed methods across varying problem dimensions.
In addition, as discussed in Appendix D, to address the bilevel problem where the lower-level problem exhibits an indefinite Hessian, the framework LancBiO (Algorithm 2) requires a minor modification. Specifically, line 13 in Algorithm 2, which solves small-size tridiagonal linear system, will be replaced by solving a www-dimensional least squares problem. We test the modified method LancBiO- MINRES the following bilevel example borrowed from Liu al. (2023a) with non-convex
38

---

Published as a conference paper at ICLR 2025
lower-level problem,
min fx,,y" == || all + 11y" ccl #ER
c|2
(71)
8. t.
V/ E arg min sin(x + - ci), for = 1, ,,
WER
where the subscript denotes the i-th component of a vector, while E R and E R' are parameters. The results, reported in Figure 13, imply the potential of extending work to non-convex scenarios.
10
LaNNEE
Optimality metric
10a
10A
108
10
1000
Outer iteration
Test LancBiO-MINNEES on the synthetic problem (71) with 100. The metric follows from the necessary conditions developed for bilevel problems without lower-level convexity in Theorem of Xiao et al. (2023).
Figure 13:
1.4 LOGISTIC REGRESSION ON 20NEWSGROUP
Consider the hyper-parameters selection task on the 20Newsgroups dataset (Grazzi et al., 2020), which contains 20 topics with around 18000 newsgroups posts represented in feature space of dimension = 130107. The goal is to simultaneously train a linear classifier w and determine the optimal regularization parameter Ç. The task formulated as follows,
min Lval(S, w*) = [Deul E(AnyEDD L(w vi)
8. t. w" arg min (er(w,))
1 d
Σ
L(wa,, v:)
2 F
Drrl
(E..)))DD
where I..) the cross-entropy loss and {c}] are the non-negative regularizers.
The experiment is implemented in the deterministic setting, where we implement all compared methods with full-batch, the training set, the validation set and the set contain 5657, 5657 and 7532 samples, respectively. For algorithms that incorporate inner iterations to approximate y" or v*, we select the inner iteration number from the set (5i i 1,2,3,4) 4). To guarantee optimality condition of the lower-level problem, we adopt decay strategy for the outer iteration step size, i.c., A = 0/10.4, for all algorithms. The constant step size 0 inner iteration is selected from the set (0.01, 0.1, 10) and the initial step size À of outer iteration is chosen from (5 x 10 i== 3, 22, 1,,0, The results are presented in Figure 6. In this setting, AmIGO- CG exhibits slightly better performance in reducing the residual norm. Nevertheless, under the same time, LancBiO implements more outer iterations update which optimizes the hyper-function more efficiently.
39


Reviews:
[{"review_id": "XaqZbm43BX", "reviewer": "Reviewer_sZwj", "summary": "This paper investigates a gradient based method for solving bilevel optimization problems with a strongly convex lower-level objective. To enhance the efficiency of hypergradient computation, the paper employs the Krylov subspace method alongside the Lanczos process to accelerate the solution of the linear systems involved. The paper presents non-asymptotic convergence results for the proposed method. Numerical experiments on data hyper-cleaning, synthetic problems, and logistic regression demonstrate its effectiveness.", "questions": "Is the proposed algorithm applicable in a stochastic setting?", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. This work leverages the Krylov subspace method and the Lanczos process to solve the linear system involved in hypergradient computation, providing a new approach for efficiently approximating the hypergradient.\n\n2. A non-asymptotic convergence rate of $O(\\epsilon^{-1})$ is established.", "weaknesses": "1. The description of Algorithm 2 is somewhat unclear, particularly in relation to the theoretical convergence result in Theorem 3.12. Theorem 3.12 requires the step size for $x$ to be zero during the initial few steps in each epoch, implying that before updating $x$, the algorithm must obtain a $y_k$ sufficiently close to $y^*(x_k)$, as also reflected in the theoretical analysis. This indicates that the epoch $h$ functions primarily as an outer loop index, where at each iteration $h$, the algorithm first solves the lower level problem to obtain a $y_k$ that is sufficiently close to the lower level solution before proceeding to update $x$. However, in other gradient based algorithms, such as SOBA, this requirement for $y_k$ to be close to the lower-level solution is not imposed.\n\n2. In the convergence analysis, the proof of Theorem 3.12 , the step size for $x$ must be very small to ensure stability of the Krylov subspace method and the dynamic Lanczos process in LancBiO as $x_k$ is updated. It appears that the inclusion of the dynamic Lanczos process constrains the step size for $x$, potentially making it smaller than in algorithms without this process, such as SOBA. However, the paper lacks a discussion of this issue.\n\n3. The parameter $m$, which defines the dimension of the Krylov subspace in LancBiO, should play a crucial role in the efficiency of the Krylov subspace method and hence the performance of LancBiO. However, the convergence analysis does not show how the choice of $m$ impacts the performance of LancBiO.\n\n4. The convergence analysis requires an unusual assumption, Assumption 3.10, which is not present in other works, such as SOBA. Why is Assumption 3.10 necessary? Is it because $y_k$ generated by LancBiO cannot be guaranteed to converge to $y^*_k$ ?", "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": "qBdc55UGpv", "reviewer": "Reviewer_H14N", "summary": "The paper addresses the complexities of bilevel optimization, a framework widely used in machine learning applications. It introduces a novel approach that utilizes the Lanczos process to construct low-dimensional Krylov subspaces, aiming to alleviate the computational challenges associated with hypergradient calculations. By avoiding the direct computation of the Hessian inverse, the proposed method demonstrates improved efficiency and achieves a convergence rate of $ O(\\epsilon^{-1}) $. The authors present a theoretical foundation for their approach and validate it through experiments on synthetic problems and deep learning tasks.", "questions": "1. How does the proposed method perform in scenarios where the lower-level problem is not strongly convex? \n2. Can the authors elaborate on how the method scales with increasing dimensions in the bilevel optimization problems?\n3. Are there any limitations observed when applying the proposed method to non-standard bilevel optimization problems, such as those with noisy or sparse data?\n4. Could the authors provide more detailed insights into the computational complexity of solving the small-size tridiagonal linear system mentioned?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- **Innovative Approach**: The incorporation of subspace techniques into bilevel optimization is novel and contributes significantly to the field, potentially opening new avenues for research.\n- **Theoretical Rigor**: The paper offers a solid theoretical framework that is well-justified, providing confidence in the proposed method's validity and effectiveness.\n- **Empirical Validation**: The experimental results show promising performance improvements over existing methods, suggesting practical applicability in real-world scenarios.\n- **Clarity of Presentation**: The paper is well-organized, making complex concepts accessible to readers, which enhances its impact.", "weaknesses": "- **Limited Scope of Experiments**: While the experimental results are promising, they are conducted on a limited set of problems. Broader validation across diverse benchmarks would strengthen the paper's claims.\n- **Assumptions in Theory**: The theoretical results rely on certain assumptions that may not hold in all contexts, potentially limiting the generalizability of the findings.\n- **Lack of Comparison with More Methods**: The paper could benefit from comparisons with additional state-of-the-art bilevel optimization methods to contextualize the improvements more clearly.", "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": "2wrfwbmTzk", "reviewer": "Reviewer_nGhD", "summary": "This work presents an approach for bilevel optimization using Krylov subspace methods and the Lanczos process to approximate inverse-Hessian vector products. The method constructs low-dimensional Krylov subspaces and solves tridiagonal linear systems, achieving convergence to an $\\epsilon$-stationary point with $\\mathcal{O}(\\epsilon^{-1})$ complexity. Experimental evaluations are conducted to illustrate its performance.", "questions": "Could the authors provide further details on extending LancBiO to stochastic scenarios? Additionally, could they elaborate on how the Lanczos process ensures convergence in stochastic settings?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper introduces a novel modification to the Lanczos process by re-linearising the objective functions within the iteration, specifically for bilevel optimization. This reduces large-scale subproblems to smaller tridiagonal linear systems, and the explanation provides a potential framework for applying subspace techniques in this context.", "weaknesses": "I identify two weaknesses:\n\n1.To my knowledge, many recent advances in bilevel optimization have focused on addressing problems where the function is not necessarily strongly convex, or even non-convex. The strong convexity assumption in this paper may be overly stringent, potentially limiting the method's applicability to a broader range of optimization problems.\n\n2.The proposed method is technically constrained to deterministic settings. While LancBiO can be extended to stochastic scenarios, its performance in these settings has been inconsistent. It remains uncertain whether LancBiO can be effectively adapted for stochastic environments, which is crucial for many practical applications. For example, despite including experiments in stochastic settings, SOBA seems to perform better, indicating challenges in extending LancBIO effectively.", "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": "CrkSiJ5rxZ", "reviewer": "Reviewer_CLSW", "summary": "This work proposes two novel algorithms, SubBiO and LancBiO, for solving bilevel optimization problems in which the lower-level problem is strongly convex. The main contribution lies in incorporating the Krylov subspace and the Lanczos process into bilevel optimization, achieving a more accurate hypergradient estimate by efficiently and dynamically solving the associated linear system. Under certain conditions, the authors establish non-asymptotic convergence for LancBiO and conduct an empirical study to validate the efficiency of both SubBiO and LancBiO.", "questions": "Apart from the questions raised in the Weaknesses section, some additional questions are as follows:\n\n$\\textbf{Q1:}$ In the implementation of SubBiO, how is the two-dimensional subproblem in line 4 solved? A bit more discussion on these choices would be helpful. \n\n$\\textbf{Q2:}$ In the convergence analysis of LancBiO, the hyperparameter $m_0$ plays an important role, but it is not included in the experiments. Could the authors clarify why? \n\n$\\textbf{Q3:}$ Can the authors provide more detail on why the Lanczos process in Algorithm 2 does not affect the final convergence guarantee? A brief proof sketch of Theorem 3.12 would be helpful.\n\nSuggestions for improvement that did not affect the score:\n\n$\\textbf{About $m_0$ in Lemma 3.11 and Theorem 3.12:}$ First, $m_0$ in Lemma 3.11 should satisfy $m_0 = \\Omega(1)$, meaning $m_0$ must be greater than some positive constant, as implied at the end of the proof of Lemma G.4. Second, $m_0$ in Theorem 3.12 should be set to $m_0 = \\Omega(\\log m)$, following from lines 1643–1647 and equation (63) in the proof of Lemma H.2.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "$\\textbf{S1:}$ The introduction of subspace techniques into bilevel optimization is well-motivated. The explanations of how the proposed SubBiO and LancBiO algorithms relate to existing methods are helpful. \n\n$\\textbf{S2:}$ A convergence rate guarantee is provided for LancBiO under certain conditions. \n\n$\\textbf{S3:}$ To empirically justify the improvement achieved by accurately solving the associated linear system, the experiments report the residual norm of the linear system.", "weaknesses": "For the theory, the main concerns are as follows:\n\n$\\textbf{W1: Assumption 3.10}.$ Based on (10), Assumption 3.10 depends on the step size, $\\lambda$. However, in Lemma 3.11 and Theorem 3.12, $\\lambda$ is specified to satisfy $\\lambda\\sim \\mathcal{O}(\\frac{1}{m^4})$. Under this setting, it seems that Assumption 3.10 may not hold when $m$ is sufficiently large. Therefore, in general, how can Assumption 3.10 be checked—either theoretically or empirically? A bit more discussion on these issues would be helpful. \n\n$\\textbf{W2: Theoretical results lack discussion}.$ In reading Theorem 3.12, it is unclear how the results compare to other bilevel algorithms, as some existing algorithms (e.g., deterministic SOBA) also reach an \\epsilon-stationary point within $\\mathcal{O}(\\epsilon^{-1})$ outer iterations.  \nAdditionally, there is no complexity analysis. For example, what is the number of oracle calls required to reach an $\\epsilon$-stationary point? Typically, bilevel optimization literature provides the number of gradient, Hessian-vector, and Jacobian-vector products required to reach a stationary point with precision $\\epsilon$. \n\n\n$\\textbf{W3: Lack of convergence analysis for SubBiO}.$ Given that SubBiO has a simpler structure than LancBiO, it would be beneficial to include a theoretical analysis or, if not feasible, to discuss the reasons for its absence.\n\nFor the experiments, the main concern is as follows:\n\n$\\textbf{W4: Experiments could be expanded}.$ The experiments could be more comprehensive. For example, it would be beneficial to include additional competing bilevel methods, especially Hessian-free algorithms like F2SA (Kwon et al., ICML 2023). Additionally, using more datasets in the data hyper-cleaning and logistic regression tasks would help validate the efficiency of both SubBiO and LancBiO.", "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)