Penalty-based Methods for Simple Bilevel Optimization under Hölderian Error Bounds

Authors: pengyu  chen
xu  shi
rujun  jiang
jiulin  wang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 14:35:44

Abstract: This paper investigates simple bilevel optimization problems where we minimize a convex upper-level objective over the optimal solution set of a convex lower-level objective. Existing methods for such problems either only guarantee asymptotic convergence, have slow sublinear rates, or require strong assumptions. To address these challenges, we propose a penalization framework that delineates the relationship between approximate solutions of the original problem and its reformulated counterparts. This framework accommodates varying assumptions regarding smoothness and convexity, enabling the application of specific methods with different complexity results. 
Specifically, when both upper- and lower-level objectives are composite convex functions, under an $\alpha$-Hölderian error bound condition and certain mild assumptions, our algorithm attains an $(\epsilon,\epsilon^{\beta})$-optimal solution of the original problem for any $\beta> 0$ within $\mathcal{O}\left(\sqrt{{1}/{\epsilon^{\max\\{\alpha,\beta\\}}}}\right)$ iterations. The result can be improved further if the smooth part of the upper-level objective is strongly convex. We also establish complexity results when the upper- and lower-level objectives are general nonsmooth functions. Numerical experiments demonstrate the effectiveness of our algorithms.

Main paper: 
 introduction: 
 Bilevel optimization involves embedding one optimization problem within another, creating a hierarchical structure where the upper-level problem's feasible set is influenced by the lower-level problem. This framework frequently occurs in various real-world scenarios, such as meta-learning [Bertinetto et al., 2018 , Rajeswaran et al., 2019] , hyper-parameter optimization [Chen et al., 2024 , Franceschi et al., 2018 , Shaban et al., 2019] , reinforcement learning [Mingyi et al., 2020] and adversarial learning [Bishop et al., 2020 , Wang et al., 2021 , 2022] . In this paper, we concentrate on a subset of bilevel optimization known as simple bilevel optimization (SBO), which has garnered significant interest in the machine learning community due to its relevance in dictionary learning [Beck and Sabach, 2014, Jiang et al., 2023] , lexicographic optimization [Kissel et al., 2020 , Gong et al., 2021] , lifelong learning [Malitsky, 2017 , Jiang et al., 2023] ; see more details in Appendix A.SBO aims to find an optimal solution that minimizes the upper-level objective over the solution set of the lower-level problem. In other words, we are interested in solving the following problem:min x∈R n F (x) s.t. x ∈ arg min z∈R n(P)Here F, G : R n → R ∪ {∞} are proper, convex, and lower semi-continuous functions. We also assume that the optimal solution set of the lower-level problem, denoted as X opt , is nonempty. Moreover, since G is convex and lower semi-continuous, it holds that X opt is closed and convex [Bertsekas et al., 2003, Proposition 1.2.2 and Page 49] .In this paper, we first reformulate problem (P) into the constrained form:min x∈R n F (x) s.t. G(x) -G * ≤ 0, (P Val )where G * represents the optimal value of the unconstrained lower-level problem.Based on this reformulation, we consider the following penalization of (P Val ),min x∈R n Φ γ (x) = F (x) + γp(x), (P γ )where p(x) := G(x) -G * is the so-called residual function and γ > 0 is the penalized parameter. Obviously, we have p(x) ≥ 0, and p(x) = 0 if and only if x ∈ X opt .Denote F * and G * as the optimal values of problem (P) and the lower-level problem min x∈R n G(x), respectively. We aim to find an (ϵ F , ϵ G )-optimal solution x * of problem (P), which satisfiesF (x * ) -F * ≤ ϵ F , G(x * ) -G * ≤ ϵ G . (1)Moreover, a point x * γ is said to be an ϵ-optimal solution of problem (P γ ) if Φ γ (x * γ ) -Φ * γ ≤ ϵ, where Φ * γ is the optimal value of problem (P γ ).Various approaches have been developed to solve problem (P) [Cabot, 2005 , Solodov, 2007 , Sabach and Shtern, 2017 , Dutta and Pandit, 2020 , Gong et al., 2021] . Among those, one category that is the most related to penalization formulation (P γ ) is the regularization method, which integrates the upper-and lower-level objectives through Tikhonov regularization [Tikhonov and Arsenin, 1977] min x∈R n η(x) := σF (x) + G(x), (P Reg )where σ is the so-called regularization parameter. When F is strongly convex and its domain is compact, Amini and Yousefian [2019] extended the IR-PG method from Solodov [2007] , which achieved a asymptotic convergence rate for the upper-level problem and a convergence rate of O ( 1/K 0.5-b ) for the lower-level problem, where b ∈ (0, 0.5). Malitsky [2017] studied a version of Tseng's accelerated gradient method and showed a convergence rate of O (1/K) for the lower-level problem, while the convergence rate for the upper-level objective is not explicitly provided. Kaushik and Yousefian [2021] proposed an iteratively regularized gradient (a-IRG) method which obtains a complexity of O ( 1/K 0.5-b ) and O ( 1/K b ) for the upper-and lower-level objective, respectively, where b ∈ (0, 0.5). Inspired by this research, and under a quasi-Lipschitz assumption for F , Merchav and Sabach [2023] introduced a bi-subgradient (Bi-SG) method. This method demonstrates convergence rates of O(1/K b ) and O(1/K 1-b ) for the lower-and upper-level objectives, respectively, where b ∈ (0.5, 1). In their framework, the convergence rate of the upper-level objective can be improved to be linear when F is strongly convex. Recently, under the weak-sharp minima assumption of the lower-level problem, Samadi et al. [2023] proposed a regularized accelerated proximal method (R-APM), showing a convergence rate of O(1/K 2 ) for both upper-and lower-level objectives. When the domain is compact and F, G are both smooth, Giang-Tran et al. [2023] proposed an iteratively regularized conditional gradient (IR-CG) method, which ensures convergence rates of O(1/K p ) and O(1/K 1-p ) for upper-and lower-level objectives, respectively, where p ∈ (0, 1).Despite the abundance of existing methodologies yielding non-asymptotic convergence outcomes, their efficacy is frequently contingent upon additional assumptions. Denote L f1 and L g1 as the Lipschitz constants for the gradients of the smooth components in the upper-and lower-level objectives, respectively. Specifically, when F is strongly convex and G is smooth, Beck and Sabach [2014] presented the Minimal Norm Gradient (MNG) method and provided the asymptotic convergence to the optimal solution set and a convergence rate of O ( L 2 g1 /ϵ 2 ) for the lower-level problem. When F is assumed to be smooth, Jiang et al. [2023] introduced a conditional gradient-based bilevel optimization (CG-BiO) method, which invokes at most O (max{L f1 /ϵ F , L g1 /ϵ G }) of linear optimization oracles to achieve an (ϵ F , ϵ G )-optimal solution. Shen et al. [2023] combined an online framework with the mirror descent algorithm and established a convergence rate of O(1/ϵ 3 ) for both upper-and lower-level objectives, assuming a compact domain and boundedness of the functions and gradients at both levels. Furthermore, they showed that the convergence rate can be improved to O(1/ϵ 2 ) under additional structural assumptions. For a concise overview of overall methodologies, including their assumptions and convergence outcomes, refer to Table 1 in Appendix B.For general bilevel optimization problems, there have been recent results on convergent guarantees [Shen and Chen, 2023 , Sow et al., 2022 , Chen et al., 2023 , Huang, 2023] . Among those, the one that is the most related to ours is [Shen and Chen, 2023] . It investigates the case when the upper-level objective is nonconvex and gives convergence results under additional assumptions [Shen and Chen, 2023, Theorem 3 and 4] . However, as the general bilevel optimization problem is nonconvex, the algorithms in the literature often converge to weak stationary points, while our method for SBO converges to global optimal solution.Our approach is straightforward. Firstly, we introduce a penalization framework delineating the connection between approximate solutions of problems (P) and (P γ ). This framework enables the attainment of an (ϵ F , ϵ G )-optimal solution by solving problem (P γ ) approximately. Subsequently, our focus shifts solely to resolving the unconstrained problem (P γ ). Depending on varying assumptions regarding smoothness and convexity, we can employ different methods such as the accelerated proximal gradient (APG) methods [Beck and Teboulle, 2009 , Nesterov, 2013 , Lin and Xiao, 2014] to solve problem (P γ ). We summarize our main contributions as follows.• We propose a framework that explicitly examines the relationship between an ϵ-optimal solution of penalty formulation (P γ ) and an (ϵ F , ϵ G )-optimal solution of problem (P). We also provide a lower bound for the metric F (x) -F * . • When F and G are both composite convex functions, we provide a penalty-based APG algorithm that attains an (ϵ, ϵ β )-optimal solution of problem (P) within O( √ 1/ϵ max{α,β} ) iterations. If the upper-level objective is strongly convex, the complexity can be improved to O( √ 1/ϵ max{α-1,β-1} log 1 ϵ ). We also apply our method for the scenario where both the upper-and lower-level objectives are generalized nonsmooth convex functions.• We present adaptive versions of PB-APG and PB-APG-sc with warm-start, which dynamically adjust the penalty parameters, and solve the associated penalized problem with adaptive accuracy. The adaptive ones have similar complexity results as their primal counterparts but can achieve superior performance in some experiments.Utilizing the penalization method to address the original SBO problem is a novel approach. While Tikhonov regularization may seem similar to our framework, its principles differ. Implementing Tikhonov regularization necessitates the "slow condition" (lim k→∞ σ k = 0, ∑ ∞ k=0 σ k = +∞), which requires iterative solutions for each iteration. In contrast, our method simply involves solving a single optimization problem (P γ ) for a given γ. Furthermore, we establish a relationship between the approximate solutions of the original bilevel problem and those of the reformulated single-level problem (P γ ) for a specific γ. This is the first theoretical result connecting the original bilevel problem to the penalization problem, accompanied by an optimal non-asymptotic complexity result.we replace the assumption with the condition that the upper-level objective is Lipschitz continuous (cf. Theorems 2.7 and 2.8). Assumption 2.2 (Hölderian error bound). The function p(x) := G(x) -G * satisfies the Hölderian error bound with exponent α ≥ 1 and ρ > 0. Namely, dist(x, Xopt) α ≤ ρp(x), ∀x ∈ dom(G), where dist(x, Xopt) := infy∈X opt ∥x -y∥.We remark that Hölderian error bounds are satisfied by many practical problems and widely used in optimization literature [Pang, 1997 , Bolte et al., 2017 , Zhou and So, 2017 , Roulet and d'Aspremont, 2020 , Jiang and Li, 2022] . There are two notable special cases: (i) when α = 1, we often refer to X opt as a set of weak sharp minima of G [Burke and Ferris, 1993 , Studniarski and Ward, 1999 , Burke and Deng, 2005 , Samadi et al., 2023] ; (ii) when α = 2, Assumption 2.2 is known as the quadratic growth condition [Drusvyatskiy and Lewis, 2018a] . Additional examples of functions exhibiting Hölderian error bound, along with their corresponding parameters, are presented in Appendix C.We are now ready to establish the connection between approximate solutions of problems (P) and (P γ ). The subsequent two lemmas build upon the work of Shen and Chen [2023] for (general) bilevel optimization. Compared with their work, we generalize the exponent α from 2 to α ≥ 1, providing a more general result. Furthermore, we also derive a lower bound for the penalized parameter for all α ≥ 1 and present a theoretical framework for these scenarios. Lemma 2.3. Suppose that Assumptions 2.1 and 2.2 hold with α > 1. Then, for any ϵ > 0, an optimal solution of problem (P) is an ϵ-optimal solution of problem(P γ ) when γ ≥ ρl α F (α -1) α-1 α -α ϵ 1-α .Lemma 2.3 establishes the relationship between an optimal solution of problem (P) and an ϵ-optimal solution of problem (P γ ) when α > 1. It also provides a lower bound for γ, which plays a pivotal role in the complexity results. The proofs of this paper are deferred to Appendix E.The lemma presented below yields a more favorable outcome when α = 1, which is referred to as exact penalization. Notably, this specific result is not discussed in Shen and Chen [2023] . Lemma 2.4. Suppose that Assumptions 2.1 and 2.2 hold with α = 1. Then an optimal solution of problem (P) is also an optimal solution of problem (P γ ) if γ ≥ ρl F , and vice versa if γ > ρl F . In this case, we say that there is an exact penalization between problems (P) and (P γ ).For simplicity, we defineEQUATION)Based on Lemmas 2.3 and 2.4, we give an overall relationship of approximate solutions between problems (P γ ) and (P). Theorem 2.5. Suppose that Assumptions 2.1 and 2.2 hold. For any given ϵ > 0 and β > 0, letγ = γ * + { 2l β F ϵ 1-β if α > 1, l β F ϵ 1-β if α = 1, with γ * defined in (2). If x * γ is an ϵ-optimal solution of problem (P γ ), then x * γ is an (ϵ, l -β F ϵ β )- optimal solution of problem (P).Particularly, we are also able to establish a lower bound for F (x * γ ) -F * under the same conditions outlined in Theorem 2.5. Theorem 2.6. Suppose that the conditions in Theorem 2.5 hold. Then, x * γ satisfies the following suboptimality lower bound,F (x * γ ) -F * ≥ -l F (ρl -β F ϵ β ) 1 α .By setting β = α, we obtain F (x * γ ) -F * ≥ -ρ 1 α ϵ. which along with Theorem 2.5 gives|F (x * γ ) -F * | ≤ max{ϵ, ρ 1 α ϵ}. 
analysis of non-convex upper-level: 
 Note that the upper-level objective F is required to be convex in the above context (cf. Theorem 2.5). This raises a question: while Theorem 2.5 establishes the relationship between approximate solutions of problems (P) and (P γ ), the distinction between the global or local optimal solutions of problem (P) and (P γ ) remains unclear when F is non-convex.We first establish the relationship between global optimal solutions of problems (P) and (P γ ) when F is non-convex, which is similar to Theorem 2.5. Theorem 2.7. Suppose that Assumption 2.2 holds, G is convex, and F is l-Lipschitz continuous on dom(F ). For any given ϵ > 0 and β > 0, letEQUATIONwhere γ * is given by (2). If x * γ is an ϵ-global optimal solution of problem (P γ ), then x * γ is an (ϵ, l -β ϵ β )-global optimal solution of problem (P).Theorem 2.7 provides the relationship between the global optimal solutions of problems (P γ ) and (P). However, the relationship between local optimal solutions of these problems is more intricate than those of the global ones [Shen and Chen, 2023] . Given r > 0 and z ∈ R n , define B(z, r) := {x ∈ R n : ∥x -z∥ ≤ r}. We present the following theorem, which demonstrates that local optimal solutions of problem (P γ ) can serve as approximate local optimal solutions of problem (P).Theorem 2.8. Suppose that Assumption 2.2 holds and G is convex. Let x * γ be a local optimal solution of problem(P γ ) on B(x * γ , r). Assume F is l-Lipschitz continuous on B(x * γ , r). Then x * γ is an approximate local optimal solution of problem (P) that satisfies F (x * γ )-F * B ≤ 0 and G(x * γ )-G * ≤ ϵ when α > 1 and γ ≥ ( ρl α ϵ α-1 ) 1 α , where F * B is the optimal value of problem (P) on B(x * γ , r) ∩ X opt . Furthermore, x *γ is a local optimal solution of problem (P) when α = 1 and γ > ρl.Indeed, the relationship between approximate local optimal solutions of problems (P γ ) and (P) is more intricate than the connection among global solutions presented in Theorem 2.5. These interactions will be the focus of our future work. The proofs of Theorems 2.7 and 2.8 are presented in Appendixes E.5 and E.6. 
main algorithms: 
 In this section, we concentrate on addressing problem (P), making various assumptions, and offering distinct convergence outcomes.In this scenario, we address problem (P) where F and G are both composite functions, i.e.,F = f 1 + f 2 and G = g 1 + g 2 .Assumption 3.1. F and G satisfy the following assumptions.(1) The gradient of f 1 (x), denoted as ∇f 1 , is L f1 -Lipschitz continuous on dom(F );(2) The gradient of g 1 (x), denoted as ∇g 1 , is L g1 -Lipschitz continuous on dom(G);(3) f 2 and g 2 are proper, convex, lower semicontinuous, and possibly non-smooth.We remark that Assumption 3.1(1)(3) is more general than many existing papers in the literature. Specifically, while previous works such as Beck and Sabach [2014] , Amini and Yousefian [2019] , Jiang et al. [2023] , Giang-Tran et al. [2023] require the upper-level objective to be smooth or strongly convex, we simply assume that F is a composite function composed of a smooth convex function and a possibly non-smooth convex function. For the lower-level objective, previous works such as Beck and Sabach [2014] , Amini and Yousefian [2019] , Jiang et al. [2023] , Giang-Tran et al.[2023] impose smoothness assumptions and, in some cases, convexity and compactness constraints on the domain; while our approach does not require these additional constraints, allowing for more flexibility and generality as presented in Assumption 3.1(2)(3).We are now prepared to introduce two algorithms: the penalty-based accelerated proximal gradient (PB-APG) algorithm and its adaptive counterpart, the aPB-APG to solve problem (P γ ) and, subsequently, to obtain an (ϵ F , ϵ G )-optimal solution of problem (P).To simplify notations, we omit the constant term -γG * , and rewrite problem (P γ ) as follows,min x∈R n Φ γ (x) := ϕ γ (x) + ψ γ (x), (P Φ )where ϕ γ (x) = f 1 (x) + γg 1 (x) and ψ γ (x) = f 2 (x) + γg 2 (x) represent the smooth and nonsmooth parts, respectively. Then, it follows that the gradient of ϕ γ (x) is L γ -Lipschitz continuous with L γ = L f1 + γL g1 .To implement the APG methods, we need another assumption concerning ψ γ (x). Assumption 3.2. For any γ > 0, the function ψ γ (x) is prox-friendly, i.e., the proximal mapping prox tψγ (y) := arg minx∈R n {ψ γ (x) + 1 2t ∥x -y∥ 2 },is easy to compute for any t > 0.The function ψ γ (x) represents the sum of two non-smooth functions, and proximal mapping for such function sums is widely studied and used in the literature [Yu, 2013 , Pustelnik and Condat, 2017 , Adly et al., 2019 , Boob et al., 2023 , Latafat et al., 2023] . This assumption is also a more general requirement compared to many existing algorithms [Sabach and Shtern, 2017, Giang-Tran et al., 2023] . It is important to note that our assumption is more general than existing literature.In the simple bilevel literature, when employing proximal mappings, researchers often consider the scenario where only one level contains a nonsmooth term (see, e.g., [Jiang et al., 2023 , Doron and Shtern, 2023 , Samadi et al., 2023 , Merchav and Sabach, 2023] ). In this case, the proximal mapping of the sum f 2 + γg 2 is then reduced to the proximal mapping of either f 2 or g 2 , which is a more easily satisfied condition.We apply the APG algorithm [Beck and Teboulle, 2009 , Lin and Xiao, 2014 , Nesterov, 2013] to solve problem (P Φ ), as outlined in Algorithm 1. Moreover, if the Lipschitz constant L γ is unknown or computationally infeasible, line search [Beck and Teboulle, 2009 ] can be adopted and will yield almost the same complexity bound. For brevity, we denote Algorithm 1 as x = PB-APG(ϕ γ , ψ γ , L f1 , L g1 , x 0 , ϵ), where x represents an ϵ-optimal solution of (P Φ ).Algorithm 1 Penalty-based APG (PB-APG) 1:Input: γ, L γ = L f1 + γL g1 , x -1 = x 0 ∈ R n , R > 0, t -1 = t 0 = 1, k = 0, ϵ > 0 and {t k }. 2: for k ≥ 0 do 3: y k = x k + t k ( t -1 k-1 -1 ) (x k -x k-1 )4:x k+1 = prox L -1 γ ψγ (y k -L -1 γ ∇ϕ γ (y k )) 5: end forIn Algorithm 1, we stop the loop of Line. 3 -4 if the number of iterations satisfies that:2(L f + γLg)R 2 (k + 1) 2 ≤ ϵ,where R is a constant that satisfies ∥x0 -x * ∥ ≤ R.Combining Theorem 2.5 and [Tseng, 2008, Corollary 2] , we establish the following complexity result for problem (P). Theorem 3.3. Suppose that Assumptions 2.1, 2.2, 3.1 and 3.2 hold and the sequence {t k } in Algorithm 1 satisfies1-t k+1 t 2 k+1 ≤ 1 t 2 k. Let γ be given as in Theorem 2.5. Algorithm 1 generates an (ϵ, l -β F ϵ β )-optimal solution of problem (P) after at most K iterations, whereK = O   √ L f 1 ϵ + √ l max{α,β} F Lg 1 ϵ max{α,β}   .Note that Theorem 3.3 encompasses all possible relationships between the magnitudes of ϵF and ϵG in (1), as α ≥ 1 and β > 0 are arbitrary. Specially, if α = 1 and β ≤ α, the number of iterations isK = O ( √ (L f 1 + lF Lg 1 )/ϵ) . This result matches the lower bound complexity for unconstrained smooth or convex composite optimization [Nemirovsky and Yudin, 1983, Woodworth and Srebro, 2016] . Additionally, if g1 ≡ 0, the number of iterations for obtaining an (ϵ, ϵ β )-optimal solution of problems (P) is independent of γ, which can be improved toK = O( √ L f 1 /ϵ).Remark 3.4. It is noteworthy that Theorem 1 in a previous paper Samadi et al. [2023] provides the first method that needs O( √ (Lg 1 + lF Lg 1 )/ϵ) iterations to achieve an (ϵ, ϵ) solution if α = 1 and F is smooth. Nevertheless, our methodology diverges in various respects. First, our approach is rooted in the penalization formulation of problem (P Val ), while the approach proposed by Samadi et al. [2023] is based on the Tikhonov regularization [Tikhonov and Arsenin, 1977] . Second, we provide a theoretical framework that clearly delineates the relationship between approximate solutions of problems (P) and (Pγ) for all cases of α ≥ 1 and F is non-convex, as indicated in Lemmas 2.3, 2.4 and Theorems 2.5, 2.7, 2.8. Therefore, we can first shift our focus from (P) to (Pγ) based on the penalization framework and then use various methods to solve (Pγ), not limited to using the APG methods. Besides, the association between approximate solutions of problem (P) and (Pγ) differs significantly based on whether α > 1 or α = 1. For the case of α > 1, the lower bound comprehensively integrates the accuracy parameter ϵ, which results in a more sophisticated analysis of the convergence result, while Samadi et al. [2023] did not consider the situation when α > 1. Third, our method applies to the case that F is composite, while Samadi et al. [2023] requires F to be smooth. Finally, we also propose an adaptive version of our algorithm (see Algorithm 2) that does not require an estimate of γ.In practice, the penalty parameter γ might be difficult to determine. This motivates us to propose Algorithm 2, which adaptively updates γ and invokes PB-APG with dynamic γ and solution accuracies.Algorithm 2 Adaptive PB-APG method (aPB-APG) 1: Input:x 0 ∈ R n , γ 0 = γ 1 > 0, L f1 , L g1 , ν > 1, η > 1, ϵ 0 > 0. 2: for k ≥ 0 do 3: ϕ k (x) = f 1 (x) + γ k g 1 (x) 4: ψ k (x) = f 2 (x) + γ k g 2 (x) 5: Invoke x k = PB-APG(ϕ k , ψ k , L f1 , L g1 , x k-1 , ϵ k ) 6: ϵ k+1 = ϵ k /η 7:γ k+1 = νγ k 8: end forIn Algorithm 2, we adaptively update the penalty parameter γ k , and invoke the PB-APG to generate an approximate solution for (Pγ) with accuracy ϵ = ϵ k . Meanwhile, a warm-start mechanism is employed, meaning that the initial point for each subproblem is the output of the preceding subproblem. The convergence result of Algorithm 2 is as follows.Theorem 3.5. Suppose that Assumptions 2.1, 2.2, 3.1, and 3.2 hold. Also assume that for every outcome of inner loop in Algorithm 2, ∥x k -x * k ∥ ≤ R. Let ϵ0 > 0 be given.• When α > 1, set ν > η α-1 , and define -1) .N := ⌈log η 1-α ν (ρL α F (α -1) α-1 α -α ϵ 1-α 0 /γ0)⌉+ and γ * k := ρL α F (α -1) α-1 α -α ϵ 1-α 0 η k(α• When α = 1, set ν > 1, and define N := ⌈log ν (ρlF /γ0)⌉+ and γ * k := ρLF .Then, for any k ≥ N , Algorithm 2 generates an( ϵ 0 η k , 2ϵ 0 η k (γ 0 ν k -γ * k ))-optimal solution of problem (P) after at most K iterations, where K satisfiesK = O ( √ L f 1 η k ϵ0 + √ Lg 1 γ0(ην) k ϵ0) .Theorem 3.5 shows that for any given initial accuracy ϵ0 > 0, Algorithm 2 can produce an approximate solution of problem (P) with the desired accuracy.Remark 3.6. From Theorem 3.5, one can obtain an (ϵ,ϵ γ 0 ν k -γ * k )-optimal solution of problem (P) within O( √ L f 1 /ϵ + √ Lg 1 /ϵ α ) iterations when ϵ/η ≤ ϵ0/η k ≤ ϵ,which is similar to the complexity results in Theorem 3.3.We investigate the convergence outcomes when the smooth part of the upper-level objective exhibits strong convexity.Assumption 3.7. f1(x) is µ-strongly convex on dom(F ) with µ > 0.Assumption 3.7 is another widely adopted setting in the existing SBO literature [Beck and Sabach, 2014 , Sabach and Shtern, 2017 , Amini and Yousefian, 2019 , Merchav and Sabach, 2023] . Here, we propose a variant of PB-APG that can provide better complexity results than existing methods. Our main integration is an APG-based algorithm, which has been studied in the existing literature [Nesterov, 2013 , Lin and Xiao, 2014 , Xu, 2022] . In this paper, we adopt the algorithm proposed in Lin and Xiao [2014] and modify it with a constant step-size for simplicity as in Algorithm 3. Similar to Algorithm 1, we denote Algorithm 3 byx = PB-APG-sc(ϕγ, ψγ, µ, L f 1 , Lg 1 , y0, ϵ). Algorithm 3 PB-APG method for Strong Convexity Case (PB-APG-sc) 1: Input: µ, γ, L γ = L f1 + γL g1 , x -1 , y 0 ∈ R n . 2: ỹ = y 0 -L -1 γ ∇ϕ γ (x -1 ) 3: x = prox L -1 γ ψγ (ỹ -L -1 γ ∇ϕ γ (ỹ)) 4: Initialization: Let x -1 = x 0 = x, k = 0 5: for k ≥ 0 do 6: y k = x k + √ Lγ - √ µ √ Lγ + √ µ (x k -x k-1 )7:x k+1 = prox L -1 γ ψγ (y k -L -1 γ ∇ϕ γ (y k )) 8: end forThe convergence analysis of Algorithm 3 is in the existing literature [Nesterov, 2013, Lin and Xiao, 2014] .Combining [Lin and Xiao, 2014 , Theorem 1] and Theorem 2.5, we have the following complexity result.Theorem 3.8. Suppose that Assumptions 2. 1, 2.2, 3.1, 3.2, and 3.7 hold. Algorithm 3 can produce an (ϵ, l -β F ϵ β )-optimal solution of problem (P) after at most K iterations, where K satisfiesK = O   √ L f 1 µ log 1 ϵ + √ l max{α,β} F Lg 1 ϵ max{α-1,β-1} log 1 ϵ   .Theorem 3.8 improves the complexity results of Theorem 3.3 significantly. Specifically, when 0 < β ≤ α = 1, the convergence rate can be improved to be linear, i.e., K = O( √ L f 1 /µ log 1 ϵ ). Additionally, we present an adaptive variant of PB-APG-sc, termed aPB-APG-sc, which adaptively executesx k = PB-APG-sc(ϕ k , ψ k , µ, L f 1 , Lg 1 , x k-1 , ϵ k )and enjoys the similar complexity results of Algorithm 3, as delineated in Algorithm 4 within Appendix D.1.In this section, we focus on the scenario where both the upper-and lower-level objectives are non-smooth, namely, f1 = g1 ≡ 0. Additionally, we assume that there is a point x ∈ C in the lower level problem, where C is either R n (the unconstrained case) or a nonempty closed and convex set satisfyingC ⊆ int (dom(F ) ∩ dom(G)).It is worth noting that in the case where both F and G are non-smooth, the convergence result may not be as favorable as those in the previous scenarios. This is primarily due to the limited availability of information and unfavorable properties concerning F and G. In this case, we employ a subgradient method to solve problem (Pγ), which has been extensively studied in the existing literature [Shor, 2012 , Bubeck et al., 2015 , Beck, 2017 , Nesterov, 2018] . Specifically, we updateEQUATIONwhereξ k ∈ ∂Φγ(x k ) is an subgradient of Φγ(x k ), and Proj C (x) is the projection of x onto C.Let x * γ be an optimal solution of problem (Pγ) and suppose that there exists a constant R such that ∥x0 -x * γ ∥ ≤ R. Motivated by Theorem 8.28 in Beck [2017] , we establish the subsequent complexity result for problem (P).Theorem 3.9. Suppose that Assumption 3.1(3) holds, f2 and g2 are l f 2 -and lg 2 -Lipschitz continuous, respectively. Set step-sizeη k = R lγ √k+1 in (4). Then, the subgradient method produces an (ϵ, l -β f 2 ϵ β )-optimal solution of problem (P) after at most K iterations, where K satisfiesK = O ( l 2 f 2 ϵ 2 + l max{2α,2β} f 2 l 2 g 2 ϵ max{2α,2β}) .For non-smooth SBO problems, our method has lower complexity compared to existing approaches. Specifically, under a bounded domain assumption, Helou and Simões [2017] simply proposed an ϵ-subgradient method with an asymptotic rate towards the optimal solution set. The a-IRG method in Kaushik and Yousefian [2021] achieved convergence rates of O(1/ϵ1 0.5-b ) and O(1/ϵ 1 b) for the upper-and lower-level objectives, respectively, where b ∈ (0, 0.5). Setting b = 0.25 yields the convergence rates of O(1/ϵ 4 ) for both upper-and lower-level objectives, which indicates that our complexity is more efficient than theirs when α < 2 and β ≤ α. Furthermore, the online framework proposed in Shen et al. [2023] performed a complexity of O(1/ϵ3 ) for both upper-and lower-level objectives. Similarly, our approach prevails over theirs when α < 1.5 and β ≤ α.Strongly convex upper-level objective. Based on Theorem 8.31 in Beck [2017] , we next explore the improved complexity result for problem (P) when f2 is additionally strongly convex.Theorem 3.10. Suppose that Assumption 3.1(3) holds, C ⊆ int (dom(F ) ∩ dom(G)), f2 is l f 2 -Lipschitz continuous and µ f 2 -strongly convex 3 , and g2 is lg 2 -Lipschitz continuous. Choose step-size η k = 2 µ f 2 (k+1)in (4). Then, the subgradient method produces an (ϵ, l -β f 2 ϵ β )-optimal solution of problem (P) after at most K iterations, where K satisfiesK = O ( l 2 f 2 µ f 2 ϵ + l max{2α,2β} f 2 l 2 g 2 µ f 2 ϵ max{2α-1,2β-1}) .To our knowledge, within the context of Theorem 3.10, current findings fail to exploit strong convexity to enhance results. However, our approach capitalizes on distinct structural characteristics that yield superior complexity outcomes relative to Theorem 3.9 in cases where α < 2 and β ≤ α. 
numerical experiments: 
 We apply our Algorithms 1, 2, 3 and 4 to two simple bilevel optimization problems from the motivating examples in Appendix A. The performances of our methods are compared with several existing methods: MNG [Beck and Sabach, 2014] , BiG-SAM [Sabach and Shtern, 2017] , DBGD [Gong et al., 2021] , a-IRG [Kaushik and Yousefian, 2021], CG-BiO [Jiang et al., 2023] , Bi-SG [Merchav and Sabach, 2023] and R-APM [Samadi et al., 2023] . For practical efficiency, we use the Greedy FISTA algorithm proposed in Liang et al. [2022] as the APG method in our approach. Detailed settings and additional experimental results are presented in Appendix F. The LRP readsEQUATIONwhere IC (x) is the indicator function of the set C = {x ∈ R n : ∥x∥1 ≤ θ} with θ = 10. Our goal is to find a solution to the lower-level problem with the smallest Euclidean norm. The upper-level objective only consists of the smooth part, which is 1-strongly convex and 1-smooth; meanwhile, the lower-level objective is a composite function, where the smooth part is 1 4m λmax(A T A)-smooth, and the nonsmooth part is prox-friendly [Duchi et al., 2008] .In this experiment, we compare our methods with MNG, BiG-SAM, DBGD, a-IRG, CG-BiO, and Bi-SG. We plot the values of residuals of the lower-level objective G(x k ) -G * and the upper-level objective over time in Figure 1 .As shown in Figure 1 , the PB-APG, aPB-APG, PB-APG-sc, and aPB-APG-sc algorithms exhibit significantly faster convergence performance than the other methods for both lower-and upper-level objectives, although R-APM attains similar outcomes, our PB-APG and PB-APG-sc ensure a more rapid decline than it, as shown in the first subfigure of Figure 1 . This is because our methods achieve lower optimal gaps and desired function values of the lower-and upper-level objectives with less execution time. This observation confirms the improved complexity results stated in the theorems above. Although the high exactness of our methods for the lower-level problem leads to larger upper-level objectives, Table 3 in Appendix F.1 shows that our methods are much closer to the optimal value. This is reasonable because the other methods exhibit lower accuracy at the lower-level problem, resulting in larger feasible sets compared to the lower-level optimal solution set Xopt. In addition, Figure 1 demonstrates that aPB-APG and aPB-APG-sc outperform PB-APG and PB-APG-sc in terms of convergence rate. This improvement can be attributed to the adaptiveness incorporated in Algorithms 2 and 4.The LSRP has the following form:EQUATION)where τ = 0.02 regulates the trade-off between ℓ1 and ℓ2 norms. We aim to find a sparse solution for the lowerlevel problem. The upper-level objective is formulated as a composite function, which consists of a τ -strongly convex and τ -smooth component, along with a proximal-friendly non-smooth component [Beck, 2017] . The lower-level objective is a smooth function with a smoothness parameter of 1 m λmax(A T A). In this experiment, we compare the performances of our methods with a-IRG, BiG-SAM, and Bi-SG. We plot the values of residuals of lower-level objective G(x k ) -G * and the upper-level objective over time in Figure 2 .Figure 2 shows that the proposed PB-APG, aPB-APG, PB-APG-sc, and aPB-APG-sc converge faster than the compared methods for both the lower-and upper-level objectives, as well. For the upper-level objective, our methods achieve larger function values than other methods, except BiG-SAM (δ = 0.01). This is because our methods attain higher accuracy for the lower-level objective than other methods. We have similar observations in Section 4.1. Furthermore, Figure 1 also demonstrates that the adaptive mechanism produces staircase-shaped curves for aPB-APG and aPB-APG-sc, which might prevent undesirable fluctuations in PB-APG and PB-APGsc.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: Please refer to Section 4 and Appendix F. 
conclusion: 
 This paper proposes a penalization framework that effectively addresses the challenges inherent in simple bilevel optimization problems. By delineating the relationship between approximate solutions of the original problem and its penalized reformulation, we enable the application of specific methods under varying assumptions for the original problem. Under the Hölderian error bound condition, our methods achieve superior complexity results compared to the existing methods. The performance is further improved when the smooth component of the upper-level objective is strongly convex. Additionally, we extend our framework to scenarios involving general nonsmooth objectives. Numerical experiments also validate the effectiveness of our algorithms.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: Please refer to Section 4 and Appendix F and the supplemental material.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: Please refer to Section 4 and Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [No] Justification: The error bars are not applicable in this paper.Guidelines:• The answer NA means that the paper does not include experiments. • The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: Please refer to Section 4 and Appendix F.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: Please refer to Section 4 and Appendix F.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Many machine learning applications involve a primary objective G, which usually represents the training loss, and a secondary objective F , which can be a regularization term or an auxiliary loss. A common approach for such problems is to optimize G fully and then use F to select the optimal solutions from the ones obtained for G. This is called lexicographic optimization [Kissel et al., 2020 , Gong et al., 2021] . Two classes of lexicographic optimization problems are the regularized problem, also known as the ill-posed optimization problem [Amini and Yousefian, 2019, Jiang et al., 2023] , and the over-parameterized regression [Jiang et al., 2023] , where the upper-level objectives are the regularization terms or loss functions, and the lower-level objectives are the loss functions and the constraint terms. We present some examples of these classes of problems as follows.Example A.1 (Linear Inverse Problems). Linear inverse problems aim to reconstruct a vector x ∈ R n from measurements b ∈ R m that satisfy b = Ax + ρε, where A : R n → R m is a linear mapping, ε ∈ R m is unknown noise, and ρ > 0 is its magnitude. Various optimization techniques can address these problems. We focus on the bilevel formulation, widely adopted in the literature [Beck and Sabach, 2014 , Sabach and Shtern, 2017 , Dempe et al., 2021 , Latafat et al., 2023 , Merchav and Sabach, 2023] .The lower-level objective in the bilevel formulation is given byEQUATIONwhereIC (x) is the indicator function of a set C that satifies IC (x) = 0 if x ∈ C, and IC (x) = +∞ if x / ∈ C. The set C is a closed, convex set that can be chosen as C = R n , C = {x ∈ R n : x ≥ 0}, or C = {x ∈ R n : ∥x∥1 ≤ θ} for some θ > 0.This problem may have multiple minimizer solutions. Hence, a reasonable option is to consider the minimal norm solution problem, i.e., find the optimal solution with the smallest Euclidean norm [Beck and Sabach, 2014 , Sabach and Shtern, 2017 , Latafat et al., 2023] :F (x) = 1 2 ∥x∥ 2 .We need to solve the simple bilevel optimization problem:min x∈R n 1 2 ∥x∥ 2 s.t. x ∈ arg min z∈R n 1 2m ∥Az -b∥ 2 + IC (z).Example A.2 (Sparse Solution of Linear Inverse Problems). Consider the same setting as in Example A.1, but with the additional goal of finding a sparse solution among all the minimizers of the linear inverse problem (7). This can simplify the model and improve computational efficiency. To achieve sparsity, we can use any function that encourages it. One such function is the well-known elastic net regularization [Friedlander and Tseng, 2008 , Amini and Yousefian, 2019 , Merchav and Sabach, 2023] , which is defined asF (x) = ∥x∥ 1 + τ 2 ∥x∥ 2 ,where τ > 0 regulates the trade-off between ℓ1 and ℓ2 norms.This example corresponds to our second experiment in Section 4.2.Example A.3 (Logistic Regression Problem). The logistic regression problem aims to map the feature vectors ai to the target labels bi. A standard machine learning technique for this problem is to minimize the logistic loss function over the given dataset [Amini and Yousefian, 2019 , Gong et al., 2021 , Jiang et al., 2023 , Latafat et al., 2023 , Merchav and Sabach, 2023] . We assume that the dataset consists of a feature matrix A ∈ R m×n and a label vector b ∈ R m , with bi ∈ {-1, 1} for each i. The logistic loss function is defined asEQUATIONOver-fitting is a common issue when the number of features is large compared to the number of instances m. A possible approach is to regularize the logistic objective function with a specific function or a constraint [Jiang et al., 2023, Merchav and Sabach, 2023] . For instance, we can use g2(x) = IC (x), where IC (x) is the indicator of the setC = {x ∈ R n : ∥x∥1 ≤ θ}, as in Example A.1.This problem may also have multiple optimal solutions. Hence, a natural extension is to consider the minimal norm solution problem [Gong et al., 2021 , Jiang et al., 2023 , Latafat et al., 2023] , as in Example A.1. This requires solving the following problem:min x∈R n 1 2 ∥x∥ 2 s.t. x ∈ arg min z∈R n 1 m m ∑ i=1 log(1 + exp(-a T i zbi)) + IC (z).When choosing C = {x ∈ R n : ∥x∥1 ≤ θ} for some θ > 0, it corresponds to our first experiment in Section 4.1.Example A.4 (Over-parameterized Regression Problem). The linear regression problem aims to find a parameter vector x ∈ R n that minimizes the training loss ℓtr(x) over the training dataset Dtr. Without explicit regularization, the over-parameterized regression problem has multiple minima. However, these minima may have different generalization performance. Therefore, we introduce a secondary objective, such as the validation loss over a validation set D val , to select one of the global minima of the training loss. This results in the following bilevel problem:EQUATION)For instance, we can consider the sparse linear regression problem, where the lower-level objective consists of the training error and a regularization term, namely, G(x) = 1 2 ∥Atrx -btr∥ 2 + IC (x). Here, IC (x) denotes the indicator of a convex set, as in Example A.2. The upper-level objective is the validation error, i.e.,F (x) = 1 2 ∥A val x -b val ∥ 2 .The linear regression problem is over-parameterized when the number of features n is larger than the number of data instances in the training set.Table 1: Summary of simple bilevel optimization algorithms. The abbreviations "SC," "C,", "diff", "comp", "WS" and "C3" represent "strongly convex," "convex,", "differentiable", "composite", "weak sharpness" and "Convex objective with Convex Compact constraints," respectively. The abbreviation α-HEB refers to Hölderian error bound with exponent parameter α. We only include the gradients Lipschitz constant in the complexity result when its relation to the complexity is clear; otherwise, we omit it. Notation l F is the upper bound of subdifferentials of F , L f1 and L g1 are the Lipschitz constants of ∇f 1 and ∇g 1 , respectively. ), ε ∈ (0, 0.5)IR-CG [Giang-Tran et al., 2023] C, smooth C3, smooth (ϵF , ϵG) O ( max{1/ϵ 1 1-p F , 1/ϵ 1 p G } ) p ∈ (0, 1)Tseng's method [Malitsky, Bi-SG [Merchav and Sabach, 2023 ]C, quasi-Lip/comp C, comp (ϵF , ϵG) O ( max{1/ϵ 1 1-a F , 1/ϵ 1 a G } ) , a ∈ (0.5, 1) µ-SC, comp C, comp (ϵF , ϵG) O ( max{ ( log 1/ϵF µ ) 1 1-a , 1/ϵ 1 a G } ) , a ∈ (0.5, 1) R-APM [Samadi et al., 2023] C, smooth C, comp, WS (ϵ, ϵ) O ( √ 1/ϵ ) Online Framework [Shen et al., 2023] C, Lip C3, Lip (ϵF , ϵG) O ( max{1/ϵ 3 F , 1/ϵ 3 G } ) Our method C, comp C, comp, α-HEB (ϵ, l -β F ϵ β ) O (√ Lf 1 ϵ + √ l max{α,β} F Lg 1 ϵ max{α,β} ) , α ≥ 1, β > 0 µ-SC, comp C, comp, α-HEB (ϵ, l -β F ϵ β ) O ( √ Lf 1 µ log 1 ϵ + √ l max{α,β} F Lg 1 ϵ max{α-1,β-1} log 1 ϵ ) , α ≥ 1, β > 0 nonsmooth, Lip nonsmooth, Lip, α-HEB (ϵ, l -β F ϵ β ) O ( l 2 f2 ϵ 2 + l max{2α,2β} f2 l 2 g2 ϵ max{2α,2β} ) , α ≥ 1, β > 0We present several examples of functions that satisfy the Hölderian error bound Assumption 2.2 and their corresponding exponent parameter α in Table 2 . We also provide some clarifications for Table 2 below. The abbreviations "Q ∈ S n " and "Q ≻ 0" stand for "Q is a symmetric matrix of order n and a positive definite matrix, respectively. We refer the reader to Pang [1997] G(x) Remarks Name α max i∈[m] {⟨a i , x⟩ -b i } a i ∈ R n , i ∈ [m], b ∈ R m piece-wise maximum 1 ∥x -x 0 ∥ Q = √ (x -x 0 ) T Q(x -x 0 ) Q ∈ S n , Q ≻ 0, x 0 ∈ R n Q-norm 1 ∥x -x 0 ∥ p x 0 ∈ R n , p ≥ 1 ℓ p -norm 1 ∥x∥ 1 + τ 2 ∥x∥ 2 τ > 0 Elastic net 1 or 2 4 ∥Ax -b∥ 2 A ∈ R m×n , b ∈ R m Least squares 2 1 m ∑ m i=1 log(1 + exp(-a T i xb i )) a i ∈ R n , i ∈ [m], b ∈ R m , A ∈ R m×n Logistic loss 2 η(x) + σ 2 ∥x∥ 2 η convex, σ > 0 Strongly-convex 2D Supplementary resultsAlgorithm 4 Adaptive PB-APG-sc method (aPB-APG-sc) 1: Input:x -1 = x 0 ∈ R n , γ 0 = γ 1 > 0, L f1 , L g1 , ν > 1, η > 1, ϵ 0 > 0. 2: for k ≥ 0 do 3: ϕ k (x) = f 1 (x) + γ k g 1 (x) 4: ψ k (x) = f 2 (x) + γ k g 2 (x) 5:Invokex k = PB-APG-sc(ϕ k , ψ k , µ, L f1 , L g1 , x k-1 , ϵ k ) 6: ϵ k+1 = 1 η ϵ k 7:γ k+1 = νγ k 8: end for Similar to Algorithm 2, we have the following convergence results of Algorithm 4.Theorem D.1. Suppose that Assumptions 2.1, 2.2, 3.1, 3.2, and 3.7 hold. Let ϵ0 > 0 be given.• When α > 1, set ν > η α-1 , N = ⌈log η 1-α ν (ρL α F (α -1) α-1 α -α ϵ 1-α 0 /γ0)⌉+ and γ * k = ρL α F (α - 1) α-1 α -α ϵ 1-α 0 η k(α-1) ; • When α = 1, set ν > 1, N = ⌈log ν (ρlF /γ0)⌉+ and γ * k = ρLF .Then, for any k ≥ N , Algorithm 2 generates an( ϵ 0 η k , 2ϵ 0 η k (γ 0 ν k -γ * k ))-optimal solution of problem (P) after at most K iterations, where K satisfiesK = O   √ L f 1 µ log η k ϵ0 + √ ν k l max{α,β} F Lg 1 ϵ max{α-1,β-1} log η k ϵ0   .The proof is similar to the proof of Theorem 3.5 in Appendix E.8. So we omit it here.In this section, we propose the proofs of our main convergence results in this paper.Proof. Since Xopt is closed and convex [Beck and Sabach, 2014] , the projection of any x ∈ R n onto Xopt, denoted as x, exists and is unique. Furthermore, it holds that dist(x, Xopt) = ∥x -x∥.Then, by Assumption 2.1, we haveEQUATIONEQUATIONwhere (a) follows from the Hölderian error bound assumption of p(x), and (b) is from the fact that y =-lF z + γ * ρ z α attains its minimum at z * = ( ρl F αγ * ) 1 α-1 .Since x ∈ Xopt is feasible for problem (P), we have F (x) ≥ F * . This along with (11) indicatesEQUATION)Let x * be an optimal solution of (P) so that F (x * ) = F * . In addition, since x * ∈ Xopt, we have p(x * ) = 0. Combine these results with (12), we haveEQUATION)This demonstrates that an optimal solution of (P) is an ϵ-optimal solution for (Pγ).Proof. The proof is motivated by Theorem 1 in Luo et al. [1996] . Denote x * , x * γ as optimal solutions of problem (P) and (Pγ), respectively.For any x ∈ R n , let x be the projection of x onto Xopt. Then x is a feasible solution of (P) and F (x) ≥ F (x * ) holds. Then we haveEQUATIONwhere (a) follows from (10) and the Hölderian error bound assumption of p(x), and (b) follows from γ ≥ ρlF . Therefore, we conclude that x * is an optimal solution of (Pγ).For the converse part, let x * γ be the projection of x * γ onto Xopt. Then x * γ is a feasible solution of (P). Therefore, it holds that F (x * γ ) ≥ F (x * ). Similarly, we haveEQUATIONwhere the inequality (c) follows from (10) and the Hölderian error bound assumption of p(x).Therefore, all inequalities in (15) become equalities. We deduce that ∥x * γ -x * γ ∥ = 0 if γ > ρlF , implying that x * γ is in Xopt, i.e., p(x * γ ) = 0. Furthermore, as the first inequality of ( 15) becomes an equality, we obtainF (x * ) = F (x * γ ) + γp(x * γ ) = F (x * γ ). Therefore, x *γ is also an optimal solution of (P).Proof. Denote x * , x * γ as optimal solutions of problem (P) and (Pγ), respectively.• Case of α > 1. Since x * γ is an ϵ-optimal solution of (Pγ), we haveF (x * γ ) + γp(x * γ ) ≤ F (x) + γp(x) + ϵ, ∀x ∈ R n .(16) Note that the arguments in the proof of Lemma 2.3 still hold. Substituting x = x * into (16) and utilizing p(x * ) = 0, we haveF (x * γ ) + γp(x * γ ) ≤ F (x * ) + ϵ = F (x * ) + γ * p(x * ) + ϵ ≤ F (x * γ ) + γ * p(x * γ ) + 2ϵ, where the last inequality follows from setting x = x * γ in (13). Then, it holds thatEQUATION)By setting x = x * in ( 16), we haveEQUATION)Combing ( 18) with ( 17), we conclude that x * γ is an (ϵ, l -β F ϵ β )-optimal solution of (P). • Case of α = 1. Since x * γ is an ϵ-optimal solution of (Pγ), we haveEQUATION)On the one hand, as γ = γ * + l β F ϵ 1-β > γ * , by Lemma 2.4, x * γ is an optimal solution of (P). On the other hand, since γ ≥ γ * , according to Lemma 2.4, x * is also an optimal solution of (Pγ). Therefore, p(x * ) = 0 and p(x * γ ) = 0, it holds thatEQUATIONwhere the first inequality follows from the fact that x * is an optimal solution of (Pγ), and the last inequality follows from the optimality of x * to (Pγ) when γ ≥ γ * .The second inequality of (20) and p(x * γ ) ≥ 0 imply thatEQUATION)In addition, from (20), we haveEQUATION)This result along with (21) demonstrate that x * γ is an (ϵ, l -β F ϵ β )-optimal solution of (P).Proof. Let x * γ be the projection of x * γ on Xopt, we have ∥x * γ -x * γ ∥ = dist(x * γ , Xopt). By Assumption 2.2, the following inequality holds,EQUATION)where (a) follows from (17) when α > 1 or from ( 22) when α = 1.By Assumption 2.1, we haveF (x * γ ) -F * ≥ F (x * γ ) -F (x * γ ) (10) ≥ -lF ∥x * γ -x * γ ∥ ≥ -lF ( ρl -β F ϵ β ) 1 α ,where the first inequality follows from F (x * γ ) ≥ F * and x * γ ∈ Xopt.Proof. For any x ∈ dom(F ), let x be the projection of x onto Xopt, where the existence and uniqueness of x follows from that Xopt is closed and convex. Since F is l-Lipschitz continuous, similar to (10), we haveEQUATIONTherefore, all the requirements of (10) in equations ( 11), ( 14) and ( 15) can be replaced by ( 24). This implies that Lemmas 2.3 and 2.4 also hold for the global solutions of problems (P) and (Pγ) when F is non-convex.Then, the final result follows a similar pattern to Theorem 2.5. Here we omit it.E.6 Proof of Theorem 2.8Proof. Let x * γ be the projection of x * γ onto Xopt and x * r) . Then, we haveγ = cx * γ + (1 -c)x * γ with c = min{1, 1 - r ∥x * γ -x * γ ∥ }, which implies that x * γ ∈ B(x * γ ,EQUATION) where inequality (i) follows from the convexity of p(x).Inequality (25) demonstrates that γ(1 -c)p(x * γ ) ≤ F (x * γ ) -F (x * γ ) ≤ l∥x * γ -x * γ ∥ = l(1 -c)∥x * γ -x * γ ∥ ≤ l(1 -c)(ρp(x * γ )) 1 α, where the second inequality follows from the l-Lipschitz continuity of F on B(x * γ , r). Therefore, it holds thatEQUATION) • Case of α > 1. By (26), we have p(x * γ ) ≤ ( ρl α γ α ) 1 α-1 , which demonstrates that p(x * γ ) ≤ ϵ if γ ≥ ( ρl α ϵ α-1 ) 1 α .Then, for any xγ ∈ B(x * γ , r) that also satisfies p(xγEQUATION) which implies that F (x * γ ) -F (xγ) ≤ γ(p(xγ) -p(x * γ )) ≤ 0.The desired result follows. Proof. From [Beck, 2017, Theorem 10 .34], the objective value after K iterations can be bounded byΦγ(xK ) -Φ * γ ≤ 2Lγ∥x0 -x * ∥ 2 (K + 1) 2 ,whereLγ = L f 1 + γLg 1 .Combining this with our stopping criterion, we find that after K iterations,Φγ(xK ) -Φ * γ ≤ ϵ.This indicates that we obtain an ϵ-optimal solution to problem (Pγ). The value of K satisfies:K = √ 2(L f 1 + γLg 1 ) ϵ R -1.Specifically, we analyze the value of K in various scenarios in the form of O(•).• Case of α > 1. In this case, γ = γ * + 2l β F ϵ 1-β comprises two components: γ * and 2l β F ϵ 1-β . Therefore, it is natural to discuss which of these two components plays the dominant role in the complexity results. First, we write K in the form:K = √ 2(L f 1 + (ρl α F (α -1) α-1 α -α ϵ 1-α + 2l β F ϵ 1-β )Lg 1 ) ϵ R -1. If β < α, the dominating term in γ is γ * = ρl α F (α -1) α-1 α -α ϵ 1-α .Then, the number of iterations isK = O (√ L f 1 + l α F ϵ 1-α Lg 1 ϵ ) = O ( √ L f 1 ϵ + √ l α F Lg 1 ϵ α ) . If β = α, we have γ = ( ρ(α -1) α-1 α -α + 2 ) l α F ϵ 1-α .Then, the number of iterations isK = O (√ L f 1 + l α F ϵ 1-α Lg 1 ϵ ) = O ( √ L f 1 ϵ + √ l α F Lg 1 ϵ α ) . If β > α, the dominating term in γ is 2l β F ϵ 1-β .Then, the number of iterations isK = O   √ L f 1 + 2l β F ϵ 1-β Lg 1 ϵ   = O   √ L f 1 ϵ + √ l β F Lg 1 ϵ β   .• Case of α = 1. In this case, γ = γ * + l β F ϵ 1-β , where γ * = ρlF . Similarly, we explore which of these two elements plays a more significant role.If β < 1, the dominating term in γ is γ * . Then, the number of iterations isK = O ( √ L f 1 + ρlF Lg 1 ϵ ) = O ( √ L f 1 ϵ + √ lF Lg 1 ϵ).If β = 1, we have γ = (ρ + 1)lF ϵ 1-α . Then the number of iterations isK = O (√ L f 1 + (ρ + 1)lF Lg 1 ϵ ) = O ( √ L f 1 ϵ + √ lF Lg 1 ϵ ) . If β > 1, the dominating term in γ is l β F ϵ 1-β .Then, the number of iterations isK = O   √ L f 1 + l β F ϵ 1-β Lg 1 ϵ   = O   √ L f 1 ϵ + √ l β F Lg 1 ϵ β   .Combining the above results, we conclude thatK = O   √ L f 1 ϵ + √ l max{α,β} F Lg 1 ϵ max{α,β}   .Proof. In this proof, we denote Φ * k as the optimal value of problem (Pγ) when γ = γ k , and x k as the output of PB-APG (Algorithm 1) in the k-th iteration.• Case of α > 1. Suppose that N is the smallest nonnegative integer such that γN ≥ γ * N := ρl α F (α - 1) α-1 α -α ϵ 1-α N .In this case, we haveEQUATIONEQUATION)From ( 29), after at mostN := ⌈log η 1-α ν ( ρl α F (α-1) α-1 α -α ϵ 1-α 0 γ 0 ) ⌉+ iterations, (28) holds. Since xN = PB-APG(ϕN , ψN , L f 1 , Lg 1 , xN-1, ϵN ), we have ΦN (xN ) -Φ * N ≤ ϵN , γN ≥ γ * N ,which shows that xN is an ϵN -optimal solution of (Pγ) with γ = γN . From the proof in Theorem 2.5 (see inequalities ( 17) and ( 18) in Appendix E.3), xN is also an( ϵ 0 η N , 2ϵ 0 η N (γ 0 ν N -γ * N ))-optimal solution of problem (P). Furthermore, note that for any iteration k ≥ N , inequality (29) always holds, which means that the following statement holds for any k ≥ N :EQUATION)Let I k be the number of iterations of PB-APG required to satisfy (30) at the k-th iteration of aPB-APG. Then, for any k ≥ N , the total number of iterations isK = I0 + I1 + • • • + I k .From [Beck, 2017, Theorem 10.34 ], the number of iterations in i-th inner loop satisfies:Ii = √ 2(L f 1 + γiLg 1 ) ϵi ∥xi-1 -x * i ∥ -1,where x * i is the optimal solution in i-th inner loop. Then we have thatK = k ∑ i=0 √ 2(L f 1 + γiLg 1 ) ϵi ∥xi-1 -x * i ∥ -k ≤ k ∑ i=0 √ 2(L f 1 + γ k Lg 1 ) ϵi R -k = η k 2 -1 η 1 2 -1 √ 2(L f 1 + γ0ν k Lg 1 ) ϵ0 -k.For simplicity, we can also use O(•) to show the value of K.K = O ( √ L f 1 + γ0Lg 1 ϵ0 ) + • • • + O ( √ L f 1 + γ k Lg 1 ϵ k ) ≤ O ( √ L f 1 + γ k Lg 1 ϵ0 ) + • • • + O ( √ L f 1 + γ k Lg 1 ϵ k ) = O ( √ L f 1 + γ k Lg 1 ϵ k ( 1 + √ 1/η + √ 1/η 2 + • • • + √ 1/η k ) ) = O ( √ L f 1 + γ k Lg 1 ϵ k ) = O ( √ L f 1 η k ϵ0 + √ Lg 1 γ0(ην) k ϵ0) .• Case of α = 1. Suppose that after N updates, we have γN ≥ ρlF , i.e.,EQUATION)This demonstrates that after for all k ≥ N := log ν( ρl F γ 0) , (31) always holds.Similar to the case of α > 1, the total iteration number is:K = O ( √ L f 1 + γ0Lg 1 ϵ0 ) + • • • + O ( √ L f 1 + γ k Lg 1 ϵ k ) = O ( √ L f 1 + γ k Lg 1 ϵ k ) = O ( √ L f 1 η k ϵ0 + √ Lg 1 γ0(ην) k ϵ0). F ϵ 1-β . Then, we haveK = O   √ L f 1 + l β F ϵ 1-β Lg 1 µ log 1 ϵ   = O   √ L f 1 µ log 1 ϵ + √ l β F Lg 1 ϵ β-1 log 1 ϵ   .Combining the above results, we conclude thatK = O   √ L f 1 µ log 1 ϵ + √ l max{α,β} F Lg 1 ϵ max{α-1,β-1} log 1 ϵ   .E.10 Proof of Theorem 3.9Proof. Denote lγ = l f 2 + γlg 2 . Define Φ K γ,best = min i=0,...,KΦγ(xi) and ΦK,j γ,best = min i=j,...,K Φγ(xi) for all 0 ≤ j ≤ K. We claim that the sequence generated by the subgradient method satisfiesEQUATION)Specifically, from Lemma 8.24 in Beck [2017] , for all 0 ≤ j ≤ K, we haveEQUATION)Define ⌊•⌋ and ⌈•⌉ as rounding up and rounding down, respectively. LetEQUATION)where the second inequality follows from that∑ K k=⌊ K 2 ⌋ 1 k+1 ≤ ∫ K ⌈ K 2 ⌉-1 1 s+1 ds ≤ 2 log 2 and ∑ K k=⌊ K 2 ⌋ 1 √ k+1 ≥ ∫ K+1 ⌈ K 2 ⌉ 1 √ s+1 ds ≥ 1 2 √ K + 2.From the fact that Φ K γ,best ≤ ΦK,j γ,best , The desired result of (35) follows. Then, inequality (35) demonstrates that the number of iterations to obtain an ϵ-optimal solution for problem (Pγ) isK = O ( l f 2 + γlg 2 ϵ ) 2 . • Case of α > 1. we have γ = γ * + 2l β f 2 ϵ 1-β and γ * = ρl α f 2 (α -1) α-1 α -α ϵ 1-α . If β < α, the dominating term in γ is γ * . Then, the number of iterations is K = O ( l f 2 + l α f 2 ϵ 1-α lg 2 ϵ ) 2 = O ( l 2 f 2 ϵ 2 + l 2α f 2 l 2 g 2 ϵ 2α).If β = α, we have γ = ( ρ(α -1) α-1 α -α + 2 ) l α F ϵ 1-α .Then, the number of iterations isK = O ( l f 2 + l α f 2 ϵ 1-α lg 2 ϵ ) 2 = O ( l 2 f 2 ϵ 2 + l 2α f 2 l 2 g 2 ϵ 2α) .If β > α, the dominating term in γ is 2l β F ϵ 1-β . Then, the number of iterations isK = O ( l f 2 + 2l β f 2 ϵ 1-β lg 2 ϵ ) 2 = O ( l 2 f 2 ϵ 2 + l 2β f 2 l 2 g 2 ϵ 2β) .• Case of α = 1. we have γ = γ * + l β f 2 ϵ 1-β and γ * = ρl f 2 . If β < 1, the dominating term in γ is γ * . Then, the number of iterations isK = O ( l f 2 + ρl f 2 lg 2 ϵ ) 2 = O ( l 2 f 2 ϵ 2 + l 2 f 2 l 2 g 2 ϵ 2) .If β = 1, we have γ = (ρ + 1)lF ϵ 1-α . Then, the number of iterations isK = O ( l f 2 + ρl f 2 lg 2 ϵ ) 2 = O ( l 2 f 2 ϵ 2 + l 2 f 2 l 2 g 2 ϵ 2).If β > 1, the dominating term in γ is l β f 2 ϵ 1-β .Then, the number of iterations isK = O ( l f 2 + l β f 2 lg 2 ϵ 1-β ϵ ) 2 = O ( l 2 f 2 ϵ 2 + l 2β f 2 l 2 g 2 ϵ 2β) .Combining the above results, we conclude thatK = O ( l 2 f 2 ϵ 2 + l max{2α,2β} f 2 l 2 g 2 ϵ max{2α,2β}) .E.11 Proof of Theorem 3.10Proof. Denote lγ = l f 2 + γlg 2 , define Φ K γ,best = min i=0,...,KΦγ(xi). From Theorem 8.31 in Beck [2017] , the sequence generated by the subgradient method satisfiesΦ K γ,best -Φ * γ ≤ 2l 2 γ µ f 2 (K + 1) .This demonstrates that the number of iterations to obtain an ϵ-optimal solution for problem (Pγ) isK = O ( (l f 2 + γlg 2 ) 2 µ f 2 ϵ) .• Case of α > 1. we have γ = γ * + 2l β f 2 ϵ 1-β and γ * = ρl α f 2 (α -1) α-1 α -α ϵ 1-α . If β < α, the dominating term in γ is γ * . Then, the number of iterations is K = O ( (l f 2 + l α f 2 ϵ 1-α lg 2 ) 2 µ f 2 ϵ ) = O ( l 2 f 2 µ f 2 ϵ + l 2α f 2 l 2 g 2 µ f 2 ϵ 2α-1).If β = α, we have γ = ( ρ(α -1) α-1 α -α + 2 ) l α F ϵ 1-α .Then, the number of iterations isK = O ( (l f 2 + l α f 2 ϵ 1-α lg 2 ) 2 µ f 2 ϵ ) = O ( l 2 f 2 µ f 2 ϵ + l 2α f 2 l 2 g 2 µ f 2 ϵ 2α-1) .If β > α, the dominating term in γ is 2l β F ϵ 1-β . Then, the number of iterations isK = O ( (l f 2 + 2l β f 2 ϵ 1-β lg 2 ) 2 µ f 2 ϵ ) = O ( l 2 f 2 µ f 2 ϵ + l 2β f 2 l 2 g 2 µ f 2 ϵ 2β-1) .• Case of α > 1. we have γ = γ * + l β f 2 ϵ 1-β and γ * = ρl f 2 . If β < 1, the dominating term in γ is γ * . Then, the number of iterations isK = O ( (l f 2 + ρl f 2 lg 2 ) 2 µ f 2 ϵ ) = O ( l 2 f 2 µ f 2 ϵ + l 2 f 2 l 2 g 2 µ f 2 ϵ) .If β = 1, we have γ = (ρ + 1)lF ϵ 1-α . Then, the number of iterations isK = O ( (l f 2 + ρl f 2 lg 2 ) 2 µ f 2 ϵ ) = O ( l 2 f 2 µ f 2 ϵ + l 2 f 2 l 2 g 2 µ f 2 ϵ ) . If β > 1, the dominating term in γ is l β f 2 ϵ 1-β .Then, the number of iterations isK = O ( (l f 2 + l β f 2 lg 2 ϵ 1-β ) 2 µ f 2 ϵ ) = O ( l 2 f 2 µ f 2 ϵ + l 2β f 2 l 2 g 2 µ f 2 ϵ 2β-1) .Combining the above results, we conclude thatK = O ( l 2 f 2 µ f 2 ϵ + l max{2α,2β} f 2 l 2 g 2 µ f 2 ϵ max{2α-1,2β-1}) .In this section, we provide supplementary experiment settings and results. Specifically, in Appendix F.1, we present the detailed experimental settings, and in Appendix F.2, we provide the detailed experimental results. Additionally, in Appendix F.3 and F.4, we conduct experiments with different values of penalty parameter γ and solution accuracy ϵ, respectively.All simulations are implemented using MATLAB R2023a on a PC running Windows 11 with an AMD (R) Ryzen (TM) R7-7840H CPU (3.80GHz) and 16GB RAM.We conduct the first experiment using the a1a.t data from LIBSVM datasets5 . This data consists of 30, 956 instances, each with n = 123 features. For this experiment, a sample of 1, 000 instances is taken from the data, denoted as A. The corresponding labels for these instances are denoted as b, where each label bi is either -1 or 1, corresponding to the i-th instance ai.The Greedy FISTA algorithm [Liang et al., 2022] is used as a benchmark to compute G * . To compute the proximal mapping of f2(x) + γg2(x) in problem (Pγ), i.e, projection onto a 1-norm ball, we utilize the method proposed in Duchi et al. [2008] , which performs exact projection in O(n) expected time, where n is the dimension of x.For the PB-APG and PB-APG-sc algorithms, we set the value of γ = 10 5 , and we terminate the algorithms when ∥x k+1 -x k ∥ ≤ 10 -10 . For the aPB-APG and aPB-APG-sc algorithms, we set γ0 = 1 2 5 , ν = 20, η = 10, and ϵ0 = 10 -6 . The iterations of these two algorithms continue until ϵ k reaches 10 -10 (meanwhile, γ = 10 5 ).We compare our methods with MNG, BiG-SAM, DBGD, a-IRG, CG-BiO, Bi-SG, and R-APM in this experiment. Specifically, for R-APM [Samadi et al., 2023] , the regularization parameter η is set to η = 1/γ, reflecting the equivalence of the penalty formulation (Pγ) to (PReg), with σ = 1/γ, as previously discussed.We note that the termination criterion ∥x k+1 -x k ∥ ≤ 10 -10 used in our experiments is different from the one proposed in our algorithms since the parameters required for the latter are not easily measurable. Nevertheless, this termination criterion is also widely used in the literature, as it corresponds to a gradient mapping [Beck, 2017 , Nesterov, 2018 , Davis and Drusvyatskiy, 2019] . Furthermore, Theorem 3.5 of Drusvyatskiy and Lewis [2018b] implies that ∥x k+1 -x k ∥ also measures the distance to the optimal solution set.In the second experiment, we address the problem of least squares regression using the YearPredictionMSD data from the UCI Machine Learning Repository 6 . This data consists of 515, 345 songs with release years ranging from 1992 to 2011. Each song has 90 features, and the corresponding release year is used as the label.For this experiment, a sample of m = 1, 000 songs is taken from the data, and the feature matrix and release years vector are denoted as A and b, respectively. Following Section 5.2 in Merchav and Sabach [2023] , we apply the min-max scaling technique to normalize the feature matrix A. Additionally, we add an intercept term and 90 collinear features to A such that the resulting matrix A T A becomes positive semi-definite, which implies that the feasible set Xopt is not a singleton.We compare our methods with a-IRG, BiG-SAM, and Bi-SG in this experiment. Specifically, for BiG-SAM [Sabach and Shtern, 2017] , we consider the accuracy parameter δ for the Moreau envelope with two values, namely δ = 1 and δ = 0.01.To benchmark the performance, we utilize the MATLAB function lsqminnorm to compute G * . Moreover, we follow the parameter settings outlined in Section 4.1.To approximate the optimal value F * , we use the MATLAB function fmincon to solve a relaxed version of the function-value-based reformulations in equation (P Val ). In this relaxed version, we replace the constraint in (P Val ) with G(x) -G * ≤ ε, where ε = 10 -10 . This allows us to obtain an approximation of the optimal value while allowing for a small deviation from the true optimal value G * .We gather the total number of iterations for our methods, as well as the lower-and upper-level objective values and the optimal gaps for all the methods, in Table 3 . Subsequently, we compare the optimal gaps of all methods, which are defined as G(x) -G * and F (x) -F * for the lower-and upper-level optimal gaps, respectively. Table 3 reveals that for the logistic regression problem (5), our PB-APG, aPB-APG, PB-APG-sc, and aPB-APGsc exhibit almost identical function values for both objectives, surpassing other methods in terms of optimal gaps for the lower-and upper-level objectives (measured by the numerical value of the upper-level objective).In the case of the least squares regression problem (6), aPB-APG achieves the smallest optimal gaps for both objectives, followed by PB-APG and PB-APG-sc. These results demonstrate that our methods, despite yielding larger upper-level function values, generate solutions that are significantly closer to the optimal solution, as depicted in Figure 1 . Additionally, for the problem in (5), both aPB-APG and aPB-APG-sc require fewer iterations than PB-APG and PB-APG-sc, respectively. This can be attributed to the warm-start mechanism employed in aPB-APG and aPB-APG-sc. Moreover, for the problem in (6), both aPB-APG and aPB-APG-sc require more iterations than PB-APG and PB-APG-sc, respectively. However, they exhibit staircase-shaped curves, which avoid the unwanted oscillations in PB-APG and PB-APG-sc, we have a similar observation in Figure 2 .In this section, we investigate the impact of different values of penalty parameter γ on the experimental results of problems ( 5) and ( 6). We set γ to be either 2 × 10 4 or 5 × 10 5 for PB-APG and PB-APG-sc, and choose the corresponding γ0 values as 0.2 2 5 or 5 2 5 for aPB-APG and aPB-APG-sc, respectively. The remaining settings are the same as in Section 4.We plot the values of the residuals of the lower-level objective G(x k ) -G * and the upper-level objective over time in Figures 3 and 4 . Additionally, we also collect the total number of iterations, the lower-and upper-level objective values, and the optimal gaps of our methods in Table 4 for problems ( 5) and ( 6) with different values of γ. As Figures 3 and 4 show, our methods consistently outperform the other methods for both the lower-and upperlevel objectives, irrespective of the penalty parameter γ, since our methods achieve lower optimal gaps and desired function values for the lower-and upper-level objectives, respectively. The only exception is problem (6) with γ = 5 × 10 5 , as the third subfigure of Figure 4 shows, since we do not set the solution accuracy of BiG-SAM (δ = 0.01), it attains a lower optimal gap than our PB-APG-sc and aPB-APG-sc for the lower-level objective. However, BiG-SAM (δ = 0.01) produces significantly worse upper-level objective values, which are much larger than the objective values of our methods. , 4 , and 5 reveal that the number of iterations for our methods increases as penalty parameter γ increases. However, it is worth noting that the accuracy of the obtained solutions also increases, as indicated by the decreasing optimal gaps of the lower-and upper-level objectives. This observation confirms that the complexity results and solution accuracies of our methods are indeed dependent on the choice of penalty parameters, specifically, Lγ, as demonstrated in corresponding Theorem 3.3 and other related theorems. In this section, we investigate the impact of different solution accuracies on the experimental results of problems ( 5) and ( 6). We set ϵ to be either 10 -4 or 10 -7 and terminate the algorithms for PB-APG and PB-APG-sc when ∥x k+1 -x k ∥ ≤ ϵ. For aPB-APG and aPB-APG-sc, we choose the corresponding ϵ0 values as 1 or 10 -3 . The remaining settings are the same as in Section 4.We also plot the values of the residuals of the lower-level objective G(x k ) -G * and the upper-level objective over time in Figures 5 and 6 . Additionally, we also collect the total number of iterations, the lower-and upperlevel objective values, and the optimal gaps of our methods in Table 6 for problems ( 5) and ( 6) with different solution accuracies. From Figures 5 and 6 , it is evident that in most cases, our methods outperform the other methods in terms of both the lower-and upper-level objectives. However, there is an exception in the case of the upper-level objective for problem (6) when ϵ = 10 -4 . As illustrated in the second subfigure in Figure 6 , our methods exhibit larger function values for the upper-level objective compared to the other methods (except BiG-SAM (δ = 0.01)), despite still achieving smaller optimal gaps for the lower-level objective. This discrepancy actually indicates that our methods have not yet achieved the desired accuracy when ϵ = 10 -4 , and it is important to note that ∥x k+1 -x k ∥ ≤ ϵ is not the termination criterion in our proposed algorithms, as explained in Appendix F.1. Therefore, the larger optimality gaps for the upper-level objective in this case may be attributed to the termination criterion. 3, 6 , and 7 demonstrate that the number of iterations for our methods also increases with the solution accuracy, while the optimal gaps of the lower-and upper-level objectives decrease correspondingly. This finding confirms that the number of iterations and the optimal gaps are influenced by the solution accuracy, as illustrated in the expressions for the number of iterations provided by Theorem 3.3 and other related theorems.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.Justification: There is no societal impact of the work performed.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).In this case, we must have C bounded, as f2 is both strongly convex and Lipschitz continuous.According to Table2ofDoron and Shtern [2023], the parameter α can take values of either 1 or 2. Particularly, when α = 1, we have ρ = 1; when α = 2, we have ρ = 2/τ .https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/a1a.thttps://archive.ics.uci.edu/ml/datasets/YearPredictionMSD 


Decision: I would like to thank both the authors and reviewers for the good discussion. The reviewers agree that the paper is interesting although two reviewers are a little hesitant only recommending "borderline accept". Having looked at the paper and read the reviews, I recommend acceptance, however I urge the authors to take the reviewers feedback into account when preparing the camera-ready version.

Reviews: review: paper_summary: paper_summary: The authors propose algorithms to solve simple bilevel optimization problems of the form $\min_x F(x)$ s.t. $x$ minimizes $G$, in the case where $F$ and $G$ are composite convex functions (i.e. some of 2 convex functions, one of which is also smooth), and assuming Lipschitz continuity of $F$ on the set of minimizers of $G$ and $G$ verifying a Holder error bound.

The proposed algorithms actually tackle a simple minimization problem, solving a relaxed version of the original problem that authors previously proved to be equivalent.
strength_weakness: strength_weakness: I am not very familiar with this specific piece of literature, so I cannot pretend to be sure of the novelty of the results.
In particular, in line 47, the authors say « Tikhonov regularization most related to $P_\gamma$ », I do not see any difference between the two formulations. Can the authors develop on this? What is exactly done in Tikhonov's paper? And what is new here?


Minors:
- Usually, a proximal-based algorithm can be applied to the particular case where the non-smooth part of the problem is null, then the prox map is simply the identity map, and all the guarantees hold. But in this paper's setting, the non-smooth part is essential and cannot be considered null as a particular case. Indeed, otherwise, $G$ would be smooth, therefore upper bounded by a quadratic function, and $G$ is also assumed to grow faster than $||x-\bar{x}||^{alpha}$ which cannot happen (if alpha<2).
- A natural oracle assumption in the case where $F$ and $G$ contain each a non-smooth part would be to assume access to the proximal map associated with each. However, authors assume access to the proximal of all the linear combinations of the 2 non-smooth parts, and this choice seems to be only motivated by the reformulation of the problem. Doing this makes the derivation of PB-APG trivial once Thm 2.5. has been proven. But how could we tackle the same problem with the sole knowledge of the 2 prox maps independently?
questions: questions: - l.32: Why exclude the particular case where $X_{opt}$ is a singleton? This assumption is not used, therefore if we do not know in advance whether $X_{opt}$ is a singleton or not, one can still run this paper's algorithms, and guarantees will follow.
- l.39: Even if what $G^*$ is seems clear, I think it is worth mentioning it before using it in line 35.
- l.50: What about the upper problem?
- l.50: « where b in (0, 0.5) ». The authors should be more specific here. Is it true for one specific b in the interval, depending on some parameters of the problem? Is it true for any (one proposed algorithm for each value of b)?
- l.55: same remark
- l.61: It would be more convenient if the authors keep a consistent way of displaying complexity.
- l.68: What is g_2? Not defined yet at this stage. Moreover, note that the authors could not look at what happens in their case when $g_2=0$, as discussed in the "weakness" section.
- l.101: What is f_2? Again not defined. And I think here that the authors meant $F$.
- l.193: typo: « an » is repeated.
- l.205: « while samadi … » If they consider the smooth case only, studying $\alpha<2$ would be pointless because it is empty.
- It would be great if authors could add explicit stopping criteria in their algorithms. Indeed, if $F^*$ and $G^*$ are not known, I guess we stop PB-APG after K iterations where K is set according to Th3.3. But then, PB-APG is not only a function of $\phi$, $\psi$, $L_{f_1}$, $L_{g_1}, x_0, \epsilon$, but also $\ell_F$, $\alpha$ and $\beta$. Which is important to actually understand the iterates of aPB-APG.
limitations: limitations: No limitation
paper_summary: paper_summary: This work proposes a penalty based algorithm for simple bilevel optimization problems. The paper studies the relationship between the solutions of the penalized problem and the original bilevel problem. It extends the existing results on general bilevel optimization problem, which are established under PL condition, to the more generic Hölderian error bounds. The results are also extended to non-smooth upper and lower-level functions.
strength_weakness: strength_weakness: 1. Missing literature comparison/review:  The method proposed in this paper is a penalty-based method. There is a large body of work for penalty-based bilevel optimization methods. While in the related work section in the main text, there is no detailed review/comparison of the previous penalty-based bilevel optimization algorithms. It might be beneficial if the authors could discuss about the advantage of the proposed method compared with the existing ones.

2. In the current state of the paper, the experiments seem to be purely synthetic and might be a little too toy.  Though the paper is majorly a theoretical one, it would benefit to have more practical experiments. It is also of great importance to find meaningful applications for the simple bilevel optimization problem.

3. The method proposed in this work falls into the penalty-based methods. It might be interesting if the authors include a comparison with existing penalty-based bilevel optimization methods in, e.g.,
questions: questions: 1. The proposed method is for convex simple bilevel optimization problems, which might be restrictive. Can the author explain how the algorithm can be extended to the non-convex case?

2. How does this algorithm compare to the penalty-based algorithms for general bilevel optimization problems? Are there any advantages demonstrated by this algorithm as compared to the generic ones?
limitations: limitations: no potential social impact
paper_summary: paper_summary: This paper deals with a simple bilevel (the lower-level problem has no dependence on the upper-level variable) optimization problem, where both the upper- and the lower-level objectives are convex and potentially non-smooth. To simplify the original bilevel problem the authors consider a penalty-based single-level reformulation and establish the relationship between the solutions of the two problems. An accelerated proximal gradient method (PB-APG)  is proposed to solve the penalty problem. Under a Holderian error bound condition and certain additional assumptions about the convexity and smoothness of the objectives, convergence of the PB-APG is shown to approximate global optimal solutions of the original bilevel problem.
strength_weakness: strength_weakness: * From an applications perspective, the problem class considered here does not seem very interesting. The bilevel problem is simple, the objectives are (strongly) convex, a Holderian error bound holds and the non-smooth term is required to be proximal-friendly. This significantly restricts the number and complexity of potential applications. Indeed, the applications presented in the appendix seem to revolve around (regularized) linear least-squares problems and simple learning models like logistic regression.
* The proposed algorithm solves the reformulated minimization problem $P_{\Phi}$ rather than the original bilevel. This minimization problem looks fairly standard, and it seems that known algorithms are applied directly for its solution. Thus, it appears that there is no significant novelty in terms of algorithm design.
questions: questions: * Are Holderian error bounds satisfied in problems involving training of complex learning models, like neural networks.
* In table 1 some of the algorithms depend on parameters $\alpha$ and b, but the property to which these parameters correspond to is not specified. Please clarify this.
limitations: limitations: No negative societal impact
paper_summary: paper_summary: This article presents a nice extension to the broad class of penalty method problems, particularly for the case when the objective functions have a Hölderian error bound. This assumption for this class of problems appears new, and a thorough coverage of interesting results are claimed.
strength_weakness: strength_weakness: - The article does not distinguish between what type of subdifferential is used. Since the authors also make a mention of the applicability of their results to nonconvex objectives, this is crucial point. Many classical results in optimization and convex analysis are proven using the *convex* subdifferential; however, even for differentiable nonconvex functions (e.g., $-x^2$), the convex subdifferential is empty, while other notions like the Clarke subdifferential is nonempty.
- Assumption 3.2 is a very strong one; while I agree with the author(s) that the prox of a sum is widely studied, the list of settings in which one can compute this is still quite limited. However, this is not a major drawback for me, since handling HEB bilevel optimization even in the L-smooth setting appears to be new.
- In my opinion, leaving lines 140--142 to stake claim to results for non-convex which are exclusively in the appendix seems improper. The appendix is technically "unreviewed" material; if the formal statement of the result cannot fit into this article, then I do not think it is worth including. I suggest that the authors either try to find a way to incorporate these results into the article (since they are indeed quite interesting), or save the results for another work.
- The graphs of the numeric are too small for me to read, especially in-print.

Minor comments:
- The fact that the argmin of  a convex lsc function is closed and convex is a classical result which has been known for decades, so it seems odd to cite [Beck and Sabach, 2014] instead of a classical book, or the article where it was first proven.
- The article makes mention to functions $f_2$, $g_2$ far before they are introduced, which is confusiung to the reader.
- Names "Pock, Tseng, etc." are not capitalized in the citations.
- Citation on lines 391-392 is missing volume number.
questions: questions: The classical article on penalty methods,

Dmitri P. Bertsekas, "Necessary and Sufficient conditions for a penalty method to be exact", Math. Program., vol. 9, pp. 87 - 99, 1975.

basically states that, if the penalty function is smooth near the solution set, then the penalty must go to infinity in order to achieve exact penalization; and, if the penalty function is nonsmooth near the solution set, then a finite penalty will suffice. My question is: are your results consistent with this article? In particular,
  - In the setting where $G$ is smooth, can you confirm that $\lim_{\varepsilon\to 0}\gamma^*\to\infty$?
  - On the other hand, if $G$ is nonsmooth at its solution set, does $\lim_{\varepsilon \to 0} \gamma^*$ approach a finite value? Violating this second requirement would not contradict the aforementioned article, but it would significantly strengthen your results.
limitations: limitations: The authors do a good job of explaining the required assumptions for their results.

Comments: title: Response to authors' rebuttal
Comment from Reviewer fmaf
References
References
References
References
comment: We deeply appreciate your acknowledgement and decision to raise the score.
I thank the authors very much for their clear and thorough responses to both my questions and concerns. I am quite happy to hear that the revised version will include the precise subdifferential definition used, as well as the nonconvex results. I have decided to increase my score.
Thank you for acknowledging our response, and for raising the score.
Thank you for acknowledging our paper and response, and for the positive score.
I thank the authors for responding to my minor concerns.
Overall, the paper is well-written, sound and, I think, interesting for the Neurips community.
If the authors process the discussed modifications/clarifications, then I am in favor of the acceptance of this paper.
I would like to thank the authors for their responses. I am raising my score to 5.
[1] Bertsekas, Dimitri, Angelia Nedic, and Asuman Ozdaglar. Convex analysis and optimization. Vol. 1. Athena Scientific (2003).

[2] Cui Y, Liu J, Pang J S. Nonconvex and nonsmooth approaches for aﬀine chance-constrained stochastic programs[J]. Set-Valued and Variational Analysis, 2022, 30(3): 1149-1211.

[3] Doron, Lior, and Shimrit Shtern. Methodology and first-order algorithms for solving nonsmooth and non-strongly convex bilevel optimization problems. Mathematical Programming 201.1 (2023): 521-558.
    
[4] Jiang R, Abolfazli N, Mokhtari A, Hamedani EY. A conditional gradient-based method for simple bilevel optimization with convex lower-level problem. In International Conference on Artificial Intelligence and Statistics. PMLR (2023): pp. 10305-10323.
    
[5] Merchav, Roey, and Shoham Sabach. Convex Bi-level Optimization Problems with Nonsmooth Outer Objective Function. SIAM Journal on Optimization 33.4 (2023): 3114-3142.

[6] Yu, Yao-Liang. On decomposing the proximal map. Advances in neural information processing systems 26 (2013).

[7] Pustelnik, Nelly, and Laurent Condat. Proximity operator of a sum of functions; application to depth map estimation. IEEE Signal Processing Letters 24.12 (2017): 1827-1831.

[8] Liu, Hongying, Hao Wang, and Mengmeng Song. Projections onto the intersection of a one-norm ball or sphere and a two-norm ball or sphere. Journal of Optimization Theory and Applications 187 (2020): 520-534.

[9] Bertsekas, Dimitri P. Necessary and sufficient conditions for a penalty method to be exact. Mathematical programming 9.1 (1975): 87-99.
[1] Beck, Amir, and Shoham Sabach. A first order method for finding minimal norm-like solutions of convex optimization problems. Mathematical Programming 147.1 (2014): 25-46.

[2] Gong, Chengyue, and Xingchao Liu. Bi-objective trade-off with dynamic barrier gradient descent. NeurIPS 2021 (2021).

[3] Jiang R, Abolfazli N, Mokhtari A, Hamedani EY. A conditional gradient-based method for simple bilevel optimization with convex lower-level problem. In International Conference on Artificial Intelligence and Statistics. PMLR (2023): pp. 10305-10323.

[4] Doron, Lior, and Shimrit Shtern. Methodology and first-order algorithms for solving nonsmooth and non-strongly convex bilevel optimization problems. Mathematical Programming 201.1 (2023): 521-558.

[5] Sepideh Samadi, Daniel Burbano, and Farzad Yousefian. Achieving optimal complexity guarantees for a class of bilevel convex optimization problems. arXiv preprint arXiv:2310.12247, 2023
    
[6] Merchav, Roey, and Shoham Sabach. Convex Bi-level Optimization Problems with Nonsmooth Outer Objective Function. SIAM Journal on Optimization 33.4 (2023): 3114-3142.

[7] Yu, Yao-Liang. On decomposing the proximal map. Advances in neural information processing systems 26 (2013).

[8] Pustelnik, Nelly, and Laurent Condat. Proximity operator of a sum of functions; application to depth map estimation. IEEE Signal Processing Letters 24.12 (2017): 1827-1831.

[9] Liu, Hongying, Hao Wang, and Mengmeng Song. Projections onto the intersection of a one-norm ball or sphere and a two-norm ball or sphere. Journal of Optimization Theory and Applications 187 (2020): 520-534.


[10] Mostafa Amini and Farzad Yousefian. An iterative regularized incremental projected subgradient method for a class of bilevel optimization problems. In 2019 American Control Conference (ACC), pages 4069–4074. IEEE, 2019.

[11] Bolte, Jerome, et al. From error bounds to the complexity of first-order descent methods for convex functions. Mathematical Programming 165 (2017): 471-507.

[12] Zeng, Jinshan, et al. Global convergence of block coordinate descent in deep learning. International conference on machine learning. PMLR, 2019.

[13] Khanh-Hung Giang-Tran, Nam Ho-Nguyen, and Dabeen Lee. Projection-free methods for solving convex bilevel optimization problems. arXiv 2023.
[1] Han Shen and Tianyi Chen. On penalty-based bilevel gradient descent method. ICML 2023.

[2] Lu, Zhaosong, and Sanyou Mei. First-order penalty methods for bilevel optimization. SIAM JOPT 2024.

[3] Kwon, Jeongyeol, Dohyun Kwon, Steve Wright, and Robert Nowak. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. arxiv 2023.

[4] Phillips, David L. A technique for the numerical solution of certain integral equations of the first kind. JACM 1962.

[5] Doron, Lior, and Shimrit Shtern. Methodology and first-order algorithms for solving nonsmooth and non-strongly convex bilevel optimization problems. MP 2023.

[6] Jiang, R., Abolfazli, N., Mokhtari, A. and Hamedani, E.Y. A conditional gradient-based method for simple bilevel optimization with convex lower-level problem. AISTATS 2023.

[7] Gong C, Liu X. Bi-objective trade-off with dynamic barrier gradient descent[J]. NeurIPS 2021, 2021.

[8] Cao J, Jiang R, Abolfazli N, et al. Projection-free methods for stochastic simple bilevel optimization with convex lower-level problem[J]. Advances in Neural Information Processing Systems, 2024, 36.

[9] Cui Y, Liu J, Pang J S. Nonconvex and nonsmooth approaches for affine chance-constrained stochastic programs[J]. Set-Valued and Variational Analysis, 2022, 30(3): 1149-1211.

[10] Mostafa Amini and Farzad Yousefian. An iterative regularized incremental projected subgradient method for a class of bilevel optimization problems. In 2019 American Control Conference (ACC), pages 4069–4074. IEEE, 2019.

[11] Ruichen Jiang, Nazanin Abolfazli, Aryan Mokhtari, and Erfan Yazdandoost Hamedani. A condi- tional gradient-based method for simple bilevel optimization with convex lower-level problem. In International Conference on Artificial Intelligence and Statistics, pages 10305–10323. PMLR, 2023.

[12] Sepideh Samadi, Daniel Burbano, and Farzad Yousefian. Achieving optimal complexity guarantees for a class of bilevel convex optimization problems. arXiv preprint arXiv:2310.12247, 2023.

[13] Lingqing Shen, Nam Ho-Nguyen, and Fatma Kılınç-Karzan. An online convex optimization-based framework for convex bilevel optimization. Mathematical Programming, 198(2):1519–1582, 2023.

[14] Risheng Liu, Yaohua Liu, Shangzhi Zeng, and Jin Zhang. Towards gradient-based bilevel opti- mization with non-convex followers and beyond. Advances in Neural Information Processing Systems, 34:8662–8675, 2021.

[15] Daouda Sow, Kaiyi Ji, Ziwei Guan, and Yingbin Liang. A primal-dual approach to bilevel optimization with multiple inner minima. arXiv preprint arXiv:2203.01123, 2022.

[16] Lesi Chen, Jing Xu, and Jingzhao Zhang. On bilevel optimization without lower-level strong convexity. arXiv preprint arXiv:2301.00712, 2023.

[17] Feihu Huang. On momentum-based gradient methods for bilevel optimization with nonconvex lower-level. arXiv preprint arXiv:2303.03944, 2023.
[1] Mikhail Solodov. An explicit descent method for bilevel convex optimization. Journal of Convex Analysis, 14(2):227–237, 2007.

[2] Andre Nikolaevich Tikhonov and V. I. A. K. Arsenin. Solutions of ill-posed problems. Wiley, 1977.

[3] Doron, Lior, and Shimrit Shtern. Methodology and first-order algorithms for solving nonsmooth and non-strongly convex bilevel optimization problems. Mathematical Programming 201.1 (2023): 521-558.

[4] Jiang R, Abolfazli N, Mokhtari A, Hamedani EY. A conditional gradient-based method for simple bilevel optimization with convex lower-level problem. In International Conference on Artificial Intelligence and Statistics. PMLR (2023): pp. 10305-10323.

[5] Khanh-Hung Giang-Tran, Nam Ho-Nguyen, and Dabeen Lee. Projection-free methods for solving convex bilevel optimization problems. arXiv preprint arXiv:2311.09738, 2023.

[6] Sabach, Shoham, and Shimrit Shtern. A first order method for solving convex bilevel optimization problems. SIAM Journal on Optimization 27.2 (2017): 640-660.

[7] Merchav, Roey, and Shoham Sabach. Convex Bi-level Optimization Problems with Nonsmooth Outer Objective Function. SIAM Journal on Optimization 33.4 (2023): 3114-3142.

[8] Yu, Yao-Liang. On decomposing the proximal map. Advances in neural information processing systems 26 (2013).

[9] Pustelnik, Nelly, and Laurent Condat. Proximity operator of a sum of functions; application to depth map estimation. IEEE Signal Processing Letters 24.12 (2017): 1827-1831.

[10] Liu, Hongying, Hao Wang, and Mengmeng Song. Projections onto the intersection of a one-norm ball or sphere and a two-norm ball or sphere. Journal of Optimization Theory and Applications 187 (2020): 520-534.

[11] Mostafa Amini and Farzad Yousefian. An iterative regularized incremental projected subgradient method for a class of bilevel optimization problems. In 2019 American Control Conference (ACC), pages 4069–4074. IEEE, 2019.

[12] Kaushik, Harshal D., and Farzad Yousefian. A method with convergence rates for optimization problems with variational inequality constraints. SIAM Journal on Optimization 31.3 (2021): 2171-2198.