Newton-Based Methods for Bilevel Optimization Beyond Lower-Level Strong Convexity

Newton-Based Methods for Bilevel Optimization Beyond Lower-Level Strong Convexity Youran Dong∗ Junfeng Yang† Wei Yao‡ Jin Zhang §

arXiv:2610.26839v1 [math.OC] 19 Oct 2026

Abstract Bilevel optimization provides a framework for machine learning tasks such as hyperparameter optimization, data reweighting, and meta-learning. Existing Newton-based bilevel methods typically rely on the classical hypergradient formula, which is generally unavailable when lower-level solutions are nonunique and the lower-level Hessian is singular. We develop two Newton-based methods beyond lower-level strong convexity by exploiting curvature shared by the lower-level and multiplier updates. Newton-based GALET handles potentially singular lower-level curvature under a positive spectral-gap condition, exploiting its range structure to avoid the squared-Hessian multiplier system used in GALET. Newton-based BAMM uses a strongly convex aggregation model when the lower-level objective is convex and the upper-level objective is strongly convex in the lower-level variable. We establish convergence guarantees and improved complexity bounds relative to GALET and sl-BAMM under the corresponding assumptions. We further introduce GNBO, a generalized Newton-based framework that captures the shared-curvature structure of the proposed methods and several existing Newton and quasi-Newton bilevel methods. Experiments on synthetic problems, data hyper-cleaning, and GPT-2 data reweighting evaluate the proposed methods across different bilevel learning tasks.

Keywords: bilevel optimization, Newton methods, lower-level nonuniqueness, KKT conditions, oracle complexity

1

Introduction

Bilevel optimization models hierarchical learning and decision-making problems in which one optimization problem depends on the solution of another. It arises in hyperparameter optimization [Pedregosa, 2016, Franceschi et al., 2018], data cleaning and reweighting [Ren et al., 2018, Pan et al., 2025], meta-learning [Franceschi et al., 2018, Rajeswaran et al., 2019, Ji et al., 2020], and adversarial learning [Bishop et al., 2020, Wang et al., 2021, Zhang et al., 2022]. We consider the optimistic bilevel problem min

x∈Rm , y∈Y ∗ (x)

f (x, y),

Y ∗ (x) := arg minn g(x, y), y∈R

(1)

where f and g are the upper-level (UL) and lower-level (LL) objectives, respectively, and the LL problem may admit multiple solutions. One of the main challenges in solving (1) stems from its nested optimization structure. When g(x, ·) is strongly convex, the LL problem admits a unique solution y ∗ (x), and the bilevel problem reduces to minimizing the hyper-objective Φ(x) := f (x, y ∗ (x)). Under standard smoothness assumptions, the implicit function theorem yields the hypergradient


−1
∇Φ(x) = ∇1 f x, y ∗ (x) − ∇212 g x, y ∗ (x) ∇222 g x, y ∗ (x) ∇2 f x, y ∗ (x) . (2) ∗ School of Mathematics, Nanjing University, Nanjing, P. R. China. (yrdong@smail.nju.edu.cn). † School of Mathematics, Nanjing University, Nanjing, P. R. China. (jfyang@nju.edu.cn).

‡ School of Artificial Intelligence, Southern University of Science and Technology, Shenzhen, P. R. China. (yaow@sustech.edu.cn). § Department of Mathematics and National Center for Applied Mathematics Shenzhen, Southern University of Science and Technology, Shenzhen, P. R. China. (zhangj9@sustech.edu.cn).

1

A large body of bilevel optimization methods approximates the hypergradient in (2) and applies a gradient-type update to the UL variable. A standard approach is approximate implicit differentiation (AID) [Pedregosa, 2016, Ghadimi and Wang, 2018, Grazzi et al., 2020, Ji et al., 2022, Chen et al., 2021, Hong et al., 2023, Dagréou et al., 2022, Chu et al., 2024, Ramzi et al., 2022, Fang et al., 2025, Dong et al., 2025]. AID involves two key approximation tasks: estimating the LL solution y ∗ (x) and 
−1
approximating the Hessian-inverse vector product ∇222 g x, y ∗ (x) ∇2 f x, y ∗ (x) . Importantly, the Hessian ∇222 g in the second task is the curvature of the same LL objective whose minimizer is estimated in the first, revealing a shared-curvature structure between the two tasks. Motivated by this observation, several recent Newton-based bilevel methods exploit LL curvature to couple these two approximation tasks in different ways. SHINE [Ramzi et al., 2022] applies a quasi-Newton method to the LL problem and reuses the resulting quasi-Newton approximation to estimate the Hessian inverse in the hypergradient. qNBO [Fang et al., 2025] instead uses separate quasi-Newton procedures for the LL solution and the Hessian-inverse vector product, with nonasymptotic guarantees based on properties of BFGS updates. More recently, NBO [Dong et al., 2025] uses two inexact Newton directions involving the same LL Hessian to approximate the LL solution and the Hessian-inverse vector product, leading to improved nonasymptotic complexity bounds over existing AID methods. Together, these methods highlight the benefits of exploiting LL curvature in both LL optimization and hypergradient approximation. The Newton-based methods discussed above rely on LL strong convexity and the resulting hypergradient formula (2). When the LL problem admits multiple solutions, however, the solution mapping becomes set-valued and the implicit differentiation underlying (2) is generally unavailable. An alternative is to characterize LL optimality through first-order stationarity and consider min

x∈Rm , y∈Rn

f (x, y)

s.t.

∇2 g(x, y) = 0.

(3)

When g(x, ·) is convex or satisfies the Polyak–Lojasiewicz (PL) condition, its stationary points are global minimizers, making (3) equivalent to the original bilevel problem (1). This reformulation underlies several methods for bilevel optimization beyond LL strong convexity [Mehra and Hamm, 2021, Shen et al., 2025, Xiao et al., 2023, Liu et al., 2023, Hallak and Suissa, 2026]. The constrained viewpoint has led to several methods beyond LL strong convexity. Xiao et al. [2023] derive KKT-type necessary conditions for bilevel problems with PL LL objectives and develop GALET, whose convergence analysis assumes a uniform positive gap in the nonzero singular values of the LL Hessian. Liu et al. [2023] introduce an iteration-dependent aggregation model and develop the single-loop method sl-BAMM to approximate the resulting KKT-type residuals. Fully firstorder alternatives include BOME [Liu et al., 2022], which uses a value-function-based constrained reformulation and a modified dynamic barrier gradient method, and HJFBiO [Huang, 2024], which combines finite-difference estimators with projection-based updates for Hessian/Jacobian-free nonconvex-PL bilevel optimization under an additional nonsingularity condition on the LL Hessian at LL solutions.

1.1

Motivation and Contributions

The preceding discussion motivates extending Newton-type methods beyond LL strong convexity. When LL solutions are nonunique, the LL Hessian may be singular, making the classical hypergradient formula generally unavailable. This raises a natural question: Can Newton-type methods effectively address bilevel optimization with possibly nonunique lower-level solutions? Our main contributions are summarized as follows. (i) We develop Newton-based GALET for bilevel problems with possibly nonunique LL solutions under a positive spectral-gap condition, as also assumed in GALET [Xiao et al., 2023]. Starting from the KKT residuals of (3), we establish a range-inclusion property of the LL Hessian and exploit the curvature shared by the LL and multiplier residuals to construct Newton-type updates. This structure allows the multiplier update to use the LL Hessian directly, rather than the squared-Hessian system employed by GALET. 2

(ii) We develop Newton-based BAMM for problems with a convex LL objective and an UL objective strongly convex in the LL variable. Building on the aggregation model of sl-BAMM [Liu et al., 2023], we construct Newton-type LL and multiplier updates using the same positive-definite Hessian. (iii) We establish convergence and complexity guarantees for both methods under assumptions comparable to those of GALET [Xiao et al., 2023] and sl-BAMM [Liu et al., 2023], respectively. The resulting bounds improve the corresponding complexity results, with gradient/Jacobianvector-product (JVP) and Hessian-vector-product costs analyzed separately. We further introduce GNBO, a generalized Newton-based framework that captures the shared-curvature structure of the proposed methods and several existing Newton and quasi-Newton approaches [Ramzi et al., 2022, Fang et al., 2025, Dong et al., 2025]. (iv) We evaluate the proposed methods on two synthetic examples, data hyper-cleaning with linear and neural models, and GPT-2 data reweighting, examining their empirical performance across different bilevel learning tasks.

1.2

Notation

For a differentiable function g, we denote its partial gradients with respect to x and y by ∇1 g(x, y) and ∇2 g(x, y), respectively. For twice differentiable g, ∇212 g(x, y) ∈ Rm×n denotes the Jacobian of ∇1 g(x, y) with respect to y, and ∇222 g(x, y) ∈ Rn×n denotes the Hessian with respect to y. Unless otherwise specified, ∥·∥ denotes the Euclidean norm for vectors and the Frobenius norm for matrices, while ∥ · ∥op denotes the spectral norm.

1.3

Organization

The remainder of the paper is organized as follows. Section 2 develops Newton-based GALET and Newton-based BAMM. Sections 3 and 4 present the convergence analysis and numerical experiments, respectively. Section 5 introduces GNBO as a unifying framework, and Section 6 concludes the paper. Detailed proofs and experimental details are provided in the appendices.

2

Newton-Based Bilevel Optimization Beyond Strong Convexity

2.1

From Hypergradients to KKT Residuals

We first revisit the shared-curvature structure in the strongly convex LL setting. When g(x, ·) is strongly convex, the LL solution y ∗ (x) is unique, and the hypergradient error induced by an inexact LL solution y is O(∥y − y ∗ (x)∥) under standard smoothness assumptions [Ghadimi and Wang, 2018]. Locally, a Newton step y + = y − [∇222 g(x, y)]−1 ∇2 g(x, y) reduces the LL solution error quadratically: ∥y + − y ∗ (x)∥ = O(∥y − y ∗ (x)∥2 ). The Hessian-inverse vector product in (2) involves the same LL Hessian. Exploiting this shared-curvature structure, NBO [Dong et al., 2025] uses the following two-step procedure: (i) Newton step. Approximately solve ∇222 g(x, y)(v, u) = (∇2 g(x, y), ∇2 f (x, y)). (ii) Hypergradient step. Set dx := ∇1 f (x, y − v) − ∇212 g(x, y − v)u and update the UL variable. Here, −v is an inexact Newton direction for the LL variable, while u approximates the Hessian-inverse vector product. Both involve the same LL Hessian. When the LL solution may be nonunique, the classical hypergradient formula (2) is generally unavailable. We instead consider (3), with Lagrangian L(x, y, u) = f (x, y) − u⊤ ∇2 g(x, y) and KKT residuals ry (x, y) = ∇2 g(x, y), ru (x, y, u) = ∇2 f (x, y) − ∇222 g(x, y)u, rx (x, y, u) = ∇1 f (x, y) − ∇212 g(x, y)u. 3

(4) (5) (6)

Under the PL condition and suitable regularity assumptions, ry = ru = rx = 0 is a necessary optimality condition [Xiao et al., 2023]. Under strong convexity, these conditions recover the LL solution, the Hessian-inverse vector product, and hypergradient stationarity, respectively [Liu et al., 2023]. The shared-curvature structure extends naturally to these residuals. For fixed x, Taylor expansion gives ry (x, y − v) = ry (x, y) − ∇222 g(x, y)v + O(∥v∥2 ). (7) Hence, if ∇222 g(x, y)v = ∇2 g(x, y) and ∥v∥ = O(∥∇2 g(x, y)∥), then ∥∇2 g(x, y−v)∥ = O(∥∇2 g(x, y)∥2 ). For the multiplier residual, ru (x, y, u) = 0 corresponds to ∇222 g(x, y)u = ∇2 f (x, y). Thus, the LL and multiplier residuals give rise to two systems involving the same LL Hessian. When the LL solution is nonunique, this Hessian may be singular, so the two systems need not have unique solutions. Rather than relying on a Hessian inverse, we therefore work directly with their residuals: (i) Residual step. Given (x, y), construct v and u by reducing r1 := ∇2 g(x, y) − ∇222 g(x, y)v and r2 := ∇2 f (x, y) − ∇222 g(x, y)u. (ii) UL step. Construct dx from rx (x, y − v, u) and use it to update x. Here, r1 is the linearized residual associated with ry (x, y−v), so −v is a Newton-type direction, while r2 = ru (x, y, u) is the multiplier residual. This residual viewpoint preserves the shared-curvature structure without requiring an invertible LL Hessian and allows the two systems to be treated according to the problem structure. We next instantiate this principle in two settings, leading to Newton-based GALET and Newton-based BAMM.

2.2

Newton-Based GALET

We first consider the positive spectral-gap setting. Assumption 1. There exists µ > 0 such that every nonzero singular value of ∇222 g(x, y) is uniformly bounded below by µ, i.e.,
+ inf σmin ∇222 g(x, y) ≥ µ, x,y

+ where σmin denotes the smallest positive singular value.

Assumption 1, also used in GALET [Xiao et al., 2023], allows the LL Hessian to be singular and is therefore compatible with nonunique LL solutions. For example, g(x, y) = 12 ∥A(x − y) − b∥2 satisfies this assumption for a fixed matrix A ∈ Rn×n with Rank (A) = r < n. It also implies ∥[∇222 g(x, y)]† ∥op ≤ 1/µ. The following lemma characterizes the resulting LL structure. Its proof is given in Appendix B. Lemma 1. Suppose Assumption 1 holds, g ∗ (x) := miny g(x, y) > −∞, and ∇222 g(x, ·) is Lg,2 Lipschitz. Then g(x, ·) is convex, ∇2 g(x, y) ∈ Range (∇222 g(x, y)), and ∥∇2 g(x, y)∥ ≥ µ dist(y, Y ∗ (x)). If, in addition, ∇2 g(·, y) is Lg,1 -Lipschitz, then Range (∇221 g(x, y)) ⊆ Range (∇222 g(x, y)). Remark 1. By convexity and the error bound in Lemma 1, g(x, y) − g ∗ (x) ≤ ∥∇2 g(x, y)∥ dist(y, Y ∗ (x)) ≤ µ−1 ∥∇2 g(x, y)∥2 . Thus, Assumption 1, together with the conditions of Lemma 1, implies the PL condition with parameter µ/2, which is assumed separately in GALET [Xiao et al., 2023].

4

Algorithm 1 GD(H, q; γ, T ) for Hp = q, with H(p) = Hp 1: Initialize p0 = 0. 2: for t = 0, . . . , T − 1 do 3: pt+1 = pt − γ(H(pt ) − q). 4: end for 5: Return pT .

Algorithm Design

To reduce the LL residual ry , we approximately solve ∇222 g(xk , y k )v = ∇2 g(xk , y k )

(8)

for v k and set y k+1 = y k − v k . Although ∇222 g(xk , y k ) may be singular, Lemma 1 ensures ∇2 g(xk , y k ) ∈ Range (∇222 g(xk , y k )), so (8) is consistent. We approximately solve the linear systems using Algorithm 1, where H denotes a linear operator. Throughout, Hessian-vector products are computed via automatic differentiation [Pearlmutter, 1994, Dagréou et al., 2024], without explicitly forming or storing the Hessian matrices. For the multiplier update, let qk := ∇2 f (xk , y k ), so r2 = qk − ∇222 g(xk , y k )u. Since qk need not lie in Range (∇222 g(xk , y k )), r2 = 0 may be infeasible. Let Pk := ∇222 g(xk , y k )[∇222 g(xk , y k )]† be the orthogonal projector onto Range (∇222 g(xk , y k )). Then r2 = (Pk qk − ∇222 g(xk , y k )u) + (I − Pk )qk , where the two components are orthogonal and the second is independent of u. We therefore reduce the controllable component through the consistent system ∇222 g(xk , y k )u = Pk qk .

(9)

The projection need not be computed explicitly. We justify this for both the multiplier residual and the UL direction. First, applying Algorithm 1 directly to (∇222 g(xk , y k ), qk ) gives pt+1 = pt − ηk (∇222 g(xk , y k )pt − qk ). Since ∇222 g(xk , y k ) is symmetric, Pk ∇222 g(xk , y k ) = ∇222 g(xk , y k )Pk = ∇222 g(xk , y k ), and hence Pk pt+1 = Pk pt − ηk (∇222 g(xk , y k )Pk pt − Pk qk ). Thus, {Pk pt } is exactly the GD sequence for the consistent system (9). Second, Lemma 1 implies ∇212 g(xk , y k )(I − Pk ) = 0. Consequently, ∇212 g(xk , y k )Pk pt = ∇212 g(xk , y k )pt ,

∇222 g(xk , y k )Pk pt = ∇222 g(xk , y k )pt .

Hence, the component of pt in ker ∇222 g(xk , y k ) affects neither the UL direction nor the Hessianweighted multiplier residual. We can therefore set uk+1 = GD(∇222 g(xk , y k ), qk ; ηk , T2 ) without explicitly forming Pk . We then set dkx = ∇1 f (xk , y k ) − ∇212 g(xk , y k )uk+1 , and update xk+1 = xk − αk dkx . Algorithm 2 summarizes Newton-based GALET. Unlike GALET [Xiao et al., 2023], which solves a squared-Hessian system, Newton-based GALET solves the LL Hessian system directly, requiring only one Hessian-vector product per inner GD step instead of two.

2.3

Newton-Based BAMM

We next consider problems with a convex LL objective and an UL objective strongly convex in the LL variable. 5

Algorithm 2 Newton-Based GALET 1: Input: x0 , y 0 ; K, T1 , T2 ; step sizes {αk , γk , ηk }. 2: for k = 0, . . . , K − 1 do 3: Define Hk (p) := ∇222 g(xk , y k )p. 4: v k = GD(Hk , ∇2 g(xk , y k ); γk , T1 ). 5: uk+1 = GD(Hk , ∇2 f (xk , y k ); ηk , T2 ). 6: Update y k+1 = y k − v k . 7: Set dkx = ∇1 f (xk , y k ) − ∇212 g(xk , y k )uk+1 . 8: Update xk+1 = xk − αk dkx . 9: end for

Assumption 2. The lower-level function g(x, y) is convex in y, and the upper-level function f (x, y) is µf -strongly convex in y. sl-BAMM [Liu et al., 2023] uses the aggregation function ψσk (x, y) := σk f (x, y) + (1 − σk )g(x, y), where 0 < σk < 1 decreases to zero. By Assumption 2, ψσk (x, ·) is (σk µf )-strongly convex, with unique minimizer yσ∗k (x) := arg miny ψσk (x, y). The corresponding residuals are ry (x, y) = ∇2 ψσk (x, y), ru (x, y, u) = ∇2 f (x, y) − ∇222 ψσk (x, y)u, rx (x, y, u) = ∇1 f (x, y) − ∇212 ψσk (x, y)u. As σk → 0, these residuals approach those in (4)–(6). The LL and multiplier residuals share the aggregation Hessian ∇222 ψσk . Algorithm Design

At iteration k, we approximately solve ∇222 ψσk (xk , y k )v = ∇2 ψσk (xk , y k )

(10)

for v k and set y k+1 = y k − v k . For the multiplier, we use the warm-start update uk+1 = uk − wk , where wk approximately solves ∇222 ψσk (xk , y k )w = ∇222 ψσk (xk , y k )uk − ∇2 f (xk , y k ).

(11)

An exact solution of (11) yields ru (xk , y k , uk+1 ) = 0. The two systems share the same aggregation Hessian. Unlike Newton-based GALET, no rangespace projection is needed because ψσk (xk , ·) is (σk µf )-strongly convex, so its Hessian is positive definite for σk > 0. As σk → 0, this modulus vanishes, and its decay must be controlled in the convergence analysis; see Theorem 3.2. Following the UL residual, we set dkx = ∇1 f (xk , y k ) − ∇212 ψσk (xk , y k )uk and update xk+1 = k x − αk dkx . Algorithm 3 summarizes Newton-based BAMM. Algorithm 3 Newton-Based BAMM K−1

1: Input: x0 , y 0 , u0 ; K, {Tk }k=0 ; step sizes {αk , γk }; aggregation parameters {σk }. 2: for k = 0, . . . , K − 1 do 3: Define Hk (p) := ∇222 ψσk (xk , y k )p. 4: v k = GD(Hk , ∇2 ψσk (xk , y k ); γk , Tk ). 5: wk = GD(Hk , Hk (uk ) − ∇2 f (xk , y k ); γk , Tk ). 6: Update y k+1 = y k − v k and uk+1 = uk − wk . 7: Set dkx = ∇1 f (xk , y k ) − ∇212 ψσk (xk , y k )uk . 8: Update xk+1 = xk − αk dkx . 9: end for

6

3

Convergence Analysis

We establish convergence and complexity guarantees for Newton-based GALET and Newton-based BAMM; proofs are given in Appendix B. Both analyses use a Lyapunov argument. For Newton-based GALET, we combine the Lagrangian-type potential Pk = f (xk , y k )−u∗ (xk , y k )⊤ ∇2 g(xk , y k ), where u∗ (xk , y k ) = [∇222 g(xk , y k )]† ∇2 f (xk , y k ), with the LL residual ∥∇2 g(xk , y k )∥2 . For Newton-based BAMM, we combine Pk = f (xk , yσ∗k (xk )) with the LL and multiplier tracking errors. The key ingredient is the local contraction induced by the inexact Newton LL update: the LL residual for Newton-based GALET and the LL tracking error for Newton-based BAMM, up to the UL perturbation and, for BAMM, the variation of σk . Suitable initialization and parameter choices keep the iterates in the required local regions.

3.1

Convergence Analysis of Newton-Based GALET

In addition to Assumption 1, we impose the following regularity conditions. Assumption 3.

(a) The UL function f (x, y) has an Lf,1 -Lipschitz continuous gradient.

(b) The LL function g(x, y) has an Lg,1 -Lipschitz continuous gradient and an Lg,2 -Lipschitz continuous Hessian. (c) There exists Cf,0 > 0 such that ∥∇1 f (x, y)∥ ≤ Cf,0 for all x, y. (d) There exists Lf,0 > 0 such that ∥∇2 f (x, y)∥ ≤ Lf,0 for all x, y. Remark 2. Except for part (c), these conditions are also assumed in Xiao et al. [2023]. Part (c) bounds the UL direction and helps keep the LL iterates in the local region required for the Newtontype error reduction. It holds, for example, in data hyper-cleaning when the validation objective is independent of the UL variable, so ∇1 f (x, y) = 0. Following the necessary optimality conditions in [Xiao et al., 2023, Theorem 1], we adopt the stationarity measure used in GALET. Definition 3.1. A point (x̄, ȳ) is ϵ-stationary if there exists ū such that Rx (x̄, ȳ, ū) ≤ ϵ, Ru (x̄, ȳ, ū) ≤ ϵ, and Ry (x̄, ȳ) ≤ ϵ, where Rx (x, y, u) := ∥∇1 f (x, y) − ∇212 g(x, y)u∥2 , Ru (x, y, u) := ∥∇222 g(x, y)(∇2 f (x, y) − ∇222 g(x, y)u)∥2 , Ry (x, y) := g(x, y) − min g(x, z). z

Here, Ru is a Hessian-projected multiplier residual, accommodating a possibly singular LL Hessian. Let L := max{Lf,1 , Lg,1 , Lg,2 } and κ := L/µ. Without loss of generality, we assume 0 < µ ≤ 1. We next establish the complexity of Newton-based GALET. Theorem 3.1. Suppose Assumptions 1 and 3 hold, and f ∗ := inf x, y∈Y ∗ (x) f (x, y) > −∞. For Newton-based GALET, choose T1 = Θ(κ), T2 = Θ(κ log(κK)), γk , ηk = Θ(1/Lg,1 ) with 0 < γk , ηk ≤ 1/Lg,1 , and αk = α = cα κ−4 for a sufficiently small constant cα > 0. The constants in T1 and T2 are chosen sufficiently large. Suppose further that ∥∇2 g(x0 , y 0 )∥ ≤ min{µ2 /(2Lg,2 ), µ2 /L}. Then, K−1

K−1

1 X Rx (xk , y k , uk+1 ) = O(κ4 /K), K

1 X Ru (xk , y k , uk+1 ) = O(1/K), K

k=0

1 K

K−1 X

Ry (xk , y k ) = O(1/(κ2 K)),

k=0

∀ K ≥ 1.

k=0

Consequently, Newton-based GALET finds an ϵ-stationary point
in O(κ4 /ϵ) outer iterations, requiring O(κ4 /ϵ) gradient/JVP computations and O κ5 ϵ−1 log(κ/ϵ) Hessian-vector products.

7

Remark 3. The initialization condition in Theorem 3.1 can be satisfied by gradient descent on g(x0 , ·), requiring O(κ log κ) LL gradient evaluations, with dependence on the initial optimality gap and other problem constants omitted. This is a one-time cost. Each subsequent outer iteration requires O(1) gradient/JVP computations and O(T1 + T2 ) Hessian-vector products, yielding the stated oracle-complexity bounds. Under comparable assumptions, Newton-based GALET improves the gradient/JVP complexity of GALET from O(κ6 /ϵ) to O(κ4 /ϵ) and the Hessian-vector-product complexity from O(κ8 ϵ−1 log(κ/ϵ)) to O(κ5 ϵ−1 log(κ/ϵ)). Two features contribute to these improvements. First, the inexact Newton LL update yields a local quadratic reduction of the LL residual up to the linear-system error (Lemma 4). Second, the range inclusion Range (∇221 g(x, y)) ⊆ Range (∇222 g(x, y)) allows the multiplier update to use the LL Hessian directly, rather than the squared-Hessian system used in GALET. Compared with the GALET analysis [Xiao et al., 2023], two technical differences are worth noting: (a) We establish a global lower bound for the Lyapunov function using the LL error bound and the lower boundedness of the bilevel objective, enabling the telescoping argument underlying our complexity guarantees. (b) Our stationarity bounds are evaluated at (xk , y k , uk+1 ) rather than (xk , y k , uk ), as the cold-start multiplier update directly controls the Hessian-weighted multiplier residual at uk+1 without requiring a uniform bound on {uk }.

3.2

Convergence Analysis of Newton-Based BAMM

In addition to Assumption 2, we impose the following regularity conditions. Assumption 4. (a) The UL function f (x, y) has an Lf,1 -Lipschitz continuous gradient and an Lf,2 -Lipschitz continuous Hessian, and is bounded below by f0 . (b) The LL function g(x, y) has an Lg,1 -Lipschitz continuous gradient and an Lg,2 -Lipschitz continuous Hessian. (c) There exists Cf,0 > 0 such that ∥∇1 f (x, y)∥ ≤ Cf,0 for all x, y. Except for part (c), these conditions are also assumed in Liu et al. [2023]. As in the Newton-based GALET analysis, part (c) is used to preserve the local region required for the Newton-type error reduction; see Remark 2. Under Assumptions 2 and 4, ψσk (x, ·) is (σk µf )-strongly convex, with Lψ,1 -Lipschitz gradient and Lψ,2 -Lipschitz Hessian, where Lψ,1 and Lψ,2 are defined in (32). Define yσ∗ (x) := arg miny ψσ (x, y) and u∗σ (x) := [∇222 ψσ (x, yσ∗ (x))]−1 ∇2 f (x, yσ∗ (x)). Following Liu et al. [2023], we call (x, y, u) ϵ-stationary if KKT(x, y, u) ≤ ϵ, where KKT(x, y, u) := ∥∇1 f (x, y) − ∇212 g(x, y)u∥2 + ∥∇2 f (x, y) − ∇222 g(x, y)u∥2 + ∥∇2 g(x, y)∥2 . As in sl-BAMM [Liu et al., 2023], we analyze three boundedness settings. Throughout this subsection, O(·) and Θ(·) suppress dependence on µf , Lf,1 , Lf,2 , Lg,1 , and Lg,2 . Theorem 3.2. Suppose Assumptions 2 and 4 hold, and let 0 < τ < 1/3. For Newton-based BAMM, choose σk = σ̄(k + k0 )−p with σ̄ > 0 and 0 < σk < 1, γk = Θ(1/Lψ,1 ) with 0 < γk ≤ 1/Lψ,1 , and Tk = Θ(1/σk ) sufficiently large so that Tk ≥ log(1/4)/ log(1 − µf σk γk ). Let ᾱ > 0 be fixed, choose k0 sufficiently large, and suppose ∥y 0 − yσ∗0 (x0 )∥ ≤ µf σ0 /(2Lψ,2 ). Then the following results hold. • Case 1: Suppose {xk }, {y k }, and {uk } are bounded. For 0 < p < (1 − 3τ )/6, choose αk = ᾱ(k + k0 )−4p−2τ . Then   min KKT(xk , y k , uk ) = O K −(1−6p−3τ ) + K −2p . 0≤k≤K−1

Choosing p = (1 − 3τ )/8 yields an ϵ-stationary iterate using O(ϵ−4/(1−3τ ) ) gradient/JVP computations and O(ϵ−(9−3τ )/(2−6τ ) ) Hessian-vector products.

8

Setting

Condition

PL with positive spectral gap

— xk , y k , uk bounded

Method

Gradient/JVP computations

Hessian-vector products

GALET [Xiao et al., 2023]

O(κ6 ϵ−1 )

O(κ8 ϵ−1 log(κ/ϵ))

Newton-based GALET

O(κ4 ϵ−1 )

O(κ5 ϵ−1 log(κ/ϵ))

sl-BAMM [Liu et al., 2023] O(ϵ−11/(2−6τ ) ) O(ϵ−4/(1−3τ ) )

O(ϵ−(9−3τ )/(2−6τ ) )

O(ϵ−7/(2−6τ ) )

O(ϵ−7/(2−6τ ) )

O(ϵ−3/(1−3τ ) )

O(ϵ−(7−3τ )/(2−6τ ) )

sl-BAMM [Liu et al., 2023]

O(ϵ−5/(2−6τ ) )

O(ϵ−5/(2−6τ ) )

Newton-based BAMM

O(ϵ−2/(1−3τ ) )

O(ϵ−(5−3τ )/(2−6τ ) )

Newton-based BAMM g(x, ·) convex, ∗ k f (x, ·) strongly convex ∇2 fσk , u sl-BAMM [Liu et al., 2023] bounded Newton-based BAMM u∗σk , uk bounded

O(ϵ−11/(2−6τ ) )

Table 1: Complexity comparison with the corresponding baselines beyond LL strong convexity. For BAMM, O(·) suppresses dependence on µf , Lf,1 , Lf,2 , Lg,1 , Lg,2 . • Case 2: Suppose {∇2 f (xk , yσ∗k (xk ))} and {uk } are bounded. For 0 < p < (1 − 3τ )/4, choose αk = ᾱ(k + k0 )−3p−2τ . Then   min KKT(xk , y k , uk ) = O K −(1−4p−3τ ) + K −2p . 0≤k≤K−1

Choosing p = (1 − 3τ )/6 yields an ϵ-stationary iterate using O(ϵ−3/(1−3τ ) ) gradient/JVP computations and O(ϵ−(7−3τ )/(2−6τ ) ) Hessian-vector products. • Case 3: Suppose {u∗σk (xk )} and {uk } are bounded. For 0 < p < (1 − 3τ )/2, choose αk = ᾱ(k + k0 )−2p−2τ . Then   min KKT(xk , y k , uk ) = O K −(1−2p−3τ ) + K −2p . 0≤k≤K−1

Choosing p = (1 − 3τ )/4 yields an ϵ-stationary iterate using O(ϵ−2/(1−3τ ) ) gradient/JVP computations and O(ϵ−(5−3τ )/(2−6τ ) ) Hessian-vector products. Remark 4. The initialization condition in Theorem 3.2 can be satisfied by applying gradient descent to ψσ0 (x0 , ·) with x0 and σ0 fixed. By strong convexity and smoothness, this requires O(σ0−1 log(1/σ0 )) LL gradient evaluations, with dependence on the initial distance and problem constants omitted. Since k0 can be chosen independently of K, this one-time cost does not change the stated asymptotic oracle-complexity bounds. We conclude this section by comparing the complexity bounds of the proposed methods with their corresponding baselines. Table 1 summarizes the results. Outer-iteration complexities coincide with gradient/JVP complexities for all methods shown and are therefore omitted. Here, ∇2 fσ∗k := ∇2 f (xk , yσ∗k (xk )). The GALET bounds are derived from the analysis in Appendix F of Xiao et al. [2023], as detailed in Appendix B.3. The sl-BAMM bounds follow from Theorems 3.4–3.6 of Liu et al. [2023].

4

Numerical Experiments

We evaluate the proposed methods on two toy examples, data hyper-cleaning, and GPT-2 data reweighting.

9

4.1

Toy Examples

We first consider the positive spectral-gap setting: f (x, y) =

m X

log cosh(xi − 2) +

i=1

m X

log cosh([y1 ]i + 3[y2 ]i − 2[x]i ),

i=1

(12)

m X ζ g(x, y) = ∥y1 − y2 − x∥2 + log cosh([y1 ]i − [y2 ]i − [x]i ), 2 i=1

where x ∈ Rm , y = (y1 , y2 ) ∈ R2m , and cosh(t) = (et + e−t )/2. The LL objective is not strongly convex in y and has solution set {y : y1 = y2 + x}. Its positive Hessian eigenvalues are uniformly bounded below by 2ζ, so Assumption 1 holds. Figures 1 and 2 compare Newton-based GALET with GALET and GALET-Q, a variant that replaces Line 7 of Algorithm 1 in Xiao et al. [2023] by
wk,t+1 = wk,t − ρ ∇2 f (xk , y k+1 ) + ∇222 g(xk , y k+1 )wk,t . GALET-Q thus applies GD directly to the multiplier linear system rather than the squared-residual problem used in GALET, isolating the effect of this modification. We set N = 1 and T = 5 for GALET and GALET-Q, and T1 = 2, T2 = 6 for Newton-based GALET. Since the GD subroutine starts from p0 = 0, the first update p1 = γq is curvatureindependent, so T iterations contain T − 1 curvature-dependent refinements. Although the analysis uses T1 = Θ(κ), we find that T1 = 2 suffices in this example, consistent with the empirical observation in Dong et al. [2025]. We use additional iterations for the cold-start multiplier solve. Figure 1 shows faster residual reduction for GALET-Q than GALET and a further improvement for Newton-based GALET. Figure 2 shows that increasing T2 reduces Ru , consistent with Lemma 5. We also vary the condition number in (12), where κ = Θ(1/ζ). Figure 3 shows that Newton-based GALET degrades less as κ increases than GALET in this experiment. 104

GALET GALET-Q Newton-based GALET

102 100

101 10 1

R y (x, y)

R x (x, y, u)

10 2

10 3

10 4

10 5

10 6 10 8

10 7

10 10

10 9

0.0

GALET GALET-Q Newton-based GALET

103

0.5

1.0

1.5 2.0 Time(s)

2.5

3.0

0.0

0.5

1.0

1.5 2.0 Time(s)

2.5

3.0

Figure 1: GALET, GALET-Q, and Newton-based GALET on (12) with ζ = 1. Left: Rx . Right: Ry . 105

GALET GALET-Q Newton-based GALET

104

101

R u (x, y, u)

R u (x, y, u)

103

10 1

102

10 3

101 0.0

T2 = 6 T2 = 11 T2 = 21

103

10 5 0.5

1.0

1.5 2.0 Time(s)

2.5

3.0

0.0

0.5

1.0

1.5 2.0 Time(s)

2.5

3.0

Figure 2: Results on (12) with ζ = 1. Left: Ru for GALET, GALET-Q, and Newton-based GALET. Right: Ru for Newton-based GALET with T2 = 6, 11, 21. 10

104

GALET (ζ = 0.5) Newton-based GALET (ζ = 0.5) GALET (ζ = 0.1) Newton-based GALET (ζ = 0.1)

102 100

GALET (ζ = 0.5) Newton-based GALET (ζ = 0.5) GALET (ζ = 0.1) Newton-based GALET (ζ = 0.1)

103 101 10 1

R y (x, y)

R x (x, y, u)

10 2

10 3

10 4

10 5

10 6 10 8

10 7

10 10

10 9

0.0

0.5

1.0

1.5 2.0 Time(s)

2.5

3.0

0.0

0.5

1.0

1.5 2.0 Time(s)

2.5

3.0

Figure 3: GALET and Newton-based GALET on (12) with ζ = 0.5, 0.1. Left: Rx . Right: Ry . For the second setting, with a convex LL objective and an UL objective strongly convex in the LL variable, consider f (x, y) =

m X

1 1 log cosh xi − [y2 ]i + ∥y1 − e∥2 + ∥y2 ∥2 , 2 2 i=1

(13)

1 g(x, y) = ∥y1 ∥2 − ⟨x, y1 ⟩, 2 where x ∈ Rm , y = (y1 , y2 ) ∈ R2m , and e is the all-ones vector. The first term of f is adapted from [Liu et al., 2023, Eq. (13)] to ensure that ∥∇1 f (x, y)∥ is uniformly bounded. The LL solution set is {y : y1 = x, y2 ∈ Rm }, while f (x, ·) is strongly convex. Figure 4 compares Newton-based BAMM with sl-BAMM. We set Tk = 2 for Newton-based BAMM. On this example, Newton-based BAMM shows faster reduction in both the LL and KKT residuals. sl-BAMM Newton-based BAMM

102

101

100 k∇ 2 g(x, y)k 2

KKT (x, y, u)

10 2 10 4 10 6 10 8

10 1 10 3 10 5 10 7

10 10 0.0

sl-BAMM Newton-based BAMM

103

0.5

1.0

1.5

Time(s)

2.0

2.5

3.0

0.0

0.5

1.0

1.5

Time(s)

2.0

2.5

3.0

Figure 4: Newton-based BAMM and sl-BAMM on (13), with stepsizes chosen according to Case 1 of Theorem 3.2.

4.2

Data Hyper-Cleaning

We next consider data hyper-cleaning, a standard bilevel optimization benchmark. The UL variable λ weights the training samples, while the LL variable ω represents the model parameters. We consider both linear and MLP models. For the linear model, we use f (λ, ω) =

1 |D′ |

X

L(ωx′e , ye′ ) + cr ∥ω∥2 ,

(x′e ,ye′ )∈D ′

g(λ, ω) =

1 |D|

X

σ(λe )L(ωxe , ye ),

(xe ,ye )∈D

where D and D′ denote the training and validation sets, respectively, L is the cross-entropy loss, and σ(λe ) is the weight assigned to the e-th training sample. The LL objective is convex in ω, while the quadratic regularizer with cr > 0 makes the UL objective strongly convex in ω, matching the setting of Newton-based BAMM. 11

We compare Newton-based BAMM with sl-BAMM on MNIST and FashionMNIST, using Tk = 2 for Newton-based BAMM. Figure 5 shows faster reduction in the reported metrics for Newton-based BAMM on both datasets. For the MLP model, we consider f (λ, ω) =

1 |D′ |

X

L(N (x′e , ω), ye′ ),

g(λ, ω) =

(x′e ,ye′ )∈D ′

1 |D|

X

σ(λe )L(N (xe , ω), ye ),

(xe ,ye )∈D

where N (·, ω) is a two-layer MLP. Unlike the linear setting, neither the convexity of the LL objective nor the strong convexity of the UL objective in ω is guaranteed. We therefore compare Newton-based GALET with GALET. Figure 6 shows that Newton-based GALET achieves higher accuracy within the reported running time on both datasets. 1.2

0.95

0.80 0.75 0.70 0.65 0.60

0

5

10

Newton-based BAMM sl-BAMM 15 20 25 30

0.6

0

Time(s)

5

10

1.2

Validation loss

0.75 0.70 0.65 0

5

10

Newton-based BAMM sl-BAMM 15 20 25 30

Time(s)

15

Time(s)

20

25

30

Newton-based BAMM sl-BAMM

1.1

0.80

Test accuracy

0.8

0.4

0.85

0.60

Newton-based BAMM sl-BAMM

1.0

0.85

Validation loss

Test accuracy

0.90

1.0 0.9 0.8 0.7 0.6 0.5 0

5

10

15

Time(s)

20

25

30

Figure 5: Data hyper-cleaning with linear models. Comparison of sl-BAMM and Newton-based BAMM on MNIST (upper) and FashionMNIST (lower).

4.3

Data Denoising with GPT-2

We further evaluate Newton-based GALET on the data-denoising task of Pan et al. [2025], using GPT-2 with 124M parameters. The validation set Dval contains 1,000 clean Alpaca samples [Taori et al., 2023], while the training set consists of 1,000 clean samples D1 and 9,000 corrupted samples D2 , whose responses are replaced by “.”. The goal is to assign a large weight to D1 and a small weight to D2 . We formulate the task as f (λ, w) =

X 1 ℓ(w; z), |Dval |

g(λ, w) =

2 X pi (λ) X i=1

z∈Dval

|Di |

ℓ(w; z),

z∈Di

where w denotes the GPT-2 parameters, ℓ(w; z) is the causal language modeling loss, and pi (λ) := P2 exp(λi )/ j=1 exp(λj ) is the sampling weight of Di . We compare Newton-based GALET with ScaleBiO [Pan et al., 2025], a single-loop first-order bilevel method. For Newton-based GALET, we set T1 = T2 = 2; both methods use AdamW for all variable updates. Figure 7 shows that both methods increase the clean-source weight and decrease the corrupted-source weight, with Newton-based GALET adjusting the weights faster initially. After three epochs, Table 2 shows similar source weights but lower validation and test losses for 12

Validation loss

Test accuracy

1.2

0.8 0.7 0.6 0.5 0.4

0

10

20

30

0.8 0.6

0.2

Time(s)

0

10

20

30

0.7 0.6 0.5 0

10

20

30

40

Time(s)

50

60

70

80

GALET Newton-based GALET

1.4

Validation loss

Test accuracy

1.0

0.4

GALET Newton-based GALET 40 50 60 70 80

0.8

0.4

GALET Newton-based GALET

1.4

0.9

1.2 1.0 0.8 0.6

GALET Newton-based GALET 40 50 60 70 80

0

Time(s)

10

20

30

40

Time(s)

50

60

70

80

Figure 6: Data hyper-cleaning with two-layer MLPs. Comparison of GALET and Newton-based GALET on MNIST (upper) and FashionMNIST (lower). 1.0

Weight

0.8 ScaleBiO (clean data) ScaleBiO (noisy data) Newton-based GALET (clean data) Newton-based GALET (noisy data)

0.6 0.4 0.2 0.0

0

100

200

300

400

500

Step

Figure 7: Clean and corrupted data weights for ScaleBiO and Newton-based GALET, averaged over five random seeds. Newton-based GALET, together with higher peak GPU memory due to the additional second-order computations. Method

Clean weight Corrupted weight Validation loss Test loss Peak memory

ScaleBiO

0.928

0.072

2.603

2.555

9.984 GiB

Newton-based GALET

0.929

0.071

2.439

2.398

18.509 GiB

Table 2: GPT-2 data denoising after three training epochs, averaged over five random seeds.

5

A Generalized Newton-Based Framework

The proposed methods share a common structure: the LL and multiplier updates exploit curvature from the same LL model. Newton-based GALET uses g with potentially singular curvature, whereas Newton-based BAMM uses the strongly convex aggregation model ψσk . Related structures also 13

appear in Newton and quasi-Newton bilevel methods [Ramzi et al., 2022, Fang et al., 2025, Dong et al., 2025]. We capture these common features in a generalized Newton-based framework for bilevel optimization (GNBO). Lower-level model

At iteration k, GNBO selects an LL model hk (x, y) with residuals ry (x, y) = ∇2 hk (x, y), ru (x, y, u) = ∇2 f (x, y) − ∇222 hk (x, y)u, rx (x, y, u) = ∇1 f (x, y) − ∇212 hk (x, y)u.

Newton-based GALET takes hk = g, while Newton-based BAMM takes hk = ψσk ; in both cases, hk supplies the curvature shared by the LL and multiplier updates. Lower-level update A Newton-type LL update targets ∇222 hk (xk , y)v = ∇2 hk (xk , y). To allow one or multiple LL corrections per outer iteration, let y 0,k = y k , . . . , y S,k = y k+1 and consider Bvs,k v = qvs,k ,

(14)

where Bvs,k approximates the curvature and qvs,k is derived from the LL residual. A subroutine Av approximately solves (14), followed by y s+1,k = y s,k − βs,k v s,k . Here, S counts LL corrections and T1k denotes the internal solver effort. The multiplier and UL updates may use either the pre- or post-update LL point: ỹ k = k δy y + (1 − δy )y k+1 , where δy ∈ {0, 1}. Thus, δy = 0 uses y k+1 , while δy = 1 uses y k , as in Newton-based GALET, Newton-based BAMM, and NBO-GD. Multiplier update At ỹ k , the ideal multiplier system is ∇222 hk (xk , ỹ k )u = ∇2 f (xk , ỹ k ). GNBO represents its treatment by Buk (δu uk − w) = quk , (15) where Buk and quk specify the curvature approximation and right-hand side. The indicator δu ∈ {0, 1} distinguishes cold and warm starts: δu = 0 gives uk+1 = −wk , as in Newton-based GALET, while δu = 1 gives uk+1 = uk − wk , as in Newton-based BAMM and NBO-GD. A subroutine Aw approximately solves (15), with T2k denoting the internal solver effort. Upper-level update

Set ũk = δu uk + (1 − δu )uk+1 and dkx = ∇1 f (xk , ỹ k ) − ∇212 hk (xk , ỹ k )ũk .

Then xk+1 = xk − αk dkx . Algorithm 4 summarizes the framework. Algorithm 4 Generalized Newton-Based Bilevel Optimization (GNBO) 1: Input: x0 , y 0 , u0 ; K, S, {T1k , T2k }; step sizes {αk , βs,k }; δy , δu ∈ {0, 1}. 2: for k = 0, . . . , K − 1 do 3: Choose an LL model hk and set y 0,k = y k . ▷ Lower-Level Update 4: for s = 0, . . . , S − 1 do

v s,k = Av (hk ; xk , y s,k ; T1k ) based on (14). y s+1,k = y s,k − βs,k v s,k . end for Set y k+1 = y S,k and ỹ k = δy y k + (1 − δy )y k+1 . wk = Aw (hk ; xk , ỹ k , uk ; T2k ) based on (15). 10: Set uk+1 = δu uk − wk and ũk = δu uk + (1 − δu )uk+1 . 11: Set dkx = ∇1 f (xk , ỹ k ) − ∇212 hk (xk , ỹ k )ũk . 12: Set xk+1 = xk − αk dkx . 13: end for 5: 6: 7: 8: 9:

▷ Multiplier Update ▷ Upper-Level Update

The number of LL corrections S is distinct from the internal solver effort: S = 1 means one LL update per outer iteration, even if computing its direction requires multiple solver iterations. Table 3 summarizes representative instances of GNBO. 14

Method

LL model LL corrections

SHINE [Ramzi et al., 2022] qNBO [Fang et al., 2025] NBO-GD [Dong et al., 2025]

Solver / curvature

(δy , δu ) (0, 0)

g

Multiple

Quasi-Newton, reused

g

Multiple

Separate quasi-Newton (0, 0)

g

One

Newton-based GALET

g

One

Newton-based BAMM

ψσk

One

T GD steps, common Hessian T1 /T2 GD steps, common Hessian Tk GD steps, common Hessian

(1, 1) (1, 0) (1, 1)

Table 3: Representative instances of GNBO.

6

Concluding Remarks

We developed Newton-based methods for bilevel optimization beyond LL strong convexity and captured their shared-curvature structure through GNBO. For the settings studied, exploiting the same curvature in the LL and multiplier updates leads to improved complexity bounds over the corresponding baselines. Several directions remain open. Hessian-vector products introduce additional memory and computational costs, particularly for large-scale neural networks. Improving scalability through parameter-efficient strategies such as LoRA [Hu et al., 2022] and LISA [Pan et al., 2024], or limitedmemory quasi-Newton approximations [Liu and Nocedal, 1989], is therefore of interest. Moreover, our methods use Newton-type LL and multiplier updates but retain a first-order UL update. Extending this approach to Newton-type UL updates under suitable higher-order smoothness conditions, as explored in recent work [Bi et al., 2026, Yang et al., 2026], is another direction for future study.

15

A

Experimental Details

All experiments are conducted on a machine with an Intel Xeon Gold 5218R CPU (2.10 GHz) and an NVIDIA H100 GPU (80 GB). The toy and data hyper-cleaning experiments are repeated over 10 runs, while the GPT-2 experiments use five random seeds. Toy examples For GALET, GALET-Q, and Newton-based GALET, we set αk = ηk = 0.1 and γk = 0.01. For sl-BAMM, we use p = 37/440 and τ = 1/40; for Newton-based BAMM, we use p = 37/320 and τ = 1/40. Data hyper-cleaning For GALET and Newton-based GALET, we jointly tune (αk , βk , ηk ) by grid search, with αk ∈ {103 , 104 , 105 } and γk , ηk ∈ {10−3/2 , 10−1 , 10−1/2 } for MNIST or {10−5/2 , 10−2 , 10−3/2 } for FashionMNIST. For sl-BAMM and Newton-based BAMM, we use (p, τ ) = (37/440, 1/40) and (37/320, 1/40), respectively. Data denoising with GPT-2 ScaleBiO uses AdamW with a constant learning rate of 10−5 for both model-variable updates. Newton-based GALET uses AdamW with the same learning rate for both linear-system solves. In both methods, the source-weight logits are updated using AdamW with a constant learning rate of 10−2 . GPT-2 parameters are stored in BF16, with LISA disabled.

B

Proofs of Main Results and Auxiliary Lemmas

This appendix provides detailed proofs of the main convergence results, together with the statements and proofs of the auxiliary lemmas used in the analysis.

B.1

Convergence Analysis of Newton-Based GALET

Proof of Lemma 1. Fix x. By Weyl’s inequality and the Lg,2 -Lipschitz continuity of ∇222 g(x, ·), its eigenvalues vary continuously with y. Assumption 1 prevents any nonzero eigenvalue from crossing zero, so the inertia and rank of ∇222 g(x, y) are independent of y. At any y ∗ ∈ Y ∗ (x), the second-order necessary condition gives ∇222 g(x, y ∗ ) ⪰ 0. Hence ∇222 g(x, y) ⪰ 0 for all y, and g(x, ·) is convex. Denote the common rank by r(x). Apply [Hartman and Nirenberg, 1959, Lemma 2] to p(y) := ∇2 g(x, y), whose Jacobian has constant rank r(x). Since p(y ∗ ) = 0, its zero set is locally an affine plane y ∗ + L(x) of dimension n − r(x). Convexity of Y ∗ (x) = {y : p(y) = 0} extends this local characterization globally: Y ∗ (x) = y ∗ + L(x). Moreover, differentiating p(y ∗ + td) = 0 at t = 0 for d ∈ L(x) gives L(x) ⊆ ker(∇222 g(x, y ∗ )). Since both spaces have dimension n − r(x), Y ∗ (x) = y ∗ + K(x),

K(x) := ker(∇222 g(x, y ∗ )).

(16)

Let S(x) := K(x)⊥ . For d ∈ K(x), convexity and (16) imply ⟨∇2 g(x, y), y − y ∗ − td⟩ ≥ 0 for every t ∈ R. Hence ∇2 g(x, y) ⊥ K(x), so ∇2 g(x, y) ∈ S(x). Differentiating this relation with respect to y gives Range (∇222 g(x, y)) ⊆ S(x). Both spaces have dimension r(x); therefore Range (∇222 g(x, y)) = S(x),

∇2 g(x, y) ∈ S(x).

(17)

Let ŷ be the projection of y onto Y ∗ (x). By (16), y − ŷ ∈ S(x). Since the positive eigenvalues of the Hessian are at least µ, (17) gives Z 1 ⟨∇2 g(x, y), y − ŷ⟩ =

(y − ŷ)⊤ ∇222 g x, ŷ + s(y − ŷ) (y − ŷ) ds ≥ µ∥y − ŷ∥2 .

0

Thus, by Cauchy–Schwarz, ∥∇2 g(x, y)∥ ≥ µ∥y − ŷ∥ = µ dist (y, Y ∗ (x)). Finally, suppose ∇2 g(·, y) is Lg,1 -Lipschitz. For d ∈ K(x), (16) and the error bound give, for any x′ and t ∈ R, Lg,1 ′ dist (y ∗ (x) + td, Y ∗ (x′ )) ≤ ∥x − x∥. µ 16

Since Y ∗ (x′ ) = y ∗ (x′ ) + K(x′ ), dividing by t > 0 and letting t → ∞ yields d ∈ K(x′ ). Interchanging x and x′ shows that K(x) = K(x′ ) =: K, so S := K ⊥ is independent of x. Since ∇2 g(x, y) ∈ S for all x, y, differentiating with respect to x gives Range (∇221 g(x, y)) ⊆ S. By (17), S = Range (∇222 g(x, y)), and hence Range (∇221 g(x, y)) ⊆ Range (∇222 g(x, y)). Following in part Xiao et al. [2023], define F (x, y; u) := f (x, y) − u⊤ ∇2 g(x, y), with u treated as fixed, and u∗ (x, y) := [∇222 g(x, y)]† ∇2 f (x, y). For the subsequent analysis, let √ Lf,0 Lg,2 Lf,1 (1 + 5)Lf,0 Lg,2 LF := Lf,1 + + , Lu∗ := . (18) µ 2µ2 µ Lemma 2. Suppose Assumptions 1 and 3 hold. Then
∥u∗ (x, y) − u∗ (x′ , y ′ )∥ ≤ Lu∗ ∥x − x′ ∥ + ∥y − y ′ ∥ . Proof. By Lemma 1, ∇222 g(x, y) has a common range, and hence constant rank, over (x, y). The result follows from the pseudoinverse perturbation bound [Stewart, 1977, Theorem 3.3] and the Lipschitz continuity of ∇222 g and ∇2 f ; see also [Xiao et al., 2023, Lemma 9]. Lemma 3. Suppose Assumptions 1 and 3 hold. Then Newton-based GALET satisfies F (xk+1 , y k+1 ; u∗ (xk+1 , y k+1 )) − F (xk , y k ; u∗ (xk , y k )) ! 3L2g,1 αk2 αk2 LF αk2 3 2 2 2 ≤ −αk + + + µ Lu∗ αk + ∥dkx ∥2 + 2 ∥∇2 g(xk , y k )∥2 2 2 2µ2 2µ ! 3L2g,1 L2g,1 LF 2 2 k 2 + µ Lu∗ + ∥v ∥ + ∥∇222 g(xk , y k )uk+1 − ∇222 g(xk , y k )u∗ (xk , y k )∥2 . + 2 2µ2 2µ2 Proof. Let gk := ∇2 g(xk , y k ) and u∗k := u∗ (xk , y k ). For fixed u, F (·, ·; u) has an (Lf,1 + Lg,2 ∥u∥)Lipschitz continuous gradient. Since ∥u∗k ∥ ≤ Lf,0 /µ, F (·, ·; u∗k ) has an LF -Lipschitz continuous gradient. Decompose the potential difference as ∆1 + ∆2 , where ∆1 := F (xk+1 , y k+1 ; u∗k ) − F (xk , y k ; u∗k ), ∆2 := F (xk+1 , y k+1 ; u∗ (xk+1 , y k+1 )) − F (xk+1 , y k+1 ; u∗k ). By LF -smoothness, ∆1 ≤ ⟨∇1 F (xk , y k ; u∗k ), xk+1 − xk ⟩ + ⟨∇2 F (xk , y k ; u∗k ), y k+1 − y k ⟩
LF + ∥xk+1 − xk ∥2 + ∥y k+1 − y k ∥2 . 2

(19)

By Lemma 1, Range (∇221 g(xk , y k )) ⊆ Range (∇222 g(xk , y k )). Since ∇222 g(xk , y k ) is symmetric, this implies ∇212 g(xk , y k ) = ∇212 g(xk , y k )[∇222 g(xk , y k )]† ∇222 g(xk , y k ). Using this identity and xk+1 − xk = −αk dkx , ⟨∇1 F (xk , y k ; u∗k ), xk+1 − xk ⟩ = −αk ∥dkx ∥2
+ ∇212 g(xk , y k )[∇222 g(xk , y k )]† ∇222 g(xk , y k )uk+1 − ∇222 g(xk , y k )u∗k , xk+1 − xk   L2g,1 α2 ≤ − αk − k ∥dkx ∥2 + ∥∇222 g(xk , y k )uk+1 − ∇222 g(xk , y k )u∗k ∥2 , 2 2µ2 where we used ∥∇212 g(xk , y k )∥op ≤ Lg,1 and ∥[∇222 g(xk , y k )]† ∥op ≤ 1/µ. For the y-term, Algorithm 1 gives TX 1 −1 v k = γk (I − γk ∇222 g(xk , y k ))i gk ∈ Range (∇222 g(xk , y k )). (20) i=0

17

Let Pk := ∇222 g(xk , y k )[∇222 g(xk , y k )]† be the orthogonal projector onto Range (∇222 g(xk , y k )). Since ∇222 g(xk , y k )u∗k = Pk ∇2 f (xk , y k ), ∇2 F (xk , y k ; u∗k ) = (I − Pk )∇2 f (xk , y k ) ∈ (Range (∇222 g(xk , y k )))⊥ . As y k+1 − y k = −v k ∈ Range (∇222 g(xk , y k )), the corresponding inner product in (19) vanishes. Hence   αk2 LF αk2 LF k 2 ∆1 ≤ −αk + + ∥dkx ∥2 + ∥v ∥ 2 2 2 L2g,1 ∥∇222 g(xk , y k )uk+1 − ∇222 g(xk , y k )u∗k ∥2 . (21) + 2µ2 For ∆2 , linearity of F in u and Young’s inequality give ∆2 = u∗k − u∗ (xk+1 , y k+1 ), ∇2 g(xk+1 , y k+1 ) ≤

µ2 ∗ k+1 k+1 1 ∥u (x ,y ) − u∗k ∥2 + 2 ∥∇2 g(xk+1 , y k+1 )∥2 . 2 2µ

By Lemma 2 and the updates of x and y,
∥u∗ (xk+1 , y k+1 ) − u∗k ∥2 ≤ 2L2u∗ αk2 ∥dkx ∥2 + ∥v k ∥2 . Moreover, the Lg,1 -Lipschitz continuity of ∇2 g gives ∥∇2 g(xk+1 , y k+1 )∥2 ≤ 3∥gk ∥2 + 3L2g,1 αk2 ∥dkx ∥2 + 3L2g,1 ∥v k ∥2 . Therefore, ∆2 ≤

2

µ L2u∗ αk2 +

3L2g,1 αk2 2µ2

! ∥dkx ∥2 +

3 ∥gk ∥2 + 2µ2

2

µ L2u∗ +

3L2g,1 2µ2

! ∥v k ∥2 .

(22)

Combining (21) and (22) proves the result. Lemma 4. Suppose Assumptions 1 and 3 hold. In Newton-based GALET, let 0 ≤ γk ≤ 1/Lg,1 and T1 ≥ log(1/4)/ log(1 − µγk ). If ∥∇2 g(xk , y k )∥ ≤ µ2 /(2Lg,2 ), then ∥∇2 g(xk+1 , y k+1 )∥ ≤

1 ∥∇2 g(xk , y k )∥ + Lg,1 ∥xk+1 − xk ∥. 2

Proof. Let gk := ∇2 g(xk , y k ). Since y k+1 = y k − v k and ∇222 g(xk , ·) is Lg,2 -Lipschitz continuous, Taylor’s theorem gives ∥∇2 g(xk , y k+1 ) − gk + ∇222 g(xk , y k )v k ∥ ≤

Lg,2 k 2 ∥v ∥ . 2

Define rk := gk − ∇222 g(xk , y k )v k . Then ∥∇2 g(xk , y k+1 )∥ ≤

Lg,2 k 2 ∥v ∥ + ∥rk ∥. 2

(23)

By Lemma 1, gk ∈ Range (∇222 g(xk , y k )). Using (20) and the fact that the positive eigenvalues of ∇222 g(xk , y k ) lie in [µ, Lg,1 ], we have ∥v k ∥ ≤

1 ∥gk ∥. µ

(24)

To bound rk , let v t,k be the t-th iterate of GD(∇222 g(xk , y k ), gk ; γk , T1 ) and set rt,k := gk − ∇222 g(xk , y k )v t,k . Then rt+1,k = (I − γk ∇222 g(xk , y k ))rt,k . Since rt,k ∈ Range (∇222 g(xk , y k )), Assumption 1 and γk ≤ 1/Lg,1 imply ∥rt+1,k ∥ ≤ (1 − µγk )∥rt,k ∥. 18

As v 0,k = 0, after T1 iterations, ∥rk ∥ ≤ (1 − µγk )T1 ∥gk ∥ ≤

1 ∥gk ∥. 4

(25)

Combining (23), (24), and (25), ∥∇2 g(xk , y k+1 )∥ ≤

Lg,2 1 1 ∥gk ∥2 + ∥gk ∥ ≤ ∥gk ∥, 2 2µ 4 2

where the last inequality uses ∥gk ∥ ≤ µ2 /(2Lg,2 ). Finally, the Lg,1 -Lipschitz continuity of ∇2 g gives ∥∇2 g(xk+1 , y k+1 )∥ ≤ ∥∇2 g(xk , y k+1 )∥ + Lg,1 ∥xk+1 − xk ∥ ≤

1 ∥gk ∥ + Lg,1 ∥xk+1 − xk ∥, 2

which proves the result. Lemma 5. Suppose Assumptions 1 and 3 hold. In Newton-based GALET, if 0 ≤ ηk ≤ 1/Lg,1 , then ∇222 g(xk , y k )uk+1 − ∇222 g(xk , y k )[∇222 g(xk , y k )]† ∇2 f (xk , y k ) ≤ Lf,0 (1 − µηk )T2 . Proof. Let ut,k be the t-th iterate of GD(∇222 g(xk , y k ), ∇2 f (xk , y k ); ηk , T2 ) and define et,k := ∇222 g(xk , y k )ut,k − ∇222 g(xk , y k )[∇222 g(xk , y k )]† ∇2 f (xk , y k ). The GD iteration gives et+1,k = (I − ηk ∇222 g(xk , y k ))et,k . Since et,k ∈ Range (∇222 g(xk , y k )), the nonzero eigenvalues of ∇222 g(xk , y k ) lie in [µ, Lg,1 ], and ηk ≤ 1/Lg,1 , ∥et+1,k ∥ ≤ (1 − µηk )∥et,k ∥. Moreover, u0,k = 0 and Pk = ∇222 g(xk , y k )[∇222 g(xk , y k )]† is the orthogonal projector onto Range (∇222 g(xk , y k )), so ∥e0,k ∥ = ∥Pk ∇2 f (xk , y k )∥ ≤ ∥∇2 f (xk , y k )∥ ≤ Lf,0 . Therefore, ∥eT2 ,k ∥ ≤ Lf,0 (1 − µηk )T2 . Since uk+1 = uT2 ,k , the result follows. Proof of Theorem 3.1 . We first establish by induction the local-region bound ∥∇2 g(xk , y k )∥ ≤

µ2 , 2Lg,2

k = 0, . . . , K.

(26)

By Lemma 1, Range (∇221 g(xk , y k )) ⊆ Range (∇222 g(xk , y k )). Using ∥u∗ (xk , y k )∥ ≤ Lf,0 /µ, ∥[∇222 g(xk , y k )]† ∥op ≤ 1/µ, and Lemma 5, we obtain ∥dkx ∥ ≤ Cf,0 +

Lg,1 Lf,0 2Lg,1 Lf,0 1 + (1 − µηk )T2 ≤ Cf,0 + =: D. µ µ

If (26) holds at iteration k, Lemma 4 and xk+1 − xk = −αdkx give ∥∇2 g(xk+1 , y k+1 )∥ ≤

µ2 + Lg,1 αD. 4Lg,2

Thus, α ≤ µ2 /(4Lg,1 Lg,2 D) preserves (26). Together with the initialization condition, induction proves (26). This requirement is satisfied by α = Θ(κ−4 ) with a sufficiently small constant factor, under the stated parameter-dependence convention. Let Fk := F (xk , y k ; u∗ (xk , y k )) and define the Lyapunov function Vk := Fk − f ∗ + c∥∇2 g(xk , y k )∥2 +

L2f,0 , cµ2

where c > 0 will be specified below. We first show that Vk ≥ 0. Let ŷ k be the Euclidean projection of y k onto Y ∗ (xk ). By Assumption 3(d) and Lemma 1, f (xk , y k ) ≥ f (xk , ŷ k ) −

Lf,0 Lf,0 ∥∇2 g(xk , y k )∥ ≥ f ∗ − ∥∇2 g(xk , y k )∥. µ µ 19

Moreover, since ∥u∗ (xk , y k )∥ ≤ Lf,0 /µ, Fk ≥ f (xk , y k ) −

Lf,0 2Lf,0 ∥∇2 g(xk , y k )∥ ≥ f ∗ − ∥∇2 g(xk , y k )∥. µ µ

Consequently, completing the square yields  2 Lf,0 k k ≥ 0. Vk ≥ c ∥∇2 g(x , y )∥ − cµ This lower bound holds globally, without requiring (26), which is needed for the Newton-type contraction analysis. Lemma 4 gives ∥∇2 g(xk+1 , y k+1 )∥2 ≤

1 ∥∇2 g(xk , y k )∥2 + 2L2g,1 α2 ∥dkx ∥2 . 2

(27)

Using ∥v k ∥ ≤ µ−1 ∥∇2 g(xk , y k )∥ from (24), Lemmas 3 and 5 together with (27) yield !! 3L2g,1 1 + LF 2 2 2 2 + 2cLg,1 + µ Lu∗ + ∥dkx ∥2 Vk+1 − Vk ≤ −α + α 2 2µ2 ! 2 L2g,1 L2f,0 3L c 3 + LF g,1 2 k k 2 + − + + L (1 − µηk )2T2 . ∥∇ g(x , y )∥ + ∗ + 2 u 2 2µ2 2µ4 2µ2 The second coefficient is at most −c/4 if c≥

6L2g,1 2(3 + LF ) + 4L2u∗ + , 2 µ µ4

while the first is at most −α/2 if α≤

2

1 + LF + 2µ L2u∗ +

3L2g,1 + 4cL2g,1 µ2

!−1 .

The bounds on LF and Lu∗ , together with the parameter-dependence convention, allow us to choose c = Θ(L2 /µ4 ) with a sufficiently large constant and α = Θ(κ−4 ) with a sufficiently small constant. These choices satisfy both descent conditions and preserve (26). Since ηk = Θ(1/Lg,1 ) and Lg,1 /µ ≤ κ, there exists cη > 0, independent of k and κ, such that µηk ≥ cη /κ. Choosing T2 ≥ κ log(κ2 K)/(2cη ) gives (1 − µηk )2T2 ≤ e−2cη T2 /κ ≤

1 . κ2 K

(28)

Thus, T2 = Θ(κ log(κK)) with a sufficiently large constant factor yields L2f,0 α c Vk+1 − Vk ≤ − ∥dkx ∥2 − ∥∇2 g(xk , y k )∥2 + , 2 4 2K

(29)

where we used L2g,1 /µ2 ≤ κ2 . It remains to bound V0 . By its definition, ∥u∗ (x0 , y 0 )∥ ≤ Lf,0 /µ, and the initialization condition, L2f,0 Lf,0 ∥∇2 g(x0 , y 0 )∥ + c∥∇2 g(x0 , y 0 )∥2 + µ cµ2 2 4 Lf,0 Lf,0 cµ ≤ f (x0 , y 0 ) − f ∗ + + 2 + . κ L cµ2

V0 ≤ f (x0 , y 0 ) − f ∗ +

Since c = Θ(L2 /µ4 ) and κ = L/µ ≥ 1, we have V0 = O(1) with respect to κ, provided the initial objective gap and other problem constants are uniformly bounded. 20

Summing (29) over k = 0, . . . , K − 1 and using VK ≥ 0 gives K−1

 4 κ , K k=0   K−1 1 X 1 k k 2 ∥∇2 g(x , y )∥ = O . K cK 1 X k 2 ∥dx ∥ = O K

(30)

(31)

k=0

Since Rx (xk , y k , uk+1 ) = ∥dkx ∥2 , (30) gives the claimed Rx bound. By the PL inequality following Lemma 1, Ry (xk , y k ) ≤ µ−1 ∥∇2 g(xk , y k )∥2 . Thus, (31) and c = Θ(L2 /µ4 ) yield K−1

 µ  1 X =O Ry (xk , y k ) = O 2 K κ K



k=0

1 2 κ K

 ,

where the last step uses 0 < µ ≤ 1. For the multiplier residual, ∇222 g(xk , y k )(I − Pk ) = 0 and Lemma 5 give Ru (xk , y k , uk+1 ) = ∥∇222 g(xk , y k )(Pk ∇2 f (xk , y k ) − ∇222 g(xk , y k )uk+1 )∥2 ≤ L2g,1 L2f,0 (1 − µηk )2T2 . By (28) and Lg,1 /µ ≤ κ, we obtain K−1

1 X Ru (xk , y k , uk+1 ) = O(K −1 ). K k=0

Combining the three bounds, there exists k ∗ ∈ {0, . . . , K − 1} such that all three residuals at ∗ ∗ (x , y k , uk +1 ) are O(κ4 /K). Thus, K = O(κ4 /ϵ) outer iterations suffice for ϵ-stationarity. Each outer iteration requires O(1) gradient/JVP computations and T1 + T2 = O(κ log(κK)) Hessianvector products. With K = O(κ4 /ϵ), the total complexities are O(κ4 /ϵ) gradient/JVP computations and O(κ5 ϵ−1 log(κ/ϵ)) Hessian-vector products. k∗

B.2

Convergence Analysis for Newton-Based BAMM

For convenience, define Lψ,1 := Lg,1 + Lf,1 sup σk ,

Lψ,2 := Lg,2 + Lf,2 sup σk ,

k≥0

k≥0

CΦ1 := L2ψ,1 Lψ,2 + 2Lψ,1 Lψ,2 µf sup σk + Lψ,2 µ2f sup σk2 , k≥0 k≥0

(32)

CΦ2 := L2ψ,1 Lf,1 + 2Lψ,1 Lf,1 µf sup σk + Lf,1 µ2f sup σk2 . k≥0

k≥0

The following lemma gives one-step tracking bounds for the LL and multiplier variables. Lemma 6. Suppose Assumptions 2 and 4 hold. In Newton-based BAMM, let 0 ≤ γk ≤ 1/Lψ,1 and Tk ≥ log(1/4)/ log(1 − µf σk γk ). Then ∥y k+1 − yσ∗k+1 (xk+1 )∥ ≤

1 k Lψ,2 ∥y − yσ∗k (xk )∥ + ∥y k − yσ∗k (xk )∥2 4 2µf σk k+1

+

2∥∇2 f (x Lψ,1 k+1 ∥x − xk ∥ + µf σk

21

(33)

, yσ∗k+1 (xk+1 ))∥ σk − σk+1 µf

σk

,

and ∥uk+1 − u∗σk+1 (xk+1 )∥
5 Lf,1 + Lψ,2 ∥u∗σk (xk )∥ 1 k ∗ k ≤ ∥u − uσk (x )∥ + ∥y k − yσ∗k (xk )∥ 4 4µf σk   Lψ,1 Lf,1 + Lψ,2 ∥u∗σk (xk )∥ k+1 + 1+ ∥x − xk ∥ (34) µf σk µf σk ! (2Lf,2 + 2Lg,2 )∥u∗σk+1 (xk+1 )∥ 4Lf,1 + 2Lg,1 σk − σk+1 + ∥∇2 f (xk+1 , yσ∗k+1 (xk+1 ))∥ . + µ2f µ2f σk2 Proof. Let yk∗ := yσ∗k (xk ) and Hkψ := ∇222 ψσk (xk , y k ). Let v t,k be the t-th iterate of GD(Hkψ , ∇2 ψσk (xk , y k ); γk , Tk ). As in [Dong et al., 2025, Lemma E.3], ∥y k − v t+1,k − yk∗ ∥ ≤ (1 − µf σk γk )∥y k − v t,k − yk∗ ∥ +

Lψ,2 γk k ∥y − yk∗ ∥2 . 2

After Tk iterations, ∥y k+1 − yk∗ ∥ ≤ (1 − µf σk γk )Tk ∥y k − yk∗ ∥ + ≤

Lψ,2 ∥y k − yk∗ ∥2 2µf σk

Lψ,2 1 k ∥y − yk∗ ∥ + ∥y k − yk∗ ∥2 . 4 2µf σk

(35)

Moreover, [Liu et al., 2023, Lemma D.8] gives k+1

∥yk∗ − yσ∗k+1 (xk+1 )∥ ≤

2∥∇2 f (x Lψ,1 k+1 ∥x − xk ∥ + µf σk

, yσ∗k+1 (xk+1 ))∥ σk − σk+1 µf

σk

.

(36)

Combining (35) and (36) proves (33). For the multiplier update, define u∗σk (xk , y k ) := (Hkψ )−1 ∇2 f (xk , y k ). Let wt,k be the t-th iterate
of GD Hkψ , Hkψ uk − ∇2 f (xk , y k ); γk , Tk . As in [Dong et al., 2025, Lemma E.3], ∥uk − wt+1,k − u∗σk (xk , y k )∥ ≤ (1 − µf σk γk )∥uk − wt,k − u∗σk (xk , y k )∥. Since w0,k = 0, after Tk iterations, ∥uk+1 − u∗σk (xk , y k )∥ ≤

1 1 k ∥u − u∗σk (xk )∥ + ∥u∗σk (xk , y k ) − u∗σk (xk )∥. 4 4

∥uk+1 − u∗σk+1 (xk+1 )∥ ≤

1 k 5 ∥u − u∗σk (xk )∥ + ∥u∗σk (xk , y k ) − u∗σk (xk )∥ 4 4 + ∥u∗σk (xk ) − u∗σk+1 (xk+1 )∥.

Consequently,

(37)

Using the (µf σk )-strong convexity of ψσk (xk , ·) and an argument similar to that for [Dong et al., 2025, Eq. (21)], we obtain ∥u∗σk (xk , y k ) − u∗σk (xk )∥ ≤

Lf,1 + Lψ,2 ∥u∗σk (xk )∥ k ∥y − yk∗ ∥. µf σk

(38)

Finally, [Liu et al., 2023, Lemma D.8] gives ∥u∗σk+1 (xk+1 ) − u∗σk (xk )∥   Lψ,1 Lf,1 + Lψ,2 ∥u∗σk (xk )∥ k+1 ≤ 1+ ∥x − xk ∥ µf σk µf σk +

(2Lf,2 + 2Lg,2 )∥u∗σk+1 (xk+1 )∥ µ2f

+

4Lf,1 + 2Lg,1 µ2f

Combining (37), (38), and (39) proves (34). 22

(39) ! ∥∇2 f (xk+1 , yσ∗k+1 (xk+1 ))∥

σk − σk+1 . σk2

Proof of Theorem 3.2. Let Φσ (x) := f (x, yσ∗ (x)), eky := y k − yσ∗k (xk ), and eku := uk − u∗σk (xk ). Define
Vk := ak Φσk (xk ) − f0 + bk ∥eky ∥2 + ck ∥eku ∥2 , where ak , bk , ck > 0 are non-increasing. Since Φσk (xk ) ≥ f0 , we have Vk ≥ 0. By [Liu et al., 2023, Lemma C.1], αk Φσk+1 (xk+1 ) − Φσk (xk ) ≤ − ∥∇Φσk (xk )∥2 2 ! CΦ1 ∥u∗σk (xk )∥ + CΦ2 1 1 − − ∥xk+1 − xk ∥2 + αk L2ψ,1 ∥eku ∥2 2 αk µ2f σk2 2∥∇2 f (xk+1 , yσ∗k+1 (xk+1 ))∥2 σk − σk+1
2 + αk Lψ,2 ∥u∗σk (xk )∥ + Lf,1 ∥eky ∥2 + . µf σk

(40)

Let Uk := ∥u∗σk (xk )∥ and Gk+1 := ∥∇2 f (xk+1 , yσ∗k+1 (xk+1 ))∥. Whenever ∥eky ∥ ≤ µf σk /(2Lψ,2 ), combining (40) with (33)–(34) and the monotonicity of the Lyapunov weights gives ak αk Vk+1 − Vk ≤ − ∥∇Φσk (xk )∥2 2  2 ! ak (CΦ1 Uk + CΦ2 ) + 8bk L2ψ,1 6ck (Lf,1 + Lψ,2 Uk )2 Lψ,1 ak − − 1+ ∥xk+1 − xk ∥2 − 2αk 2µ2f σk2 µ2f σk2 µf σk ! bk 75ck (Lf,1 + Lψ,2 Uk )2 2 − − ak αk (Lψ,2 Uk + Lf,1 ) − ∥eky ∥2 2 8µ2f σk2   7ck 2 − ak αk Lψ,1 ∥eku ∥2 − 8   2ak G2k+1 σk − σk+1 16bk G2k+1 σk − σk+1 2 + + µf σk µ2f σk   2
2 σk − σk+1 2 6ck Gk+1 + (2Lf,2 + 2Lg,2 )Uk+1 + 4Lf,1 + 2Lg,1 . (41) µ4f σk2 Define Πk := ∥∇Φσk (xk )∥2 + ∥xk+1 − xk ∥2 + ∥eky ∥2 + ∥eku ∥2 . We next verify the local tracking conditions and Lyapunov descent for each boundedness setting, with k0 chosen independently of K. Suppose {xk }, {y k }, and {uk } are bounded. By [Liu et al., 2023, Lemma D.5],

Case 1

∥∇2 f (xk , yσ∗k (xk ))∥ = O(σk−1 ),

∥u∗σk (xk )∥ = O(σk−2 ).

Choose σk =

σ̄ , (k + k0 )p

αk =

ᾱ , (k + k0 )4p+2τ

ak =

1 , (k + k0 )2p+τ

bk = 1,

ck =

1 . (k + k0 )6p+τ

We first verify the local condition in (41). By Assumption 4, ∥dkx ∥ ≤ Cf,0 + Lψ,1 (∥eku ∥ + Uk ). Since Uk = O((k + k0 )2p ), ∥xk+1 − xk ∥ = O((k + k0 )−4p−2τ )∥eku ∥ + O((k + k0 )−2p−2τ ). Moreover, (σk − σk+1 )/σk = O((k + k0 )−1 ) and (σk − σk+1 )/σk2 = O((k + k0 )−1+p ). Whenever ∥eky ∥ ≤ µf σk /(2Lψ,2 ), (33)–(34) give 1 k ∥e ∥ + O((k + k0 )−3p−2τ )∥eku ∥ + O((k + k0 )−p−2τ ) + O((k + k0 )−1+p ), 2 y 1 k 3p k −2τ ∥ek+1 )∥eku ∥ + O((k + k0 )2p−2τ ) u ∥ ≤ ∥eu ∥ + O((k + k0 ) )∥ey ∥ + O((k + k0 ) 4 + O((k + k0 )−1+4p ).

∥ek+1 ∥≤ y

23

(42)

(43)

We prove by induction that ∥eky ∥ ≤

µf σk , 2Lψ,2

∥eku ∥ ≤ (k + k0 )2p+τ .

(44)

Under (44), the perturbations in (42) are O((k + k0 )−p−2τ ), O((k + k0 )−p−τ ), and O((k + k0 )−1+p ), all o(σk ) for p < 1/2. Thus, for sufficiently large k0 , ∥ek+1 ∥ ≤ µf σk+1 /(2Lψ,2 ). y Similarly, (43) gives ∥ek+1 u ∥≤

1 (k + k0 )2p+τ + O((k + k0 )2p ) + O((k + k0 )2p−2τ ) + O((k + k0 )2p−τ ) + O((k + k0 )−1+4p ), 4

which is at most (k + k0 + 1)2p+τ for sufficiently large k0 and p < 1/2. The initial LL bound holds by assumption, while ∥e0u ∥ ≤ ∥u0 ∥ + O(k02p ) ≤ k02p+τ for sufficiently large k0 . Hence (44) holds for all k. Substituting these bounds into (41) and choosing k0 sufficiently large gives ak αk ak 1 7ck k 2 ∥∇Φσk (xk )∥2 − ∥xk+1 − xk ∥2 − ∥eky ∥2 − ∥e ∥ 2 4αk 4 16 u + O((k + k0 )−1−τ ) + O((k + k0 )−2+2p ) + O((k + k0 )−2−τ +2p ).

Vk+1 − Vk ≤ −

(45)

For p < 1/2 and τ > 0, the remainder terms are summable. Since ak αk = Θ((k + k0 )−6p−3τ ) is, up to constants, the smallest descent weight, telescoping (45) and using VK ≥ 0 give K−1

1 X Πk = O K k=0





1 KaK αK

  = O K −(1−6p−3τ ) .

(46)

By the boundedness of {uk } and [Liu et al., 2023, Eq. (35)],   min KKT(xk , y k , uk ) = O K −(1−6p−3τ ) + K −2p . 0≤k≤K−1

Balancing the exponents gives p = (1 − 3τ )/8, and hence an ϵ-stationary iterate within K = O(ϵ−4/(1−3τ ) ) outer iterations. Each outer iteration requires O(1) gradient/JVP computations and Tk = Θ((k + k0 )p ) HessianPK−1 vector products. Since k0 is independent of K, k=0 Tk = O(K 1+p ). With p = (1−3τ )/8, the total complexities are O(ϵ−4/(1−3τ ) ) gradient/JVP computations and O(ϵ−(9−3τ )/(2−6τ ) ) Hessian-vector products. Case 2 Suppose {∇2 f (xk , yσ∗k (xk ))} and {uk } are bounded. By the (µf σk )-strong convexity of ψσk (xk , ·), ∥∇2 f (xk , yσ∗k (xk ))∥ ∥u∗σk (xk )∥ ≤ = O(σk−1 ). µf σk Choose σk =

σ̄ , (k + k0 )p

αk =

ᾱ , (k + k0 )3p+2τ

ak =

1 , (k + k0 )p+τ

bk = 1,

ck =

1 . (k + k0 )4p+τ

By Assumption 4 and the choice of αk , ∥xk+1 − xk ∥ = O((k + k0 )−3p−2τ )∥eku ∥ + O((k + k0 )−2p−2τ ). Whenever ∥eky ∥ ≤ µf σk /(2Lψ,2 ), (33)–(34) give 1 k ∥e ∥ + O((k + k0 )−2p−2τ )∥eku ∥ + O((k + k0 )−p−2τ ) + O((k + k0 )−1 ), 2 y 1 k 2p k −2τ ∥ek+1 )∥eku ∥ u ∥ ≤ ∥eu ∥ + O((k + k0 ) )∥ey ∥ + O((k + k0 ) 4 + O((k + k0 )p−2τ ) + O((k + k0 )−1+2p ). ∥ek+1 ∥≤ y

24

(47)

(48)

We prove by induction that ∥eky ∥ ≤

µf σk , 2Lψ,2

∥eku ∥ ≤ (k + k0 )p+τ .

(49)

Under (49), the perturbations in (47) are o(σk ) for 0 < p < 1, so the LL bound holds at iteration k + 1 for sufficiently large k0 . Similarly, ∥ek+1 u ∥≤

1 (k + k0 )p+τ + O((k + k0 )p ) + O((k + k0 )p−τ ) + O((k + k0 )−1+2p ) ≤ (k + k0 + 1)p+τ . 4

The initial LL bound holds by assumption, while ∥e0u ∥ ≤ ∥u0 ∥ + O(k0p ) ≤ k0p+τ for sufficiently large k0 . Hence (49) holds for all k. Substituting these bounds into (41) and choosing k0 sufficiently large gives ak αk ak 1 7ck k 2 ∥∇Φσk (xk )∥2 − ∥xk+1 − xk ∥2 − ∥eky ∥2 − ∥e ∥ 2 4αk 4 16 u + O((k + k0 )−1−p−τ ) + O((k + k0 )−2 ) + O((k + k0 )−2−τ ).

Vk+1 − Vk ≤ −

(50)

The remainder terms are summable. Since ak αk = Θ((k + k0 )−4p−3τ ) is, up to constants, the smallest descent weight, telescoping (50) and using VK ≥ 0 give K−1

1 X Πk = O K k=0





1 KaK αK

  = O K −(1−4p−3τ ) .

(51)

By the boundedness of {uk } and [Liu et al., 2023, Eq. (35)],   min KKT(xk , y k , uk ) = O K −(1−4p−3τ ) + K −2p . 0≤k≤K−1

Balancing the exponents gives p = (1 − 3τ )/6, and hence an ϵ-stationary iterate within K = O(ϵ−3/(1−3τ ) ) outer iterations. Each outer iteration requires O(1) gradient/JVP computations and Tk = Θ((k + k0 )p ) HessianPK−1 vector products. Since k0 is independent of K, k=0 Tk = O(K 1+p ). With p = (1−3τ )/6, the total complexities are O(ϵ−3/(1−3τ ) ) gradient/JVP computations and O(ϵ−(7−3τ )/(2−6τ ) ) Hessian-vector products. Case 3

Suppose {u∗σk (xk )} and {uk } are bounded. Since ∇2 f (xk , yσ∗k (xk )) = ∇222 ψσk (xk , yσ∗k (xk ))u∗σk (xk ),

we have ∥∇2 f (xk , yσ∗k (xk ))∥ ≤ Lψ,1 ∥u∗σk (xk )∥ = O(1). Choose σk =

σ̄ , (k + k0 )p

αk =

ᾱ , (k + k0 )2p+2τ

ak =

1 , (k + k0 )τ

bk = 1,

ck =

1 . (k + k0 )2p+τ

By Assumption 4 and the choice of αk , ∥xk+1 − xk ∥ = O((k + k0 )−2p−2τ )(1 + ∥eku ∥). Whenever ∥eky ∥ ≤ µf σk /(2Lψ,2 ), (33)–(34) give 1 k ∥e ∥ + O((k + k0 )−p−2τ )(1 + ∥eku ∥) + O((k + k0 )−1 ), 2 y 1 k p k −2τ ∥ek+1 )(1 + ∥eku ∥) + O((k + k0 )−1+p ). u ∥ ≤ ∥eu ∥ + O((k + k0 ) )∥ey ∥ + O((k + k0 ) 4

∥ek+1 ∥≤ y

(52) (53)

We prove by induction that ∥eky ∥ ≤

µf σk , 2Lψ,2

∥eku ∥ ≤ (k + k0 )τ . 25

(54)

Under (54), the perturbations in (52) are o(σk ) for 0 < p < 1, so the LL bound holds at iteration k + 1 for sufficiently large k0 . Similarly, ∥ek+1 u ∥≤

1 (k + k0 )τ + O(1) + O((k + k0 )−τ ) + O((k + k0 )−1+p ) ≤ (k + k0 + 1)τ . 4

The initial LL bound holds by assumption, while ∥e0u ∥ = O(1) ≤ k0τ for sufficiently large k0 . Hence (54) holds for all k. Substituting these bounds into (41) and choosing k0 sufficiently large gives ak αk ak 1 7ck k 2 ∥∇Φσk (xk )∥2 − ∥xk+1 − xk ∥2 − ∥eky ∥2 − ∥e ∥ 2 4αk 4 16 u + O((k + k0 )−1−τ ) + O((k + k0 )−2 ) + O((k + k0 )−2−τ ).

Vk+1 − Vk ≤ −

(55)

The remainder terms are summable. Since ak αk = Θ((k + k0 )−2p−3τ ) is, up to constants, the smallest descent weight, telescoping (55) and using VK ≥ 0 give K−1

1 X Πk = O K



k=0



1 KaK αK

  = O K −(1−2p−3τ ) .

(56)

By the boundedness of {uk } and [Liu et al., 2023, Eq. (35)],   min KKT(xk , y k , uk ) = O K −(1−2p−3τ ) + K −2p . 0≤k≤K−1

Balancing the exponents gives p = (1 − 3τ )/4, and hence an ϵ-stationary iterate within K = O(ϵ−2/(1−3τ ) ) outer iterations. Each outer iteration requires O(1) gradient/JVP computations and Tk = Θ((k + k0 )p ) HessianPK−1 vector products. Since k0 is independent of K, k=0 Tk = O(K 1+p ). With p = (1−3τ )/4, the total complexities are O(ϵ−2/(1−3τ ) ) gradient/JVP computations and O(ϵ−(5−3τ )/(2−6τ ) ) Hessian-vector products.

B.3

Derivation of the GALET Baseline Complexity Bounds

We derive condition-number-dependent complexity bounds for GALET from the analysis in Appendix F of Xiao et al. [2023]. Following Xiao et al. [2023], we take N = Θ(1). Under our parameter-dependence convention, the constants Lw , LF , and Lg in their analysis satisfy Lw = O(L/µ2 ), LF = O(L/µ), and Lg = O(L2 /µ). The parameter restrictions following Eq. (42) then permit c = Θ(L3 /µ4 ) and β = Θ(1/L). In particular, the UL stepsize restriction α ≤ µ2g β/(8cL2g,1 ), with µg = µ/2, permits  2  µ β α=Θ = Θ(κ−6 ). cL2 GALET uses the Lyapunov function V k = F (xk , y k ; u∗ (xk , y k )) + cRy (xk , y k ). We first establish a uniform lower bound. By the LL error bound and ∥u∗ (x, y)∥ ≤ Lf,0 /µ, as in the proof of Theorem 3.1, 2Lf,0 F (x, y; u∗ (x, y)) ≥ f ∗ − ∥∇2 g(x, y)∥. µ By smoothness of g(x, ·), Ry ≥ ∥∇2 g∥2 /(2Lg,1 ), completing the square yields V k ≥ f∗ −

2L2f,0 Lg,1 = f ∗ − O(κ−2 ). cµ2

Thus, we may choose a lower bound V = f ∗ − C, where C = O(1) with respect to κ. Summing [Xiao et al., 2023, Eq. (42)] and using V K ≥ V gives K−1

K−1

k=0

k=0

16(V 0 − V ) 8L2g,1 X 2 1 X k 2 ∥dx ∥ ≤ + bk , K αK K 26

where bk measures the multiplier approximation error. Following Xiao et al. [2023], this error is controlled by choosing T = Θ(κ2 log(κK)) inner iterations. The first term therefore depends on the initial Lyapunov gap V 0 − V . To control this gap, we additionally impose Ry (x0 , y 0 ) ≤

µ4 , L3

which ensures cRy (x0 , y 0 ) = O(1). Moreover, by the smoothness of g and ∥u∗ (x0 , y 0 )∥ ≤ Lf,0 /µ, |F (x0 , y 0 ; u∗ (x0 , y 0 )) − f (x0 , y 0 )| ≤

Lf,0 µ

q

2Lg,1 Ry (x0 , y 0 ) = O(κ−1 ).

Together with a uniformly bounded initial objective gap and the lower bound above, this yields V 0 − V = O(1). This initialization condition is not imposed in the original GALET analysis. Without it, an initialization satisfying Ry (x0 , y 0 ) = Θ(1) contributes Θ(c) = Θ(κ3 /µ) to V 0 , potentially increasing the resulting complexity bound. As discussed in Remark 3, this condition can be satisfied by a one-time LL initialization whose cost is independent of the target accuracy. Substituting α = Θ(κ−6 ) and V 0 − V = O(1) into the descent bound, and using the residual estimates in Appendix F of Xiao et al. [2023], yields O(κ6 /ϵ) gradient/JVP computations and O(κ8 ϵ−1 log(κ/ϵ)) Hessian-vector products.

References Jiawen Bi, Jiaxiang Li, Mingyi Hong, and Shuzhong Zhang. On second-order methods for bilevel optimization. arXiv:2606.20534, 2026. 15 Nicholas Bishop, Long Tran-Thanh, and Enrico Gerding. Optimal learning from verified training data. In Advances in Neural Information Processing Systems, volume 33, pages 9520–9529, 2020. 1 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, volume 34, pages 25294–25307, 2021. 2 Tianshu Chu, Dachuan Xu, Wei Yao, and Jin Zhang. SPABA: A single-loop and probabilistic stochastic bilevel algorithm achieving optimal sample complexity. In International Conference on Machine Learning, volume 235, pages 8848–8903, 2024. 2 Mathieu Dagréou, Pierre Ablin, Samuel Vaiter, and Thomas Moreau. A framework for bilevel optimization that enables stochastic and global variance reduction algorithms. In Advances in Neural Information Processing Systems, volume 35, pages 26698–26710, 2022. 2 Mathieu Dagréou, Pierre Ablin, Samuel Vaiter, and Thomas Moreau. How to compute hessianvector products? ICLR Blogposts, 2024. URL https://iclr-blogposts.github.io/2024/ blog/bench-hvp/. 5 Youran Dong, Junfeng Yang, Wei Yao, and Jin Zhang. Efficient curvature-aware hypergradient approximation for bilevel optimization. In International Conference on Machine Learning, volume 267, pages 14231–14259, 2025. 2, 3, 10, 14, 15, 22 Sheng Fang, Yongjin Liu, Wei Yao, Chengming Yu, and Jin Zhang. qNBO: quasi-Newton meets bilevel optimization. In International Conference on Learning Representations, 2025. 2, 3, 14, 15 Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In International Conference on Machine Learning, volume 80, pages 1568–1577, 2018. 1

27

Saeed Ghadimi and Mengdi Wang. arXiv:1802.02246, 2018. 2, 3

Approximation methods for bilevel programming.

Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration complexity of hypergradient computation. In International Conference on Machine Learning, volume 119, pages 3748–3758, 2020. 2 Nadav Hallak and Nitay Suissa. An augmented lagrangian approach to bi-level optimization via a smooth equilibrium constrained problem. Journal of Optimization Theory and Applications, 208 (3):105, 2026. 2 Philip Hartman and Louis Nirenberg. On spherical image maps whose jacobians do not change sign. American Journal of Mathematics, 81(4):901–920, 1959. 16 Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A two-timescale stochastic algorithm framework for bilevel optimization: Complexity analysis and application to actor-critic. SIAM Journal on Optimization, 33(1):147–180, 2023. 2 Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen. LoRA: Low-rank adaptation of large language models. In International Conference on Learning Representations, 2022. 15 Feihu Huang. Optimal Hessian/Jacobian-free nonconvex-PL bilevel optimization. In International Conference on Machine Learning, volume 235, pages 19598–19621, 2024. 2 Kaiyi Ji, Jason D Lee, Yingbin Liang, and H. Vincent Poor. Convergence of meta-learning with task-specific adaptation over partial parameters. In Advances in Neural Information Processing Systems, volume 33, pages 11490–11500, 2020. 1 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–3023, 2022. 2 Bo Liu, Mao Ye, Stephen Wright, Peter Stone, and Qiang Liu. BOME! bilevel optimization made easy: A simple first-order approach. In Advances in Neural Information Processing Systems, volume 35, pages 17248–17262, 2022. 2 Dong C Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1):503–528, 1989. 15 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, pages 21839–21866, 2023. 2, 3, 4, 6, 8, 9, 11, 22, 23, 24, 25, 26 Akshay Mehra and Jihun Hamm. Penalty method for inversion-free deep bilevel optimization. In Asian Conference on Machine Learning, pages 347–362, 2021. 2 Rui Pan, Xiang Liu, Shizhe Diao, Renjie Pi, Jipeng Zhang, Chi Han, and Tong Zhang. LISA: Layerwise importance sampling for memory-efficient large language model fine-tuning. In Advances in Neural Information Processing Systems, volume 37, pages 57018–57049, 2024. 15 Rui Pan, Dylan Zhang, Hanning Zhang, Xingyuan Pan, Minrui Xu, Jipeng Zhang, Renjie Pi, Xiaoyu Wang, and Tong Zhang. Scalebio: Scalable bilevel optimization for llm data reweighting. In Proceedings of the 63rd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 31959–31982, 2025. 1, 12 Barak A Pearlmutter. Fast exact multiplication by the hessian. Neural Computation, 6(1):147–160, 1994. 5 Fabian Pedregosa. Hyperparameter optimization with approximate gradient. In International Conference on Machine Learning, volume 48, pages 737–746, 2016. 1, 2

28

Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine. Meta-learning with implicit gradients. In Advances in Neural Information Processing Systems, volume 32, pages 113–124, 2019. 1 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. 2, 3, 14, 15 Mengye Ren, Wenyuan Zeng, Bin Yang, and Raquel Urtasun. Learning to reweight examples for robust deep learning. In International Conference on Machine Learning, pages 4334–4343, 2018. 1 Han Shen, Quan Xiao, and Tianyi Chen. On penalty-based bilevel gradient descent method. Mathematical Programming, 214:539–589, 2025. 2 Gilbert W Stewart. On the perturbation of pseudo-inverses, projections and linear least squares problems. SIAM Review, 19(4):634–662, 1977. 17 Rohan Taori, Ishaan Gulrajani, Tianyi Zhang, Yann Dubois, Xuechen Li, Carlos Guestrin, Percy Liang, and Tatsunori B Hashimoto. Stanford alpaca: An instruction-following llama model, 2023. 12 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 Machine Learning, volume 139, pages 10708–10716, 2021. 1 Quan Xiao, Songtao Lu, and Tianyi Chen. A generalized alternating method for bilevel optimization under the polyak-lojasiewicz condition. In Advances in Neural Information Processing Systems, volume 36, pages 63847–63873, 2023. 2, 3, 4, 5, 7, 8, 9, 10, 17, 26, 27 Sheng Yang, Chengchang Liu, Lesi Chen, and John Lui. Second-order bilevel optimization with accelerated convergence rates. arXiv:2605.06431, 2026. 15 Yihua Zhang, Guanhua Zhang, Prashant Khanduri, Mingyi Hong, Shiyu Chang, and Sijia Liu. Revisiting and advancing fast adversarial training through the lens of bi-level optimization. In International Conference on Machine Learning, volume 162, pages 26693–26712, 2022. 1

29

Related documents

• A Review of Bilevel Optimization: Methods, Emerging Applications, and Recent Advancements#1117272
• To Solve Bilevel Optimization with Nonconvex Lower Levels, We Need Second-Order Stationarity#1104700
• Bi-ZOL: Bilevel Zeroth-Order Learning with Nonsmooth Responses#985422
• A Hybrid Subgradient Method for Nonsmooth Nonconvex Bilevel Optimization#160891
• A true single-level reformulation for pessimistic bilevel optimization#1142315