We thank the reviewer for appreciating our presentation and recognizing our work. Our response to your comments follows.

Q1. Bounded optimal Lagrange Multiplier assumption for all $x\in\mathcal{X}$.

We 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.

We restate in the following of Theorem 1 in [65].

1. The set $\Lambda(f)$ is closed and convex for all $f\in \mathcal{F}$.
2. The set $\Lambda(f)$ is non-empty for all $f\in \mathcal{F}$ if and only if (GCQ) issatisfied.
3. Let (MFCQ) be satisfied. Then, the set $\Lambda(f)$ is compact for all $f\in \mathcal{F}$.
4. If there exists $f\in \mathcal{F}$, such that $\Lambda(f)$ is compact, then (MFCQ) is satisfied.
5. Let (LICQ) be satisfied. Then, the set $\Lambda(f)$ is a singleton for all $f\in \mathcal{F}$.

where $\mathcal{F}$ is the set of all differentiable function, $\Lambda(f)$ is the set of Lagrange multipliers with respect to $f$.

This 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.

Thank you again for your careful review and valuable feedback.

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).

We would like to clarify this from the following two aspects.

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.

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. However, 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.

Q3. It is unclear how restrictive Assumption 5 is.

Assumption 5 is mild and the mildness is justified for the following reasons:

• 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.

• 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.

Q4 The meaning of Assumption 3 of $g^c$ satisfies the LICQ in $y$.

Asssumption 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.

Q5 Whether the experiments satisfy the assumptions Yes, the experiments satisfy the assumptions.

In 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.

In 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).

Q6. Typos. In equation (7): the + symbol is misplaced 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$

Thank you for pointing out. We will correct it.

With the above clarifications, we hope the reviewer can kindly reconsider the evaluation of our work!

We thank the reviewer for appreciating our contributions and recognizing our algorithm on large scale experiments "particularly appealing".

Q1 More clarification for large-scale problems and textualizaing in a broader BLO literature view. e.g. as in value function reformulation [1]

We 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.

We 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.

Q2 More analysis of the comparison on LV-HBA and GAM for the large-scale problem in 4.3.

Thank 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.

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)

Thank 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$.

Q4 How easy/robust is the tuning of the penalty parameters?

We can fairly conclude that the tunning of the penalty parameter is not very hard and it is very robust.

In 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.

The practical steps for tuning $\gamma$ and the step size $\eta$ are as follows:

1. 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$.
2. 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.
3. Try larger $\eta$ that can keep the evolution of the objective function across iterations being smooth and monotonic, therefore fasten the iteration.

Regarding 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.

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].

We thank the reviewer for this question. Indeed, the constraints of BLO problems can be expressed as in equations (7)-(8) of [1]. According 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.

Decoupling 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.

Q6. Which parameters are difficult to tune.

We 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.

Here, 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)).

In 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.

We hope this can help in addressing your problem.

We thank the reviewer for appreciating our work and the constructive review. Our response to your comments follows.

Q1.1 Stating strongly convexity in introduction

Thank you for your suggestion. We have noted this in line 43 and will clarify it further in the introduction, such as around line 26.

Q1.2 Whether the problem is difficult assuming:

• $g$: stongly convex in y
• $g^c$: convex & LICQ in y

We thank the reviewer for this intersting question. We answer this affirmatively to say even with these assumptions, the constrained BLO remains challenging.

Firstly, 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.

Secondly, 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.

Q1.3 Whether the experiments satisfy the assumptions.

G3 Whether the experiments satisfies the assumptions

Yes, the experiments satisfy the assumptions.

In 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.

In 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).

Q2 Computational cost of inner minimax oracle

Thank you for raising this question. We provide analysis on computational cost in General Response G1.

Q3. Theoretical results on stepsizes.

Thank you for your question. We will explicitly present the results to avoid confusion. Below are the detailed theoretical results for stepsize choices.

For BLOCC (Algorithm 1), the stepsize $\eta \leq (l_{F,1})^{-1}$ is presented in Theorem 1 (line 202).

For 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.

Q4. Experimental choice of hyper-parameters.

Thank you for your question. We are sorry for not including this explicitly in appendix. We will clarify it.

For the SVM experiment, we used:

• BLOCC: $\gamma = 12$, stepsize $\eta = 0.01$
• LV-HBA [69]: $\alpha = 0.01$, $\gamma_1 = 0.1$, $\gamma_2 = 0.1$, $\eta = 0.001$
• GAM [67]: $\alpha = 0.05$, $\epsilon = 0.005$

The parameters for LV-HBA and GAM are consistent with those used in their respective papers for SVM experiments.

For 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.

Q5. Another minimax solver in the SVM experiment.

Thank 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.

Furthermore, 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.

Q6 Typos

Thank you for pointing out, we will correct it.

Q7Iteration number for the two inner minimax solvers?

Thank 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:

• $T_g, T_F \approx \gamma\epsilon^{-0.5}$ as the number of iteration the two maxmin solvers, and
• $T_y \approx \ln(\epsilon^{-1})$ as the inner iterate number for $y$

for a target level of accuracy $\epsilon$, e.g. $\epsilon = 1e-2$.

If $g^c$ is affine in $y$, we use the non-accelerated version (Theorem 4) with:

• $T_g , T_F \approx \ln(\epsilon^{-1})$, and
• $T_y = 1$, making the solver being single-loop.

In 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}$.

We thank the reviewer for appreciating our contribution and the constructive review. Our response to your comments follows.

Q1. Computational complexity and scalability of BLOCC

Thank you for raising this question. We answer this in General Response G1.

Q2. Limited baselines in the experiments.

We 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.

Q3 How can BLOCC be integrated with other solvers?

We 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.