A Primal-Dual-Assisted Penalty Approach to Bilevel Optimization with Coupled Constraints
Authors: Liuyuan Jiang
Quan Xiao
Victor M. Tenorio
Fernando Real-Rojas
Antonio Marques
Tianyi Chen
Venue: NeurIPS 2024 poster (2024)
Keywords: Bilevel Optimization
Constrained Optimization
Hessian-free
Convergence Analysis
Penalty Based
Primal Dual
Primary area: optimization

Abstract: Interest in bilevel optimization has grown in recent years, partially due to its relevance for challenging machine-learning problems. Several exciting recent works have been centered around developing efficient gradient-based algorithms that can solve bilevel optimization problems with provable guarantees. However, the existing literature mainly focuses on bilevel problems either without constraints, or featuring only simple constraints that do not couple variables across the upper and lower levels, excluding a range of complex applications. Our paper studies this challenging but less explored scenario and develops a (fully) first-order algorithm, which we term BLOCC, to tackle BiLevel Optimization problems with Coupled  Constraints.  We establish rigorous convergence theory for the proposed algorithm and demonstrate its effectiveness on two well-known real-world applications - support vector machine (SVM) - based model training and infrastructure planning in transportation networks.

A Primal-Dual-Assisted Penalty Approach to Bilevel Optimization with Coupled Constraints
Liuyuan Jiang', Quan Xiao*, Victor M. Tenorio', Fernando Real-Rojas'
Antonio G. Marques", Tianyi Chen
Rensselaer Polytechnic Institute, Troy, NY, United States
Kiigg Juan Carlos University, Madrid, Spain
[jiang17, xiaoq5, chent18 @rpii edu
(victor. antonio.. aattonio.gariiaa (aaariiaaaruuusaa....
marques @urje. es; f.real 2012200111alaaa cummmossuujee. urjc.es
Abstract
Interest in bilevel optimization has grown in recent years, partially due its appli- cations to tackle challenging machine-learning problems. Several exciting recent works have been centered around developing efficient gradient-based algorithms that can solve bilevel optimization problems with provable guarantees. However, the existing literature mainly focuses bilevel problems either without constraints, or featuring only simple constraints that do not couple variables across the upper and lower-levels, excluding a range of complex applications. Our paper studies this challenging but less explored scenario and develops (fully) first-order algorithm, which we term BLOCC, to tackle BiLevel Optimization problems with Coupled Constraints. We establish rigorous convergence theory for the proposed algorithm and demonstrate effectiveness on two well-known real-world applications hyperparameter selection in support vector machine (SVM) and infrastructure planning transportation networks using the real data from the city of Seville.
Introduction
Bilevel optimization (BLO) approaches are pertinent in various machine learning problems, including
hyperparameter optimization 44 24], meta-learning [22], and reinforcement learning 667 63)
Moreover, the ability to handle BLO with constraints is particularly important, these constraints
appear in applications such as pricing [5], transportation [. and kernelized SVM [30]. Although there is extensive research
BLO problems without constraints or with uncoupled constraints 42 622 solutions for BLO problems with coupled constraints (CCs) remain limited;
31 12
39
details in Table T However, is of particular interest investigate BLO with lower-level CCs.
Taking infrastructure planning in a transportation network as an example, the lower-level seeks to optimize utility constrained by the upper-level parameter, network configuration.
Motivated by this, we consider the coupled-constraing BLO problem the following form
min [a, (x)) 2E
(1a)
s.t. V (x) == arg min g(z, v) with V(x) = {v E y 9xx v) S 0}} VE)(x>
(1b)
where ʃ: Rd x Rdd R the upper-level objective, g Rd Rdd
the lower-level
The work was supported by NSF project 2412486, NSF SCAL project 2401297, NSF CAREER project 2047177, Cisco Research Award, Amazon Research Award, the IBM-Rensselaer Future of Computing Research Collaboration, and the Comunidad de Madrid (via URIC grants F840 and F978).
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

objective which is strongly convex, ga y) R x R& + Rd defines the lower-level CCs, and X' C TR , C Rd are the domain of x and y that are easy to project, such as the Fuclidean ball.
The challenge in solving (1) arises from the coupling of the upper and lower-level problems. Prior work addressed this by starting with the unconstrained BLO problem, using implicit gradient descent
(IGD) methods 27 34 12 13 40 61 43
66 and penalty-based methods [45]
44
61
41 42
To solve BLO with CCs, AiPOD 71 and GAM 722 investigated the IGD method under different constraint settings. However, AIPOD only considered equality constraints, and GAM lacked finite- time convergence guarantees. Leveraging a penalty reformulation, (75) developed a Hessian-fres method with finite-time convergence. However, key algorithm step 77 is joint projection of the current iterate (s, y) onto the coupled constraint set. This projection, required at each iteration, becomes partieularly challenging when g ( y) is not jointly convex and can be computationally expensive for large-scale problems with a high number of variables (dz, dy) or constraints (do).
To this end, this paper aims to address the following question
Can we develop an efficient algorithm that bypasses joint projections g
(x. y) and quickly solves the BLO problem with coupled inequality constraints in !
We address this question affirmatively, focusing on the setting where the lower-level objective, g(x, y). strongly convex in y, and the constraints g² y) are convex in y. To avoid implementing a joint projection, we put forth novel single-level primal-dual-assiste r and y. Specifically, with μ E Id denoting the Lagrange multiplier of (Ib). we propose solving
penalty reformulation that decouples
min #E
Fx(x) = max min (xx v) + Hlo(,, v) v(x)) + (1, (x, (x, v) MER4 yEy
(2a)
penalty term
Lagrangian term
where v(x) = min g(z, y)
(2b)
y
where the penalty constant controls the distance between y and vg (z) by penalizing g(I, y) to its value function v(r), and the Lagrangian term penalizes the constraint violation of (I, v).
However, recognizing the max-min subproblems involved in (2a). it becomes computationally costly to evaluate the penalty function (2) and gradient. To this end, we pose the following question
Can we develop efficient algorithms to solve the max-min subproblem and evaluate V FF)
We answer this question by proving that this reformulation exhibits several favorable properties including smoothness. These properties are critical for designing gradient-based algorithms and characterizing their performance. However, the presence of the CCs renders the calculation of the gradient FF, more challenging than for its unconstrained counterpart [62]. Building upon this, we design a primal-dual gradient method with rigorous convergence guarantees for the BLO with general inequality CCs, and provide an improved result for the case of of being affine in y.
1.1 Main contributions
In a nutshell, our main contributions are outlined below.
C1) In Section 2 leveraging the Lagrangian duality theorem, introduce the function F, (z) in (2a) as penalty-based reformulation of (I), establish the continuity and smoothness of F, (r). and develop a novel way compute its gradient.
C2) In Section 3 we develop BLOCC, fully first-order algorithm to tackle BLO problems
with CCs. With being the target error for the generalized gradient norm square of F. Ff(z). we establish that the itcration complexity under the generic constraint g (z, y) in Tb is
Ö(c-2.5). We establish, for the first time, the linear convergence of a strongly convex- concave max-min problem with linear interaction and a constrained maximization parameter, reducing BLOCC's complexity to Ō(c-1.5) when the constraint gf (z., y) affine y.
C3) In Section 4 we apply our BLOCC algorithm to two real-world applications: SVM model training and transportation network planning. By comparison with LV-HBA (75) and GAM 722. we demonstrate the algorithm's effectiveness and its robustness to large-scale problems.

---

LL constraint
First Order
Complexity O(e-2,);
g0 v) < convex in y and LICQ holds;
BLOCC
Special case: ( y) affine in
O
LV-HBA
(a y) & 0 convex x y
Olc
GAM
9 (2,9)
9) S convex in x y and LICQ holds
BVESM
9 g
x satislying other requirements
/ O(e -1.)
AiPOD
9 [z 4) Ay- b(x)
X
Table I: Comparison of work with LV-HBA 775,, BVFSM 46), AiPOD [71], and Gradient Approximation method (GAM) [72]. LL convergence on metric the squared distance of y to its optimal solution, and UL convergence on squared (generalized) gradient norm.
1.2 Related works
BLO has a long history, dating back to the seminal work of [8. It has inspired rich body of literature, e.g. [76, 68 14 [65] The recent focus on BI.O is centered on developing efficient gradient-based approaches with provable finite-time guarantees.
Methods for BLO without constraints. A cluster of BLO gradient-based approaches gravitates around the implicit gradient descent (IGD) method 55 where the key idea is to approximate the hypergradient by the implicit function theorem. The finite-time convergence
IGD was first lower-level problems. Subsequent works improved the convergence rates and/or relaxed the assumptions under various settings; see [31 34
established in [27] for unconstrained strongly-conve
12
13
40
61 43
66
73 35 77] Another cluster of works is based on iterative differentiation (ITD)
60 which estimates the hypergradient by differentiating the entire iterative algorithm used to solve the lower-level problem with respect to the upper-level variables. The finite-time
49
23
53
guarantee was first established in 228 47 33 6). Viewing the lower-level problem as a constraint such in [64]. penalty- methods have also emerged as promising approach for BLO. Dated back to 777. this line of works 45
62 25 48) reformulated the original BLO as the single-level problems with various penalty terms and leveraged first-order methods to solve them.
51
BLO with constraints. While substantial progress has been made for unconstrained BLO, the analysis for constrained BLO is more limited. Upper-level constraints of the form x E X' were considered in [2].2] For the lower-level uncoupled constraint, SIGD
3g considered the uncoupled employed penalty reformulation and considered both upper and lower uncoupled constraints. However, the literature on BL. with CCs is scarce. BVFSM (46) conducted penalty-based method to avoid the calculation of the Hessian, as IGD methods do. However, only asymptotic convergence was achieved. GAM (72) investigated the IGD method under inequality constraints while failing to provide finite-time convergence results as well. AiPOD (71) also applied IGD and successfully achieved finite-time convergence, but only considered equality constraints. LV-HBA (75] considered inequality constraints and constructed a penalty- based reformulation. However, it employed joint projection of (z, y) onto {X x x gx y) < 0} which is computationally inefficient when there many constraints when (x, v) is not jointly convex. After our initial submission, we found a concurrent work 77 posted on ArXiv, which used Lagrange duality theory differently from ours, applying it construct new smoothed penalty term. However, it does not quantify the relationship between the relaxation of this penalty and the relaxation of lower-level optimality, and it does not guarantee lower-level feasibility. We summarized prior works on BLO with lower-level CCs in Table [
constraint Ay & b and achieved asymptotic convergence, 42 62
2 Primal-dual Penalty-based Reformulation
In this section, our goal to construct a caanaa-dal--ssssstee penalty reformulation for our BLOCC problem. The technical challenge comes from finding suitable penalty function for BLO with CCs.
2.1 The challenges in BLO with coupled constraints
Here, we will elaborate on the two technical challenges of BLO with CCs.
The first challenge associated with the presence of CCs is the difficulty to find the descent di- rection for z, which involves finding the closed-form expression of the gradient Vu(z). The ex- pression without CCs, Ve(x) = Vag(z, vy (x)) provided in (62 42 44), is not applicable.
3

---

For example, when v) = (y 2z)² and (z. y) = 3x y. the optimal lower-level solution is Vž (2) 3r and thus, v(x) x² with Ve(x) = 2r. However, (r. (x))
4r + (v()). The closed form gradient for BLO with CCs should be ( which
2
considers a Lagrange term that will be illustrated later in this paper. In Figure we present the gradient Vv(x)) without the Lagrange multiplier in [62] [2 41 [44] and ours with the Lagrange term. this example, the gradient without the Lagrange term leads to the opposite direction to the true gradient.
wish agle) wishouu M9(a)
-1
0
2
3
Calculation of Vu(z). The blue line is v(z). the the yellow dashed line is calculated by the formulation given in (62 41), while red dashed line is derived by our BLOCC. It can be seen that Ve(x) without the Lagrange multiplier is very biased from the true gradient.
Figure
The second challenge associated with the presence of CCs is the difficulty of performing the joint projection. If directly extend the penalty reformulation in 62 41, we could treat the coupling constraint set (z)) as joint constraint {(x,y) 9 gx v) 0), and employ a joint projection of (x.y) to ensure the feasibility. However, this can be computationally inefficient when the problem is large-scale, i.e. the number of variables or constraints is large. The detailed analysis of computational cost can seen Appendix H Moreover, when (x, y) is complex, the projection may not have a closed form and may not even be well-defined if gf (æ, y) is not jointly convex.
2.2 The Lagrangian duality-based penalty reformulation
Before proceeding, we summarize the assumptions considered as follows.
Assumption 1 (Lipschitz Continuity). Assume that f. f. Vg 9e and Vgg and are respectively 15.0- 1,,1 1g,2, so, and 1g-, -Lipschit jointly over (z. y) E X V, and g is L-ipsshht in x E X. Assumption 2 (Convexity in y). For any given x E X, g(,, v) and gx, y) are Otg-strongly convex and convex in y E y, respectively.
Assumption 3 (Domain Feasibility). Domain X C Rd and ) C TR.dv are non-empty, closed, and convex. For any given E X, y(x) == (y E ) g (x, v) S 0} is non-empty.
Assumption 4 (Constraint Qualification). For any given E X, of satisfies the Linear Independence Constraint Qualification (LICQ) condition for every y in a neighborhood of (z)
The Lipschitz continuity in Assumption I and the strong convexity of g over in Assumption 2 are conventional 31 41 7 33
Moreover, we only require g" ga( y) be convex in y rather than: i) jointly convex in (I, y) as in 75 or ii) linear as in [39 71]. Assumption! 3 pertains to the convexity and closeness of the domain, which is also conventional, and Assumption a standard constraint qualification condition.
To build penalty reformulation for BLO with lower-level CCs, the first challenge is find penalty that regulates |yy - vy (z)|]- In the following lemma, we show that g(z, y) - v(x) is a good choice.
Lemma 1. Suppose that Assumptions
2 hold and v(x) ) is defined as in (bb). Then, it holds that
cl) g( y) - v(x) 2 +1yy vn (1) for some d' > O. for all y E Max); and
c2) g(x,)) v) = v(x) if and only if vz). for all y E V(x).
The technical challenge of proving this lemma lies in showing the quadratic growth property in cl) of Lemma For BLO without the CCs, this naturally holds as the lower-level objective g(z, y) is strongly convex in y. For BLO with CCs, one needs to the Lagrangian duality theorem. The full proof of Lemma given in Appendix B. where the key challenge is to show the invariance of the modulus of strong convexity in the Lagrangian reformulated lower-level objective.
With y denoting a penalty constant, we consider the following penalty reformulation:
min
f(r y) + ((g(,, y) - v(x))
where the penalty term confines the squared Euclidean distance from y to vg (x). Furthermore, to avoid projecting (z, y) onto the (X x xY: g v) < 0}. we propose the dualaal-dl----ssttee penalty

---

reformulation, which was defined in ② The approximate equivalence of the reformulation to the original BLO problem ①
established in the following theorem.
Theorem 1 (Equivalence), Suppose that I is Lipschitz in y and V is 11,I-LLipcchit in (2r, y) in 2 hold. Then, solving the e-approximation problem of
Assumption hold, and Assumptions
min
f(z, y) s.t. l vr((2)] s,
is equivalent to solve the primal-dual penalty reformulation in with 7 (0-5)) and 7
aa
The detailed proof is provided in Appendix B By setting O (())) we effectively state the equivalence between the penalty reformulation (3) and the approximated original problem 4 We can then decouple the joint minimization on (a, y) to min-min problem on x and & constrained y. For the inner minimization problem in y, we choose y in way that yag 11,1 > 0 holds. This ensures the objective in (2a) being strongly convex in y. y -msoothhessss ensures lower bound for negative curvature fzz y). Furthermore, the convexity of gx y) in y validates the strong
duality theorem in 59 32 thereby enabling the max-min primal-dual reformulation in
2.3 Smoothness of the penalty reformulation
To evaluate (x) defined in 2a we can find the solution to the inner max-min problem:
() == arg max min [(x, v) + 7(((x,)) - v(x)) + (1, g (1. v))
(5)
pCR yE)
L, (1.2a)
The uniqueness of (x()) and u (x) is guaranteed under Assumptions and 4 (see Lemma 5 in Appendix A). Therefore, F-(x) in (2a) can be evaluated using the unique optimal solutions by
(x) = Lr (x)) yf (.)
(6)
Similarly, we can evaluate vzx) = Lo (1g ()); where
(1g(1), (x)) = arg max min g(r, v) + (M, (x, y)) MERT yEy
= =.g (M,y:))
The penalty reformulation (x) can hardly be convex, as (x) may be concave, even when g(2r, y) g(y) and f y) A y x. see Lemma 4.24 in (59] Instead, in the subsequent lemma, we not only show that v(x) differentiable, but also provide a closed-form expression of Vv(x). Lemma 2 (Danskin-like, theorem for v(x)). Suppose that Assumptions 7 hold, and let By co be a constant such that [1)](=))] < for all x E X. Then, it holds that
(x) and (x) defined 7 are L.g-Lipschitz some finite constant T.g 2 0.
) 2. v(x) defined in (2b) is lmmoooh where 1,,1 I (lg.1 + Bylge, Lg) gg and Vv(z) V.g(x. (z)) + (z). Vag (z. v7 (z))).
(8)
The assumption of the existence of the upper bound for the Lagrange multiplier a consequence of the LICQ condition (69, Theorem 1]. This assumption is mild and traditional 75) Finding Vv(ar) is crucial for the design gradient V FF,(z can be obtained by the next lemma.
gradient- based method to solve mingex (x).). Leveraging Lemmal 2 the
Lemma 3 (Danskin-like theorem for F, (z)). Suppose that the conditions in Lemma 2 hold. Moreover, ay" and there exist B < such that (1 Br, Vz E X. Then, it holds that
assume that y > 11.1
1. yo (x) and (x) defined in 5 are L Lipiscchhtt some constant LF 2
(x) is -moothh with la.1 I (l + 7g + Bplg Lg) + No,1 + LL and ( (x)) Ve(x)) + (x), (x,
(x)) V (z))
(9)
The proof of Lemmas 2 and is provided in Appendix Similar to 62 the Danskin-like theorems in Lemmas 2 and 3 rely on the Lipschitzness of the solutions 5 and Different from BLO
without CCs 62 the results here also hinge on the Lagrange multipliers ut, (x) and μ (z).
5

---

3
Main Results
We will first introduce algorithm tailored for BiLevel Optimization problems with inequality Coupled Constraints, and present its convergence analysis. In Section 3.2 we will propose a primal-dual solver for the inner max-min problems and characterize the overall convergence.
3.1 BLOCC algorithm
As F- (z) features differentiability and smoothness, we can apply a projected gradient descent (PGD)- based method to solve min (z).). At iteration update
Algorithm 1 Meta algorithm: BLOCC
I: inputs: initial points Tu. Ig.0s Mg.0s FF.. MFo stepsize Ng. nes counters Tg. TF.
J++1 Projx ngF,;)
(10)
2: for t = 1......... T do
with stepsize 1) and gF as estimate of (x))
3:
(Ugh1 P9g.)
To To
MaxMin Ty)
In the previous section, we have obtained the closed- form expressions of Vu(x) in (8) and VF, (z) in (9) Evaluating the closed-form expression requires find- ing (vğ (az)) ) (*x)) and (aa)) (1) the solu- tions to the max-min problem Lg (H, v: in and respectively.
4:
Tp Te MaxMin T,f)
(F.t)
5:
update I++1 via (10)
with gr. in (3)
6: end for
7: outputs: (2T, VFTT
LF y; 2)
Given 2t E X, we can use some minmax optimization solvers with Ty iterations on 7 and Tr To iterations on (5) to find an -s-ouutionn (Vg,t (g..)) and ep-solution (WP TE) satisfying
To (ggh) - (G (z)))) O I(0F, 1 for target estimation accuracy Eg. EF > Such an effective minmax optimization solver will be
(f))
O(EF) (11)
introduced in the following Section 3.2 In this way, Ve(x;) in (9) can be estimated as
To 9u.t Veg(ZE. wg.t) (ught Veg"(ar (ot, Vg.)).
(12)
Leveraging Iu.r, the gradient (F, can be estimated via (8)
9F., V. + Vig(Ie, TP) gv.t) (4F Veg (z.8F,))
(13)
We summarize the oracle for finding min (xr as in Algorithm| I which we term BLOCC, an algorithm designed for BiLevel Optimization with Coupled Constraints. Notably, our BLOCC algorithm can be seamlessly integrated with any MaxMin Solver (or min-max solver) that converges to the optimal solutions the max-min subproblems by achieving [II.. In the following, we present the convergence result of it allowing estimation eπror Eg: EF > 0 from the MaxMin Solver.
Theorem 2. Suppose that the assumptions Lemma 3 Run Algorithm 7 with some effective To and (yEt) ( respectively ((eg) and ((er) (optimal in inner MaxMin solver to find (Ughr a.t))
squared distance in Set m I TFA1
I. with some 1p,1 defined in Lemma 3 It then holds that
Σ IGG,
- =:)==-=))1==0111 + Per + 22gg)
(14)
Tn
The proof is available in Appendix D.1 Using the projected gradient G., = 17 (- xx) as the convergence metric for constrained optimization problems is standard [26]. The term O(> 7 arises from estimation errors using specific MaxMin Solver. Although the inner oracle can any one that achieves (1I, we value the computational effectiveness and therefore present particular efficient solver Algorithm 2 in the following section.
3.2 MaxMin Solver for the BLO with inequality CCs
In this section, we specify the MaxMin solver in Algorithm I By viewing L,(H, y; in 7 and LF yi x) in as L(u, y) for fixed given E X, we consider the following max-min problem
max min T( y)
(15)
R yey

---

which is concave (linear) in z and strongly convex in y. We can evaluate the dual function of [1, v) defined below by finding U (m)- D())= == min T(u,)) v) = T.,, ((())) (16)
Algorithm Subroutine on MaxMin Tτ)
I: inputs: initial points you Ho; stepsizes fh, 725 counters T,, TyTy Blue part is the version with acceleration; and red part is without. 2: for t ...... T do
JEY
where (u( = arg min yEy
migmin (H, y)-
According to Danskin's theorem, we have DD(p) V,LL(,, (f)) Taking either [1g(U, y; x) or L (p as L(u,y) y) for given E X, D(J) defined as (T6) exhibits favorable properties namely smoothness and concavity, with details illustrated in Lemma
update M++ via (17)
M++ Me
for by = ,,..., Ty do update yr,,,+1 via (18)
5:
D set 1/+,0 Yr
end for
7: update Me+1 via 19
D set y1+1 = ye,T,
end for
10 in Appendix D22 We can, therefore, ap- ply accelerated gradient methods designed for smooth and convex functions such as [52
9: outputs: (yT, PT)
At each iteration t. we first perform momentum-basee update step to update Me+ as
+ (14 - 14-1), with μl1 Mo- 1-
P44 = 14 +
(17)
To evaluate (D(4++4) = V vpllet (14++)))
4). with arbitrary small target accuracy >> 0, can run T, = O(I((-))))))))) PGD steps on L(M in y via
yt,t, = Projy (1.., m (y (..1,,)))
(18)
Defining the output after Ty iterations ye+1, since strongly convexity of [(p,,) ensures that PGD) converges linearly 19] Theorem 3.] it implies that [| - 97, (1r4444) (c)) We can conduct
M++1 = Projg (M++3 + 72 V ( 1t+1)).
(19)
We summarize this oracle in Algorithm@bastt based on an accelerated method 52
When we skip the
momentum update, i.e. setting M++ f4,. it is simple PGD method on D(u).
However, for convex function - (()), the standard results only provide convergence of the function
value, i.e. pER. D(u) D(4,) - 0. To establish the convergence of [1 *g(a))) and
|| ux we define the dual functions associated with inner problems 77) and (5) as
Dy(m) = min (1, x) and Dr(m) min Le(u, 21 yEY
(20)
yey
and make the following curvature assumption near the optimum. Assumption 5. There exist 6g. 6r and Cs., Cop 0 such that (x)())), p- 17(()) 20, 10, p1)33)1 Vp E B(rG(z)::g) VDr(H) + VDz (MT
(2la)
(=(()) 2 20rl1 11
Au (())
(21b)
It is worth noting that (- (-VDg(M) + V (1,,(()))) 1- -=)),,,;;(() [1]xx))))) and in a neighborhood of the optimal, the equality only happens at the optimal due the uniqueness of (z). The same argument applies to DF due to the concavity of the dual functions. Therefore, Assumption Jessenttallyy asserts a positive lower bound on the curvature the left-hand side term, which is mild as it only applies to the neighborhood of the optima (5(x) and (x)) It is also weaker than the local strong concavity or global restricted secant inequality (RSI) conditions.
2 0 holds for all due to concavity
By choosing Algorithm 2 with acceleration as the MaxMin solver. we provide the convergence analysis of Algorithm next, the proof of which can be found in Appendix D.2 D.2
Theorem 3. Suppose that Assumptions
i and the conditions in Theorem hold. Let > Lي
Cg S Cag by and < Caf dy. If we choose Algorithm 2 with acceleration as the inner loop and $2 2g imput T -((,)05), TF O(e and Tу (1n((==)) T. ((n((7) with proper constant stepsizes in Remark 3 then the iterates generated by the Algorithm Isatisfy IT.
1 T. - (, (.)))) O(Eg): II IF) (vF(2))) (G.))) OleF).
7

---

Theorem concludes the 00gg0.9) complexity for achieving Eg op-tima solutions for the constrained
concave-stronglly-oonee, problem max " minyey Lif, y), and so as for eF. For solving an e- approximation problem of BLO defined in 4. we solve min (x) with y (c--0.)) according to Theorem In this way, to achieve Et= 1- (z)]l= O + -2eF )2,)) by
14). we need Cg> Cp = O(e²), i.e. complexity (( for the MaxMin solver in Algorithm 2 and
= One ') = O(e 1.5) for the number of iteration in BLOCC (Algorithm Therefore, the overall complexity is (e)-) where O omits the In terms.
3.3 Special case of the MaxMin Solver: g((2,,) 3) being affine in and y = R
In this section, investigate special case of BLO with CCs where Assumption 5 aauommatically holds. Specifically, we focus on the case where gf is affine in y and y R i.c.
gx y)
f y g6(x).
(22)
In this case, fixing æ, taking Lg (1, y1 ) in and Lp(u, yi x) as L(u, v). (66 gives (g(H) = min (99((=),)) p) + (v. (((())))) + g( v) and
(23a)
yER""
(r(u) = min
- (0S(()) p) (u. (f(m)) fzz v) + +ol(,, v(z)).
(23b)
yER-y
When (I) is of full column rank, both (1) and Dr(u) are strongly concave according to Lemma
13 in Appendix so that Assumption 5 holds globally. Moreover, when applying PGD on -D (1)
and (I) strongly convexity guarantees linear convergence (Theorem 3.10 in (9). Therefore, even without acceleration, Algorithm with T,, sufficiently large performs PGD on ((x)) and converges linearly up to inner loop accuracy. Moreover, PGD on L(u, y) in also converges linearly, This motivates us to implement single-loop version (T, = of Algorithm 2
When both y and j are unconstrained, the analysis has been established in [19]. However, as constrained to Te in the Lagrangian formulation for inequality constraints, the convergence analysis in [19] is not applicable, and the extension is nontrivial due the non-differentiability of the projection. We address this technical challenge by treating RC as an inequality constraint and reformulating it as an unconstrained problem using Lagrange duality theory. This approach demonstrates that single-loop Ty = 1 update in Algorithm without acceleration, achieves linear convergence for the max-min problem (15). The detailed proof is provided in Appendix [D.3]
Thus, by selecting the single-loop version (Ty 1) of Algorithm 2 MaxMin solver in Algorithm [ we establish the following theorem with proof available in Appendix
without acceleration as the
D33
Theorem 4 (Inner linear convergence). Consider BLO with y = Rdy and gz, y) defined in (2). Suppose Assumptions 7 and the conditions in Theorem 2 hold. Suppose for any E X, there exist constants 9mn and Smax such that < Smin I m m (=(())) S (max(9((2)) S Smax oo.
Let 7 > 4111 If we choose Algorithm 2 without acceleration as the inner loop and impul Tg = ((m(c,))
Og ")) TF = (In(c=')) and T TF = 1 with proper constant stepsizes in Remark then the iterates generated by Algorithm 7
satissf |1(,,, To (y) (1x). p2 (xx) (.
O(eg); ( (TE) (1x)
O(e)).
Theorem 4 esstablishes for the first time, the linear convergence of a strongly convex-concave max- min problem with constrained maximization parameter. Similar to the previous analysis, we choose
y = (( 0.5) to solve the equivalent (2a) to e-approximation problem of BLO 4 To achieve 1 Et T-
(14). As the inner MaxMin solver convergences linearly, the overall complexity is Õ -155 where Ô omits the In terms. We summarized the overall iteration complexity of our BLOCC algorithm in different
(x)())) < the number of iteration in BLOCC T = O(yc-1) = O(c
settings in Table I
4 Numerical Experiments
This section reports the results of numerical experiments for three different problems: a toy example used to validate our method, an SVM training application, and network design problem in both
8

---

synthetic and real-world transportation scenarios. We provide sensitivity analysis and insights for pyper-parameter choices in Appendix G.I] In the two real-world experiments, we compare the
proposed algorithm with two baselines, LV-HBA 75 and GAM (22). The code is available at
https//gfttu...ommm///uuuu//gegn Lagranannnaa Bilevel
4.1 Toy example
Consider the BLO problem with an inequality coupled con- straint V(x) = {y V: x S 0}, given by
f(z, (z)) = é (x)
1
Rx, Myy x (ST, Vo. Fum, fu.. 12:1)
min =E[0,3]
+ 2 In ((44 2)) 2)2+10 +1)
+ cos(62) 6z)
with (1(z)) E arg min g(x,v) v) = (u - 2x)²
(24)
y))(x)
Problem 24 satisfies all assumptions for Theorem 2 and Theorem The lower-level problem is strongly convex and (z(2)) = x. Therefore, the BLO problem with incquality con- straint in (24) reduces to min (z, 1) In Figure 2 we plot the dashed line as the intersected line of the surface f(z and the plane f(r, vy (z)), and the red points as the converged points by running BLOCC with y 5 and with 200 different initialization values. It can be seen that BLOCC consistently finds the local minima, verifying the effectiveness.
y
Figure 2: 3-D plot of the upper- level objective f(x, y) of the toy example, with the line [xx v)ly shown in dashed red and the con- vergence points marked as red dots.
4.2 Hyperparameter optimization for SVM
We test the performance of our algorithm BLOCC with y = 12 when training linear SVM model on the diabetes (20) and four- class datasets
Methods
diabetes
fourclass
BLOCC
0.767 0.039
0.761 0.014
(Ours)
(1.729 1 0.599)
(1.922 0.108)
[29]. The model is trained via the BLO formulation (I) with inequal- ity CCs: see the details Appendix E We compared performance with two base- lines, LV-HBA 75 and GAM 72,, as they are the only existing algorithms address- ing coupled constraints and experimentally evaluated on SVM problems in their respec- tive papers.
0.766 0.039
0.745 1 0.060
LV HBA
(2.899 + 1.378)
(2.404 4 0.795)
0.721 + 0.047
0.715 + 0.056
GAM
(8.752 + 4.736)
(13.481 + 0.970)
Table 2: Numerical results the training outcome of our BLOCC in comparison with LV-HBA (75] and GAM (72). The first row represents accuracy mean + standard deviation, and the second row between brack- ets represents the running time until the upper-level objective's 's update is smaller than 1e-5
Table 2 shows that the model trained by our BLOCC algorithm outperforms that of GAM 722 significantly and is of similar level as that of LV-HBA [75]. We present some of the performance plots for the diabetes dataset as in Figure Looking into the test accuracy (left), our algorithm achieves more than 0.76 accuracy in the first 2 iterations, which is significantly better than the other ones. For the upper-level objective (middle), in the first few iterations, the loss decreases significantly under all algorithms, in which our BLOCC achieves the lowest results. Moreover, in Figure 3 (right), the lower-level optimum for LV-HBA was not attained until the very end. This indicates that the decrease of upper loss between 0-40 iterations in the middle figure may be due to the suboptimality of lower-level variables. For our BLOCC and GAM, the lower-level minimum is attained as lower-level objective g 0 in this case.
4.3 Transportation network design problem
BLO is particularly relevant transportation, where network planning must consider different time horizons and actors. Those problems are large-scale with a large number of upper- and lower-level variables and CCs, challenging the use of traditional BLO techniques, which involve expensive calculations of second order information. As a result, our final experiment considers network design problem, where we act as an operator whose profit is modeled as the upper-level objective that determined by the passengers' behavior, modeled in the lower-level. We considered three networks:

---

Test Accuracy
Validation Loss (Upper Level)
Lower Level Objective
0.80
BLOCC
BLOCC
3.0
LV-HBA
LV-HBA
Accarncy
0.75
GAM
GAM
2.5-
Loss
BLOCC
0.65
LV-HBA
2.0
GAM
0.60
1.5
75
50 Iteration Count
75
Ireration Count
Iteration Count
Figure 3: Test accuracy (left), upper loss / (ar, y) (middle), and lower loss g( y) (right) for the SVM on the diabetes dataset. The experiments are executed for 50 different random vain-validatiin splits, with the bold line representing the mean, and the shaded regions being the standard deviation.
N. 48
NN= =C= NZ 6,222
NN =6,, NV= 49,216 NC= 13,336, NZ= 129,256
=C=24, 24, NZ 126
Methods
Runtime (8) UL utility
Runtime (8) UL utility
Runtime (s)
UL utility
LV-HBA
3.51e2
1.53
/
BLOCC (Ours)-y = 2
2.02el
1.07
7.27e2
8,09
6.82e4
98,40
BLOCC (Ours)-7 = 3
2.00e1
1.69
8.50e2
10.37
6.42e4
111.39
BLOCC (Ours) = 4
2.01e1
1.71
8.68e2
11.04
6.70e4
138.78
Table 3: Results of the transportation experiment, both terms running time (Runtime) and convergenced upper-level objective value (UL utility, larger there better), with stepsize r) 1.66-4. We use "/" for algorithms that cannot converge within 24 hours of execution. NN (Number of Nodes) the number of stations in the network. Analogously, NV (Number of Variables), NC (Number of Constraints), and NNZ (Number of non-zero elements) are the number of optimization variables, constraints, and non-zero elements of the constraints matrix, respectively.
two synthetic networks of 3 and 9 nodes, respectively, and a real-world network of 26 nodes in the city of Seville, Spain. Further details about the formulation and the experiment be found in Appendix
In this experiment, only compare our BLOCC with LV-HBA 75 as the sole baseline, which is the only existing algorithm that addresses both coupled inequality of (z. y) S 0 and domain constraints V. GAM [72) and BVFSM [46] cannot handle lower-level domain constraints they rely on the hypergradient of lower-level and require LL stationarity in an unconstrained space. From Table 3 we that LV-HBA failed work efficiently, especially for large networks. This is mainly because the increased constraints render the projection step impracticable. Our BLOCC, in contrast, is much faster, and it successfully converges even with large real-world networks. We provided computational complexity analysis Appendix and we can conclude that our BLOCC is robust to large-scale problems.
5 Conclusions and Future Work
This paper proposed novel primal-dual lal-dddduls--seeeee penalty reformulation for BLO problems with coupled lower-level constraints, and developed new first-order method BLOCC to solve the resultant problem. The non-asymptotic convergence rate of our algorithm is (-5-2.), tightening to ((c-1.5) when the lower-level constraints are affine in y without any other constraints. Our method achieves the best-known convergence rate and projection-free, making more favorable for large-scale, high- dimensional, constrained BLO problems. Experiments on SVM model training and transportation network planning showcased the effectiveness of our algorithm.
10

---

References
[1] Seyed Mchdi Alizadeh, Patrice Marcotte, and Gilles Savard. Two-stage stochastic bilevel programming over transportation network. Transportation Research Part B: Methodological, 58:92-105, 2013.
[2] Aram V Arutyunov, Evgeniy R Avakov, and Alexey F Izmailov. Directional regularity and metric regularity. SIAM Journal on Optimization, (8(3):80-8335 2007.
[3] Dominique Azé and Jean-Paul Penot. Uniformly convex and uniformly smooth convex functions. In Annales de la Faculté des sciences de Toulouse: Mathématiques, 1995.
[4] Amir Beck and Marc Teboulle. A fast iterative shrinkage-threshoolinng algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183-202, 2009.
[5] Moshe E Aenva and Steven R Lerman. Discrete choice analysis: theory and application to travel demand, volume 9. MIT press, 1985.
[6] J. Bolte, E. Pauwels, and S. Vaiter. Automatic differentiation of nonsmooth iterative algorithms. In Advances in Neural Information Processing Systems, 2022.
[7] J Frédéric Bonnans and Alexander Shapiro. Perturbation analysis of optimization problems. Springer Science& & Business Media, 2013.
[8] Jerome Bracken and James McGill. Mathematical programs with optimization problems in the constraints. Operations Research, 21(():37-44, 1973.
[9] Sébastien Bubeck. Convex Optimization: Algorithms and Complexity. Foundations and Trends in Machine Learning, 2015.
[10] Luis Cadarso and Angel Marin. Combining robustness and recovery in rapid transit network design. Transportmetricca A: Transport Science, 12:1-26, 11 2015.
[11] Ennio Cascetta. Transportation systems analysis: models and applications, volume 29. Springer Science & Business Media, 2009.
[12] Tianyi Chen, Yuejiao Sun, Quan Xiao, and Wotao Yin. A single-timescal method for stochastic bilevel optimization. International Conference on Artificial Intelligence and Statistics, virtual, 2022.
[13] Tianyi Chen, Yuejiao Sun, and Wotao Yin. Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems. In Advances in Neural Information Processing Systems, Virtual, 2021.
[14] Benoît Colson, Patrice Marcotte, and Gilles Savard. An overview of bilevel optimization. Annals of operations research, 153(1)223--256, 2007.
[15] dean-Philippe Côté, Patrice Marcotte, and Gilles Savard. A bilevel modelling approach to pricing and fare optimisation in the airline industry. Journal of Revenue and Pricing Management, 2:23-36, 2003.
[16] Olivier Devolder, François Glineur, and Yurii Nesterov. First-order methods of smooth convex optimization with inexact oracle. Mathematical Programming, 146:37-75, 2014.
[17] Asen L Dontchev and R Tyrrell Rockafellar. Implicit functions and solution mappings, volume 543. Springer, 2009.
[18] Dmitriy Drusvyatskiy and Adrian S Lewis. Error bounds, quadratic growth, and linear conver- gence of proximal methods. Mathematics of Operations Research, 43(3):919-948, 2018.
[19] Simon S Du and Wei Hu. Linear convergence of the primal-dual gradient method for convex- concave saddle point problems without strong convexity. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 196-205 2019.
[20] Dheeru Dua and Casey Graff. Uci machine learning repository, 2017. Accessed: 2024-05-21.
11

---

[21] Laureano Escudero and Susana Muñoz. An approach for solving a modification of the extended rapid transit network design problem. Top, 17:320-334, 334, 12 2009.
[22] Chelsea Finn, Pieter Abbeel, and Sergey Levine. aodel-agnostic meta-learning for fast adapta- tion of deep networks. In International Conference on Machine Learning, pages 1 1126-1155 115,, Sydney, Australia, 2017.
[23] Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil. Forward and reverse gradient-based hyperparameter optimization. In International Conference on Machine Learning, pages 1165-11733 2017.
[24] Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning, In International Conference on Machine Learning. Stockholm, Sweden, 2018.
[25] Lucy L Gao, Jane Ye, Haian Yin, Shangzhi Zeng, and Jin Zhang. Value function based difference-of-cone algorithm for bilevel hyperparameter selection problems. In International Conference on Machine Learning, pages 7164-7182, Baltimore, MD, 2022.
[26] Saeed Ghadimi, Guanghui Lan, and Hongchao Zhang. Mini-batch stochastic approxima- tion methods for nonconvex stochastic composite optimization. Mathematical Programming, 155(1):267-305 2016.
[27] Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. arXiv preprint 1802v:880000224: 2018.
[28] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo, On the iteration complexity of hypergradient computation. In International Conference on Machine Learning, pages 3748 -3558, virtual, 2020.
[29] Tin Kam Ho and Eugene Kleinberg. Building projectable classifiers of arbitrary complexity. In Proceedings of 13th International Conference on Pattern Recognition, volume 2, pages 880-885, 1996.
[30] Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola. Kernel methods in machine learning. 2008
[31] Mingyi Hong. Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic. SIAM Journal on Optimization, 33(1), 2023.
[32] Kazufumi Ito and Karl Kunisch. Lagrange Multiplier Approach to Variational Problems and Applications, volume 15 of Advances in Design and Control. Not specified in the provided information, 2008. Includes bibliographical references and index.
[33] Kaiyi Ji, Mingrui Liu, Yinghin Liang, and Lei Ying. Will bilevel optimizers benefit from loops. arXiv preprint arXiv:2205. /4224, 2022.
[34] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Provably faster algorithms for bilevel optimization and applications to meta-learning. In Advances in Neural Information Processing Systems, 2020.
[35] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning, virtual, 2021.
[36] Sham Kakade, Shai Shalev-Shwartz Ambuj Tewari, et al. On the duality of strong convexity and strong smoothness: Learning applications and matrix regularization. Unpublished Manuscript, http://tic. uchicago.
pdf. 2 2(1):35, 2009.
[37] Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear convergence of gradient and proximal- gradient methods under the polyak- lojasewicz condition. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16, pages 795-811, 2016.
12

---

[38] Narendra Karmarkar. A new polynomial-timg time algorithm for linear programmingg In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 302-311, 1984.
[39] Prashant Khanduri, Ioannis Tsaknakis, Yihua Zhang, Jia Liu, Sijia Liu, Jiawei Zhang, and Mingyi Hong. Linearly constrained bilevel optimization: A smoothed implicit gradient approach. In International Conference on Machine Learning, pages 16291-16325, 2023.
[40] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang. A near- near-optimal algorithm for stochastic bilevel optimization via nouble-momentum. In Advances in Neural Information Processing Systems, Virtual, 2021.
[41] Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert Nowak. A fully first-order method for stochastic bilevel optimization. In International Conference on Machine Learning, pages 18083-18113, 2023.
[42] Jeongyeol Kwon, Dohyun Kwon, Steve Wright, and Robert Nowak. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. In International Conference on Learning Representations, Vienna, Austria, 2024.
[43] Junyi Li, Bin Gu, and Heng Huang. A fully single loop algorithm for bilevel optimization without hessian inverse. In Association for the Advancement of Artificial Intelligence, pages 7426-7434, virtual, 2022.
[44] Bo Liu, Mao Ye, Stephen Wright, Peter Stone, and Qiang Liu. Bome! bilevel optimization made easy: A simple first-order approach. Advances in neural information processing systems, 15:17248-17662 2022.
[45] Risheng Liu, Xuan Liu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. A value-function--asee interior-point method for non-convex bi-level optimization. In International conference on machine learning, pages 6882-6892, 2021.
[46] Risheng Liu, Xuan Liu, Shangzhi Zeng, Jin Zhang, and Yixuan Zhang. Value-function-basee sequential minimization for bi-level optimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023.
[47] Risheng Liu, Yaohua Liu, Shangzhi Zeng. and Jin Zhang. Towards gradient- bilevel optimization with non-convex followers and beyond. In Advances in Neural Information Processing Systems, volume 34, pages 8662-8675, 2021.
[48] Songtao Lu. Slm: A smoothed first-order lagrangian method for structured constrained noncon- vex optimization. Advances in Neural Information Processing Systems, 36, 2023.
[49] Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient-based hyperparameter opti- mization through reversible learning. In International Conference on Machine Learning, pages 2113-2122, Lille, France, 2015.
[50] Patrice Marcotte. Network design problem with congestion effects: A case of bilevel program- ming. Mathematical programming, 34(2): 142-162, 1986.
[51] Akshay Mehra and Jihun Hamm. Penalty method for inversion-free deep bilevel optimization. In Asian conference on machine learning, pages 347-362, 2021.
[52] Yu Nesterov. Smooth minimization of non-smooth functions. Mathematical programming, 103:127-152, 2005.
[53] Alex Nichol, Joshua Achiam, and John Schulman. On first-order meta-learning algorithms. arXiv preprint arXiv: 1803. 02999, 2018
[54] Norbert Oppenheim. Equilibrium trip distribution/aasiigmment with variable destination costs. Transportation Research Part B: Methodological, 27(3):20-2177, 1993.
[55] Fabian Pedregosa. Hyperparameter optimization with approximate gradient, In International conference machine learning, pages 737-746, 2016.
13

---

[56] Fernando Real Real-Rojas. Diseño de redes de transporte: integrando la competencia modal con técnicas de optimización convexa, BSc Thesis, King Juan Carlos University. June 2024.
[57] Fernando Real-Rojas, Victor M. Tenorio, and Antonio G. Marques. A sparse nonlinear approach for designing hub-and-spoke air transportation networks. In Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, 2024.
[58] R. Tyrrell Rockafellar. Convex Analysis. Princeton University Press, Princeton, NJ, USA,
1970s.
[59] Andrzej Ruszezyúski. Nonlinear Optimization. Princeton University Press, 2006.
[60] Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots. Truncated back- propagation for bilevel optimization. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1723-1732, 2019.
[61] Han Shen and Tianyi Chen. A single-timescale analysis for stochastic approximation with multiple coupled sequences. Advances in Neural Information Processing Systems, 35:17415- 17429, 2022.
[62] Han Shen, Quan Xiao, and Tianyi Chen. On penalty- bilevel gradient descent method. In International Conference on Machine Learning, Honolulu, HI, 2023.
[63] Han Shen, Zhuoran Yang. and Tianyi Chen. Principled penalty- based methods for bilevel reinforcement learning and RLHF. In International Conference on Machine Learning, Vienna, Austria, 2024.
[64] Ankur Sinha, Samish Bedi, and Kalyanmoy Deb. Bilevel optimization based on kriging approximations of lower level optimal value function. In 2018 IEEE congress on evolutionary computation (CEC), pages 1-8. IEEE, 2018.
[65] Ankur Sinha, Pekka Malo, and Kalyanmoy Deb. A review on bilevel optimization: from classical to evolutionary approaches and applications. IEEE Transactions on Evolutionary Computation, 22(2))276-295, 2017.
[66] Daouda Sow, Kaiyi Ji, and Yinghin Liang. On the convergence theory for hessian-free bilevel algorithms. volume 35, pages 4136-4149, 2022.
[67] Bradly Stadie, Lunjun Zhang, and Jimmy Ba. Learning intrinsic rewards as a bi-level optimiza- tion problem. In Conference on Uncertainty Artificial Intelligence, virtual, 2020.
[68] Luis Vicente and Paul H Calamai. Bilevel and multilevel programming: bibliography review. Journal of Global optimization, 5(3):291-306, 1994.
[69] Gerd Wachsmuth. On licq and the uniqueness of lagrange multipliers. Operations Research Letters, 41(1):78 80, 2013.
[70] Quan Xiao, Songtao Lu, and Tianyi Chen. A generalized alternating method for bilevel learning under the polyal-tojasiewicz condition. In Advances in Neural Information Processing Systems, 2023.
[71] Quan Xiao, Han Shen, Wotao Yin, and Tianyi Chen. Alternating implicit projected sgd and its efficient variants for Artificial Intelligence and Statistics, 2023.
bilevel optimization. In International Conference on
[72] Siyuan Xu and Minghui Zhu. Efficient gradient approximation method for constrained bilevel optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, 2023.
[73] Haikuo Yang, Luo Luo, Chris Junchi Li, Michael Jordan, and Maryam Fazel. Accelerating inexact hypergradient descent for bilevel optimization. OPT 2023: Optimization for Machine Learning, 2023.
[74] Wei Yao, Haian Yin, Shangzhi Zeng, and Jin Zhang. Overcoming lower-level constraints in bilevel optimization: novel approach with regularized gap functions. arXiv preprint M4Xiv:240600199 2024.
14

---

[75] Wei Yao, Chengming Yu, Shangzhi Zeng, and Jin Zhang. Constrained bi-level optimization: Proximal lagrangian value function approach and hessian-free algorithm. arXiv preprint arXiv:2401. 16164, 2024.
[76] Jane J Ye and Daoli Zhu. Optimality conditions for bilevel programming problems. Optimization, 33(1):9-25 27, 1995.
[77] Jane Ye, Daoli Zhu, and Qiji Jim Zhu. Exact penalization and necessary optimality conditions for generalized bilevel programming problems. SIAM Journal optimization, 7(2):88-5507 1997.
15

---

Supplementary Material for "A Trimal-Dual-Assissee Penalty Approach to Bilevel Optimization with Coupled Constraints"
Table of Contents
Preliminaries
16
Analysis of the Penalty-Basee Lagrangian Reformulation
B. Proof Lemmal
17
B.2 Proof of Theorem
17
C Analysis of the Differentiability of Value Functions
18
C.I Proof Lemma
19
Lemma
D Convergence Analysis of the Main Result
20
D. Proof of Theorem
20
D 2 Proof of Theorem
................................
21
D.3 Proof of Theorem
24
E
Applications to Hyperparameter Optimization for SVM
29
F.1
Problem formulation
29
E.2 Experiment details
30
F
A Applications to Transportation Network Planning
31
F1
Problem formulation
................................
31
Numerical results for 3-node network
................................
34
?
Numerical results for the 9-node network
................................
35
44
Numerical results for the Seville network
36
G Sensitivity Analysis
37
Sensitivity analysis for the toy example
37
Sensitivity Analysis for the 3-node network
37
H Analy Analysis of the Computational Complexity
38
H.1
Complexity comparison
38
Complexity analysis of BLOCC
39
A Preliminaries
This section will provide some preliminaries for our subsequent analysis. Definition 5. For a convex function h. : Rdd -> IR whose domain is Q C Rd,, the Legendre conjugate of Q* R. is defined as:
h" h²())= == sup ((''g) - h(g')) inf (4'' + h(g')}
I'EQ
9'00
4q E Q- (q E Rd sup (g) .q) h(g')} < co).
y'EQ
Remark 1. When is strongly convex in Rdd, it is lower bounded and therefore Q = TR...
Lemma 4. Suppose h : Rdd, + R is 14.1- -mooth and on -trrongyy convex and its domain Q S Rd is convex, closed and non non-empty.
1. If Rd, the gradient mappings Vh and Vh are inverse of each other (158); and
ht R + R is 1 --moothh and -strongly convex (Proposition 2.6 (3)
04
T,
16

---

2. If C CR, h² is 1 -mooth ((36)) and and convex (Theorem 4.43 (32)). on
Lemma 5. Suppose Q S Rd is convex, closed and non-empty, h R24 + R. is strongly convex on
2 he RR gd is convex
and is finite, and {g Q ht S 0} non-empty.
1. The problem li h(g) has unique feasible solution.
2. When linear independence constraint qualification (LICQ) condition holds for (4))), the Lagrange multiplier for the problem m (-1e(a) h(q), i.e. solution to the problem max NER mingge h(q) + (n, (g)) is unique (69).
Lemma (Lemma 3.1 in 19. Suppose Q S Rd is convex, closed, and nonempty. For any q1 R& and any 92 Q it follows that
(P)oja((1) 92, (roje(()) -a) S 0.
(25)
In this way, take q1 = q3 ui for any q3 Q and denote q3 gradient update in direction g with stepsize 1), we have,
Proje ng)) as projected
1 n-9 q2, 43 7- (g3 1)
(g, 13 - Q2) S
43), V4z, Q3 E EQ
(26)
Lemma (Theorem 3.10 9.. Suppose a differentiable function 14,1- and ah, strongly convex. Consider the constrained problem mingeR h(g) where Q is non-empty, closed and convex. Projected Gradient Descent with S 1g.1
converges linearty to the unique q" arg minge h(g):
Projg(q v)Vh())) - q0 s (1 12 901l S (1 -
Ag Q. (27)
B Analysis of the Penalty-Base Lagrangian Reformulation
B.1 Proof of Lemma
According Lemma for any fixed there exists unique ug Therefore, according Lagrange duality theorem, for any fixed x, the primal problem
min g( v) s.t. (x SO yey
min g(z, v) + (z). yEy
As g(z, y) is ag -trongly convex in y and gx y) is convex in y. we know g y) + ((g g (x, y)) ag -trongly convex in y. where the modulus a, is independent of æ. Therefore, according to Appendix F and G 337. and Theorem 3.3 in [88, the quadratic growth in statement can be concluded.
As g(x, y) strongly convex y. and (*)) is non-empty, closed and convex under assumption 3. there exists unique solution (1r) such that g(2, v, (x)) = v(x) by Lemma 5 In this way, if * Vi (x)
and y E Vxr), we have g(r,, y) >> v(x). This completes the proof of statement 2.
B.2 Proof of Theorem
We know from Lemma that g(2, y) v(x) 2 2 1y (z) and g(x, y) = v(x) ) if and only
y = (z). This squared-distance bound follows Definition in [62]. Under Lipschitzness of f(x y) with respect to y, finding the solutions to the e-approximate problem in 4 is equivalent to
finding the solutions to its penalty reformulation
min (x. (Xxx X3'gg (a.w))))
f(z, + +igla, v) v(x))
(28)
with = Occ-0.5) following Theorems and in 62
Moreover, jointly finding solution for (z, y) in 28 is in equivalence to finding solutions in
min min f(r, v) + (0(a, v) v(x)).
(29)
xEA (E)(x)
17

---

The proof of this equivalence is as follows. Suppose (zo. vo) (X'x) gz(z.)) - 0) being v(x)). We
solution to (28). Suppose for any 2 E X. ff(() E arg mig) know that for any a E X. y E (xx), it follows that
[wo, y0) + 7lg(20 vo) - r(zo)) SI(,, + + u(x)
<ʃ ((((,))) v(x))
This means any solution (28) is solution (29 On the other hand, suppose xo X, (z))
(z0) is a solution to 29) We know that for any (r, v) E (x x of 5 0}.
f(Fo, (xo)) + rl( (( v(x))) <fx (x)) + Hal,, - v(x))
s {(x,))) + 70((,, v(x))
This means any solution to 29 is a solution to (88
Besides, we know [(x, y) is f.1 -mooohh g( y) is ag -trongly convex in y. By the definitions of strongly convexity and smoothness, we know for fixed x, for any y1,92 E V.
f(z,81) + 7 - v(x) - f(x, - +(o(1, v(xx))
=f (1) f(x, + (0( y1) 92x,
41 ll v2ll2 + V v (z. U2) 11 122) +
0gg Ily1
f v1 12)
y2
2
2
700g - 4111
=(Vyf(z,,)) v2) + v 1 x. U2), V/ - y2) +
13 32112
(30)
2
This proves that (x, y) +) y) -(xx)) is 170g-11 >-srongly convex in y. Moreover, accorring
to Assumption 2 the constraint (2, w) is convex in y, and minyer(c) ((r ((g(x,))) v(x)) is equivalent to its equivalent Lagrangian Dual Form [59]
max min f(x v) + 700(,, u) - v(x)) + (p, (x, v)).
(31)
HERR yEyr
Therefore, 28 can be recovered to 2a and this completes the proof.
C Analysis of the Differentiability of Value Functions
Lemma (Theorem 2.16 in 32 Suppose h(a, y) is strongly convex in y E y and is Lipschitz with respect to E X. (x, y) is convex in and is Lipschitz with respect to and both y and { E y: ( v) S 0} are non-empty. closed, and convex. For the problem minye (vEVa (a.))))) h(x. v). the unique solution vi, (a)) and unique Lagrange multiplier 15, (x). defined as
(x), pT (x)) = arg max min h1x,, v) + (1,16 (x) HERt yEy
(32)
is Lipschitz in In other words, there exist La. 2 0 that, for all æ1, E X,
(1))
(z s T x2ll-
Before proving Lemmas 2 and would like to introduce more general form.
Lemma 9. Suppose y and {y y 10 v) 0) are both non-empty, closed and convex, h(, y) is jointly smooth in (x. y) and is strongly convex in y, 1e y) is convex in y, and both h( y) and he (a y) are Lipschitz with respect to x. Then we have
un(ar (2)
min (xx, y) s.l. (x S wEY
differentiable with its gradient as
Vel(z) = (z, af (x)) + (4(()) (2 (x))).
(33)
where (x), 11, (z)) defined in 32 are unique.
18

---

Proof. We prove this using Theorem 4.24 in 77.
As h(r, y) being strongly convex in y, Condition in Theorem 4.24 in 7 is satisfied and the solution sets are of singleton value (v% (1), (2)) according to Lemma
ii) Moreover, the smoothness of lt(z, y) guarantees Robinson's 's constraint qualification [2]. which implies the directional regularity condition (Definition 4.8 in 7 for any direction d (Theorem 4.9. (ii) in (7).. This guarantees Condition 2 in Theorem 4.24 in can be satisfied for all directions d.
iii) Additionally, under the Lipschitzness of h(z, y) and (x, y) with respect to x, y (z). .4, (z) are Lipschitz according to Lemma This implies condition 3 in Theorem 4.24 in [7] holds.
In this way, all conditions in Theorem 4.24 in 7 hold and it gives the gradient in 33 for unique (x), 1% (I)). This completes the proof.
C.1 Proof of Lemma
Proof. The problem minyey g(z, v) s.t. g v) < 0 fits in the setting of Lemma by taking Wx,, y) = g(,, y) and (x. y) = g (x Therefore the derivative (8) can be obtained accordingly. Moreover, for any L2 E X,
[o((11) ve(
Va
(g((22) (22
(
([(g((z)) Vag - (1](21), Veg
(x1) (22 - (1)(22), v (z2 (22, (x2))))
5 (1.11 + () r2ll + (0))(()) 1 (x2) + b (lll 16 (12))
63 (06.1 11 + T.y) + 1g L |x1 x211.
where (a) follows triangle incquality; (b) leverage on the Lipschitzness of Vg, de and V in and
the upper bound for (I)(())] and (c) uses the Lipschitzness of vg (x) and 1 (2) from Lemma 8 As the bound loose due to the use of triangle inequality, we can conclude that v(x) is -smoohhh where 4,,1 ((1 + (1) + 19 lg (gg))
C.2 Proof of Lemma
By assumption, f(x, y) is 11.1- and g(z, y) is ag -trongly convex in y. We know f(x, y) r(g(,, y) - v(z)) is (ya,- (f,1))ssnnnn -stronlyyy convex when > Im as discussed in (30). Moreover, constraint of v) is convex in y by Assumption
ag 2 In this way, the problem
min
f(x, v) + 7(g(, v) - -xx)) DEy
s.t. g. 0
features strong convexity according Chapter 4 in 55 and it equals
Fx(x) = max min f(r.8) + ((((2 y) v(x)) + (m. (x,
"ER yey
Considering the smoothness of u(x)) presented in Lemma 2 all assumptions in Lemma are satisfied. Therefore the derivative ( can be obtained. In addition, for any 3 B2 E X,,
VF(z)) V
= / (ad ||
(1)(11) [olf - (1F
19

---

+ (21), (22, (1)()2),
(1 + 71.11 + B 22ll (m)) (f((2)))) + ,,, x2l
-
+ + - (x2))
2 (( + N4,1 B 1) LF) + N4,1 (g =121,,
where (a) follows triangle inequality; (6) leverage the Lipschitzness of f.. Vg. 9 and Vg in
and the upper bound for [ |: and (c) uses the Lipschitzness of (z) and (z( from Lemma
8 As the bound is loose due the use of triangle equality, we can conclude that F F(x)) is smooooo
where 1r,1 S (4,, + 74.11 + Brlge + Lp) + yo, + I...
D Convergence Analysis of the Main Result
D.1 Proof of Theorem 2
Define the bias term b(zz) of the gradient VF,, (x,) as
b(z)) =VVF gF.t
Val((r (x:)) + (Vog(VV
(z(1)) + (Mz(at), v ))
+ (u)(),) Vag (f (zx)))
-(vaf(an ( + Vag(zeNBF1) (Veg(zr, + (ugie (x.. ));))))
(u Vag (xL )
In this way, we have
(0)
[((z,]] 2 V.f(zr, vi:)) (x)(
+ 7 | 1 (
Vag(ze. 32 (x,))
(1] (2,, (
(M (z:), (ag m
Ngght (
+ 11 V
- (6 (z:),
(x..
+ 11 (x:), V (z:, o7 ( V (Z+, (+((x))))
(
2/1211uf, - (x + 1 1,.11177 vx(())) 2.1 1 (xx))
+ + To (1((x) + Bybg .lll v (()0))
1oralll - (x:)))) Brl 1 11 18 (2:))
(1 + N6,1 + Brlg a) I v + 1g oll "f(xx))
+ (1.11 + Bylgrggggg.o 1) 11 To - (x,)) + lgg.ollugg - - ((())))))
where (a) uses triangle inequality, (b) relies on the Lipschitzness of V.. Vy, gf, and Vg in Jr, the upper bounds for [,)](())] and (x)], and Cauchy-Schward inequality, and (c) is by rearrangement.
20

---

Furthermore, according to Young's inequality. it follows that
[b<x))]1 <2 (1111 + Ma + Bplg 0) oll To - VF, ll + low.ooollll..... Tr. -
27² (10 + Bg
V (xeel + fg olllig 16 (
=0 + kg)
According to Lemma (z) is Ismothooth in X. In this way, by the smoothness, we have
F(a++1) SF(x) + (VF(T)), 2x+1 - zz) + [F1
Tt+1 xel?
2
<F(x)) + (gF 2(41 xz) + 2n | - xz + (b(x)), 2141 - xz),
(34)
where the second inequality is by 1 S T..1 and (z((z)) = 9r. - b(xz).
The projection guarantees that Tt-1 and It are in X. Following Lemma we know that
(9F4, - nn) S - -| - =x,11²
Plugging this back to 34). follows
1
F(I++1) (z()))
| - 2el1? + (b(xz), 21+1 F(+11 -x1)
2m
<F(x)
|1111 - xel1² n)|b(())))' -
-|10,-111 - xe11
21j
4n
1
=F(x))
|| - zell nlbb(z))11
4n
where the second inequality is from Young's inequality. Telescoping therefore gives
T-1
4 - (F(zo)
4
Σ
|G,, (x,))
<
l1b(,;)11
7/
T
=0(n 1 + 0 72egg
=0 + 22 p20g
where last equality comes from 0(>-1) 17 S ب and 1F $ (1 + Ng.1 + Bplge (1
LF) + y0,,1 + L = O(y). This completes the proof.
D.2 Proof of Theorem 3
To restate, are viewing Ly (M, yi x) in 7 and LF x)) in ( respectively as L v) and considering the following max-min problem
max min L(u,v) v). MER4E yE
In 16 we defined the minimization part
P(4) = min L(p,y) = ((,,((()) where (1](1) = arg min [(,, y)-
yEY
yEy
To evaluate D(u) for (g v5 in 7 and Le(u,r in 5 L(u, y), we define the following mappings for fixed E X:
(p) = arg min Ly y x)x,,
(35)
yEy ( arg min Lp 2p) ye
(36)
21

---

In this way. for (1, yi v: in 7 and (u, v; in 5 respectively as L(,,, y), D(p) defined in T6] equals Do and Dr(u) respectively. where
D²(u) = min (M, vx Lo(u, (u, Gin (u);z), PEy
(37)
Dr(1) = min F y: Lr( D (()))
(38)
yED
In the following lemma, we show that D(μ) exhibits concavity and smoothness, which are favorable properties for conducting gradient maddeent-based algorithm.
Lemma 10 (Smoothness and Concavity of D(u)). Suppose all the assumptions in Theorem 3 For fixed x E X, the following holds
hold.
1. Up.g (1) in 335 and Vg, in 36 are respectively J and - Lipschitz to p.
ag
2. Dz(n) in 37 is concave and le -moothh and Dy (p(p) in 38 is concave and L0 7ag -lpa
0g
smooth.
Proof. To restate,
(,(,,,,)) =9(2, v) + gz,,0), Lr(,, y; [f y) + (g((, y) v(x)) (u. (f v)).
Under Assumption 1 and for fixed x and given μr, Lg (1, y: x) is -g--trongly convex and (1g.1
(ualllys.))smmotee in y, and (p x; is (1aa --st-stoyy convex and + + || ()) smooth in y. Therefore, know (p (n) in (35) and V (4) in (36) respectively 1 and
0g
you-lya Lipschitz to by directly quoting Theorem F.10 in 17 or Theorem 4.47 in 332 This proves the first part of the Lemma.
1
For the second part, the concavity of (n) and DF(p) can be directly obtained by Lemma 2.58 in [59 as (1, yı x) and (1, yi are both convex in y.
Moreover, following Theorem 4.24 in 7. have
(D) = vV Lhglu = (2, (1)), VDr(H) = Vubr(u = g ( (1)))
As (x, 3) is byr o-Lipschitz by Assumption for any fl: E B::::
[[VD)(m)) - VD2(12) =11g u (1z)) (x, ynn (12))
Slo (WI) Vn (H2)) S 19:00 1000 p1 -
Hell.
(39)
0g
and similarly,
-VDp(H))-DDr(a)))) = -19 Vz. (x, (..
10:00
Slgeoolb ( S
|| -
H2ll-
(40)
70g - 11,1
We conclude that (1) and (F are respectively
4..0 and 4"0 smooth. yo,
og
1,11
In Algorithm we are implementing an accelerated projected gradient descent on -D))) where the gradient bias is controlled by T... Following (16), the following lemma presents the convergence analysis of the accelerated method on smooth and convex functions.
Lemma 11 (Section 5 and in (16). Suppose D(μ) is concave and Ismooohh Consider the constrained problem max "E D(u), At iteration i = O,..., T - 1, perform accelerated projected
gradient update with stepsize S T. and initial value Ho p
1-1 (Mu 14-1)
M4+ =14 + ++2
(41a)
22

---

Me+1 = =rrojed (M+= + 79++
(41b)
where 9++3 is an c..5 cpprrpppoee to (D((++)) ) satisfying 1196+3 - V D(4++))] = O(e1.)), for
given accuracy > 0. Denote D max HERR D(u), performing T ((E-0.5) iterations leads to
D(uT) D(e))
Remark 2. The domain for μ E R4 can be replaced by any closed, convex, non-empty domain.
In this way, we are ready to proceed to the proof of Theorem
proof of Theorem 3 Algorithm 2 solves and (5) by taking Lg(1, y1 2) and LF v; respec- tively as [(,,)) y) and run respectively iterations 7 equals Ta and TF
We begin our proof with analysis on 7 The accelerated version of Algorithm solves this by taking Tig (1,, y3 x) in (7) with fixed x E as (u,.y)
Fixing P : steps 6 are -tep projected gradient descent in y. As Lg x) is (l,,1 + ())
smooth and 0g 0g-- strongly convex in y, taking the inner loop stepsize !n as fl..1 S 1..1 tlg,, ensures
linear convergence according Lemma
7 Choosing T2 (In(((e = (1n((()) leads to
|| - Vi (e+++ )1l = O(e),) 1.5)
for target accuracy > 0 where (1) is defined in 335
In step 7. we update μ with a projected gradient descent step which takes (pL 3c+11) (x, V. The estimation bias is bounded by
gx 11+1) as an estimate of V (
V [p Bt+1) - VD (14 f )1 91 1c+1) (z, Viu ( ))1l
5/g 0111-11 Vi (Me+ f)ll = O(eg.5).
By Lemma 11 we can conclude the complexity is ((g=0.6) for conducting the accelerated version of Algorithm on Lg (1, x)) in 7 as L(J,,) v) to achieve
(o(n((=)) (y(TT,) = O(e)).
(42)
Similarly, the complexity of the accelerated version of Algorithm 2 for (5) is Ö(e i.e.
Dr (=)p - D=(ur>) = O(ep).
(43)
In the following, we are going to show that 42 and 43 and respectively sufficient to bound ||RT, - -1g(2) (x)] and || (x))l considering Assumption is satisfied.
As Dy(H) and DF (M) are both concave in and E 1R is equivalent to e 2 0, the problems
max (n(u) and max Dr(H) MERS
MER+
are respectively equivalent to the unconstrained problems max
= D=(x))=D,(m)) + 2/5 and max DP(1) = DF(N) + ATF dE
nERdc
with the Lagrange multipliers Ag, Ar being non-negative and finite in all dimension, i.e. 0 < >g oc. 0 < App < 00 according to Lagrange duality Theorem. Moreover, it is well known (Chapter 4 in
59 that the Largrangian terms respectively equal zero when the problems attain the optimals. i.e.
A M (z) 0 and (z =0
(44)
Moreoever, the first-order stationary condition requires (Dy(H(())) VDG(G(()) A 0 and (D)( = VD (1))()) + Af = and therefore
(D)(H(())) -> and VDF(H) -AP.
(45)
23

---

In this way, for for all (11; n R7.
Do(1)(t)) D2(u)= = (VD -(((()))) 1)dT
1 - =((((z)
((u())(((()) -+)), T(u - Ms (z))))
(VPg( (2)), -1((2))
Cs |1 - (g(()) 1 (VDg(rg()))
C5 || - (z)) + (g - -x;(x))
(c) C5
|| 12 (x))²
2
5 and the fact that the p. (G(()) (((((z):&,) n R4 implies f [pg (z) μ³ (()G()); n R4; (6) solves the integral and uses Ag = - (z)) in
where (a) uses (1a in Assumption
and (c) follows from the fact that (A, ug
0 in (44) and μ,, Ag 2
Analogously, for all E (z) n 184, it follows that
Cór 2
Dz(μf(z)) DF(H) 2
1l (x)))
Cāg In this way, for arbitrary Eg 2 8g. the complexity of Algorithm 2 to solve 7 is Ö(EF 0.5, i.c.
||T,, (1)x1) O(Eg),
and llyr, (x))l Sllvr, vy (pr, + 1|r,,- - (=)=1|? <(1/ag + 1))l
At each iteration / in Algorithm
x It, and the output (IT,, is chosen as (Ug.3 Hg,1)
Ct ig applying Algorithm with complexity Ō 0.5 to
Similarly, to solve (). for arbitrary < achieve (33) can achieve
2
UPTE (x)) O(er),
and [1yT (x)1 SlvtT (fTr: IIIt - (z)1
<(1/ar + 1)1lITe - (=)))) O(er).
At each iteration in Algorithm Fr, and the output (UTT (TT) chosen (.))
This completes the proof.
Remark 3. Under the same assumptions as in Theorem we cam choose 17
< (1g,1 + Lge 1)-
19. 10 as stepsizes for running the accelerated version of Algorithm 2 to solve
and
(F,1 < (4 + 19.1 + (gr,)) 1F,22 S 3ag=1e2 as the ones for 5
-go00 Ly-.0
D.3 Proof of Theorem
In this section, we consider
g0 v) = gf(x) -
(46)
being affine y, and y IRdy
Therefore, for fixed x, taking either (, yı x) in (7) or LF(u, ys x) in (5) as L(p, y) fits into a special case of strongl-ccomex-concon saddle point problems in the following form:
max min L(u, v) - (p) + v Au + h2(2)
(47)
PERdo yERdy
24

---

where h1 (p) is smooth and linear (concave) in and h2( is smooth and strongly convex in y. Specifically, for (1, v) as L(,, y), hn (14 (z)(z), p)) O-smooth and linear (concave),
A = (xr), and h2(v) = g(z,)) v) is sm.o-shooth and -trrnyyygly convex.
gradient descent step since } = _d The following lemma summarizes the error of ||2/ yz [1m)] where (z) is defined in (16) and the update μe+1 and yt+1 is the non-acee con-accelerated version of Algorithm with Tg =
In this context, performing PGD L(,,))
y is equivalent
Lemma 12 (Update of 11
Consider the problem (47 where of column
rank and (2( is Oh2 -strongly convex and f,, -smooth. (1) defined in
16 satisfies
va(n)) Vhz(
(48)
where (y) is the conjugate function of hy(y)) by Definition 5 At iteration t, pt: Yt are known, and conduct Vt+1 Yt 1 V ( Vt), a gradient descent step for L(ps, y) in y. This gives
1y++1 Vh2(- (1 (11-110 11/00 /2000//22)llrr
4h2( Apx)
(49)
when Tn. S Tag.
Additionally, given another 1,++ we know
mmax | - Mell.
|1u41 Vhz( 103(-A44+1] S(1 7mah//2)llls V (44)l
(50)
Proof. Recall the definition (1) = arg miny L(p y) following 16. The first-order stationary optimality condition requires that for any given μ4, holds
V L(,, v(A) Ap + (h2(v,) =
Vh2(k[)) AAu.
As the mapping Wh2 and Vhh are the inverse of each other according to Lemma 4 for any
Vhg(-Ap
At iteration 4, conducting gradient descent step on L( y) gives yr+1 = - 111 V, ((ee We). As
[(,,3) y) is -Sronnly convex and 1 -smooth in following Lemmal take 1/1 S Ing.1 we have
|| Vhi(-Apu
(olmmman,/)lell Vhi(-Ap))l-
Following triangle inequality, we also have
Bt+1 (r][-Apc+1]l
Mt+1 -
+
(hi(-AAt+1))
Tmax mmax(A) || - uell
<(1 - m,,h/ -
(51)
where the second term comes from the smoothness of the conjugate function (see Lemma
In 50 the update behavior of llye Vhi((-Appp depends on || pell Therefore, we are interested in the the update behavior of || - [eg where #4+1 = Projg (Hr V I(4, Mc+1)) as in the acc-acceled version of Algorithm 2 with T,, 1. Before proceeding, we would like to look into the properties of ((4) defined in (16) under the setting (47)
((4) = min (1) - (p) + h2(2)
(52)
where h1 (a) is smooth and linear (concave) in μ., h2(y) is smooth and strongly convex in y, and is of full rank in column. The next lemma shows that D(#) features strong concavity and smoothness.
Lemma 13 (Smoothness and strongly concavity of D(n)). Suppose h is concave and L, 1-smooth, h Ohe -trongly convex and 1hg,1 -smooth, and is full column rank. Then D(1)) in (52) satisfies
D(u) -b1(1)) hxp(-Ap,
and is (A) -strongly concave and (la,,1 + (msx(A) smooth with respect to p. or2
25

---

Proof. Following Definition 5 we have
D(u) h1 (1² - (-(-Au)
where (1) is - -tronggl convex and For all 11,
-mmoohh according Lemma 4
Ing.1
aing
-D(1) - (-D())) =h2( Ap1) 1 + h1(11)) 11 (12)
0hg( (-1u)) A11 + Ap2) + 1/1 11A
2 0 Ap2
Ap2ll
+ (hh (u2), ,1 -22)
(VD(22)) 11 -22) + o (A)/1/a.l 11 juz 11?
where the first inequality follows the strong convexity of (y) and the fact that h1 (() is convex as hy)) is concave. and the second inequality follows the chain rule to formulate (D(2)). Therefore,
-(( is (A) (
-strongly convex, and D(u) is (A) 1nga1
-strongly concave.
lig.1
an (m) -smooth as
Moreover (()) is (lh,,1 +
D(m)) (22)
= -11(-A41) - (-h2(-App)) - h2(11) + h11(22))
0 hī (A11)) App1 (-Ap22)
1/00/2
<
Ap1
(-A
Ap2
2
+ (- Vh1 M (2)) + U || 1221
2
Ouma (A)
(VD()) 11 m2)) + 1, +
|
142 11
2
The first inequality holds as h (3) and h (u) are smooth. The second follows the chain rule,
Note mmax (A) Tmin (A) > 0 as A is full column rank. This completes the proof.
Knowing (1)) has such favorable properties, we next analyze the update of || - Mell 1l, where {Fe} is the sequence generated in the non-accelerated version of Algorithm 2 with T2 1.
Lemma 14 (Update of || - pell). Consider the problem in (47 where hy is concave and 1 -soothh, h2 is Ohgy -strongly convex and 1 -smoohh,, and A is full column rank. Running the non-accelerated version of Algorithm with Ty and 1/2
S 1 -1 gives
1
||pt11 pell S
[hh,11 +
T mmax(A)
|ike 1=ll
122
ah2
+ Tmux (A)(1 1/1 1/ 12)l Vhū(-Api) + ||||I
(53)
where the constant X satisfies
and n² arg max HER D(1) with D(1)) defined in (52)
Proof. According to Lemma
13
((n) min hy (n) + n4 Ap h2(2) (n) - h5(-Ap)
yERdy
(4) smooth with respect μ. Moreover, the problem max pERde D(1) is equivalent to the unconstrained problem with the Lagrange multiplier
2 (A)
Ing.1 strongly concave and (I,,,
(A)
max D(H) D(v) + AT ME
where unique À is non-negative and finite in all dimension, i.c. 0 X oc, as D(u) is strongly
convex and μ E Rdc is equivalent to p. 2 0 satisfying the LICQ condition (Lemma 5). In this way,
VD(μ) VD(M) + + (n) A Vha(-Ap) A.
(54)
26

---

We can see that ((p) is smooth and strongly concave with the same modulus as D(H). The first-order stationary condition requires
= - (/a
(55)
In this way,
1
11411 14ell _
Projac (1t + -gt + (=(++t)) 14r
122
1/2
(a)
Vhi(pe) + A + yo-
ATM44 -| (p)) + A (pp) + + AT 3/+11 A Vh5(-Apu) All
V VD())l + Tmax (A)ll30+1 V - + MAII
WV W
VD(m) VD(n")I (mun)) (A)(1 - 1/20,, /2) Vhī(-Ap4) + ||l||
Imax
1n..1 +
|| - 1a 111l max (A) - -1000, -
+ ||l||
Othy
Inequality (a) comes from the non-expansiveness (1- (1-hiicchsssssss) of the projection operation, (b)
follows triangle inequality and uses (54), (c) uses (53) and (49) in Lemma 12, and (d) comes from the smoothness of D(u), which is of the same modulus D(u). This completes the proof.
In 53 the update behavior of || - [4+ depends on l|ke -='l|- We further look into the update of 144 14"1l and summarize in the following lemma the bound of the update of l| ||=-==
Lemma 15 (Update of I| - ===l)- Consider the problem in (47) where bn is concave and 1 smooth, l22 is Ofly -hhg-srooggl convex and 1 --mmooth and A is full column rank. Conduct the
non-accelerated version of Algorithm with Ty = 1, gives
l| 101l 5 (1 1/2
mmin (A)
l1fe 1²11 1l + 772 max (A)(1 710h//2)llu. Vhū(-Apr))
21
(56)
when 1)1 S 1h2.1 and 172
(In
Here μ² arg min min MERO. D(u) where D(u) is
defined in (52)
Proof. Define an auxiliary update as
Pt+1 == Projee: 4 (M + 72 D(H))) = Projed (pt +
+ A AA (57)
This projected gradient descent on strongly convex -D(N). As RR is closed and convex,
following Lemma 7 for 7/2
[1h,,1
o ) where (Ir.,1 + omac(A) is the modulus for
and
anz
smoothness of (²r) by Lemma 13
we have
li 10ll (1 (1-12 172
o minn
11ke p=ll.
21
As the real update M4+1 = Projge ((44 + (1 (I-Lipschitzness) of projection operation, have
+ A y.)). by the non-expansiveness
|| M4+1l| < 7// (y++1 - Vhb(-Apt))ll S 122 720 (A)] Vhū(-Au)]
By triangle inequality and (49). we have
(A) ||Me - 111l + m2x (A)[ Vhī(-Au:)l
|| -1'11 < 1
21,z,11
< <1-112 22 (A)
11:e 11 + max (p)mma((A)()) m1ang/2)llvx -aa-----,j--- Vh3(-Apt))
(58)
21
This completes the proof.
27

---

We are ready to proceed with the convergence analysis for the single-loop algorithm (Algorithm 2 without acceleration and T. = L. on the problems (47), which general form to and (5). In this way, Theorem 4 follows directly from the following theorem.
Theorem 6. Suppose L(,,)) u) is in the form of (47 where A is full column rank, hn is concave and -smoohh, b2 O1h2 -strongly convex and 1 smoothh satisfying 1,,11 (1), 1b2.1 Laz Z
(1). and 1ngl1 (1). Conduct the non non-accelerated version of Algorithm with T3 1. For
Oh2
Alng.15maxl (A) 0t2g o (A) (lbg,1 Emaa when 1n1 = O Tng.1 S Ing.1
arbitrary small positive S
and
72 = O(e) S 1 max (A)/asz"
the algorithm yields outpur (Pr, yr) such that
1IT 1011² < and -21 <e
with complexity T = ((I((c)))) Here, (1",, ,,","""== = arg maX min [(, v)-
Proof. For some positive constant p > 0, denote
Pr = plpt 14 +1 1lue -
(59)
Plugging 50 in Lemma 12 (53) in Lemma 14 and (56) in Lemma 15 to 59 we know
P2+ = p Mc+1 - 1 || + |1 (+
Juin (A)
Sp
1/2 21
11:ee 1 1²1 1 + 1/20 (A) - (10001/0////ll3 - AWill
+ (1 mmng/ /2)l)/ - Thg(-Auih + mma
1/2
Oh2
(Ih1.1 + o (max((A)
)lke 10"1 + mmax (A)(x(AA11)))) - 117 m/a /2)n///22lllv/// vhg(-Aua)) + [|AI]
ohin (A)
1
(max(())
112 1 + Tmax (A)) pll
+
pllite -=1ll
p
Orioz
Orhiz
oa (A)
maxx (A)
(1 7/1 a//2/22 /2)
1 prp20 max (A)
7/2
11ye Vhū(-
7a
|A|]-
To construct Pa-1 < (1 - Pt +
mmax A) nellA|| for some constant 0 < c < 1, it sufficient to find
an
< Tng,1' (1n
and 0 such that
0 (1 1/2 2(wrg1
oiin(() (m)r(l)
1/2 (1
)mme(d))
5 1/2 alsg,1
oiin (A) <
anz
ab2
0 <<<11= 7/1 0/// /2)
1 0 (A) +
(d)
1/2
<1 0/2 / (1 + 71 71hh, (2) <
ab2
This can be obtained when
p 1max (A)
(Uhg,1
ong tiin (1)
Ohh2
(60)
1/2
po (A)+ ((1ax(1))
Conditions in 160 41 (0mme(()
can be satisfied when > O is sufficiently small such that p = c- >
(I
22 (A) ), 11 = ( and 7/2 = O
Gig
r-1) = O(c1). In this way.
an2
[ny.1
PAI S(1 - c)P, 0(a
where e 0 is of the order O(e). Iteration gives
P; <1 - c + 0(a
(61)
28

---

Notice O(a),)) < ((1) and Po O as p e-l In this way, there exist Ti = (In(c=))) such that for all f Ti, (1 c) Po 0(1) and ( - 0(1) Accordingly, we can achieve
Pt = 01), A > Ti.
Moreover, as Pt P:=6-1114 1= I| + 11ye v 1 A44))l
l 1= & cP. = =O(c, Vt > Ti.
(62)
- for t Ti, Furthermore, choose 7)1 0 Tng. satisfying 711 S Ing.1
1l V )I <(1 1/1 //2/ /2) |1r: Vhū(-App:) O(c).
This is an iteration outcome using 55 in Lemma 12 62 and the fact that 0/2/2 /2 = 0 (11g.11
(1). In this way, for another T22 (In((=))) steps, we have
11 7 I1 (c),
and
|/y.
|1 S 11 - Vhi(-Ape)) + (1V Ay4)) Vii(-AAn)]
o max (A)
- 11y Vhi(-Apu +
1114 -1011
4hh
o (e + / - 0
Vt > 7i +22
We can see that the algorithm converges linearly with complexity O + T2) = (lI((c))). In this way, obtaining
l1or - y ==12==00 and 1Pr - 1=1 ((e),
requires complexity T ((I(((V)) ((n((=)))) This completes the proof.
Remark 4. Under the same assumptions as in Theorem 4
we cam choose 77.22 S (lg,1)=1, ng
(28 = -1n--nn ng as stepsizes for running the single-loop version of Algorithm 10
solve (7), and nF S (1g,)) 1J.22 S (709-
(700g11111 m
as the ones for
E Applications to Hyperparameter Optimization for SVM
In this section, we provide additional details about the SVM model training experiment for the linear SVM model, including the problem formulation and analysis of the results.
E.1 Problem formulation
SVMs train a machine learning model by finding the optimal hyperplane that separates data points of different classes with the maximum margin Misclassification is not tolerated for hard-margin SVMs. contrast, some samples are allowed be misclassified for soft-margin SVMs. Specifically, one first introduces variables ξ. which measure the violation associated with the classification of sample i, and then augments the original SVM objective with the norm of &, which vector collecting the values of & for all the samples in the training set.
BLO can be applied to the hyperparameter selection task of SVM. For example, it can be used to choose the value of the regularization parameters during soft-margin linear SVM training. Let us consider a classification problem and define Der (Zarr iy (u.:) lur.))) 1D ! as the training set, with and being the input feature vector for sample and lrr.i being its associated binary label. Similarly, let
D = ((vval, (vat,) [Pual be the validation set. For linear classification problem, the parameters of the SVM are w, a vector of coefficients with the same size as Str,i and the intercept b.
In short, we interested in the following constrained BLO problem
min (w. =
Σ
exp (1 Lval (zTal"* w" + +b*) + ||c|
(63a)
29

---

1
with w²,b"," ,,--" = arg min |
(63b)
wibe 2
s.1. ler (r,,(Znr,))) T
(1r, + 6) 2 1 1 Vie (1,..., Durll (63c)
& SsG
ViE f1....., (Der }. (63d)
The first term of upper-level objective (63a is a validation loss, evaluated on the validation set Daal+ the second is a regulation term on the upper-level variable c; the lower-level problem is designed to train the parameters of the hyperplane on the training set Dtr. with the soft margin violation & being upper bounded by the hyperparameter o; [ef. 63d). which are the CCs]. The lower-level objective (63b) focuses on maximizing the margin by minimizing ||l|| while allowing for violations to the separating hyperplane, which are regulated by the hyperparameter The BLO formulation aims to adjust the hyperparameter e using the validation loss (upper-level objective), ensuring that the model parameters (lower-level variables) are optimal for the training dataset.
E.2 Experiment details
In this section, we present detailed experimental results for the training of the SVM model in (63) using our BLOCC algorithm. We compare our algorithm to two baselines, LV-HBA [75] and GAM 72.. both designed for BLO problems with inequality CCs. We use 7 12 and rj = 0.01 for running our BLOCC and we apply a 0.01, "1 0.1, 72 = 0.1, = 0.001 for LV-HBA and a 0.05, c = 0.005 for GAM respectively. The hyper-parameters for LV-HBA and GAM are the ones used in the SVM experiments in their paper. The GAM algorithm often encounters issues with matrix inversion, as the matrix required for the GAM algorithm to solve the problem in (63 can be singular, preventing it from finding solution. Consequently, only BLOCC and LV-HBA solve the problem in 63). while GAM uses a different formulation as introduced in [72]. We compare the algorithms based validation loss (the first term in the upper-level objective) and classification accuracy. These metrics are computed for both validation and test datasets. We evaluate the algorithms on two datasets: diabetes 20 and fourclass 22
The detailed results are shown in Figure and Ten realizations of the problem are run, each with different training and validation sets. The plots show both the mean (line) and standard deviation (shaded region). The results reveal that our algorithm converges faster in terms of accuracy and loss, achieving lower loss value than the alternatives for both datasets, in validation and test sets. Furthermore, in terms of accuracy, we observe that for the diabetes dataset, our algorithm reaches a higher accuracy value faster than the alternatives in both validation and test sets.
Test Accuracy
4.0 Validation Loss (Upper Level)
Lower Level Objective
0.80
BLOCC
BLOCC
3.5
LV-HBA
LV-HBA
Acguracy
0.75
GAM
GAM
Loss
0.70
BLOCC LV-HBA
2.5-
2
0.65
GAM
2.0-
0.60
1.5
50
75
Running time
Running Time is
Iteration Count
(a) Test accuracy.
(b) Validation loss.
(c) Lower-level objective.
Figure 4: Plots for diabetes dataset
Test Accuracy
Validation Loss (Upper Level)
Lower Loss
4.0
0.8
BLOCC
BLOCC
3.5
LV-HBB
LV-HBA
Acouacy
0.7
g
3.0
GAM
Loxs
0.2
GAM
BLOCC
2.5
0.1
0.6
LV-HBA
2.0
GAM
1.50
25
50
75
Running time
Running Time is
Iteration Count
(a) Test accuracy.
(b) Validation loss.
(c) Lower-level objective.
Figure 5: Plots for fourclass dataset
30

---

F
Applications to Transportation Network Planning
This section applies the proposed BLOCC algorithm to a transportation network design problem, comparing it with the baselines from 72 and [75]
E.1 Problem formulation
In transportation network planning, the planning operator is to construct a new transportation network (upper level) connecting a collection of stations S, and the passengers will decide whether use the new network (lower-level), considering the options given by the new network and existing alternatives (constraints). The goal to design network that maximizes the operator's benefit, knowing the passengers will make rational choices depending on the given design.
The operator considers building a new network on collection of links A S S x S. For any link (i.j) E A connecting stations E S to j E S, the operator needs to design the capacity 24j. the link's capacity is set zero, then the link is not constructed. The larger the capacity, the larger the number of travelers, thus the larger the revenue while the higher the construction cost Cij-
The passengers are in demand to travel in the network K. S S x S. For every origin-destinatie pair (o, d). the traffic demand wdd and the traffic time on the existing route tod known. We assume that there is only one existing network. The proportion of passengers choosing the new network ynd and vl the fraction of passengers using link ji to travel from o to d will be determined following passengers' rational choice, modeled by logit choice model [T that will be explained later. If vi = 0, then the link does not belong to the route that passengers follow to travel from to d.
In a) of Figure 6, the station S are presented the dots, and links A
the dashed lines (input
topology): b) of Figure shows heatmap for the demand matrix for each (o, d) E KC (input to the design); c) is example of the constructed network connecting some of the stations (output of the design); and d) illustrates the number of passengers using the constructed network (design output).
To summarize, in this experiment, we assume that
A1) Only one existing alternative network exists;
A2) Passengers decide whether to use our network rationally, which is modeled by the logit choice model considering utility depends on travel attributes (travel time) [1; and,
A3) The demand per market, the trip prices, and the travel times per link are known to operators.
We summarize the optimization variables as follows
Fig E IR,,, the capacity constructed for the link (i,j) E A. The number of capacity variables is [A]. The link is not constructed if mig O.
yd E [0, 1]. the fraction (proportion) of passengers from market (o. d) E K. choosing the new network for their travel. The number of flow variables is [AC]. Since we consider only competitor, 1 ydd represents the fraction of incumbent network passengers. If y then the passengers of market (o, d) do not use the new network.
v [0, 1]. the fraction of passengers from market (o, d) E K. that, when choosing the new network to travel from to d, use the link (i,j) E A in their route to the destination. The number of link flow variables is [||][C]. With this definition, it holds that v(j < ydd, and, if
Vij 0, then link (i, does not belong to the route when traveling from o to d.
To simplify the notation and to make it consistent with the notations used in the paper, we denote
x = (I, (vis)ii represents the upper-level variables to be optimized.
X = | represents the domain of x.
y = ay (ug
represents the lower-level variables to be optimized
y = [e, 1 UK x [e, |||I where e is small positive number set by the designer of the network, represents the domain of y.
31

---

a) Input topology
Demand matrix
28000
21000
Es
lax
- 0.8
14000
7
1.5
fap=
destination
to 1.
2.2 F 1.6
6
7000
1.4
1og
6
1.8
top
1-
5
6
7
8
9
origin
Constructed network
d)
Passengers served by the constructed network
28000
21000
2
77.6
14000
1112
23.4
destination
Z45
28.0
6
7000
54.5
951.6
6
133.7
7
27.
7 8
9
origin
Figure 6: Example of transportation network design. (a) represents the of stations S 21, 4,5, 5, 6, 7, 8, 9} and the set of links A C S x S, where the number of links is [AA] = 30 (15 segments with two orientations each) and the value represents the demand [w"d} between all (o, d) pairs, where |C| = 9 x 8 72 values are provided (those in the off-diagonal elements of the heatmap). (c) represents the constructed network, where, to facilitate visualization, we have assumed that the capacity is symmetric and used the width of
each edge represents the travel time. (b)
the edge represent the capacity of every link. (d) represents (udd (V(.d)EE)) the number of passengers served by the constructed network.
Besides the optimization variables, our objective and constraints include the following parameters
po the total estimated demand (number of passengers) for the market (o, d) K.
m od the revenue obtained by the operator from a passenger in the market (o, d) E K.
Cij. the construction cost per passenger associated with link (i.j) E A.
tij. the travel time for link (i,j) E A.
pod travel time on the alternative network for passengers in the market (o, d) K.
We the coefficient associated with the travel time in passengers' utility function.
In transportation network design, a bilevel formulation is essential due to the interaction between two players with different levels of influence: the operator, who constructs the network, and the passengers, who choose their routes based on the network's 's characteristics. The operator's goal is maximize their benefit by minimizing construction costs and maximizing attracted demand, with link capacity as the optimization variable. Conversely, passengers aim to maximize their trip utility, determining the proportion demand using each link. This dual optimization requires that passenger choices comply with link capacity constraints set by the operator, coupling the variables at both levels. This necessitates a bilevel optimization approach, specifically using BLOCC, to address the interdependent decisions and constraints effectively.
32

---

Now, we are ready to introduce the objective formulations of our BLO problem. For the upper level, the network operator aims to maximize profits and minimize costs, and therefore, its interest is
min flz Vz =cX
Σ
m m (x))
CijTaj
(64)
H(o,d)ER
((JJ)EA
profit
cost
where yod (a) are optimal lower-level passenger flows associated with the network design æ.
For the lower-level, we model the passenger's behavior by finding the flow variables that maximize utility and minimizes flow entropy cost.
min ( 3) = yE)
Σ
w
wrtjihy +
Σ
wol yod)
(65)
(o,d)EK (A.))EA
(o.d)E
passengers ucility
yodyod (on(1 1)
Σ (o,
wod 31-100) - 1)
(o.d)EK
flow entropy cost
The passengers' utility considers the time cost of choosing the network and the existing network. The users set the flow variables y so that the transportation network yielding a higher utility is preferred. Here, the probability of chosing new network modeled by logistic (softmax) model. However, setting the objective as simple linear utility maximization would lead to an all-or-nothing policy, which is not the behavior observed in practice. Hence, the second term is brought consider. The approach considered here is to formulate an objective given by the Legendre transform (cf. Definition of the softmax sharing (see [54] 57 (56) for additional details). Intuitively, this means that rather than imposing the softmax sharing a fortiori, we formulate convex problem whose KKT conditions lead to the softmax sharing. Using the fact that the Legendre transform of an exponential e" is the negative entropy function y(in(y) 1), the lower-level objective is traditional as in [57]
Having introduced the optimization variables, parameters, and objective functions, next formulate our BLO problem, where we also incorporate the network constraints:
min xEX
Σ
m od odv
+
Σ
CijTij
(66a)
V(o,)EK
V(i.)EA
s.t. odx
Vij
oddt)
arg min yEy
Σ
Σ
w wellsyee
Σ od (1 -1- yydd) (a,))EK
(66b)
(o.d)EX
(i.))A
Σ
wod, ad
- 1)+
Σ
w wo (1 1od yod)
(o.))E
(o.d)EK
yod
if
Σ
Vij
Σ
yod
if d otherwise
(i, S
(66c)
Vi(()))
0
Σ
od Usj
S Saj
v(ij) EA
(66d)
Y(o,d)EK:
where (66c) are the flow-conservation constraints, (66d) are the capacity constraints that involve both upper and lower-level variables, coupling the optimization and motivating the of BLOCC. Note that, for networks with n stations (nodes): i) the number of CCs is approximately n and ii) each of the constraints involves approximately ni variables. Hence, even for a moderate-size network (say 30-50 nodes), we may have thousands of CCs involving millions of variables.
Experiment roadmap. To provide numerical results illustrating the behavior of our algorithm, solve the optimization in (66 for three scenarios:
S1) The design of a a 3-node simple synthetic network;
33

---

2.5
LV-HBA
gix, (x..)).)))))))) (x)))
-0.5
BLOCC y y=
2.0
BLOCC y
Rx: vd
-1.0
BLOCC y=
1.5
LV-HBA
BLOCC
-1.5
BLOCC y= BLOCC
1.0
-2.0
0.5
2.5
-3.0
0.0
20
60
Running time /s
Running
Figure 7: Upper-level objective (ax, v/) for 3-node network design problem; Solid lines show mean value of f( My. with the shaded region standard deviation; Dashed lines show the mean value of fIt, TF) with the shaded region as a standard deviation; Three different y values (red, purple, blue) and fixed stepsize 7) 1.6 x 10:4; The orange color represents the result of the LV-HBA algorithm.
Figure 8: Optimality gap of the lower-level problem for a 3-node network design prob- lem g ye) g(2,, Vī Solid lines repre- sent the mean value of the 10 realizations of the upper-level variables, dashed lines represent g(vr, ( g(r,, vi). and the shaded region is the standard deviation. Three different y values are represented in our algorithm, and fixed step- size 1 1.6 10-4
S2) The design of a 9-node synthetic network from the prior transportation literature; and S3) The design of a (real-world) subway network for the city of Seville, Spain, with 24 nodes
where S1) involves CCs and variables, and S3) involves around 100 CCs and 50,000 variables.
For the 3-node network scenario, we will conduct a comparative analysis against other algorithms to evaluate the efficacy of our approach. the other two scenarios, the baseline algorithms cannot find a solution; hence, for the 9-node and Seville networks, we will focus on providing insights into the performance and behavior and adaptability of our approach to real-world transportation networks.
algorithm under varying parameters, shedding light on the versatility
Before delving into the presentation and analysis of the results, two additional remarks are in order:
R1) While one of the goals of these experiments was to compare our BLOCC algorithm against LV-HBA 75 and GAM [72), for the scenario at hand, the GAM algorithm cannot be implemented, since the inverse of matrix at each iteration for the problem in (66) is not tractable. In this way, we only conducted the experiments using our BLOCC and LV-HBA
Tr and y While Theorems R2) The BLOCC algorithm produces two lower-level variables: VFT Tr and guarantee that (æT, FF) is an e-approximate solution, the pair (IT, yg.r) is strictly feasible, meaning that v (IT). Since strict fcasibility is important in the Tp transportation network, this section will report the results using both VFT and To
To
F22 Numerical results for the 3-node network
In this section, address the problem defined in equation (66) for network comprising 3 nodes (stations). We assume that the graph of potential links A is complete, leading to a total of link capacities that need be determined in the upper level. Additionally, as in the rest of the manuscript. we assume that there is demand for all markets, meaning that the set I. is complete and includes all possible origin-destination pairs. The specific values the key parameters can found in Table with further details about the simulated scenario available in the online code repository.
For BLOCC, set the stepsize three different values of 7: = 2, 3, and 4. The upper-level objective values are computed using fær () and fa TT). The results are presented in Figures and
1.6 10-4 and analyze the algorithm's convergence for
34

---

Market (o. d)
Parameter
(1,2)
(1.3) (2,1) (2.3) (3.1)
(3,2)
demand w
revenue mh"
6
2
6
1
travel time incumbent Fext
3
3
3
3
Link (t.j)
Parameter
(1,2)
(1.3)
(2,1) (2,2) (3.1)
(3.2)
link construction cost Ct link travel time Tup
10
3
10
2
Table Value of the parameters for scenario 1 (3-node network). The value wt is set to 0.1.
BLOCC y= 2
BLOCC y=
-5
BLOCC 4
-100
-150
fixe
-10
-20
BLOCC Y-2
-15
BLOCC
300
BLOCC y
-20
200
400
600
800
350
10000
20000
30000
Running time is
Running time
Figure 10: Upper-level objective f yo for metro network design problem in Seville, Spain, for 2 random initializations of the upper-level variables. Solid lines represent the mean of [æx, wgg3). and the shaded region is the standard deviation. The dashed lines represent the mean f(z,, vF,)). and the shaded region is the standard deviation. Three different y values are tested with fixed stepsize 17 1.6 x 10
Figure 9: Upper-level objective (ax, y1) for a 9-node network design problem; Solid lines show the mean value of f(,,, vg.3) with the shaded re- gion standard deviation; Dashed lines show
the mean value of f(2,, v with the shaded region standard deviation; Three different
y values (red, purple, blue) and fixed stepsize 7 1.6 10-4, The orange color represents the
result of the LV-HBA algorithm.
Figure illustratee Zillustaals the evolution of (xr ye) for the LV-HBA algorithm, while the other six lines represent different implementations of BLOCC. Each color represents a different value of 7; solid lines represent f( v1.t) and dashed lines represent f(r, v) Each simulation is conducted 10 times with 10 different random initializations of the upper-level variables, and both the mean and the standard deviation values are displayed. The results indicate that all versions our BLOCC algorithm converge in less than 10 seconds, while LV-HBA shows slight fluctuations even after running for more than 50 seconds. This may attributed to the fact that the LV-HBA algorithm requires a joint projection into (X ) of u) < 0) at each iteration, involving 42 variables and 6 CCs.
the performance our BLOCC algorithm over time. The orange line represents
Figure shows the lower-level optimality gap, namely g y1) g( y (z:)) for LV-HBA and BLOCC with different values of 7. In the case of LV-HBA, lower-level optimality is not attained within 70 seconds of running time. Conversely, for BLOCC, lower-level optimality is achieved by construction when is set to Mg resulting in a zero gap. Additionally, when y: is set to Ur.r. we observe that: i) lower-level optimality is accomplished for r 2 3, and ii) when y = 2, an optimality gap exists, but is one order of magnitude smaller than that for LV-HBA.
F33 Numerical results for the 9-node network
In this case, we consider the network in 21 see also Figure 6 which has |S] nodes and
[A| = 30 potential links. As before, we consider that all markets exist, so that [K] 9 72. The remaining parameters are described in Figure o and the code repository.
35

---

Figure is the counterpart of Figure 7 for the 9-node scenario, showing the behavior of our BLOCC algorithm for 7) = 1.6 10-4 and 7 E {2, 3. 4... Each simulation repeated 10 times (using 10 different random initializations of the upper-level variables), and both the mean and the standard deviation values are shown. Since the number of variables and constraints almost one order of magnitude larger, the algorithm requires more time to converge. However, convergence takes place in a reasonable amount of time (20-40 times longer than in the previous 3-node test case), Regarding the optimal value, we observe that: i) the sensitivity of the (steady-state) optimal value with respect to is not too large; ii) the solutions based on VFT yield better upper-level values than those based on yg,ri and iij) the gap between fIT 1 r) and f V decreases as y increases. Observation ii) is due to the fact that VF Ty is not feasible (meaning that it violates the optimality of the lower-level); hence, it is able to achieve better upper-level objective. In addition, the behavior observed in iii) is consistent with the discussion in Section 2.2 with the value of y having an impact on the suboptimality of VF., Tr at the lewer-level. Specifically, higher values of push V closer to Vg. lower-level optimality gap. Finally, we must note that for the solid lines (associated with vg it holds that (xre / = fx,, y" (az)). This implies that if we need a solution that feasible at the lower-level, then better objective values are associated with higher values of 7.
To and, hence, decrease the
F.4 Numerical results for the Seville network
In this section, we demonstrate the practical use of BLOCC in a real transportation network de- sign problem. Specifically, we address the de- sign of metro network in the city Seville, which has approximately one million inhabitants and is located in the south of Spain, with the bus system as its competitor. The data for the de- mand, number of stations, and locations have been taken from (11]. Information about the construction costs, capacity, and travel time has been obtained from Spanish rapid transit oper- ators (see references in [10 (56) for full details). The city authorities considered [S] 24 poten-
17
18
16
23
Figure 11: Topology of the Seville network. tial station locations. Regarding the links, the following assumptions are made:
A1) The link between nodes (i,j) E S S only exists if node j is one of the three closest neighbors to i, vice versa. The distance here is measured in terms of travel time. This assumption limits the number of lines in any station to be at most 3, which a very mild assumption for network with 24 stations.
A2) The link between nodes (i.j) E S x S only exists if the travel time tay is less than 7 minutes. This is also a mild assumption, since it enables all the stations to be connected, including those that are further away from the city center (the airport and the university campus (21)
Under these two conditions, the set A of potential links contains |A| stations, along with their actual locations, are shown in Figure markets between nodes, so AC| = 24 23- = 552.
88 links. The sets of links and Finally, we consider all possible
Following the narrative in Sections F.2 and F3
Figure 10 presents the evolution of the upper-level objective function with time for three values of the parameter 7 E (2,3, 4). The stepsize value has been set to 1 1.6 10-1, and two different initializations have been considered. Regarding the behavior of the algorithm with respect to the value of 7 and the particular output chosen (yT', vs. Te the findings are very similar to those in Figure Namely, smaller gaps are found for larger values of 7, and if feasibility at the lower-level (i.e., consistency with the user preference level) must be preserved, better objective values are achieved for larger values of 7-
The most important observation, however, is related to the running time. Specifically, Figure TO
reveals that, for this real-world scenario, our BLOCC algorithm converges in 10-20 hours. While this is more than 1,000 times larger than the convergence interval for the 3-node network, the number
36

---

of variables and constraints here is 100 times larger. More importantly, in the context of network transportation design. optimization times of 100 hours are widely accepted even for single-level formulations. Overall, we believe that the numerical results demonstrate that the BLOCC algorithm proposed in this work can solve problems with large number of variables and CCs, which can have practical value in real-world applications, such as the one studied in this section.
G Sensitivity Analysis
Regarding the selection and impact of seyper--arameteer on the performance of BLOCC, we conducted ablation study various values of the two critical parameters, y and 7), and measured their effects on the optimal value and computational time. this section, we present the sensitivity analysis on the toy example that we introduced in Section 4. and the 3-node network example for the transportation network planning problem in Section 4.3
G.1 Sensitivity analysis for the toy example
We present in Table the results of lower-level optimality (1v]
yerTl using different ? and r)
to conduct BLOCC (Al Algorithm d
on the toy example in Section 4.1
(xT)
VETT
7
Y == 0.001
Y= =
=.1
==
0.028 1 0.042
0.035 0.076
0.027 0.064
0.000 0.000
0.001
(3.314 0.042)
(3.186 0.076)
(2.049 + 0.064)
(1.967 + 0.000)
0.020 + 0.031
0.020
0.009 0.019
0.000 0.000
(2.242 + 0.031)
(1.66
(1.118 0.019)
(0.585 + 0.000)
0.006 # 0.010
0.027
0.011 = 0.033
0.000 0.000
(1.616 + 0.010)
(1.475 0.027)
(1.257 + 0.033)
(0.383 + 0.000)
0.023 + 0.019
0.030 + 0.022
0.020 + 0.045
0.000 - 0.000
1.0
(7.118 + 0.019)
(4.570 -: 0.022)
(3.477 + 0.045)
(1.645 + 0.000)
0.034 + 0.123
0.019 = 0.014
0.025 + 0.070
0.000 + 0.000
10.0
(6.565 1 0.123)
(4.365 1 0.014)
(2.808 1 0.070)
(1.750 1 0.000)
Table Sensitivity analysis for the hyperparameters in Section 4.1
Top line in each cell represents the optimality gap [luz yrrll. while the bottom line represents the time required for the algorithm to converge. Both the mean and the standard deviation are for 40 simulations.
From the table, draw the following empirical observations:
01) Larger values oy n bring the upper-level objectives closer to optimal. This is consistent with Theorem [ which illustrated that larger improves the accuracy of the lower-level optimality. Since the obtained solution yF is closer to y (IT) for larger y values, the distance between the f(rr, Vrr)) and f(r Vo (FT) will be closer well.
02) For a sufficiently small fixed 1, larger n lead to faster convergence. This implies that smaller ") S Te.1 choice due larger ? will not significantly dampen the convergence time. This is because a large value of increases [F sharpening the function (z(x)) and its gradient will be larger for most points. Thus, the gradient update V/V (x will not be very small and thus won't 't make the convergence slower.
G.2 Sensitivity Analysis for the 3-node network
We present Table 6 the results of the upper-level objective value fxx,, v) in network planning problem for 3-node network, whose detailed framework is introduced in Section
37

---

7-44
2
32x 30
$1.83
15.01 1:31
LIARE 202
6.79
36.43 8.36
1027
7.29
1
0.56
OLIAMT 1880
11.35 1360
3000
7.58 1964
K,08
30,30 19.85
3008 9.50
13.00
59.53
9.90
19,37
12.04
40.31 15.95
1.47
12:45 4.27
32x 10
855.36
218.22
923 00
1000 113
258.39
200.00
1 321122
48.41
200
173.00 28
1.9071 10.0
168.97
133000 3:34
85.02
16x
405.09 2006
109000-5:000
3TT.0 74.15
15.16
196.31
199
199.30 04.90
018 35 07.43
390.74 05
356.00 L 09.70
3
50A 36 61.03
1.415
305.70 57.00
Table 6: Sensitivity analysis for hyperparameters 7/ and in the experiment of the 3-node network transportation design described in Section 4.3 and Appendix [ Top line in each cell represents the upper-level objective value [(rr, ((FT.,,.T)) Hg,r), while the bottom line represents convergence time. Both the mean and the standard deviation of 10 simulations are provided.
We know from Section F that the goal is minimize fxx y), and negative values of f (az y) are expected. Therefore, (ær, yr) 0 in the table indicates a failure to converge and we can see that
03) Smaller y is needed to satisfy v S TF.1 to ensure convergence y is chosen larger. The algorithm tends to converge for different "y values when vj is small enough. This is because the Lipschitz smoothness constant 1F.1 of (x) increases with y (Lemma 3), and the BLOCC algorithm is guaranteed to converge y) S TF.1 (Theorem
In summary, the experimental results align with the theoretical analysis, demonstrating that the penalty constant 7 is robust and can be effectively set around 10, and the stepsize should be adjusted to ensure smooth and monotonic convergence.
H Analysis of the Computational Complexity
H.1
Complexity comparison
We present in the following the computational complexity of our BLOCC algorithm in comparison
with LV-HBA 75 and GAM 72 baselines.
Method
Iteration Cost
Computational Cost per Iteration
BLOCC
Outer loop (Algorithm
General setting (Theorem 3): O(ded + d (dedy dy));
T O(E
Special case (Theorem 4 (d + d' + delds + dy))
LV-HBA
O(d,d, de (d, dyy
GAM
Asymptotic
More than O(d) + di + de(dz + d,))
Table 7: Complexity comparison of our work with LV-HBA (75 and GAM [72].
The iteration costs are detailed in the earlier sections, and in 75 and 72 Hence, we discuss next the per-iteration cost. In the following discussion, we assume
A1) the complexity of calculating a function is proportional to the dimension of the inputs. For example, The complexity of finding Vag(z, v), V fx. y) are Odzz), and the one for (x, y)) is (dzddy)
A2) Projection cost on X' and y is respectively no more than O and Odd)). As discussed in the introduction, X' and Y are assumed to be easy-to-project domains, i.c. projection can be done using projection matrix or a simple formula, and the projection costs are no more than (dl) and (dl)). respectively.
In this way, BLOCC is first-order method with gradient calculation costs of Odg,de) and (dydd)) and the projection cost of ((d)) and (dl)). We present the detailed analysis of achieving the computational cost in the table in the following Section H.2
In the SVM model training and network planning experiments in Section 4 projections are simpler (c.g., truncation), resulting in ((dz) and (d,y) costs. In this scenario, for the general setting and it is -15-1.5..dd d,)). Therefore, our BLOCC is especially robust large-scale problems.
(Theorem 3). the complexity is 1 d.d, +e + -.ddddd and for the of affine in setting (Theorem
38

---

LV-HBA [75]. also a first-order method, has similar gradient calculation costs. However. its projection onto {X yy; (x, y) < 0} is expensive. with a complexity of ((d,, + d.,) of using interior point method to find the projected point [8)] and evaluating g(( y) I 0 adds (dgzddde) as it
requires de inequality judgment on functions taking input dimension dz, dy.
GAM (72) lacks an explicit algorithm for lower-level optimality and Lagrange multipliers. Even if we omit this, calculating Vig and its inverse incurs a cost of O(d)). Additionally, it has cost Odl + (d2 + dly)) for projection onto X' and calculating gradients.
We can see that BLOCC stands out with the lowest computational cost in both iterational and overall complexity thanks to its first-order and joint-projection-fre nature. Consequently, BLOCC is particularly well-suited for large-scale applications with high-dimensional parameters.
H.2
Complexity analysis of BLOCC
At each iteration our proposed BLOCC algorithm in Algorithm involves:
Step 1:
Solve two max-min problems.
Step 2: Calculate OF., in 13 and update F4+1 Projx (22 ngr;;))
where Step can be achieved by our proposed max-min solver in Algorithm 2 In the following, we provide analysis for using the accelerated version of the Algorithm 2 in the general case (Theorem 3) and the single-loop version for the special case (Theorem 4 as discussed in Section 3.3
In the general case, we use the accelerated version of Algorithm too solve both and For solving (7) at each inner loop iteration line of 3 the algorithm of computational cost (d The
update of y in line 4 involves calculating V Vv9(z, y) with a cost of (ddy)), finding (p, V Vyg (z, y)
for fixed a of Oddyde) following assumption A1), and the projection Projy of complexity (de)) according to A2). Moreover, converges linearly as L(p, y) is strongly convex in (Lemma 7. this
inner update for gives an iteration complexity of ((n((=1) ()). Therefore, the cost of the update for totals up to (In((=)) (dydo + d.)). The update of ş in line involves O(dydo) for calculating (pL((1+1/22 = (z, 10+1) with x being fixed, and (dd)) for projection onto P4:
Moreover, to achieve target accuracy e on the metric + Et= [[,,( S e, with y = O(e-1) using the BLOCC algorithm, the iteration complexity of the max-min solver is O according to Theorem 3 Therefore, the complexity for solving 7 in the general case using the accelerated version of Algorithm 2 (Theorem 3
OlE (d + d + ded, + de)) =
+ d.)).
(67)
Solving 5 is of the same order as the only difference is in calculating V v v). In this way, we can conclude that Step is at a cost of (ddyde + df))) in the general case.
In the special case (Theorem 4) the non-accelcrated (single-lopp) version of Algorithm involves 1 = (rojy(yt 7/, V ( 3t)) with complexity Oddy + dyde), and Mc+1 = Projege (Me 7/2 p. L U++1)) with complexity Old. + dyde) following similar analysis as in the general case.
(ee,
Moreover, the algorithm converges linearly in the special case of gf being affine in y and the computational complexity is (ln((=))), accordinng to Theorem Therefore, the complexity for the max-min Step in the special case is
OdIn(E))((,dd + di + dd)) Old,d. + dy + d.)).
(68)
Step 2 involves calculating V.g((z)) y) with with fixed y. and (/.. V v)) for fixed μr, y. which are of complexity (dz) and (d,dd) respectively, and conducting a projection on X of cost O (df) according to assumption A2).
In this way, the computational complexity of Step 2 is
Od + de)).
(69)
We can therefore conclude the complexity of BLOCC in Table 7
39

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction clearly state the claims made, including the contributions made in the paper and important assumptions and limitations.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: yes, we have discussed in the conclusions.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: All have been clearly listed.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: The code used to run the experiments published on GitHub.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The data files considered and the code used to run the experiments are published on GitHub.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All the necessary details to run the experiments are available in the code.
7. 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: [Yes]
Justification: The figures in the paper include confidence intervals representing the standard deviation for different realizations of the experiments. Also, the tables include the standard deviation of the simulations.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time exccution) needed to reproduce the experiments?
Answer: [NA]
40

---

Justification: The paper is theory paper, and does not include much computation-heayy experiments.
9. Code Of Ethies
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
Justification: We follow the NeurIPS Code of Ethics.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The work is foundational research, and there is no societal impact of the work performed.
11. Safeguards
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 (c.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses no such risks.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: The paper does not use existing assets.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: the paper does not release new assets.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects. 15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects. and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]-
Justification: The paper does not involve crowdsourcing nor research with human subjects.
41


Reviews:
[{"review_id": "znbw2q1srC", "reviewer": "Reviewer_SFco", "summary": "This paper presents an algorithm to solve bilevel optimization problems with coupled constraints using a primal-dual-assisted penalty reformulation. The study establishes rigorous convergence theory and demonstrates the algorithm's effectiveness through real-world applications, including SVM model training and transportation network planning.", "questions": "How can the proposed algorithm be seamlessly integrated with other solvers?", "limitations": "NA", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Strength:\n- The paper is well-written with clear introduction of the problem, assumptions and technical presentation for their methodology.\n- The primal-dual-assisted penalty reformulation and the main optimization algorithm effectively addresses the challenges of this class of problems.\n- Comprehensive theoretical analysis are provided to show the theoretical convergency. \n- Application examples including SVM and transportation network design, are provided to demonstrate the effectiveness of the proposed algorithm.", "weaknesses": "Weakness:\n- It would be great if more discussion on the computational complexity and scalability of the proposed method can be provided.\n- In the experiments part, the compared baselines are very limited.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating our contribution and the constructive review. Our response to your comments follows.\n\n**Q1. Computational complexity and scalability of BLOCC**\n\nThank you for raising this question. We answer this in **General Response G1**.\n\n**Q2. Limited baselines in the experiments.**\n\nWe did not include many baselines because bilevel optimization with coupled constraints is a less explored area with limited existing comparative work. The considered domain constraints impede the use of implicit gradient methods, while coupled constraints complicate finding descent directions for $x.$ Table 1 summarizes existing works addressing coupled constraints. Among them, only LV-HBA [69], GAM [67], and BVFSM [44] handle inequality constraints. We select LV-HBA and GAM as baselines for the SVM training problem, as these algorithms have conducted experimental analysis on this issue in their papers. For the network planning problem, LV-HBA is chosen as the sole baseline, being the only algorithm that addresses both coupled inequality and domain constraints. GAM [67] and BVFSM [44] cannot handle lower-level domain constraints as they require LL stationarity in an unconstrained space.\n\n**Q3 How can BLOCC be integrated with other solvers?**\n\nWe thank the reviewer for posing this question. BLOCC (Algorithm 1) is a general framework with convergence established (Theorem 2) by calling **any inner max-min** oracle with a given accuracy. Therefore, a) Algorithm 1 can be seamlessly integrated with any max-min or min-max algorithm that converges to the optimal solutions of the max-min subproblems of BLOCC and b) the analysis in Theorem 2 is valid regardless of the max-min solver used. For instance, for the original simulations carried out in our SVM training example, we used CVXpy, an effective solver for convex programming, to obtain the optimal $y$ and $\\mu$. In constrast, Theorems 3 and 4 analyze the computational complexity of the proposed scheme when a particular type of  min-max solver (the one described in Algorithm 2, which is one of the most advanced in the current state of the art) is used.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "J5qfqxbzSO", "reviewer": "Reviewer_YCHS", "summary": "The paper proposed a fully first-order algorithm named BLOCC for the bilevel optimization problems with lower level coupled constraints. They provide convergence theory for the algorithm and demonstrate its effectiveness through numerical experiments in SVM-based model training and infrastructure planning in transportation networks.", "questions": "1. There are some typos in the paper. For example, in line 165, it should be \\mathbb{R}^{d_x}_{+}, in line 125, there is an extra “and”. \n2. How did the authors choose the iteration number T for the two inner minimax solvers in the experiments?", "limitations": "Relies on the strong convexity of the lower level problem.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper designs a first-order algorithm for bilevel programming problems with coupling constraints in the lower level, and it does not require joint projection operations. The effectiveness of the algorithm is validated through extensive experiments. The writing is clear and well-structured.", "weaknesses": "1. This paper relies on some restrictive assumptions compared to other works, such as the strong convexity of both lower level objective $g$ and constraint $g^c$ with respect to $y$ (Assumption 2) and the Linear Independence Constraint Qualification (LICQ) condition in lower level (Assumption 4). These assumptions reduce the challenge of the bilevel problem considered. It is unclear whether the application problems tested in the numerical experiments satisfy these assumptions. The authors should clearly state in the introduction that their work depends on the strong convexity of the lower level problem.\n\n2. The proposed algorithm requires solving two minimax problems at each iteration, which can be computationally expensive for practical applications.\n\n3. The theoretical results do not clearly state how to choose the step size \\eta.\n\n4. The experiments in the paper lack details on the chosen hyperparameters of the algorithm.\n\n5. I have checked the code provided in the Supplementary Material. For the SVM problem, the implementation does not match the proposed algorithm; the minimax problems at each iteration are solved by different solvers than those stated in the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating our work and the constructive review. Our response to your comments follows.\n\n**Q1.1 Stating strongly convexity in introduction**\n\nThank you for your suggestion. We have noted this in line 43 and will clarify it further in the introduction, such as around line 26.\n\n**Q1.2 Whether the problem is difficult assuming:**\n- **$g$: stongly convex in y**\n- **$g^c$: convex & LICQ in y**\n\nWe thank the reviewer for this intersting question. We answer this affirmatively to say even with these assumptions, the constrained BLO remains challenging. \n\nFirstly, the problem remains difficult under the assumption of $g(x,\\cdot)$ being strongly convex. Although $y^*_g(x)$ is a singleton, it remains non-differentiable [36], thus impeding using implicit gradient methods. \n\nSecondly, the challenge remains under the assumption on $g^c(x,\\cdot)$ because finding $\\mu^*_g(x)$ effectively is difficult while it is important for the descent direction for $x$ under $g^c(x,y)\\leq 0$. Most existing solvers do not provide convergence results on the distance to optimal points, thus impeding the convergence analysis. Moreover, linear convergence of Algorithm 2 on special case (Theorem 4) is an innovative, first-time established result, and its theoretical analysis is challenging.\n\n**Q1.3 Whether the experiments satisfy the assumptions.**\n\n**G3 Whether the experiments satisfies the assumptions**\n\nYes, the experiments satisfy the assumptions.\n\nIn the **SVM** problem, the LL objective (57b) is strongly convex in LL variable $w$. The coupled constraint $1 - l_{{\\rm tr}, i}(z_{{\\rm tr}, i}^\\top w + b) \\leq c_i$ is linear and satisfies LICQ, given full row rank training data (**Assumptions 2, 3, & 4** hold). The smooth and Lipschitz conditions for the LL hold, and for the UL,  $c$ is effectively in a compact set since the training data is finite. Thus, Assumption 1 holds. **Assumption 5** is satisfied as the constraint function is affine in $y$, a sufficient condition as shown in Section 3.2.1.\n\nIn the **network design** problem, the Lipschitz and smooth conditions hold for (61a), (61b), and (61c) (**Assumption 1**). The LL problem (61b) is strongly convex on $y \\in \\mathcal{Y}$. Constraint (61d) uncoupled while the coupled constraint (61c) is linear and satisfies LICQ since $w^{od}$ does not have repeating elements in real-world problems, including this experiment (**Assumptions 2, 3, & 4** hold).\n\n**Q2 Computational cost of inner minimax oracle**\n\nThank you for raising this question. We provide analysis on computational cost in **General Response G1**. \n\n**Q3. Theoretical results on stepsizes.**\n\nThank you for your question. We will explicitly present the results to avoid confusion. Below are the detailed theoretical results for stepsize choices.\n\nFor BLOCC (Algorithm 1), the stepsize $\\eta \\leq (l_{F,1})^{-1}$ is presented in Theorem 1 (line 202). \n\nFor Algorithm 2, we use $(\\eta_{g,1}, \\eta_{g,2})$ for solving (7), and $(\\eta_{F,1}, \\eta_{F,2})$ for (5). In the accelerated version (Theorem 3), we specify $\\eta_{g,1} \\leq (l_{g,1} + l_{g^c,1})^{-1}$ and $\\eta_{F,1} \\leq (l_{f,1} + \\gamma l_{g,1} + l_{g^c,1})^{-1}$ in line 692 and 693. The choices of $\\eta_{g,2}$ and $\\eta_{F,2}$ depend on the rule in Lemma 11 (line 667). This, combined with the smoothness parameter for the duality (line 690), leads to $\\eta_{g,2} \\leq \\frac{\\alpha_g}{l_{g^c,0}}$ and $\\eta_{F,2} \\leq \\frac{\\gamma \\alpha_g - l_{f,1}}{l_{g^c,0}}$. In the single-loop version (Theorem 4), the stepsize choices follow Theorem 7. Specifically, $\\eta_{g,2} \\leq (l_{g,1})^{-1}$, $\\eta_{g,1} \\leq \\alpha_g^2 \\epsilon (2 s_{\\max} \\alpha_g + s_{\\max}^2 \\epsilon^{-1})^{-1} \\eta_{g,1}$, and $\\eta_{F,2} = (l_{f,1} + \\gamma l_{g,1})^{-1}$, $\\eta_{F,2} \\leq (\\gamma \\alpha_g - l_{f,1})^2 \\epsilon (2 s_{\\max} (\\gamma \\alpha_g - l_{f,1}) + s_{\\max}^2 \\epsilon^{-1})^{-1} \\eta_{F,1}$, where $\\epsilon$ is the target accuracy.\n\n**Q4. Experimental choice of hyper-parameters.**\n\nThank you for your question. We are sorry for not including this explicitly in appendix. We will clarify it.\n\nFor the SVM experiment, we used:\n- **BLOCC:** $\\gamma = 12$, stepsize $\\eta = 0.01$\n- **LV-HBA [69]:** $\\alpha = 0.01$, $\\gamma_1 = 0.1$, $\\gamma_2 = 0.1$, $\\eta = 0.001$\n- **GAM [67]:** $\\alpha = 0.05$, $\\epsilon = 0.005$\n\nThe parameters for LV-HBA and GAM are consistent with those used in their respective papers for SVM experiments.\n\nFor the network planning, we experimented with BLOCC using $\\gamma = 2, 3, 4$ as  in Table 3, and stepsize $\\eta = 1.6 \\times 10^{-6}$ as in lines 910, 928, and 953 in the Appendix.\n\n\n**Q5. Another minimax solver in the SVM experiment.**\n\nThank you for this question. BLOCC is a general framework with convergence ensured (Theorem 2) by calling **any inner max-min oracle** with specific accuracy, such as proposed Algorithm 2 and CVXpy solver as in the SVM problem. \n\nFurthermore, we recognized potential confusion and conducted additional experiments on the SVM example using Algorithm 2 (details in **the Figure in the PDF file**.) The results is similar to the presented one in our paper using CVXpy.\n\n**Q6 Typos**\n\nThank you for pointing out, we will correct it. \n\n**Q7Iteration number for the two inner minimax solvers?**\n\nThank you for this question. If $g^c$ is non-affine in y, we use the accelerated version of Algorithm 2 (Theorem 3) and we choose:\n- $T_g, T_F \\approx \\gamma\\epsilon^{-0.5}$ as the number of iteration the two maxmin solvers, and\n- $T_y \\approx \\ln(\\epsilon^{-1})$ as the inner iterate number for $y$\n\nfor a target level of accuracy $\\epsilon$, e.g. $\\epsilon = 1e-2$.\n\nIf $g^c$ is affine in $y$, we use the non-accelerated version (Theorem 4) with:\n- $T_g , T_F \\approx \\ln(\\epsilon^{-1})$, and\n- $T_y = 1$, making the solver being single-loop.\n\nIn this way, we can achieve $\\frac{1}{T}\\sum_{t=0}^{T-1}\\| G_{\\eta}(x_t)\\|^2 = {\\cal O}(\\epsilon),\\quad \\text{with}~\\gamma \\approx \\epsilon^{-0.5}$.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Uyyid9m9D6", "reviewer": "Reviewer_BSgg", "summary": "This paper proposes an algorithm named BLOCC for solving a particular class of bilevel optimization problems where so-called coupled constraints are present. Coupled constraints are constraints w.r.t. the upper and lower level variables which are applied in addition to simple constraints like bound/ball constraints which are easy to solve for. Bilevel problems with coupled constraints have not been addressed in prior works in this literature, which rely on solving the inner problem and using the implicit function theorem (IFT) to pass gradients through the inner problem.\n\nBLOCC is particular suited for extremely large-scale optimization problems where other existing approaches which use IFT are intractable to solve due to the size of the KKT systems of the inner problems. Notably, BLOCC avoids solving the inner problem directly, instead replacing the inner problem with a simpler problem which only involves the simple constraints (and can be solved with a fast projection instead of requiring an iterative method). The coupled constraint is moved to the objective function of the outer problem, which is also modified as a result to be a saddle point problem. A custom first-order min-max solver is proposed to solve the outer problem and is an additional contribution of the paper. \n\nThe convergence properties of the algorithm is established, along with equivalence of the penalty reformulation to the original bilevel problem and BLOCC is evaluated on some interesting bilevel optimization problems against related methods. Particularly appealing are the extremely large-scale transportation network design experiments.", "questions": "How easy/robust is the tuning of the penalty parameters?\n\nTypically you just have a set of constraints for a problem. What is the significance of decoupling the inner problem into two sets of constraints, i.e., g^c and \\mathcal{Y}(x)? This appears to be a quirk of this literature and not say something like [1].\n\n[1] Sinha et al. Bilevel Optimization based on Kriging Approximations of Lower Level Optimal Value Function. (2018 IEEE Congress on Evolutionary Computation)", "limitations": "I believe there is a general lack of discussion of limitations in this paper. It would be nice to understand at the least which parameters are difficult to tune.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The proposed method handles a broader range of bilevel problems compared to prior works in this literature, i.e., large-scale problems with coupled constraints in the inner problem. The characteristics of the algorithm itself (projection free, first-order only) are appealing.\n\nFurthermore, a reasonably detailed convergence analysis for all components of the algorithm are provided. The technical writing overall is of a reasonably high quality.\n\nFinally, the large-scale experiments for the network planning setting help demonstrate that BLOCC is useful to solve large-scale problems in practice which are currently unattainable by methods which use implicit function theorem.", "weaknesses": "The main weakness in the paper in my mind is the statement of contribution, it was not immediately obvious why BLOCC is necessary on the first read-through. Only after reading the experimental section did I realize the importance of the algorithm design for large-scale problems. \n\nAs a more concrete suggestion, I think contextualising this work in the broader bi-level optimization literature would be beneficial. For instance, the value function reformulation (see [1] for example) is another way to solve an arguably broader class of optimization problems. However, I can see why directly solving the resultant single level problem may be computationally intractable for very large-scale problems using standard solvers, which motivates BLOCC. I would suggest referencing that literature better and emphasizing why the algorithmic properties of BLOCC enable application to large-scale bi-level problems earlier in the paper.\n\nI also think more analysis around the comparison methods LV-HBA and GAM for the large-scale problem in 4.3 is required. For instance, line 330-331 states that for GAM a Hessian must be inverted (presumably as part of using the IFT and differentiating through the inner problem). However, the actual linear system of equations may be solved using something more scalable such as conjugate gradient, which may permit the IFT methods to perform reasonably on this task. Some discussion or a direct comparison would greatly improve the experimental section in the paper.\n\nFinally, a more detailed description of GAM and LV-HBA would be helpful to better understand the differences to BLOCC. For instance, it is claimed that LV-HBA projection steps are expensive, however looking at their paper it seems that the projection is a fairly simple operation (Proj_Z defined in section 2.3 in their paper)? Adding this information will make the paper more self-contained.\n\n[1] Sinha et al. Bilevel Optimization based on Kriging Approximations of Lower Level Optimal Value Function. (2018 IEEE Congress on Evolutionary Computation)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating our contributions and recognizing our algorithm on large scale experiments \"particularly appealing\".\n\n**Q1 More clarification for large-scale problems and textualizaing in a broader BLO literature view. e.g. as in value function reformulation [1]**\n\nWe appreciate the reviewer's suggestions. We will highlight the significance of BLOCC for large-scale problems earlier in the paper, and demonstrate its robustness through computational cost analysis as in **General Response G1**.\n\nWe appreciate the reviewer for directing us to the relevant literature. We will incorporate the referenced literature to show how BLO problems relate to the value function, as outlined in eq. (5)-(8) in the literature. We will emphasize that decoupling $x$ and $y$ allows us to access the value function and its gradient information, which has been crucial for establishing smoothness properties and our max-min oracle.\n\n**Q2 More analysis of the comparison on LV-HBA and GAM for the large-scale problem in 4.3.**\n\nThank you for the suggestion. We will include the comparison of the computational complexity of the algorithms in Section 4.3. The detailed analysis is provided in **General Response G1**.\n\n**Q3 More detailed description of GAM and LV-HBA. e.g. why projection is hard for LV-HBA (Proj_Z in sec. 2.3 in their paper)**\n\nThank you for this suggestion. We will include more details of the baselines and include the comparison of computational cost as in **General Response G1**. Regarding the difficulty of projection in LV-HBA, the algorithm includes $\\text{Proj}_C$ where $C:=\\{\\mathcal{X}\\times \\mathcal{y}: g^c(x,y)\\leq 0\\}$ is defined in section 2.1 in their paper [69]. This projection is of high cost of ${\\cal O}((d_x+d_y)^{3.5})$ as is equivalent to a convex programming problem and it requires additional ${\\cal O}(d_x d_y d_c)$ cost on analyzing $g^c(x,y)\\leq 0$.\n\n**Q4 How easy/robust is the tuning of the penalty parameters?**\n\nWe can fairly conclude that the tunning of the penalty parameter is not very hard and it is very robust.\n\nIn practice, $\\gamma \\approx 10$ can achieves satisfactory lower-level optimality. For simpler cases, such as the toy example in our study, $\\gamma = 1$ suffices, as shown in **General Response G2**. For more complex applications, such as our SVM and network planning experiments, we use $\\gamma = 12$ and $\\gamma = 3, 4, 5$, respectively.\n\nThe practical steps for tuning $\\gamma$ and the step size $\\eta$ are as follows:\n\n1) If the upper-level objective is non-convex, choose $\\gamma \\geq l_{f,1}/\\alpha_g$ for strong convexity with respect to $y$. If the upper-level problem is convex or if $l_{f,1}$ and $\\alpha_g$ are unknown, start with $\\gamma \\approx 5$. Avoid excessively large values of $\\gamma$.\n2) Apply a small step size $\\eta \\leq l_{F,1}^{-1}$ and run the algorithm. Increase $\\gamma$ if the LL optimality gap is large or if the objective function’s evolution is not smooth and monotonic, which indicates that $\\gamma < l_{f,1}/\\alpha_g$. Stop when satisfactory lower-level optimality and smooth, decreasing iterates are achieved.\n3) Try larger $\\eta$ that can keep the evolution of the objective function across iterations being smooth and monotonic, therefore fasten the iteration.\n\nRegarding the robustness and stability, we provide sensitivity analysis in **General Response G2**. We can see that the result is robust. Larger $\\gamma$ ensures better lower level optimality while it is of no necessity to bring $\\gamma$ too large for desired accuracy. \n\n**Q5. What is the significance of decoupling the inner problem into two sets of constraints, i.e., $g^c$ and $\\mathcal{Y}$? This appears to be a quirk of this literature and not say something like [1].**\n\nWe thank the reviewer for this question. Indeed, the constraints of BLO problems can be expressed as in equations (7)-(8) of [1]. \nAccording to our notation, we use $g^c(x,y)\\leq 0$ to denote the LL coupled constraints, $y\\in \\mathcal{Y}$ to denote uncoupled domain constraints and $\\mathcal{Y}(x):=\\{y\\in \\mathcal{Y}:g^c(x,y)\\leq 0\\}$ to denote the combination of the previous two. \n\nDecoupling these constraints makes it easier to conduct the theoretical analysis for the primal-dual method. Specifically, the functional constraint $g^c(x,y)\\leq 0$ can be incorporated into the objective function through a Lagrangian term, while the domain constraint cannot be handled in this way.\n\n\n**Q6. Which parameters are difficult to tune.**\n\nWe would say the parameters are generally not difficult to tune. Comparably, the step size $\\eta$ and the penalty hyper-parameter $\\gamma$ are harder. We have provide practical method and analysis on tunning the $\\eta$ and $\\gamma$ as in the **answer to Q4**.\n\nHere, we additionally provide some practical tips for choosing the number of iterations required to run BLOCC (Algorithm 1) and the min-max problem (Algorithm 2) ($T$, $T_y$, $T_g$, $T_F$), and the step sizes of the min-max problems ($(\\eta_{g,1}, \\eta_{g,2})$ for solving $(\\mu_g^*(x), y_g^*(x))$ defined in (7), and $(\\eta_{F,1}, \\eta_{F,2})$ for solving $(\\mu^*_F(x), y^*_F(x))$ in (5)).\n\nIn practice, we can choose the number of iteration of BLOCC (Algorithm 1) as $T\\approx \\gamma \\epsilon^{-1}$ where $\\epsilon$ is the target accuracy. If the constraint $g^c$ is not affine in $y$, we can use the accelerated version of Algorithm 2 as the inner maxmin solver with $T_y \\approx \\ln(\\epsilon^{-1})$ and choose $T_g,T_F\\approx \\gamma \\epsilon^{-0.5}$. For the stepsizes, small enough step size is fine to use in the practice, as we only provides theoretical upper-bound requirement on them. If the constraint $g^c$ is affine in $y$, we can use the non-accelerated version of Algorithm 2 as the inner maxmin solver ($T_y = 1$) and choose $T_g,T_F\\approx \\ln( \\epsilon^{-1})$. Similarly, we can use small enough step sizes. \n\nWe hope this can help in addressing your problem.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lulqMYDXZi", "reviewer": "Reviewer_QYWG", "summary": "This works addresses a bilevel problem where the lower-level has constraints that couple both the upper- and lower-level variables. The proposed method leverages a penalty approach and Lagrangian duality theory that transform the original bilevel problem into a single-level minimization problem where the x and y variables are decoupled. A first-order algorithm is developed which relies on the solution of two min-max subproblems and a rigorous convergence theory is developed. The effectiveness of the proposed method is validated on two real-world applications: infrastructure planning in transportation networks and support vector machine model training.", "questions": "* Can you elaborate on Assumption 3? What do you mean by saying that $g^{c}$ satisfies the LICQ in $y$? Does this mean that for a given x the LICQ holds: a) for every $y \\in Y(x)$ or b) for every KKT point $y^{*}(x)$ of the lower-level problem?\n* Do the Assumptions 1-5 and the conditions within the theorems hold in the bilevel problems of the experiments (SVM model training and Transportation network design)?\n* Typos\n  * In equation (7): the + symbol is misplaced\n  * In the proof of Lemma 10 in page 17: in the formula of $\\nabla u_{h}(x)$ we should have $\\nabla h^{c}$ instead of $h^{c}$", "limitations": "No negative societal impact.", "rating": 5, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "* This works deals with a challenging bilevel problem, as not only the lower-level is constrained, but in fact these constraints couple both the upper- and lower-level variables. \n* The proposed algorithm is first-order, even though the problem is a challenging one due to the presence of coupled constraints in the lower-level. In contrast, methods relying on implicit gradients, regardless if the problem is constrained or not, usually require access to second-order information.\n* Evaluation of the BLOCC algorithm on both toy examples and real-world applications (SVM model training and Transportation network design).", "weaknesses": "* Lemmas 2 and 3, which derive some key properties of the inner min-max problems, assume that the Lagrangian multipliers are bounded for every $x \\in X$. This is not a common assumption in literature and thus I am not sure why it is seen as mild. The authors cite Theorem 1.6 from [30] as a justification, however it is not immediately clear how the bounded property follows from that theorem. In fact, Theorem 1.6 is written for a more general (abstract) setting. As this is an important point, I think the authors should explain this in more detail. They might even include the theorem (or perhaps a different version of the theorem, tailored to the specific setting of this problem) in the Appendix and provide a clearer explanation.\n* Requires the solution of two min-max subproblems where the solution of each subproblem requires a double-loop algorithm. As a result, the proposed BLOCC method is a triple loop algorithm. On the contrary, most hypergradient methods rely only on the solution of a minimization problem (the lower problem).\n* The convergence result requires additional assumptions (assumption 5) beyond the ones concerning the objectives ($f,g$) and the constraints ($g_{c}$) of the original bilevel problem. It is unclear how restrictive assumption 5 is, as it is not expressed directly over $f,g,g_c$, but rather over some dual function. On the other hand, at least from the perspective of the dual function, the assumption is mild as the dual function is already concave.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for appreciating our presentation and recognizing our work. Our response to your comments follows.\n\n**Q1. Bounded optimal Lagrange Multiplier assumption for all $x\\in\\mathcal{X}$.**\n\nWe appreciate the reviewer for highlighting this issue. We apologize for the incorrect reference to \"Theorem 1.6 in [30]\". The accurate reference is Theorem 1 in [65], and we will make this correction. We calrify in the following how how Theorem 1 in [65] supports the boundness of the Lagrange multiplier.\n\nWe restate in the following of Theorem 1 in [65]. \n\n1) The set $\\Lambda(f)$ is closed and convex for all $f\\in \\mathcal{F}$.\n2) The set $\\Lambda(f)$ is non-empty for all $f\\in \\mathcal{F}$ if and only if (GCQ) issatisfied.\n3) Let (MFCQ) be satisfied. Then, the set $\\Lambda(f)$ is compact for all $f\\in \\mathcal{F}$.\n4) If there exists $f\\in \\mathcal{F}$, such that $\\Lambda(f)$ is compact, then (MFCQ) is satisfied.\n5) Let (LICQ) be satisfied. Then, the set $\\Lambda(f)$ is a singleton for all $f\\in \\mathcal{F}$.\n\nwhere $\\mathcal{F}$ is the set of all differentiable function, $\\Lambda(f)$ is the set of Lagrange multipliers with respect to $f$.\n\nThis theorem states in 3 that MFCQ, a weaker condition than LICQ, is sufficient for the compactness of the Lagrange multiplier. Since we assume LICQ in $y$ holds for all $x \\in \\mathcal{X}$, we believe the introduction of upper bound parameters $B_g$ and $B_F$ is reasonable.\n\nThank you again for your careful review and valuable feedback.\n\n**Q2. Due to the max-min subproblems, BLOCC is a triple loop algorithm while most hypergradient methods rely only on the solution of a minimization problem (the lower problem).**\n\nWe would like to clarify this from the following two aspects. \n\n*A1. BLOCC is an efficient triple-loop algorithm.* Firstly, having triple loops is not uncommon in solving challenging BLO problems with coupled constraints (CCs) since two of the loops are associated with the fact of BLO dealing with two nested optimization problems and the third one with the the presence of the CCs, see e.g. BVFSM [44] and GAM [67]. GAM [67], a hypergradient method, effectively is of triple-loop as it requires solving for $\\mu_g^*(x),y_g^*(x)$ as well. Secondly, a triple-loop algorithm is not necessarily costly, with the key being how many iterations are required in each loop. The loop in $y$ (steps 5-7 in Algorithm 2) has negligible computational complexity of ${\\cal O}(\\ln(\\epsilon^{-1}))$. This is insignificant (approximately ${\\cal O}(1)$), as even for $\\epsilon=10^{-9}$, $\\ln(\\epsilon^{-1}) \\approx 20$. Thirdly, and more importantly, we provide a detailed computational cost per iteration (see **General Response G1**) that shows that the max-min solvers (Algorithm 2) are highly efficient compared to the baselines. \n\n*A2. Hypergradients-based methods can hardly be applied in this case.* Finding the closed-form hypergradients in the Hypergradients-based methods is always non-trivial. In the literature, they mainly rely on the implicit function theory to derive the closed-form. \nHowever, the LL domain constraint $y \\in \\mathcal{Y}$ impedes the use of the implicit function theory as it needs a set of differentiable lower-level optimality conditions (e.g., $\\nabla_y g(x,y) + \\langle \\mu^*_g(x),\\nabla_y g^c(x,y)\\rangle = 0$ in an unconstrained domain), such as GAM [67] and BVFSM [44]. A detailed analysis of the non-differentiability can be seen in section 2 in [36]. Therefore, we do not use hyper-gradient methods as we are considering $y\\in\\mathcal{Y}$ and want to construct a fully first-order method using the penalty-based reformulation.\n\n**Q3. It is unclear how restrictive Assumption 5 is.**\n\nAssumption 5 is mild and the mildness is justified for the following reasons:\n\n- The convexity nature of the dual function and the unique solution guarantee via LICQ (lemma 5) imply that the condition $\\langle -\\nabla D_g(\\mu) + \\nabla D_g(\\mu_g^*(x)), \\mu - \\mu_g^*(x) \\rangle > 0$ holds locally for $\\mu \\neq \\mu_g^*(x)$. This ensures that Assumption 5 is particularly mild.\n\n- We provide a detailed example in Section 3.2.1 that illustrates a specific case being sufficient to this assumption, which further clarifies its applicability.\n\n\n**Q4 The meaning of Assumption 3 of  $g^c$ satisfies the LICQ in $y$.**\n\nAsssumption 3 means that LICQ holds **(b) for every KKT point  $y^*_g(x)$** of the lower-level problem, $\\nabla_yg^c(x,y^*_g(x))$ of size $d_c\\times d_y$ are linear independent in rows. Thank you for asking this question. We are sorry for the confutsion and we will clarify it.\n\n\n**Q5 Whether the experiments satisfy the assumptions**\nYes, the experiments satisfy the assumptions.\n\nIn the **SVM** problem, the lowL objective (57b) is strongly convex in LL variable $w$. The coupled constraint $1 - l_{{\\rm tr}, i}(z_{{\\rm tr}, i}^\\top w + b) \\leq c_i$ is linear and satisfies LICQ, given full row rank training data (**Assumptions 2, 3, & 4** hold). The smooth and Lipschitz conditions for the LL hold, and for the UL,  $c$ is effectively in a compact set since the training data is finite. Thus, Assumption 1 holds. **Assumption 5** is satisfied as the constraint function is affine in $y$, a sufficient condition as shown in Section 3.2.1.\n\nIn the **network design** problem, the Lipschitz and smooth conditions hold for (61a), (61b), and (61c) (**Assumption 1**). The LL problem (61b) is strongly convex on $y \\in \\mathcal{Y}$. Constraint (61d) uncoupled while the coupled constraint (61c) is linear and satisfies LICQ since $w^{od}$ does not have repeating elements in real-world problems, including this experiment (**Assumptions 2, 3, & 4** hold).\n\n**Q6. Typos.** \n*In equation (7): the + symbol is misplaced*\n*In Lemma 10 in page 17: in the formula of  $\\nabla v_h(x)$, the second term should be $\\nabla_x h^c$  instead of $h^c$*\n\nThank you for pointing out. We will correct it.\n\nWith the above clarifications, we hope the reviewer can kindly reconsider the evaluation of our work!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
# General Response

We sincerely thank the reviewers for their constructive feedback. We appreciate their recognition of our work in tackling the challenging constrained bilevel problem with an efficient algorithm, rigorous convergence theory, and demonstrating effectiveness and robustness through real-world experiments.

The rebuttal is divided into two sections: a general response and four individual responses for each reviewer. In the current section, we focus on two issues that were raised by multiple reviewers and that we feel warrant a more detailed explanation.

**G1. Computational cost of BLOCC.**

We present in the following the computational complexity of our BLOCC algorithm in comparison with LV-HBA [69] and GAM [67] as baslines.

|Method |Iterational cost|Computational cost per iteration|
|-|-|-|
| BLOCC | Outer loop (Algorithm 1) $T = {\cal O}(\epsilon^{-1.5})$ | **General setting (Theorem 3)**:$\tilde{\cal O}(d_c d_x + d_x^2 + \epsilon^{-1}(d_c d_y +d_y^2))$; **Special case (Theorem 4)**: $\tilde{\cal O}(d_x^2+d_y^2+d_c (d_x+d_y))$|
| LV-HBA |${\cal O}(\epsilon^{-3})$|  ${\cal O}(d_xd_yd_c + (d_x+d_y)^{3.5} )$ |
| GAM | Asymptotic |More than ${\cal O}(d_y^3+d_x^2+d_c(d_x+d_y))$|

**BLOCC** is a first-order method with gradient calculation costs of $\mathcal{O}(d_x d_c)$ and $\mathcal{O}(d_y d_c)$. Given the assumption (in line 28) that $\mathcal{X}$ and $\mathcal{Y}$ are "easy to project, such as the Euclidean ball," we presume that projections can be performed using projection matrices or simple functions. Thus, the projection cost is generally $\mathcal{O}(d_x^2)$ and $\mathcal{O}(d_y^2)$. In our SVM and network planning experiments, projections are simpler (e.g., truncation), resulting in $\mathcal{O}(d_x)$ and $\mathcal{O}(d_y)$ costs. Thus, the cost per iteration of BLOCC (Algorithm 1) is $\tilde{\cal O}(d_c d_x + \epsilon^{-1}d_c d_y)$ for the general setting (Theorem 3) and $\tilde{\cal O}( d_c (d_x+d_y))$ for the special case of affine coupling constraints (Theorem 4).

**LV-HBA**, also a first-order method, has similar gradient calculation costs. However, its projection onto $\{\mathcal{X} \times \mathcal{Y} : g^c(x, y) \leq 0\}$ is expensive, with a complexity of $\mathcal{O}((d_x + d_y)^{3.5})$, and evaluating $g^c(x, y) \leq 0$ adds $\mathcal{O}(d_x d_y d_c)$ as it requires $d_c$ inequality judgement on functions taking input dimension $d_x,d_y$.

**GAM** lacks an explicit algorithm for lower-level optimality and Lagrange multipliers. Even if we omit this, calculating $\nabla_{yy}^2 g$ and its inverse incurs a cost of $\mathcal{O}(d_y^3)$. Additionally, it has cost ${\cal O}(d_x^2+d_c(d_x+d_y))$ for projection onto $\mathcal{X}$ and calculating gradients.

We can see that BLOCC stands out with the lowest computational cost in both iterational and overall complexity thanks to its first-order and joint-projection-free nature. Consequently, BLOCC is particularly well-suited for large-scale applications with high-dimensional parameters.

**G2. Sensitivity analysis of the hyperparameters.**

Regarding the selection and impact of hyper-parameters on the algorithm's performance, we conducted an ablation study on various values of the two critical parameters, $\gamma$ and $\eta$, and measured their effects on the optimal value and computational time.

We present in the following the results of the toy example in section 4.1 of our paper. We also present the results on network planning problem in the *attached PDF file*. The real-world result is consistent with the one on toy example. 

|$\eta$ \ $\gamma$ |0.001|0.01|0.1 | 1.0|
|-|-|-|-|-|
| 0.001 | 0.028 ± 0.042| 0.035 ± 0.076| 0.027 ± 0.064| 0.000 ± 0.000|
| | 3.314 ± 0.042| 3.186 ± 0.076| 2.049 ± 0.064    |1.967 ± 0.000|
| 0.01| 0.020 ± 0.031| 0.020 ± 0.041| 0.009 ± 0.019| 0.000 ± 0.000|
|| 2.242 ± 0.031| 1.575 ± 0.041| 1.118 ± 0.019|0.585 ± 0.000|
| 0.1   | 0.006 ± 0.010| 0.017 ± 0.027| 0.011 ± 0.033|0.000 ± 0.000|
|| 1.616 ± 0.010| 1.475 ± 0.027|1.257 ± 0.033|0.383 ± 0.000|
| 1.0| 0.023 ± 0.019  | 0.030 ± 0.022| 0.020 ± 0.045|0.000 ± 0.000|
|| 7.118 ± 0.019    | 4.570 ± 0.022| 3.477 ± 0.045| 1.645 ± 0.000|

*Table 1: Sensitivity analysis of the toy example in our paper. (Mean ± SD on 40 simulations) **Top line** in each cell represents the optimality gap* $|y_g^* (x_T) - y_{F,T}|$*, while **bottom line** represents the time required for the algorithm to converge.*

From Table 1, we can draw the following empirical observations:

1) *Smaller $\eta$ is needed to satisfy $\eta \leq \frac{1}{l_{F,1}}$ to ensure convergence if $\gamma$ is chosen larger*. The algorithm tends to converge for different $\gamma$ values when $\eta$ is small enough. This is because the smoothness constant $l_{F,1}$ of $F_{\gamma}(x)$ increases with $\gamma$ (Lemma 2), and the algorithm is guaranteed to converge if $\eta \leq \frac{1}{l_{F,1}}$ (Theorem 1).

2) *Larger values of $\gamma$ bring the upper-level objectives closer to optimal*. Theorem 1 illustrated that larger $\gamma$ improves the accuracy of the LL optimality. Since the obtained solution $y_{F,T}$ is closer to $y_g^*(x_T)$ for larger $\gamma$ values, the distance between the $f(x_T,y_{F,T})$ and $f(x_T,y^*_g(x_T))$ will be closer as well. 

3) *For a sufficiently small fixed $\eta$, larger $\gamma$ lead to faster convergence*. This implies that *smaller $\eta\leq \frac{1}{l_{F,1}}$ choice due to larger $\gamma$ will not significantly dampen the convergence time*. This is because a large value of $\gamma$ increases $l_{F,1}$, sharpening the function $F_{\gamma}(x)$ and its gradient will be larger for most points. Thus, the gradient update $\eta \nabla F_{\gamma}(x_t)$ will not be very small and thus won't make the convergence slower.

In summary, while our theoretical results highlight the impact of $\gamma$ and $\eta$ on the optimality gap and computational complexity, we believe the experimental results summarized in the tables provide additional evidence and valuable insights.

Meta-review: This paper addresses the complex issue of bilevel optimization problems with coupled constraints, a topic that has been relatively underexplored in existing literature. The authors propose a first-order algorithm, BLOCC, specifically designed to tackle these challenging problems. The paper includes a rigorous theoretical analysis supporting the convergence of the proposed algorithm and demonstrates its effectiveness through applications in support vector machine (SVM) model training and infrastructure planning within transportation networks.

The paper is well-structured and clearly written, with a comprehensive introduction to the problem, assumptions, and methodology. The proposed primal-dual-assisted penalty reformulation and optimization algorithm effectively address the difficulties associated with bilevel problems involving coupled constraints. The theoretical analysis provided is thorough and establishes the convergence of the algorithm. 
All reviewers agree that the paper is acceptable pending revisions. The authors are advised to revise the paper in accordance with the reviewers' feedback provided in the rebuttal.
Decision: Accept (poster)