Overcoming Lower-Level Constraints in Bilevel Optimization: A Novel Approach with Regularized Gap Functions

## Abstract

Abstract Constrained bilevel optimization tackles nested structures present in constrained learning tasks like constrained meta-learning, adversarial learning, and distributed bilevel optimization.
However, existing bilevel optimization methods mostly are typically restricted to specific constraint settings, such as linear lower-level constraints.
In this work, we overcome this limitation and develop a new single-loop, Hessian-free constrained bilevel algorithm capable of handling more general lower-level constraints.
We achieve this by employing a doubly regularized gap function tailored to the constrained lower-level problem, transforming constrained bilevel optimization into an equivalent single-level optimization problem with a single smooth constraint.
We rigorously establish the non-asymptotic convergence analysis of the proposed algorithm under the convexity of lower-level problem, avoiding the need for strong convexity assumptions on the lower-level objective or coupling convexity assumptions on lower-level constraints found in existing literature.
Additionally, the generality of our method allows for its extension to bilevel optimization with minimax lower-level problem.
We evaluate the effectiveness and efficiency of our algorithm on various synthetic problems, typical hyperparameter learning tasks, and generative adversarial network.

## 1 Introduction

Bilevel optimization (BiO), which subsumes minimax optimization as a special case, is a hierarchical optimization comprising two levels, with one problem nested within another.
There is growing interest in BiO that is driven by an abundance of applications. Examples in machine learning (ML) include hyperparameter optimization ; ; , meta-learning ; ; ; , and
reinforcement learning ; .
Typically, the lower-leve (LL) problems of BiO in ML literature handle learning tasks as unconstrained optimization problems.
However, constraints are crucial for learning tasks and are becoming increasingly important in designing robust, fair, and safe ML systems ; ; .
The resulting constrained bilevel problem applies to a wider range of applications than the unconstrained one, such as adversarial learning
; ; ,
federated learning ; ; , continual learning ,
and meta-learning for few-shot learning .

Existing BiO methods mainly focus on the LL unconstrained case ; ; ; ; ; ; ; ; ; ; . Recently, a few approaches have emerged to tackle BiO problems with constraints.
However, most of these approaches are limited to linear constraints or constraints that rely solely on the LL variable. In this study, we address more general constrained BiO problems, where the LL problems involve constraints coupling both upper-level (UL) and LL variables. Explicitly, we consider the constrained bilevel optimization problem of the following form:

$$ $\min_{x\in X,y\in Y}\,F(x,y)\quad\mathrm{s.t.}\quad y\in S(x),$ (1) $$

with $S(x)$ being the set of optimal solutions of the constrainted lower-level problem,

$$ $\min_{y\in Y}\,f(x,y)\quad\mathrm{s.t.}\quad g(x,y)\leq 0,$ (2) $$

where $x\in\mathbb{R}^{n}$, $y\in\mathbb{R}^{m}$, $X$ and $Y$ are closed convex sets in $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$, respectively. The UL objective $F:X\times Y\rightarrow\mathbb{R}$, the LL objective $f:X\times Y\rightarrow\mathbb{R}$, and the LL constraint mapping $g:X\times Y\rightarrow\mathbb{R}^{p}$
are continuously differentiable functions.

The constrained BiO problem (1) is challenging due to the implicit constraint $y\in S(x)$, which requires $y$ to be a solution of an optimization problem. The LL constraints introduce extra complexity, especially when dealing with coupled LL constraints. This complexity hinders the straightforward extension of existing methods to constrained BiO problems .

In the pursuit of solving constrained BiO problem, several algorithms have been developed to address specific cases of (1).
For instance, the works ; handle constrained BiO problems with $g(x,y)=Ay-b$ and $f(x,y)$ being strong convexity with respect to (w.r.t.) $y$.
As for coupling constrained LL problem, the paper introduces an alternating projected SGD approach for a family of BiO problems with coupling LL linear equality constraints $Ay+h(x)-c=0$, where $A$ represents a matrix, $c$ denotes a vector, and $h(x)$ is a vector-valued function. Beyond the confines of a specific setting, the authors of develop a gradient-based approach for general BiO problems, where the LL problem is convex with equality and inequality constraints. Since the methods presented in these works predominantly utilize implicit gradient-based techniques, it is unexpected that they all presuppose the strong convexity of the LL objective w.r.t. $y$, alongside other regularity conditions, to ensure the uniqueness and smoothness of the LL solution mapping. More importantly, the algorithms in these studies would tolerate the heavier memory and computational cost of using second-order calculations in large-scale applications.

Therefore, it is imperative to develop first-order methods that do not necessitate explicit estimation of the implicit gradients, thereby facilitating the resolution of a broader class of constrained BiO problems.
In pursuit of this goal, the value function approach has garnered considerable attention due to its efficacy in designing first-order, single-loop numerical algorithms ; ; . However, challenges arise when the LL problem involves constraints. Specifically, the value function, $v(x):=\min_{y\in Y}\{f(x,y)\,|\,g(x,y)\leq 0\}$, is usually an implicit nonsmooth function, even when the underlying problem functions possess favorable properties. Typically, it is difficult to compute or estimate its generalized gradient.

A notable distinctive approach for addressing general constrained BiO problems has recently been presented in the work , which introduces a novel proximal Lagrangian value function to tackle constrained LL problems. By utilizing this function, they convert the constrained BiO problem into an equivalent optimization problem with smooth constraints. Notably, their reformulation preserves the coupling LL constraints, akin to the value function approach. Consequently, the proposed algorithm (LV-HBA) in involves Euclidean projection onto the coupled LL constraint set. However, such an operation typically requires the assumption of coupling convexity and can be costly when the set is complex.
Therefore, a natural question arises:
Can we develop a first-order algorithm to overcome possibly coupled lower-level constraints in bilevel optimization?

Our response to this question is affirmative. Next, we highlight the main contributions of this study. Additional related work is provided in Appendix A.1.

- •
Reducing constrained bilevel optimization into an equivalent optimization problem with only a single smooth inequality constraint.
We propose a novel single-level smoothed reformulation for constrained BiO with possibly coupled LL constraints. A key technique is the doubly regularized gap function defined in (3), which serves as an optimality metric for the LL problems and allows for straightforward Hessian-free gradient evaluation, as shown in (4). Furthermore, this type of gap function can be readily extended to tackle more complex bilevel optimization scenarios involving minimax lower-level problems.
- •
Developing a single-loop first-order algorithm without projection onto the coupled lower-level constraint set.
Building upon the newly introduced single-leve reformulation, we propose Bilevel Constrained GAp Function-based First-order Algorithm (BiC-GAFFA), a first-order algorithm that can be implemented entirely within a single-loop framework.
Furthermore, we rigorously establish the non-asymptotic convergence analysis of BiC-GAFFA under the convexity of LL problem, avoiding the necessity for either the strong convexity assumption on the LL objective or the full convexity assumption on the LL constraints.
- •
We validate the effectiveness and efficiency of BiC-GAFFA on various synthetic problems, typicial hyper-parameter learning tasks, and generative adversarial network. These experiments collectively substantiate the superior performance of BiC-GAFFA.

## 2 Regularized gap function and equivalent reformulation

Transforming a BiO problem to a single-level optimization problem is a useful strategy from both theoretical and computational perspectives.
In this section, we introduce a novel smoothed reformulation tailored for constrained BiO problems with potentially coupled LL constraints. For this purpose, we define the following doubly regularized gap function for the LL problem (2):

$$ $\mathcal{G}_{\gamma}(x,y,z):=\max_{\theta\in Y,\lambda\in\mathbb{R}^{p}_{+}} \Big{\{}\mathcal{L}(x,y,\lambda)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}- \mathcal{L}(x,\theta,z)-\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\Big{\}},$ (3) $$

where $z\in\mathbb{R}^{p}$, the Lagrangian function
$\mathcal{L}(x,y,z):=f(x,y)+z^{\text{{T}}}g(x,y)$, and $\gamma:=(\gamma_{1},\gamma_{2})>0$ is the regularization (or proximal) parameters.
This regularized gap function concept has been previously applied in various contexts, such as variational inequalities , standard Nash games , and the generalized Nash equilibrium problem . However, its application in studying BiO problems with constrained lower-level problem remains unexplored.

When the LL problem is convex, and the associated multipliers exist for any $y\in S(x)$, that is,
$\mathcal{M}(x,y):=\big{\{}\lambda\in\mathbb{R}^{p}_{+}~{}|~{}0\in\nabla_{y}f(x
,y)+\lambda^{\text{{T}}}\nabla_{y}g(x,y)+\mathcal{N}_{Y}(y),~{}\lambda^{\text{
{T}}}g(x,y)=0\big{\}}\neq\varnothing$,
the single inequality constraint $\mathcal{G}_{\gamma}(x,y,z)\leq 0$ can equivalently characterizes the solution set of LL problem.

###### Lemma 2.1 .

Assume that both $f(x,\cdot)$ and $g(x,\cdot)$ are convex for any $x\in X$. Let $\gamma_{1},\gamma_{2}>0$, we have $\mathcal{G}_{\gamma}(x,y,z)\geq 0$ for any $(x,y,z)\in X\times Y\times\mathbb{R}_{+}^{p}$. Furthermore, $\mathcal{G}_{\gamma}(x,y,z)\leq 0$
if and only if $y\in S(x)$ and $z\in\mathcal{M}(x,y)$.

Another advantageous property of $\mathcal{G}_{\gamma}$ is its continuously differentiability when both functions $f$ and $g$ exhibit continuous differentiability.

###### Lemma 2.2 .

Assume that both $f(x,y)$ and $g(x,y)$ are convex in $y$ on $Y$ for any $x\in X$ and be continuously differentiable on an open set containing $X\times Y$. Then $\mathcal{G}_{\gamma}(x,y,z)$ is continuously differentiable on $X\times Y\times\mathbb{R}^{p}_{+}$, and for any $(x,y,z)\in X\times Y\times\mathbb{R}^{p}_{+}$,

$$ $\nabla\mathcal{G}_{\gamma}(x,y,z)=\begin{pmatrix}\nabla_{x}f(x,y)+(\lambda^{*} )^{\text{{T}}}\nabla_{x}g(x,y)\\ \nabla_{y}f(x,y)+(\lambda^{*})^{\text{{T}}}\nabla_{y}g(x,y)\\ -\left(z-\lambda^{*}\right)/\gamma_{2}\end{pmatrix}-\begin{pmatrix}\nabla_{x}f (x,\theta^{*})+z^{\text{{T}}}\nabla_{x}g(x,\theta^{*})\\ \left(y-\theta^{*}\right)/\gamma_{1}\\ g(x,\theta^{*})\end{pmatrix},$ (4) $$

where $\theta^{*}$ and $\lambda^{*}$ denote $\theta^{*}(x,y,z)$ and $\lambda^{*}(x,y,z)$, respectively, defined as

$$ $\displaystyle\theta^{*}(x,y,z)$ $\displaystyle:=\underset{\theta\in Y}{\mathrm{argmin}}\left\{f(x,\theta)+z^{ \text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\},$ (5) $\displaystyle\lambda^{*}(x,y,z)$ $\displaystyle:=\underset{\lambda\in\mathbb{R}^{p}_{+}}{\mathrm{argmax}}\left\{ f(x,y)+\lambda^{\text{{T}}}g(x,y)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right \}=\mathrm{Proj}_{\mathbb{R}^{p}_{+}}\left(z+\gamma_{2}g(x,y)\right).$ $$

Now we derive a smooth single-level reformulation for the constrained BiO problem (1):

$$ $\min_{(x,y,z)\in X\times Y\times\mathbb{R}_{+}^{p}}F(x,y)\quad\mathrm{s.t.} \quad\mathcal{G}_{\gamma}(x,y,z)\leq 0.$ (6) $$

This reformulation problem (6) is equivalent to the original BiO problem (1).

###### Theorem 2.3 .

Assume that both $f(x,\cdot)$ and $g(x,\cdot)$ are convex for any $x\in X$. Let $\gamma_{1},\gamma_{2}>0$, the reformulation (6) is equivalent to the bilevel optimization problem (1), provided that for any feasible point $(x,y)$ of (1), a corresponding multiplier of the lower-level problem (2) exists at $(x,y)$, i.e., $\mathcal{M}(x,y)\neq\varnothing$.

###### Remark 2.4 .

The equivalent reformulation (6) possesses two noteworthy characteristics. First, the formulation includes only one inequality constraint. Second, the LL constraints $g(x,y)\leq 0$ are not explicitly stated in the reformulation (6), in contrast to the value function-based reformulation as well as the formulation presented in .
Note also that the proofs of Lemmas and Theorem presented in this section are available in Appendix A.3.

###### Lemma 2.1 .

###### Lemma 2.2 .

###### Theorem 2.3 .

###### Remark 2.4 .

## 3 Gap function-based first-order algorithm

Building upon the newly introduced reformulation, we proceed to the algorithm design phase.

First, to enhance the stability of the proposed numerical algorithms, we introduce an upper bound constraint on the variable $z$. Consequently, we consider the following truncated variant of the reformulation (6):

$$ $\min_{(x,y,z)\in X\times Y\times Z}\ F(x,y)\quad\mathrm{s.t.}\quad\mathcal{G}_ {\gamma}(x,y,z)\leq 0,$ (7) $$

where $Z:=[0,r]^{p}$ with $r\geq 0$ is a compact subset of $\mathbb{R}^{p}_{+}$.
It is demonstrated that for sufficiently large $r$, any optimal solution to (7) is also optimal to the original reformulation (6).

###### Proposition 3.1 .

Suppose $\gamma_{1},\gamma_{2}>0$ and that an optimal solution $(x^{*},y^{*},z^{*})$ to (6) exists with $z^{*}\in Z$,
then any optimal solution of (7) is also optimal for the reformulation (6).

Second, to develop a gradient-based algorithm for solving (7), we explore its penalty formulation:

$$ $\min_{(x,y,z)\in X\times Y\times Z}F(x,y)+c\,\mathcal{G}_{\gamma}(x,y,z),$ (8) $$

where $c>0$ is a penalty parameter.
This work focuses on first-order penalty methods for the constrained BiO problem (1). To this end, we propose algorithms to find approximate stationary solutions of (8), in alignment with several previous works ; ; ; .

Third, we study the relationship between the penalty formulation (8) and a relaxed problem of (7):

$$ $\min_{(x,y,z)\in X\times Y\times Z}\ F(x,y)\quad\mathrm{s.t.}\quad\mathcal{G}_ {\gamma}(x,y,z)\leq\varepsilon,$ (9) $$

where $\varepsilon>0$ is the relaxation parameter. The relation between the penalized and relaxed problems of BiO problems has been previously studied in across various single-level reformulations. Herein, we discuss the relationship between (8) and (9).

###### Proposition 3.2 .

Assume that $F(x,y)$ is bounded below by $\underline{F}$ on $X\times Y$. For any $\varepsilon>0$, there exists $\bar{c}>0$ such that any global solution $(x_{c},y_{c},z_{c})$ to the penalty formulation (8) with penalty parameter $c\geq\bar{c}$ is also a global solution to the relaxed problem (9) with some relaxation parameter $\varepsilon_{c}$ satisfying $\varepsilon_{c}\leq\varepsilon$. Moreover, if $(x_{c},y_{c},z_{c})$ is a local solution to the penalty formulation (8), then it is also a local solution to the relaxed problem (9) with relaxation parameter $\varepsilon_{c}:=\mathcal{G}_{\gamma}(x_{c},y_{c},z_{c})$.

If we consider the penalty formulation (8) with an increasing sequence of penalty parameter $c_{k}$ such that $c_{k}\rightarrow\infty$, any accumulation point of the sequence of solutions associated penalty formulation problem (15) with varying values of $c_{k}$ is a solution to the problem (7).

###### Theorem 3.3 .

Assume that $X$ and $Y$ are closed, and functions $F$, $f$ and $g$ are continuous on $X\times Y$. Suppose that $c_{k}\rightarrow\infty$ and let

$$ $(x_{k},y_{k},z_{k})\in\underset{(x,y)\in X\times Y\times Z}{\mathrm{argmin}}F( x,y)+c_{k}\,\mathcal{G}_{\gamma}(x,y,z).$ $$

Then, any accumulation point $(\bar{x},\bar{y},\bar{z})$ of the sequence $\{(x_{k},y_{k},z_{k})\}$ is a solution to problem (7).

The proofs of Lemmas and Theorem presented in this section are available in Appendix A.4.

###### Proposition 3.1 .

###### Proposition 3.2 .

###### Theorem 3.3 .

### 3.1 The proposed algorithm

Now we introduce a first-order gradient-based, single-loop algorithm to solve the truncated and penalized approximation problem (8) with a potentially varying penalty parameter $c_{k}$. Note that solving (8) still requires care since it involves an optimal value function $\mathcal{G}_{\gamma}(x,y,z)$ for another optimization problem.

Lemma 2.2 establishes that $\mathcal{G}_{\gamma}$ is continuously differentiable, implying that the objective function of (8) is also continuously differentiable. Consequently, we can apply a gradient descent-type method to solve (8). However, as also noted in Lemma 2.2, computing the gradient $\nabla\mathcal{G}_{\gamma}(x^{k},y^{k},z^{k})$ at the current iterate $(x^{k},y^{k},z^{k})$ requires solving the minimization problem in (5) to obtain the exact solution $\theta^{*}(x^{k},y^{k},z^{k})$, a process which can be computationally intensive.

To mitigate this computational challenge, we introduce an auxiliary sequence $\{\theta^{k}\}$ as an approximation to $\theta^{*}$. At each iteration, we employ a single projected gradient descent step to update $\theta^{k+1}$ to approximate $\theta^{*}(x^{k},y^{k},z^{k})$, as follows:

$$ $\theta^{k+1}=\mathrm{Proj}_{Y}\left(\theta^{k}-\eta_{k}d_{\theta}^{k}\right),$ (10) $$

where $\eta_{k}>0$ is the step size, and

$$ $d_{\theta}^{k}:=\nabla_{y}f(x^{k},\theta^{k})+(z^{k})^{\text{{T}}}\nabla_{y}g( x^{k},\theta^{k})+\frac{1}{\gamma_{1}}(\theta^{k}-y^{k}).$ (11) $$

Furthermore, we introduce iterate $\lambda^{k+1}$ to represent $\lambda^{*}(x^{k},y^{k},z^{k})$ as follows:

$$ $\lambda^{k+1}=\lambda^{*}(x^{k},y^{k},z^{k})=\mathrm{Proj}_{\mathbb{R}^{p}_{+} }\left(z^{k}+\gamma_{2}g(x^{k},y^{k})\right).$ (12) $$

By substituting $(\theta^{k+1},\lambda^{k+1})$ for $(\theta^{*},\lambda^{*})$ in (4), we can approximate the gradients of the objective function in
(8) to define the update directions:

$$ $\displaystyle d_{x}^{k}:=$ $\displaystyle\frac{1}{c_{k}}\nabla_{x}F(x^{k},y^{k})$ (13) $\displaystyle+\nabla_{x}f(x^{k},y^{k})+(\lambda^{k+1})^{\text{{T}}}\nabla_{x}g (x^{k},y^{k})-\nabla_{x}f(x^{k},\theta^{k+1})-(z^{k})^{\text{{T}}}\nabla_{x}g( x^{k},\theta^{k+1}),$ $\displaystyle d_{y}^{k}:=$ $\displaystyle\frac{1}{c_{k}}\nabla_{y}F(x^{k},y^{k})+\nabla_{y}f(x^{k},y^{k})+ (\lambda^{k+1})^{\text{{T}}}\nabla_{y}g(x^{k},y^{k})-(y^{k}-\theta^{k+1})/ \gamma_{1},$ $\displaystyle d_{z}^{k}:=$ $\displaystyle-(z^{k}-\lambda^{k+1})/\gamma_{2}-g(x^{k},\theta^{k+1}).$ $$

Finally, we implement an update for the variables $(x,y,z)$ using a step size $\alpha_{k}>0$:

$$ $(x^{k+1},y^{k+1},z^{k+1})=\mathrm{Proj}_{X\times Y\times Z}\left((x^{k},y^{k}, z^{k})-\alpha_{k}(d_{x}^{k},d_{y}^{k},d_{z}^{k})\right).$ (14) $$

The complete algorithm is presented in Algorithm 1. Notably, our proposed algorithm uses only the gradient information of the problem’s functions and can be easily implemented when the projections onto the sets $X$ and $Y$ are computationally simple. Furthermore, the updates of $x^{k+1},y^{k+1}$ and $z^{k+1}$ can be executed in parallel, enhancing computational efficiency.

Figure: Algorithm 1 Bilevel Constrained GAp Function-based First-order Algorithm (BiC-GAFFA)

## 4 Non-asymptotic convergence analysis

In this section, we rigorously establish the non-asymptotic analysis for BiC-GAFFA towards the truncated and penalized approximation problem

$$ $\min_{(x,y,z)\in X\times Y\times Z}\psi_{c}(x,y,z):=F(x,y)+c\,\mathcal{G}_{ \gamma}(x,y,z).$ (15) $$

The proofs of Lemmas and Theorem presented in this section are available in Appendix A.5.

### 4.1 General assumptions

The following assumptions formalize the smoothness property of the UL objective $F$, and smoothness and convexity properties of the LL objective $f$ and the LL constraints $g$.

###### Assumption 4.1 ( UL objective).

The UL objective $F$ is $L_{F}$-smooth on $X\times Y$.
Additionally,
$F$ is bounded below on $X\times Y$, i.e., $\underline{F}:=\inf_{(x,y)\in X\times Y}F(x,y)>-\infty$.

###### Assumption 4.2 ( LL objective).

Assume that the following conditions hold:

- (i)
For each $x\in X$, $f$ is convex w.r.t. LL variable $y$ on $Y$.
- (ii)
$f$ is continuously differentiable on an open set containing $X\times Y$ and is $L_{f}$-smooth on $X\times Y$.

###### Assumption 4.3 ( LL constraints).

Assume that the following conditions hold:

- (i)
For each $x\in X$, $g$ is convex w.r.t. LL variable $y$ on $Y$.
- (ii)
$g(x,y)$ is $L_{g}$-Lipschitz continuous on $X\times Y$.
- (iii)
$g(x,y)$ is continuously differentiable on an open set containing $X\times Y$, $\nabla_{x}g(x,y)$ and $\nabla_{y}g(x,y)$ are $L_{g_{1}}$ and $L_{g_{2}}$-Lipschitz continuous on $X\times Y$, respectively.

Our assumptions solely require the first-order differentiability of the problem functions, avoiding (possibly) higher-order smoothness.
The setting of this study substantially relaxes the existing requirement for second-order differentiability in constrained BiO literature.
Moreover, we do not impose either the strong convexity on the LL objective $f(x,\cdot)$ or the fully convexity of the LL constraints $g(x,y)$.

###### Assumption 4.1 ( UL objective).

###### Assumption 4.2 ( LL objective).

###### Assumption 4.3 ( LL constraints).

### 4.2 Convergence results

The proof of non-asymptotic convergence for BiC-GAFFA primarily hinges on establishing the sufficient descent property of the merit function defined as follows:

$$ $V_{k}:=\phi_{c_{k}}(x^{k},y^{k},z^{k})+C_{\theta}\left\|\theta^{k}-\theta^{*}( x^{k},y^{k},z^{k})\right\|^{2},$ $$

where

$$ $\phi_{c_{k}}(x,y,z):=\frac{1}{c_{k}}\big{(}F(x,y)-\underline{F}\big{)}+ \mathcal{G}_{\gamma}(x,y,z),$ $$

and $C_{\theta}:=L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+L_{g}.$ With appropriately chosen step sizes $\alpha_{k}$ and $\eta_{k}$, the following inequality holds

$$ $V_{k+1}-V_{k}\leq\,-\frac{1}{4\alpha_{k}}\left\|w^{k+1}-w^{k}\right\|^{2}- \frac{\eta_{k}C_{\theta}}{2\gamma_{1}}\|\theta^{k}-\theta^{*}(w^{k})\|^{2},$ $$

where $w^{k}:=(x^{k},y^{k},z^{k})$. See Lemma A.5 for a detailed description.

We define the residual function for problem (15) as follows

$$ $R(x,y,z):=\mathrm{dist}\left(0,\nabla\psi_{c}(x,y,z)+\mathcal{N}_{X\times Y \times Z}(x,y,z)\right),$ $$

noting that $R(x,y,z)=0$ if and only if $(x,y,z)$ is a stationary point to (15), i.e.,

$$ $0\in\nabla\psi_{c}(x,y,z)+\mathcal{N}_{X\times Y\times Z}(x,y,z).$ $$

###### Theorem 4.4 .

Under Assumptions 4.1, 4.2 and 4.3, assume that $X$, $Y$ are compact sets, $\gamma_{1}>0$, $\gamma_{2}>0$, $c>0$ and $\eta_{k}\in(\underline{\eta},1/(L_{f}+rL_{g_{2}}+1/\gamma_{1}))$ with $\underline{\eta}>0$. Then there exists $c_{\alpha}>0$ such that when $\alpha_{k}\in(\underline{\alpha},{c_{\alpha}})$ with $\underline{\alpha}>0$, the sequence of $(x^{k},y^{k},z^{k},\theta^{k},\lambda^{k})$ generated by Algorithm 1 satisfies

$$ $\min_{0\leq k\leq K}R(x^{k+1},y^{k+1},z^{k+1})=O\left(\frac{1}{K^{1/2}}\right).$ $$

In subsequent analysis, we consider Algorithm 1 and the penalty approximation problem (15) with an iteratively increasing penalty parameter $c_{k}$, leading to the following convergence result based on the residual function $R_{k}(x,y,z):=\mathrm{dist}\left(0,\nabla\psi_{c_{k}}(x,y,z)+\mathcal{N}_{X
\times Y\times Z}(x,y,z)\right).$

###### Theorem 4.5 .

Under Assumptions 4.1, 4.2 and 4.3, assume that $X$, $Y$ are compact sets, $\gamma_{1}>0$, $\gamma_{2}>0$, $c_{k}=c(k+1)^{\rho}$ with $c>0$, $\rho\in[0,1/2)$ and $\eta_{k}\in(\underline{\eta},1/(L_{f}+rL_{g_{2}}+1/\gamma_{1}))$ with $\underline{\eta}>0$. Then there exists $c_{\alpha}>0$ such that when $\alpha_{k}\in(\underline{\alpha},{c_{\alpha}})$ with $\underline{\alpha}>0$, the sequence of $(x^{k},y^{k},z^{k},\theta^{k},\lambda^{k})$ generated by Algorithm 1 satisfies

$$ $\min_{0\leq k\leq K}R_{k}(x^{k+1},y^{k+1},z^{k+1})=O\left(\frac{1}{K^{(1-2\rho )/2}}\right).$ $$

Furthermore, if $\rho>0$ and $\psi_{c_{k}}(x^{k},y^{k},z^{k})$ is uniformly bounded above, then the sequence of $(x^{k},y^{k},z^{k})$ satisfies

$$ $0\leq\mathcal{G}_{\gamma}(x^{K},y^{K},z^{K})=O\left(\frac{1}{K^{\rho}}\right).$ $$

###### Remark 4.6 .

The hypergradient norm is commonly employed as a stationary measure for BiO problems, when the LL problem is unconstrained and strongly convex. However, for constrained BiO problems, even if the LL objective is (strongly) convex w.r.t. the LL variable, the differentiability of variants of hypergradient, including optimistic and pessimistic ones, is not fully understood. To our best knowledge, no universally recognized stationary measure is known in this scenario. Different methods use various stationary measures. For instance, the KKT residual function is utilized in ; . The residual functions in the above theorems originate from the corresponding penalized problems, aligning with several previous prior on first-order penalty methods ; ; ; .

###### Theorem 4.4 .

###### Theorem 4.5 .

###### Remark 4.6 .

## 5 Extension to bilevel optimization with minimax lower-level problem

In this section, we explore the extension of our proposed gradient-based, single-loop, Hessian-free algorithm, originally designed for bilevel optimization problems with constrained lower-level problems, to bilevel optimization problems with minimax lower-level problem ; ,

$$ $\displaystyle\min_{x\in X,y\in Y,z\in Z}F(x,y,z)\quad\mathrm{s.t.}\quad(y,z) \in\mathcal{SP}(x),$ (16) $$

where $\mathcal{SP}(x)$ denotes the set of saddle points for the convex-concave minimax problem,

$$ $\displaystyle\min_{y\in Y}\max_{z\in Z}f(x,y,z),$ (17) $$

where $x\in\mathbb{R}^{n}$, $y\in\mathbb{R}^{m}$ and $z\in\mathbb{R}^{p}$, the sets $X$, $Y$ and $Z$ are closed convex sets in $\mathbb{R}^{n}$, $\mathbb{R}^{m}$ and $\mathbb{R}^{p}$, respectively. The UL objective $F:X\times Y\times Z\rightarrow\mathbb{R}$, and the LL objective $f:X\times Y\times Z\rightarrow\mathbb{R}$ are continuously differentiable with $f$ being convex in $y$ and concave in $z$.
Building upon the idea applied in the development of the regularized gap function (3), we introduce the doubly regularized gap function for lower-level minimax problems, defined as:

$$ $\displaystyle\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z):=\max_{\theta\in Y, \lambda\in Z}\Big{\{}f(x,y,\lambda)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}-f(x ,\theta,z)-\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\Big{\}}.$ (18) $$

This newly introduced gap function enables the following equivalent single-level reformulation of the problem (16), detailed in Theorem A.8,

$$ $\min_{(x,y,z)\in X\times Y\times Z}\ F(x,y,z)\quad\mathrm{s.t.}\quad\mathcal{G }^{\mathrm{saddle}}_{\gamma}(x,y,z)\leq 0.$ (19) $$

Utilizing this reformulation, analogous to BiC-GAFFA, we can propose a gradient-based, single-loop, Hessian-free algorithm for problem (16), detailed in Algorithm 2 in the Appendix. While the primary focus of this paper is the bilevel optimization with constrained lower-level problems, we defer the convergence analysis of this algorithm to future work.

## 6 Numerical experiments

To validate both the theoretical and practical performance of our proposed algorithm (BiC-GAFFA), we conduct experiments on both synthetic tests and real-world applications. We compare BiC-GAFFA with various related algorithms on the synthetic tests and some real-world applications. Detailed information about the experiments can be found in Appendix A.2. Our code is available at [https://github.com/SUSTech-Optimization/BiC-GAFFA](https://github.com/SUSTech-Optimization/BiC-GAFFA).

### 6.1 Synthetic experiments

Here we consider the following synthetic bilevel optimization problem:

$$ $\begin{array}[]{>{\displaystyle}c @{\,\,} >{\displaystyle}l}\min_{\begin{ subarray}{c}\mathbf{x}\in\mathcal{X},\\ (\mathbf{y}_{1},\mathbf{y}_{2})\in\mathcal{Y}\end{subarray}}&(\mathbf{y}_{1}-2 \cdot\mathds{1}_{n})^{\text{{T}}}(\mathbf{x}-\mathds{1}_{n})+\|\mathbf{y}_{2}+ 3\cdot\mathds{1}_{n}\|^{2}\\ \text{s.t.}&(\mathbf{y}_{1},\mathbf{y}_{2})\in\mathop{\mathrm{argmin}}_{( \mathbf{y}_{1},\mathbf{y}_{2})\in\mathcal{Y}}\big{\{}\frac{1}{2}\|\mathbf{y}_{ 1}\|^{2}-\mathbf{x}^{\text{{T}}}\mathbf{y}_{1}+\mathds{1}_{n}^{\text{{T}}} \mathbf{y}_{2}\ \text{ s.t. }\ \sum_{i=1}^{n}h(\mathbf{x}_{i})+\mathds{1}_{n}^ {\text{{T}}}\mathbf{y}_{1}+\mathds{1}_{n}^{\text{{T}}}\mathbf{y}_{2}=0\big{\}} ,\end{array}$ (20) $$

where $\mathcal{X}=\mathbb{R}^{n}$, $\mathcal{Y}=\mathbb{R}^{n}\times\mathbb{R}^{n}$, and $h(t):\mathbb{R}\mapsto\mathbb{R}$ is defined as $h(t)=t^{q}$. We consider the cases $q=1$ and $q=3$. Note in both cases, the optimal solution is $\mathbf{x}^{*}=\mathds{1}_{n}$, $\mathbf{y}_{1}^{*}=2\cdot\mathds{1}_{n}$, $\mathbf{y}_{2}^{*}=-3\cdot\mathds{1}_{n}$, and $S(\mathbf{x})=\{(\mathbf{x}+\mathds{1}_{n},\mathbf{y}_{2})|\sum_{i=1}^{n}h(
\mathbf{x}_{i})+\mathds{1}_{n}^{\text{{T}}}(\mathbf{x}+\mathds{1}_{n})+\mathds
{1}_{n}^{\text{{T}}}\mathbf{y}_{2}=0\}$, where $\mathds{1}_{n}$ denotes the $n$-dimensional vector with all elements equal to $1$.

Figure: Figure 1: The first two pictures are results for the problem (20) with $q=1,n=1000$, and the third and fourth pictures are results for the problem (20) with $q=3,n=1000$.
Refer to caption: x1.png

We compare our algorithm (BiC-GAFFA) with AiPOD and LV-HBA on the synthetic tests. Hyperparameters are collected in Appendix A.2.1.
The comparison results are collected in Figure 1. From this, we can see that BiC-GAFFA converges fast and correctly in both cases, while AiPOD fails to optimize either of the problems, the algorithm LV-HBA converges well in the case of $q=1$ but runs slowly in the case of $q=3$ due to the high complexity of the projection step.

Figure: Table 1: Sensitivity analysis on the problem (20) with $q=3,n=100$.
Refer to caption: x3.png

Sensitivity analysis.
We use the time when $\|\mathbf{x}_{k}-\mathbf{x}^{*}\|/\|\mathbf{x}^{*}\|<0.01$ as a metric to measure the performance of BiC-GAFFA. The comparison results for different hyperparameters are presented in Table 1, which demonstrate the robustness of BiC-GAFFA. Figure 2 shows the results based on the problem’s dimension, which is plotted on logarithmic scales for both axes, indicating that BiC-GAFFA scales well for large-scale problems. Detailed reports and additional information are provided in Appendix A.2.1.

### 6.2 Hyperparameter optimization

**Table 2: Results on the sparse group Lasso hyperparameter selection problem with nTr = 100, nVal = 100, nTest = 300.**
| Method | Time(s) | Val Err | Test Err |
| --- | --- | --- | --- |
| Grid | 17.3 $\pm$ 0.9 | 35.9 $\pm$ 7.2 | 37.7 $\pm$ 6.7 |
| Random | 17.4 $\pm$ 0.7 | 33.6 $\pm$ 6.7 | 35.7 $\pm$ 6.2 |
| TPE | 16.9 $\pm$ 0.7 | 33.9 $\pm$ 7.0 | 36.0 $\pm$ 5.6 |
| IGJO | 21.2 $\pm$ 2.2 | 19.7 $\pm$ 2.8 | 25.6 $\pm$ 4.4 |
| VF-iDCA | 12.4 $\pm$ 0.5 | 14.6 $\pm$ 2.6 | 25.4 $\pm$ 3.9 |
| BiC-GAFFA | 21.4 $\pm$ 0.7 | 7.3 $\pm$ 1.3 | 22.3 $\pm$ 3.0 |

Hyperparameter optimization (HO) is an inherent bilevel optimization.
According to Theorem 3.1 of Gao’s work , we can turn the hyperparameter optimization problem of a statistical learning model into an equivalent constrained BiO problem.
Three HO problems are considered in this section.
Detailed results and information are provided in Appendix A.2.2.

Sparse group LASSO problem. For the HO problem for the sparse group LASSO model, we compare BiC-GAFFA with some widely-used algorithms, including Grid Search, Random Search, TPE , IGJO , and VF-iDCA . Detailed settings are provided in Appendix A.2.2. From Table2, we can see that BiC-GAFFA outperforms the other algorithms in terms of both validation and test errors.

Support vector machine (SVM).
In this part, we apply BiC-GAFFA to a HO problem where the lower model (base learner) is a SVM and the upper model is a signed distance based loss.
We compare BiC-GAFFA with GAM and LV-HBA and the results are shown in Figure 4. We can see that BiC-GAFFA converges faster than GAM and achieves a better performance.

Figure: Figure 3: Results of the SVM problem on the diabetes dataset.
Refer to caption: x4.png

Data hyper-cleaning.
Data hyper-cleaning is to train a classifier based on data where part of their labels are corrupted with a probability $p_{c}$, to achieve this goal, we need to train a classifier that can recognize the wrong labeled data. Such a process can be modeled as a HO problem. Results are shown in Figure 4, from which we can see that BiC-GAFFA outperforms the other algorithms in terms of both validation and test accuracies.

### 6.3 Generative adversarial networks

Generative adversarial networks (GAN) models are also investigated in our paper. It involves a hierarchical structure of optimization with two interacting components: the generator and the discriminator. We compared different training strategies for a distribution recovery problem. Detailed information is provided in Appendix A.2.3. The numerical results in Figure 5 and Table 3 show BiC-GAFFA can train such a generative model efficiently.

Figure: Figure 5: Earth mover’s distance v.s Time.
Refer to caption: x6.png

## 7 Discussion and conclusion

This work introduces BiC-GAFFA, a single-loop first-order algorithm tailored for a broader class of bilevel optimization problems, where the LL constraints may depend on the UL variables.
A key technique employed here to address constrained LL problem is the newly introduced doubly regularized gap function, which serves as an optimality metric for the LL problems and allows for straightforward gradient evaluation. We also study the non-asymptotic performance guarantees of BiC-GAFFA. The experimental results validate the effectiveness of BiC-GAFFA.

Limitations and future work.
We acknowledge that the proposed method is currently limited to the deterministic setting and we plan to
extend our method to incorporate uncertainty and distributed settings in future work. Note that BiC-GAFFA is projection-free onto the coupled LL constraint set but requires projections onto the sets $X$ and $Y$. Consequently, it is interesting to develop fully projection-free algorithms for solving general constrained BiO problems.

## Appendix A Appendix

The appendix is organized as follows:

- •
We give a brief review of additional related work in Section A.1.
- •
Experimental details are provided in Section A.2.
- •
Proofs in Section 2 can be found in Section A.3.
- •
Proofs in Section 3 are presented in Section A.4.
- •
Proofs in Section 4 are given in Section A.5.
- •
Details regarding the extension to bilevel optimization with minimax lower-level problem are discussed in Section A.6.

### A.1 Additional related work

In the section we give a brief review of additional related work that are directly related to ours.

Gap functions.
Transforming a BiO problem into a single-level optimization problem is useful both theoretically and computationally. Gap functions are crucial in this process.
Among them, the most commonly used one is the value function, originally proposed in [, ]. When the LL problem is unconstrained and its smooth objective is strongly convex w.r.t. the LL variable, the value function of the LL problem is smooth and has a closed-form gradient with the LL solution. Benefiting from this property, fully first-order gradient-based algorithms have been developed, see, e.g., [, , , , ].
However, the value function is often nonsmooth when the LL problem is constrained or its objective is merely convex.
As a result, the value function reformulation usually leads to a nonsmooth single-level optimization problem.
Recently, the Moreau envelope-based gap functions have emerged, see, e.g., [, , ]. Among these works, by utilizing the Moreau envelope-based reformulation for unconstrained LL problem, the very recent study [] proposes a novel single-loop gradient-based algorithm for general BiO problems with nonconvex and potentially nonsmooth LL objective functions. In contrast, both [] and [] address constrained BiO problems. The key difference is that the Moreau envelope-based reformulation in [] is nonsmooth, while the reformulation in [] is smooth. The smoothness is achieved by leveraging a novel proximal Lagrangian value function to handle constrained LL problem. When the LL problem has (coupling) inequality constraints, it is important to note that all the reformulations in these works explicitly preserve these (coupling) inequality constraints.
In contrast, this work introduces a novel reformulation that integrates the LL (coupling) inequality constraints into a doubly regularized gap function. Consequently, it achieves the minimal smooth constraint, with only one smooth constraint in the single-level reformulation. This provides a significant advantage in algorithm design because it removes the need for Euclidean projection onto the (coupled) LL constraint set. This expands the range of applications and significantly reduces computational costs.

First-order algorithms. Many applications involve optimization problems with thousands or millions of variables. First-order algorithms are popular because their storage and computational costs can be kept at a tolerable level.
In bilevel optimization, the value function approach avoids recurrent second-order calculations like those involving the Hessian matrix in implicit gradient-based methods. This makes it particularly useful for developing efficient first-order, single-loop numerical algorithms, see, e.g., [, ], in cases where the LL problem is unconstrained and its smooth objective is strongly convex w.r.t. the LL variable. When the LL problem involves constraints, challenges arise. For instance, the value function is often nonsmooth and lacks an exact straightforward gradient evaluation.
To address the nonsmooth issue, the work [] introduces a sequential minimization algorithm framework using penalty and barrier functions. The recent study [] provides a sufficient condition under which the value function is (Lipschitz) smooth and proposes first-order algorithms through the lens of the penalty method for bilevel problems with constraints that rely solely on the LL variable. Also, through the lens of the penalty method, a first-order method with complexity guarantees is developed in [] using a novel minimax approach, and the recent study [] proposes a first-order stochastic bilevel optimization algorithm when the LL constraints depend only on the LL variable. Other advances in first-order algorithms for bilevel optimization include: primal-dual algorithms [];
primal nonsmooth reformulation-based algorithm []; and low-rank implicit gradient-based methods [].

### A.2 Details of experiments

All the experiments in this paper are performed on a computer with Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz and 16.00 GB memory. Except for the manual experiments, the reported data for all the other experiments (including the sparse group LASSO problems, SVM, data hyperclean and GAN) are the statistical results after repeating each experiment 20 times.

#### A.2.1 Synthetic problems

For all the experiments on this problem, we use the initial point $(\mathbf{x}_{0},\mathbf{y}_{0},\mathbf{y}_{1})=(\mathbf{0}_{n},\mathds{1}_{n},
\mathds{1}_{n})$, where $\mathbf{0}_{n}$ represents the $n$-dimensional zero vector, and $\mathds{1}_{n}$ represents the $n$-dimensional vector of all ones.
For the synthetic experiments recorded in Figure 1, we use the following settings:

- •
For AIPOD, we set $S=5$, $T=2$, $\beta=0.0001$, $\alpha=0.0005$;
- •
For LV-HBA, we set $\gamma_{1}=10,\gamma_{2}=1,\alpha=0.005,\beta=0.002,\eta=0.03$;
- •
For BiC-GAFFA, we set $\gamma_{1}=1,\gamma_{2}=0.1,\alpha=0.001,\eta=0.01$, $\rho=0.3$.

They are applied for both problems, i.e., for both $n=1000$, $q=1$ problem and $n=1000$, $q=3$ problem. Note that AIPOD requires the calculation of matrix inversion, which is computationally expensive, however, due to the special structure of the problem, the matrix inversion in each iteration returns the same result, so we only calculate it once and use the result in all iterations. For LV-HBA, the projection step is the most time-consuming part generally, it gets benefits from the special structure of the problem where $q=1$ since the projection is efficient in that case, which is the reason why it performs well in the first case. However, for the problem where $q=3$, the projection step is much more complex, which requires solving a nonlinear optimization problem, therefore, it performs poorly in the second case.

**Table 4: Sensitivity analysis on the hyperparemeters for the problems with $n=3$, where $q=1$ in the left table and $q=3$ in the right table.**
| $\gamma_{1}$ | $\gamma_{2}$ | $\alpha_{k}$ | $\eta_{k}$ | $\rho$ | Time (s) |
| --- | --- | --- | --- | --- | --- |
| 1 | 1.0 | 0.001 | 0.01 | 0.3 | 2.96 |
| 3 | 1.0 | 0.001 | 0.01 | 0.3 | 2.84 |
| 5 | 1.0 | 0.001 | 0.01 | 0.3 | 2.95 |
| 7 | 1.0 | 0.001 | 0.01 | 0.3 | 2.63 |
| 10 | 0.1 | 0.001 | 0.01 | 0.3 | 3.02 |
| 10 | 0.3 | 0.001 | 0.01 | 0.3 | 3.14 |
| 10 | 0.5 | 0.001 | 0.01 | 0.3 | 2.98 |
| 10 | 0.7 | 0.001 | 0.01 | 0.3 | 3.09 |
| 10 | 1.0 | 0.0005 | 0.005 | 0.3 | 4.50 |
| 10 | 1.0 | 0.0007 | 0.007 | 0.3 | 3.68 |
| 10 | 1.0 | 0.003 | 0.03 | 0.3 | 0.83 |
| 10 | 1.0 | 0.005 | 0.05 | 0.3 | 0.42 |
| 10 | 1.0 | 0.001 | 0.01 | 0.1 | 0.59 |
| 10 | 1.0 | 0.001 | 0.01 | 0.2 | 1.76 |
| 10 | 1.0 | 0.001 | 0.01 | 0.4 | 2.66 |
| 10 | 1.0 | 0.001 | 0.01 | 0.49 | 2.00 |
| 10 | 1.0 | 0.001 | 0.01 | 0.3 | 2.20 |

To see how parameters impact the performance of the algorithms, we conduct a series of sensitivity analyses, the results are collected in Table 5. Besides, we also conduct experiments to investigate the relationship between the problem’s dimension and the time cost, the results are shown in Figure 2. For these expriments, we set $\gamma_{1}=10$, $\gamma_{2}=1$, $p=0.3$, and choose $(\alpha_{k},\eta_{k})=(20/n,2/n)$ for problems with $q=1$, $(\alpha_{k},\eta_{k})=(10/n,1/n)$ for problems with $q=3$. Such a choice is quite rough, but it is sufficient to demonstrate the scalability of BiC-GAFFA.

**Table 5: Sensitivity analysis on the hyperparameters for the synthetic experiments.**
| Settings | $n=100$, $q=1$ | $n=100$, $q=3$ |
| --- | --- | --- |
| $\begin{aligned} \gamma_{2}=&1\\<br>\alpha_{k}=&0.001\\<br>\eta_{k}=&0.01\\<br>\rho=&0.3\end{aligned}$ |  |  |
| $\begin{aligned} \gamma_{1}=&10\\<br>\alpha_{k}=&0.001\\<br>\eta_{k}=&0.01\\<br>\rho=&0.3\end{aligned}$ |  |  |
| $\begin{aligned} \gamma_{1}=&10\\<br>\gamma_{2}=&1\\<br>\eta_{k}=&0.01\\<br>\alpha_{k}=&0.001\end{aligned}$ |  |  |
| $\begin{aligned} \gamma_{1}=&10\\<br>\gamma_{2}=&1\\<br>\rho=&0.3\end{aligned}$ |  |  |

#### A.2.2 Hyperparameter optimization

Hyperparameter optimization is an inherent bilevel optimization.
According to Theorem 3.1 of Gao’s work [], we can turn the hyperparameter optimization problem of a statistical learning model
which consists of data fitting $\mathcal{L}_{\text{data}}(\mathbf{x})$ and convex regularization terms $P_{i}(\mathbf{x})$ to the following constrained bilevel optimization problem.

$$ $\min_{\mathbf{x}\in\mathbb{R}^{n},\mathbf{r}\in\mathbb{R}^{J}}\mathcal{L}_{ \text{val}}(\mathbf{x})\ \text{ s.t. }\ \mathbf{x}\in\mathop{\mathrm{argmin}}_ {\mathbf{x}\in\mathbb{R}^{n}}\Big{\{}\mathcal{L}_{\text{val}}(\mathbf{x})\ \text{ s.t. }\ P_{i}(\mathbf{x})\leq r_{i},\ i=1,\cdots,J\Big{\}}.$ $$

Sparse Group LASSO Problem.
The sparse group LASSO problem [] is an advanced statistical learning model, which aims to find the grouped structure of the predictors and select the relevant ones. It is critical to select the weights of the regularization terms. The direct form of the problem is as follows:

$$ $\begin{array}[]{>{\displaystyle}c @{\,\,} >{\displaystyle}l}\min_{\bm{\beta} \in\mathbb{R}^{p},\bm{\lambda}\in\mathbb{R}_{+}^{M+1}}&\frac{1}{2}\sum_{i\in I _{\text{val}}}|b_{i}-\bm{\beta}^{\text{{T}}}\mathbf{a}_{i}|^{2}\\ \text{s.t. }&\bm{\beta}\in\mathop{\arg\min}_{\hat{\bm{\beta}}\in\mathbb{R}^{p} }\bigg{\{}\frac{1}{2}\sum_{i\in I_{\text{tr}}}|b_{i}-\hat{\bm{\beta}}^{\text{{ T}}}\mathbf{a}_{i}|^{2}+\sum_{m=1}^{M}\lambda_{m}\|\hat{\bm{\beta}}^{(m)}\|_{2 }+\lambda_{M+1}\|\hat{\bm{\beta}}\|_{1}\bigg{\}}.\end{array}$ (21) $$

By decoupling the data fitting term and regularization terms, we can reformulate the problem to the following bilevel optimization problem:

$$ $\begin{array}[]{>{\displaystyle}c @{\,\,} >{\displaystyle}l}\min_{\bm{\beta} \in\mathbb{R}^{p},\mathbf{r}\in\mathbb{R}_{+}^{M+1}}&\frac{1}{2}\sum_{i\in I_{ \text{val}}}|b_{i}-\bm{\beta}^{\text{{T}}}\mathbf{a}_{i}|^{2}\\ \text{s.t. }&\bm{\beta}\in\begin{array}[t]{>{\displaystyle}c @{\,\,} >{ \displaystyle}l}\mathop{\arg\min}_{\hat{\bm{\beta}}\in\mathbb{R}^{p}}&\frac{1} {2}\sum_{i\in I_{\text{tr}}}|b_{i}-\hat{\bm{\beta}}^{\text{{T}}}\mathbf{a}_{i} |^{2}\\ \text{ s.t. }&\|\hat{\bm{\beta}}^{(m)}\|_{2}\leq r_{m},m=1,\dots,M,\quad\|\hat {\bm{\beta}}\|_{1}\leq r_{M+1}.\end{array}\end{array}$ (22) $$

From our practice, we find that the proposed BiC-GAFFA algorithm works better on the squared two norms rather than two norms directly, so by introducing $\mathbf{u}\in\mathbb{R}_{+}^{M+1}$ such that $u_{m}=r_{m}^{2}$ for $m=1,\dots,M$, and $u_{M+1}=r_{M+1}$, we can reformulate the problem to the following bilevel optimization problem:

$$ $\begin{array}[]{>{\displaystyle}c @{\,\,} >{\displaystyle}l}\min_{\bm{\beta} \in\mathbb{R}^{p},\mathbf{u}\in\mathbb{R}_{+}^{M+1}}&\frac{1}{2}\sum_{i\in I_{ \text{val}}}|b_{i}-\bm{\beta}^{\text{{T}}}\mathbf{a}_{i}|^{2}\\ \text{s.t. }&\bm{\beta}\in\begin{array}[t]{>{\displaystyle}c @{\,\,} >{ \displaystyle}l}\mathop{\arg\min}_{\hat{\bm{\beta}}\in\mathbb{R}^{p}}&\frac{1} {2}\sum_{i\in I_{\text{tr}}}|b_{i}-\hat{\bm{\beta}}^{\text{{T}}}\mathbf{a}_{i} |^{2}\\ \text{ s.t. }&\|\hat{\bm{\beta}}^{(m)}\|^{2}_{2}\leq u_{m},m=1,\dots,M,\quad\| \hat{\bm{\beta}}\|_{1}\leq u_{M+1}.\end{array}\end{array}$ (23) $$

Grid Search, Random Search, Bayesian Optimization and VF-iDCA are applied to the model (22), IGJO is applied to the model (21), and BiC-GAFFA is applied to the model (23). The experimental results are collected in Table 6.

Data used in these experiments are generated by the following procedures: we random smapled $\mathbf{a}_{i}\in\mathbb{R}^{p}$ from the standard normal distribution, and set the true weights $\bm{\beta}$ to be a grouped sparse vector, specifically, it was defined as $\bm{\beta}=[\bm{\beta}^{(1)},\bm{\beta}^{(2)},\dots,\bm{\beta}^{(5)}]$, where the first 5 entries of $\bm{\beta}^{(i)}$ are 1, 2, 3, 4, 5, and the rest entries are all zeros. The responses are defined as $b_{i}=\bm{\beta}^{T}\mathbf{a}_{i}+\sigma\epsilon_{i}$, where the $\epsilon_{i}$ are generated from the standard normal distribution, and $\sigma$ is chosen such that the signal-to-noise ratio is 3. The training set, validation set, and test set are split randomly, and the size of the training set is denoted as nTr, the size of the validation set is nVal, and the size of the test set is nTest. For all the experiments in this part, $p=150$, $M=30$.

The detailed settings of each algorithm are provided as follows:

- •
For Grid Search, we set the grid size to be 20, the range of $r_{m}$ is $[1,10]$ for $m=1,\dots,M$, and the range of $r_{M+1}$ is $[1,100]$;
- •
For Random Search, we set the number of iterations to be 400, the range of $r_{m}$ is $[0,10]$ for $m=1,\dots,M$, and the range of $r_{M+1}$ is $[0,100]$, uniformly sampled method is adopted;
- •
For TPE, we set the number of iterations to be 400, and use uniform distribution for each hyperparameter, the range of $r_{m}$ is $[0,10]$ for $m=1,\dots,M$, and the range of $r_{M+1}$ is $[0,100]$;
- •
For IGJO, we set the number of iterations to be 100 for each hyperparameter, and the initial guess is $0.1\times\mathds{1}_{M+1}$;
- •
For VF-iDCA, we set the number of iterations to 50 for each hyperparameter. As to the initial guess, we firstly solve the lower level problem of (21) with hyperparameters $0.1\times\mathds{1}_{M+1}$, denote that $r_{m}=\|\hat{\bm{\beta}}^{(m)}\|_{2}$ for $m=1,\dots,M$, $r_{M+1}=\|\hat{\bm{\beta}}\|_{1}$, then take the $\mathbf{r}$ as the initial guess for the bilevel problem of (22);
- •
For BiC-GAFFA, we set the number of iterations to be 30000 for each hyperparameter. As to the initial guess, we firstly solve the lower level problem of (21) with hyperparameters $0.1\times\mathds{1}_{M+1}$, denote that $u_{m}=\|\hat{\bm{\beta}}^{(m)}\|_{2}^{2}$ for $m=1,\dots,M$, $u_{M+1}=\|\hat{\bm{\beta}}\|_{1}$, then take the $\mathbf{u}$ as the initial guess for the bilevel problem of (23). For this problem, we always take $\gamma_{1}=10$, $\gamma_{2}=1$, $\eta_{k}=0.1$, $\alpha_{k}=0.01$, and $\rho=0.3$.

**Table 6: Results on the sparse group Lasso hyperparameter selection problem.**
| Method | nTr = 100, nVal = 100, nTest = 300 | nTr = 300, nVal = 300, nTest = 300 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Time (s) | Val Err | Test Err | Time (s) | Val Err | Test Err |  |
| Grid | 17.3 $\pm$ 0.9 | 35.9 $\pm$ 7.2 | 37.7 $\pm$ 6.7 | 78.7 $\pm$ 1.9 | 18.9 $\pm$ 2.3 | 19.8 $\pm$ 1.8 |
| Random | 17.4 $\pm$ 0.7 | 33.6 $\pm$ 6.7 | 35.7 $\pm$ 6.2 | 78.6 $\pm$ 2.5 | 18.7 $\pm$ 2.4 | 19.5 $\pm$ 1.9 |
| TPE | 16.9 $\pm$ 0.7 | 33.9 $\pm$ 7.0 | 36.0 $\pm$ 5.6 | 74.7 $\pm$ 2.2 | 18.9 $\pm$ 2.3 | 19.8 $\pm$ 1.9 |
| IGJO | 21.2 $\pm$ 2.2 | 19.7 $\pm$ 2.8 | 25.6 $\pm$ 4.4 | 49.9 $\pm$ 2.6 | 16.5 $\pm$ 2.5 | 18.1 $\pm$ 1.4 |
| VF-iDCA | 12.4 $\pm$ 0.5 | 14.6 $\pm$ 2.6 | 25.4 $\pm$ 3.9 | 40.7 $\pm$ 1.7 | 14.9 $\pm$ 2.1 | 17.2 $\pm$ 1.3 |
| BiC-GAFFA | 21.4 $\pm$ 0.7 | 7.3 $\pm$ 1.3 | 22.3 $\pm$ 3.0 | 22.0 $\pm$ 1.0 | 12.8 $\pm$ 1.4 | 17.1 $\pm$ 1.3 |

From Table 6, we can see that BiC-GAFFA outperforms other algorithms in terms of the validation error and test error, and when the number of data is increased, the time cost of BiC-GAFFA does not increase significantly, which indicates that BiC-GAFFA is scalable to large-scale problems.

Support Vector Machine.
The mathematical model can be written as follows:

$$ $\begin{array}[]{>{\displaystyle}c @{\,\,} >{\displaystyle}l}\min_{\mathbf{c} \in\mathbb{R}^{N_{\text{tr}}},\mathbf{w}\in\mathbb{R}^{p},b\in\mathbb{R},\bm{ \xi}\in\mathbb{R}^{N_{\text{tr}}}}&\mathcal{L}_{\mathcal{D}_{\text{val}}}( \mathbf{w},b),\\ \text{ s.t. }&(\mathbf{w},b,\bm{\xi})\in\begin{array}[t]{>{\displaystyle}c @{ \,\,} >{\displaystyle}l}\mathop{\text{argmin}}_{\mathbf{w}\in\mathbb{R}^{p},b \in\mathbb{R},\bm{\xi}\in\mathbb{R}^{N_{\text{tr}}}}&\frac{1}{2}\|\mathbf{w}\| ^{2}+\frac{1}{2}\sum_{i=1}^{N_{tr}}e^{c_{i}}\xi_{i}^{2}\\ \text{ s.t. }&y_{i}(\mathbf{w}^{\text{{T}}}\mathbf{x}_{i}+b)\geq 1-\xi_{i}, \quad i=1,\ldots,N_{\text{tr}}.\end{array}\end{array}$ $$

Here we use $(\mathcal{D}_{\text{tr}},\mathcal{D}_{\text{val}})$ to denote the split of the training set and validation set of data, and $N_{\text{tr}}$ to denote the number of data in the training set. The upper-level objective function is defined in the following way:

$$ $\mathcal{L}_{\mathcal{D}_{\text{val}}}(\mathbf{w},b):=\frac{1}{|\mathcal{D}_{ \text{val}}|}\sum_{(\mathbf{x}_{i},y_{i})\in\mathcal{D}_{\text{val}}}\text{ Sigmoid}\left(-\frac{y_{i}(\mathbf{w}^{\text{{T}}}\mathbf{x}_{i}+b)}{\|\mathbf {w}\|^{2}}\right).$ $$

The function inner the sigmoid function gives an opposite of a signed distance between point $(\mathbf{x}_{i},y_{i})$ and the decision hyperplane $\{\mathbf{x}\in\mathbb{R}^{p}|\mathbf{w}^{\text{{T}}}\mathbf{x}+b=0\}$. Specifically, the inner part is positive when the sign of the prediction $\mathbf{w}^{\text{{T}}}\mathbf{x}_{i}+b$ is different with its label $y_{i}$, and negative otherwise. And the sigmoid function converts the distance values to some probability value. Such a composition makes the objective function to be a smooth approximation of the validation accuracy.

The dataset we used in this experiment is the scaled diabetes dataset provided in the repository [], which contains 768 data points and 8 features. In each experiment, the training set, validation set, and test set are split randomly, and the size of the training set is 400, and the size of the validation set is 150, the rest part is the test set.

The initial values of $\mathbf{c}$ in this problem are sampled from a uniform distribution on $[-6,-5]$, the initial values for other parameters are solutions of the lower level problem with hyperparameters $\mathbf{c}$. The hyperparameters for all the three algorithms are:

- •
For GAM, we keep the same setting as the original paper [] did, the maximal iteration number is set to be 80.
- •
For LV-HBA, we use $\alpha=0.001$, $\eta=0.001$, $\gamma_{1}=0.1$, $\gamma_{2}=0.1$, the maiximal iteration number is set to be 400.
- •
For BiC-GAFFA, we use $\gamma_{1}=10$, $\gamma_{2}=0.01$, $\eta_{k}=0.01$, $\alpha_{k}=0.001$, $\rho=0.3$, the maximal iteration number is set to be 5000.

Data Hyperclean.
For this problem, we use the same model as the one in the SVM experiments, but we change the dataset to the scaled breast-cancer dataset provided in the repository [], which contains 683 data points and 10 features. In each experiment, the training set, validation set, and test set are split randomly, and the size of the training set is 400, and the size of the validation set is 180, the rest part is the test set. For the training dataset, we change the label with probability $p_{c}=40\%$, and the validation set and test set are kept unchanged.
This means the data used in lower training is not reliable, the researcher need to find out which data are more reliable and give them higher weights while giving the unreliable data lower weights. Such a weighting process can be viewed as a data clean procedure. Here we do such a procedure by regarding the weights as hyperparameters and evaluating their effects on the validation set, such a process is regarded as data hyper cleaning. Therefore, such a problem is quite similar to the aforementioned SVM, the main difference occurs in the training data.

The hyperparameters for all the three algorithms are:

- •
For GAM, we keep the same setting as the original paper [] did, the maximal iteration number is set to be 30. The initial values of $\mathbf{c}$ are chosen from a uniformed distribution $(1.,2.)$.
- •
For LV-HBA, we use $\alpha=0.001$, $\eta=0.001$, $\gamma_{1}=0.1$, $\gamma_{2}=0.1$, the maiximal iteration number is set to be 400. The initial values of $\mathbf{c}$ are chosen from a uniformed distribution $(-5,-4)$;
- •
For BiC-GAFFA, we use $\gamma_{1}=10$, $\gamma_{2}=0.01$, $\eta_{k}=0.1$, $\alpha_{k}=0.001$, $p=0.3$, the maximal iteration number is set to be 2000. The initial values of $\mathbf{c}$ are choosen from a unifromed distribution $(-5,-4)$.

#### A.2.3 Generative adversarial network

We also apply our algorithm to the GAN, which is a popular model in the field of deep learning. It can be written as the following bilevel optimization problem:

$$ $\min_{G,D}\ \mathcal{L}_{\text{gen}}(G,D)\quad\text{s.t.}\quad D\in\mathop{ \text{argmin}_{D\in\mathcal{D}}}\mathcal{L}_{\text{det}}(G,D)$ $$

The main idea of GAN is to find a network that learns the distribution of the given training data.
This paper investigates the fitting of two-dimensional distributions (i.e., 8-Gaussians model and 25-Gaussians model, refer to []) using GANs.
For such simple distributions, we can approximate the earth mover’s distance (EMD) between the generated distribution and the true distribution by calculating the Sinkhorn distance[].
This allows us to observe the optimization progress of the network during the iteration process, and make comparisons between different strategies.

Figure: Figure 6: How EMD changes w.r.t time for different GAN models. The true distribution in the left picture is the 8 Gaussians mixture model, and the true distribution in the right picture is the 25 Gaussians mixture model.
Refer to caption: x15.png

For all the 8 GANs we make a comparison here, we use similar structures and similar parameters. The generator and discriminator are both three-layer neural networks with 128 hidden units, and the activation function except the last layer is ReLU. The output size of the generator is 256. For the WGAN and WGAN-GP, the last layer of the discriminator is the linear layer while other models have a sigmoid activation function after the linear layer in their discriminator.
The batch size is 512, and the number of iterations is 2701. According to the practical advice in [], we train the discriminator 5 steps and then the generator 1 step in 1 loop. We use Adam with $b1=0.5,b2=0.999$ for all the learnable parameters in all the models. Other parameters are listed below:

- •
GAN[]: we set learning rate of generator as $10^{-3}$, learning rate of discriminator as $10^{-4}$.
- •
WGAN[arjovsky2017wasserstein]: we set learning rate of generator as $10^{-3}$, learning rate of discriminator as $10^{-4}$, the clip value as 0.001;
- •
WGAN-GP[]: we set learning rate of generator as $10^{-4}$, learning rate of discriminator as $10^{-4}$, lambda of gradient penalty as 0.1;
- •
Con-GAN[]: we set learning rate of generator as $10^{-3}$, learning rate of discriminator as $10^{-4}$, lambda of constant penalty as 0.3;
- •
UGAN[]: we set the learning rate of the generator as $10^{-3}$, the learning rate of the discriminator as $10^{-4}$, the unrolled step as 5 and let it update the discriminator just 1 step in 1 loop;
- •
Bi-GAN (BiC-GAFFA): we set the learning rate of the generator as $10^{-3}$, the learning rate of the discriminator as $10^{-4}$, $\rho=0.1$, the learning rate of the auxiliary generator as $10^{-3}$;
- •
Bi-WGAN (BiC-GAFFA): we set the learning rate of the generator as $10^{-3}$, the learning rate of the discriminator as $10^{-4}$, $\rho=0.1$, the learning rate of the auxiliary generator as $10^{-3}$, the learning rate of the auxiliary dual variable as $10^{-4}$, the upper bound for the auxiliary dual variable as 0.1. Note here we turn the gradient penalty introduced in [] to one constraint by taking a maximum over the sample points, and by our practice, the original GAN structure seems more stable with this optimization method, further study is still required;
- •
Bi-ConGAN (BiC-GAFFA): we set the learning rate of the generator as $10^{-3}$, the learning rate of the discriminator as $10^{-4}$, $\rho=0.1$, the learning rate of the auxiliary generator as $10^{-3}$, the learning rate of the auxiliary dual variable as $10^{-4}$, the upper bound for the auxiliary dual variable as 0.1. We realize the constraint proposed in [] with $\varepsilon=0.1$.

Different analyses and assumptions lead to various regularization requirements for these models.
When constrained optimizers are unavailable, such a regularization method can only be achieved by penalizing them to the objective[, ], potentially compromising the interpretability of the model and complicating analysis. Our algorithms mitigate these issues without significant computational overhead. The effectiveness of our algorithms is demonstrated in the numerical results presented in Figures 6.

### A.3 Proofs in Section 2

#### A.3.1 Proof of Lemma 2.1

Lemma 2.1 can be proved using proof techniques similar to those in Theorem 3.3 from []. For completeness, we present an alternative proof of Lemma 2.1 here.

###### Lemma 2.1 .

Assume that both $f(x,\cdot)$ and $g(x,\cdot)$ are convex for any given $x\in X$. Let $\gamma_{1},\gamma_{2}>0$, we have $\mathcal{G}_{\gamma}(x,y,z)\geq 0$ for any $(x,y,z)\in X\times Y\times\mathbb{R}_{+}^{p}$. Furthermore,

$$ $\mathcal{G}_{\gamma}(x,y,z)\leq 0,$ $$

if and only if $y\in S(x)$ and $z\in\mathcal{M}(x,y)$, where $\mathcal{M}(x,y)$ denotes the set of multipliers of the lower-level problem at $(x,y)$, i.e.,

$$ $\mathcal{M}(x,y):=\left\{\lambda\in\mathbb{R}^{p}_{+}~{}|~{}0\in\nabla_{y}f(x, y)+\lambda^{\text{{T}}}\nabla_{y}g(x,y)+\mathcal{N}_{Y}(y),~{}\lambda^{\text{{ T}}}g(x,y)=0\right\}.$ $$

###### Proof.

For any $(x,y,z)\in X\times Y\times\mathbb{R}_{+}^{p}$, we have

$$ $\displaystyle\min_{\theta\in Y}\left\{f(x,\theta)+z^{\text{{T}}}g(x,\theta)+ \frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\}$ (24) $\displaystyle\leq$ $\displaystyle f(x,y)+z^{\text{{T}}}g(x,y)$ $\displaystyle\leq$ $\displaystyle\max_{\lambda\in\mathbb{R}^{p}_{+}}\left\{f(x,y)+\lambda^{\text{{ T}}}g(x,y)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right\},$ $$

which implies that

$$ $\displaystyle\mathcal{G}_{\gamma}(x,y,z)=$ $\displaystyle\max_{\lambda\in\mathbb{R}^{p}_{+}}\left\{f(x,y)+\lambda^{\text{{ T}}}g(x,y)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right\}$ $\displaystyle-\min_{\theta\in Y}\left\{f(x,\theta)+z^{\text{{T}}}g(x,\theta)+ \frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\}$ $\displaystyle\geq$ $\displaystyle 0.$ $$

Therefore, $\mathcal{G}_{\gamma}(x,y,z)=0$ if and only if

$$ $\max_{\lambda\in\mathbb{R}^{p}_{+}}\left\{f(x,y)+\lambda^{\text{{T}}}g(x,y)- \frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right\}=\min_{\theta\in Y}\left\{f(x, \theta)+z^{\text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\}.$ $$

Then, (24) yields that $\mathcal{G}_{\gamma}(x,y,z)=0$ if and only if

$$ $\displaystyle y$ $\displaystyle\in\underset{\theta\in Y}{\mathrm{argmin}}\left\{f(x,\theta)+z^{ \text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\},$ (25) $\displaystyle z$ $\displaystyle\in\underset{\lambda\in\mathbb{R}^{p}_{+}}{\mathrm{argmax}}\left \{f(x,y)+\lambda^{\text{{T}}}g(x,y)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2} \right\}.$ $$

Given that the function $f(x,\theta)+z^{\text{{T}}}g(x,\theta)$ is convex with respect to variable $\theta$, and that $\lambda^{\text{{T}}}g(x,y)$ is concave with respect to $\lambda$, (25) is equivalent to

$$ $\displaystyle y$ $\displaystyle\in\underset{\theta\in Y}{\mathrm{argmin}}\left\{f(x,\theta)+z^{ \text{{T}}}g(x,\theta)\right\},$ (26) $\displaystyle z$ $\displaystyle\in\underset{\lambda\in\mathbb{R}^{p}_{+}}{\mathrm{argmax}}\left \{f(x,y)+\lambda^{\text{{T}}}g(x,y)\right\},$ $$

which is equivalent to that $(y,z)$ is a saddle point to

$$ $\min_{\theta\in Y}\max_{\lambda\in\mathbb{R}^{p}_{+}}~{}f(x,\theta)+\lambda^{ \text{{T}}}g(x,\theta).$ $$

Consequently, applying the classical minimax theorem to this convex-concave min-max problem, we obtain that $\mathcal{G}_{\gamma}(x,y,z)=0$ if and only if $y\in S(x)$ and $z\in\mathcal{M}(x,y)$.
∎

###### Lemma 2.1 .

###### Proof.

#### A.3.2 Proof of Lemma 2.2

###### Lemma 2.2 .

Assume that both $f(x,y)$ and $g(x,y)$ are convex in $y$ on $Y$ for any given $x\in X$ and are continuously differentiable on an open set containing $X\times Y$. Then $\mathcal{G}_{\gamma}(x,y,z)$ is continuously differentiable on $X\times Y\times\mathbb{R}^{p}_{+}$, and for any $(x,y,z)\in X\times Y\times\mathbb{R}^{p}_{+}$,

$$ $\nabla\mathcal{G}_{\gamma}(x,y,z)=\begin{pmatrix}\nabla_{x}f(x,y)+(\lambda^{*} )^{\text{{T}}}\nabla_{x}g(x,y)\\ \nabla_{y}f(x,y)+(\lambda^{*})^{\text{{T}}}\nabla_{y}g(x,y)\\ -\left(z-\lambda^{*}\right)/\gamma_{2}\end{pmatrix}-\begin{pmatrix}\nabla_{x}f (x,\theta^{*})+z^{\text{{T}}}\nabla_{x}g(x,\theta^{*})\\ \left(y-\theta^{*}\right)/\gamma_{1}\\ g(x,\theta^{*})\end{pmatrix},$ $$

where $\theta^{*}$ and $\lambda^{*}$ denote $\theta^{*}(x,y,z)$ and $\lambda^{*}(x,y,z)$, respectively, defined as

$$ $\displaystyle\theta^{*}(x,y,z)$ $\displaystyle:=\underset{\theta\in Y}{\mathrm{argmin}}\left\{f(x,\theta)+z^{ \text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\},$ $\displaystyle\lambda^{*}(x,y,z)$ $\displaystyle:=\underset{\lambda\in\mathbb{R}^{p}_{+}}{\mathrm{argmax}}\left\{ f(x,y)+\lambda^{\text{{T}}}g(x,y)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right \}=\mathrm{Proj}_{\mathbb{R}^{p}_{+}}\left(z+\gamma_{2}g(x,y)\right).$ $$

###### Proof.

We first define two auxiliary functions,

$$ $\displaystyle\mathcal{G}_{1,\gamma}(x,y,z)$ $\displaystyle:=\min_{\theta\in Y}\left\{f(x,\theta)+z^{\text{{T}}}g(x,\theta)+ \frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\},$ $\displaystyle\mathcal{G}_{2,\gamma}(x,y,z)$ $\displaystyle:=\max_{\lambda\in\mathbb{R}^{p}_{+}}\left\{f(x,y)+\lambda^{\text {{T}}}g(x,y)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right\}.$ $$

Then, it follows that $\mathcal{G}_{\gamma}(x,y,z)=\mathcal{G}_{2,\gamma}(x,y,z)-\mathcal{G}_{1,
\gamma}(x,y,z)$. Since by assumptions that $f$ and $g$ are both continuous differentiable on an open set containing $X\times Y$, it can be easily shown that $f(x,\theta)+z^{\text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}$ satisfies the inf-compactness condition in Theorem 4.13 of [] with respect to $\theta\in Y$ on any point $(\bar{x},\bar{y},\bar{z})\in X\times Y\times\mathbb{R}^{p}_{+}$, i.e., for any $(\bar{x},\bar{y},\bar{z})\in X\times Y\times\mathbb{R}^{p}_{+}$, there exist $c\in\mathbb{R}$, compact set $D$ and neighborhood $W$ of $(\bar{x},\bar{y},\bar{z})$ such that the level set $\{\theta\in Y~{}|~{}f(x,\theta)+z^{\text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}
}\|\theta-y\|^{2}\leq c\}$ is nonempty and contained in $D$ for any $(x,y,z)\in W$.
And because of the convexity of $f(x,y)$ and $g(x,y)$ with respect to $y\in Y$ for any $x\in X$, $\underset{\theta\in Y}{\mathrm{argmin}}\left\{f(x,\theta)+z^{\text{{T}}}g(x,
\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\}$ is a singleton for any $(x,y,z)\in X\times Y\times\mathbb{R}^{p}_{+}$.
Then, by the differentiablility of $f$ and $g$, we can derive from Theorem 4.13, Remark 4.14 of [] that $\mathcal{G}_{1,\gamma}(x,y,z)$ is differentiable at any point on $X\times Y\times\mathbb{R}^{p}_{+}$ and for any $(x,y,z)\in X\times Y\times\mathbb{R}^{p}_{+}$,

$$ $\nabla\mathcal{G}_{1,\gamma}(x,y,z)=\left(\nabla_{x}f(x,\theta^{*})+z^{\text{{ T}}}\nabla_{x}g(x,\theta^{*}),\left(y-\theta^{*}\right)/\gamma_{1},g(x,\theta^ {*})\right).$ (27) $$

By using similar arguments, we can also demonstrate that
$\mathcal{G}_{2,\gamma}(x,y,z)$ is differentiable at any point on $X\times Y\times\mathbb{R}^{p}_{+}$ and for any $(x,y,z)\in X\times Y\times\mathbb{R}^{p}_{+}$,

$$ $\nabla\mathcal{G}_{2,\gamma}(x,y,z)=\left(\nabla_{x}f(x,y)+(\lambda^{*})^{ \text{{T}}}\nabla_{x}g(x,y),\nabla_{y}f(x,y)+(\lambda^{*})^{\text{{T}}}\nabla_ {y}g(x,y),-\left(z-\lambda^{*}\right)/\gamma_{2}\right).$ (28) $$

And then the conclusion follows from $\mathcal{G}_{\gamma}(x,y,z)=\mathcal{G}_{2,\gamma}(x,y,z)-\mathcal{G}_{1,
\gamma}(x,y,z)$.
∎

###### Lemma 2.2 .

###### Proof.

#### A.3.3 Proof of Theorem 2.3

###### Theorem 2.3 .

Assume that both $f(x,\cdot)$ and $g(x,\cdot)$ are convex for any given $x\in X$. Let $\gamma_{1},\gamma_{2}>0$, the reformulation (6) is equivalent to the bilevel optimization problem (1), provided that for any feasible point $(x,y)$ of (1), a corresponding multiplier of the lower-level problem (2) exists at $(x,y)$, i.e., $\mathcal{M}(x,y)\neq\varnothing$.

###### Proof.

Let $(x,y,z)$ be any feasible point of problem (6), then we have $(x,y)\in X\times Y$, $z\in\mathbb{R}^{p}_{+}$, $\mathcal{G}_{\gamma}(x,y,z)\leq 0$. And it follows from Lemma 2.1 that $\mathcal{G}_{\gamma}(x,y,z)=0$, $y\in S(x)$
and thus $(x,y)$ is feasible to bilevel program (1).

Now, suppose $(x,y)$ is an feasible point of bilevel program (1), then we have $(x,y)\in X\times Y$ and $y\in S(x)$.
According to the assumption that a multiplier $z\in\mathbb{R}^{p}_{+}$ of the lower-level problem (2) exists at $(x,y)$, i.e., $z\in\mathcal{M}(x,y)$. Then it follows from Lemma 2.1 that $\mathcal{G}_{\gamma}(x,y,z)=0$ and
thus $(x,y,z)$ is feasible to reformulation problem (6).
∎

###### Theorem 2.3 .

###### Proof.

### A.4 Proofs in Section 3

#### A.4.1 Proof of Proposition 3.1

###### Proposition 3.1 .

Suppose $\gamma_{1},\gamma_{2}>0$ and an optimal solution $(x^{*},y^{*},z^{*})$ to (6), with $z^{*}\in Z$, exists,
then any optimal solution of (7) is optimal to reformulation (6).

###### Proof.

For any feasible point $(x,y,z)$ of problem (7), $(x,y,z)$ is also feasible to problem (6) and thus the optimal value of problem (7) is larger or equal to that of problem (6). Let $(x^{*},y^{*},z^{*})$ be an optimal solution of reformulation problem (6) with $z^{*}$ belonging to the set $Z$, then $(x^{*},y^{*},z^{*})$ is also feasible to problem (7). Therefore, the optimal value of problem (7) is equal to that of problem (6). Then, because any feasible point $(x,y,z)$ of problem (7) is feasible to problem (6), we get the conclusion.
∎

###### Proposition 3.1 .

###### Proof.

#### A.4.2 Proof of Proposition 3.2

###### Proposition 3.2 .

Assume that $F(x,y)$ is bounded below by $\underline{F}$ on $X\times Y$. For any $\varepsilon>0$, there exists $\bar{c}>0$ such that any global solution $(x_{c},y_{c},z_{c})$ to the penalty formulation (8) with penalty parameter $c\geq\bar{c}$ is also a global solution to the relaxed problem (9) with some relaxation parameter $\varepsilon_{c}$ satisfying $\varepsilon_{c}\leq\varepsilon$. Moreover, if $(x_{c},y_{c},z_{c})$ is a local solution to the penalty formulation (8), then it is also a local solution to the relaxed problem (9) with relaxation parameter $\varepsilon_{c}:=\mathcal{G}_{\gamma}(x_{c},y_{c},z_{c})$.

###### Proof.

Let $(\bar{x},\bar{y},\bar{z})\in X\times Y\times Z$ be a feasible point to problem (7) and $(x_{c},y_{c},z_{c})\in X\times Y\times Z$ be a global solution of problem (8) with penalty parameter $c>0$. We then have

$$ $F(x_{c},y_{c})+c\,\mathcal{G}_{\gamma}(x_{c},y_{c},z_{c})\leq F(\bar{x},\bar{y })+c\,\mathcal{G}_{\gamma}(\bar{x},\bar{y},\bar{z})=F(\bar{x},\bar{y}),$ $$

implying

$$ $\mathcal{G}_{\gamma}(x_{c},y_{c},z_{c})\leq(F(\bar{x},\bar{y})-\underline{F})/c.$ $$

Thus, for any $\varepsilon>0$, there exists $\bar{c}>0$ such that for any $c\geq\bar{c}$,

$$ $\varepsilon_{c}:=\mathcal{G}_{\gamma}(x_{c},y_{c},z_{c})\leq\varepsilon.$ $$

Next, we demonstrate that $(x_{c},y_{c},z_{c})$ is a global solution to problem (9) with relaxation parameter $\varepsilon_{c}$. Assume, for the sake of contradiction, that there exists $(x,y,z)\in X\times Y\times Z$ with $\mathcal{G}_{\gamma}(x,y,z)\leq\varepsilon_{c}$ and $F(x,y)<F(x_{c},y_{c})$. This leads to

$$ $F(x,y)+c\,\mathcal{G}_{\gamma}(x,y,z)<F(x_{c},y_{c})+c\varepsilon_{c}=F(x_{c}, y_{c})+c\,\mathcal{G}_{\gamma}(x_{c},y_{c},z_{c}),$ $$

which contradicts the global optimality of $(x_{c},y_{c},z_{c})$to problem (8) with penalty parameter $c>0$.

By analogous reasoning, the assertion concerning the local optimality of $(x_{c},y_{c},z_{c})$ for problem (8) holds similarly.
∎

###### Proposition 3.2 .

###### Proof.

#### A.4.3 Proof of Theorem 3.3

###### Theorem 3.3 .

Assume that $X$ and $Y$ are closed and functions $F$, $f$ and $g$ are continuous on $X\times Y$. Suppose $c_{k}\rightarrow\infty$ and let

$$ $(x_{k},y_{k},z_{k})\in\underset{(x,y)\in X\times Y\times Z}{\mathrm{argmin}}F( x,y)+c_{k}\,\mathcal{G}_{\gamma}(x,y,z).$ $$

Then, any accumulation point $(\bar{x},\bar{y},\bar{z})$ of the sequence $\{(x_{k},y_{k},z_{k})\}$ is a solution to problem (7).

###### Proof.

Applying the proof techniques used in Lemma 2 of [], we can establish that $\mathcal{G}_{\gamma}(x,y,z)$ is lower semi-continuous on $X\times Y\times Z$. The theorem’s conclusion follows by employing the same proof techniques from Theorem A.4 of [].
∎

###### Theorem 3.3 .

###### Proof.

### A.5 Proofs in Section 4

The proof of non-asymptotic convergence for BiC-GAFFA primarily hinges on establishing the sufficient descent property of the merit function defined as follows

$$ $V_{k}:=\phi_{c_{k}}(x^{k},y^{k},z^{k})+C_{\theta}\left\|\theta^{k}-\theta^{*}( x^{k},y^{k},z^{k})\right\|^{2},$ $$

where

$$ $\phi_{c_{k}}(x,y,z):=\frac{1}{c_{k}}\big{(}F(x,y)-\underline{F}\big{)}+ \mathcal{G}_{\gamma}(x,y,z),$ $$

and

$$ $C_{\theta}:=L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+L_{g}.$ $$

To establish the sufficiently decreasing property of the merit function, we initially derive several crucial auxiliary lemmas. These lemmas establish the Lipschitz continuity of $\theta^{*}(x,y,z)$ and $\lambda^{*}(x,y,z)$ (as detailed in Lemma A.1), the Lipschitz continuity of $\nabla\mathcal{G}_{\gamma}(x,y,z)$ (as detailed in Lemma A.2) and a descent property with bounded errors for the function $\phi_{c_{k}}(x,y,z)$ at each iteration (as detailed in Lemma A.4).

#### A.5.1 Auxiliary lemmas

The Lipschitz properties of $\theta^{(}x,y,z)$, $\lambda^{(}x,y,z)$, and $\nabla\mathcal{G}_{\gamma}(x,y,z)$ are crucial for the convergence analysis. We establish these properties in the subsequent lemmas.

###### Lemma A.1 .

Under Assumptions 4.2 and 4.3, and let $\gamma_{1}>0$, $\gamma_{2}>0$, for any $({x},{y},z),(x^{\prime},y^{\prime},z^{\prime})\in X\times Y\times Z$, the following inequalities hold

$$ $\displaystyle\|\theta^{*}(x,y,z)-\theta^{*}(x^{\prime},y^{\prime},z^{\prime})\|$ $\displaystyle\leq\left(\gamma_{1}L_{f}+\gamma_{1}rL_{g_{2}}\right)\|x-x^{ \prime}\|+\|y-y^{\prime}\|+\gamma_{1}L_{g}\|z-z^{\prime}\|$ (29) $\displaystyle\leq L_{\theta}\|({x},{y},z)-(x^{\prime},y^{\prime},z^{\prime})\|,$ $\displaystyle\|\lambda^{*}(x,y,z)-\lambda^{*}(x^{\prime},y^{\prime},z^{\prime})\|$ $\displaystyle\leq\gamma_{2}L_{g}\|(x,y)-(x^{\prime},y^{\prime})\|+\|z-z^{ \prime}\|$ $\displaystyle\leq L_{\lambda}\|({x},{y},z)-(x^{\prime},y^{\prime},z^{\prime})\|,$ $$

where $L_{\theta}:=\sqrt{3}\max\{\gamma_{1}L_{f}+\gamma_{1}rL_{g_{2}},1,\gamma_{1}L_{
g}\}$ and $L_{\lambda}:=\sqrt{2}\max\{\gamma_{2}L_{g},1\}$.

###### Proof.

To simplify notation, we denote $(x,y,z),(x^{\prime},y^{\prime},z^{\prime})\in X\times Y\times Z$ by $w$ and $w^{\prime}$, respectively.
Considering that $\theta^{*}(w)$ and $\lambda^{*}(w)$ are optimal solutions to optimization problems in (5),
it follows from the stationary conditions that

$$ $\displaystyle 0$ $\displaystyle\in\nabla_{y}f(x,\theta^{*}(w))+z^{\text{{T}}}\nabla_{y}g(x, \theta^{*}(w))+(\theta^{*}(w)-y)/\gamma_{1}+\mathcal{N}_{Y}(\theta^{*}(w)),$ (30) $\displaystyle 0$ $\displaystyle\in-g(x,y)+(\lambda^{*}(w)-z)/\gamma_{2}+\mathcal{N}_{\mathbb{R}^ {p}_{+}}(\lambda^{*}(w)).$ $$

Same results apply to $\theta^{*}(w^{\prime})$ and $\lambda^{*}(w^{\prime})$

$$ $\displaystyle 0$ $\displaystyle\in\nabla_{y}f(x^{\prime},\theta^{*}(w^{\prime}))+(z^{\prime})^{ \text{{T}}}\nabla_{y}g(x^{\prime},\theta^{*}(w^{\prime}))+(\theta^{*}(w^{ \prime})-y^{\prime})/\gamma_{1}+\mathcal{N}_{Y}(\theta^{*}(w^{\prime})),$ (31) $\displaystyle 0$ $\displaystyle\in-g(x^{\prime},y^{\prime})+(\lambda^{*}(w^{\prime})-z^{\prime}) /\gamma_{2}+\mathcal{N}_{\mathbb{R}^{p}_{+}}(\lambda^{*}(w^{\prime})).$ $$

Defining

$$ $T(x,y,z,\theta):=\nabla_{\theta}\left(f(x,\theta)+z^{\text{{T}}}g(x,\theta)+ \frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right).$ $$

and exploiting the monotonicity of $\mathcal{N}_{Y}$, we have from (30) and (31) that

$$ $\displaystyle\big{\langle}-T(w,\theta^{*}(w))+T(w,\theta^{*}(w^{\prime})), \theta^{*}(w)-\theta^{*}(w^{\prime})\big{\rangle}$ (32) $\displaystyle+\big{\langle}-T(w,\theta^{*}(w^{\prime}))+T(w^{\prime},\theta^{* }(w^{\prime})),\theta^{*}(w)-\theta^{*}(w^{\prime})\big{\rangle}\geq 0$ $$

Under Assumptions 4.2 and 4.3, and given that $\gamma_{1}>0$, it holds that for any $(x,y,z)\in X\times Y\times\mathbb{R}^{p}_{+}$,

$$ $f(x,\theta)+z^{\text{{T}}}g(x,\theta)+\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}$ $$

is $\frac{1}{\gamma_{1}}$-strongly convex with respect to $\theta$. According to [], $T(x,y,z,\theta)$ is $1/\gamma_{1}$-strongly monotone. Consequently, we have that

$$ $\left\langle T(w,\theta^{*}(w))-T(w,\theta^{*}(w^{\prime})),\theta^{*}(w)- \theta^{*}(w^{\prime})\right\rangle\geq\|\theta^{*}(w)-\theta^{*}(w^{\prime}) \|^{2}/\gamma_{1}.$ (33) $$

Under Assumptions 4.2 and 4.3, we establish that

$$ $\displaystyle\|-T(w,\theta^{*}(w^{\prime}))+T(w^{\prime},\theta^{*}(w^{\prime} ))\|$ $\displaystyle\leq$ $\displaystyle\|\nabla_{y}f(x,\theta^{*}(w^{\prime}))-\nabla_{y}f(x^{\prime}, \theta^{*}(w^{\prime}))\|+\|z^{\text{{T}}}\nabla_{y}g(x,\theta^{*}(w^{\prime}) )-(z^{\prime})^{\text{{T}}}\nabla_{y}g(x^{\prime},\theta^{*}(w^{\prime}))\|+\| y-y^{\prime}\|/\gamma_{1}$ $\displaystyle\leq$ $\displaystyle L_{f}\|x-x^{\prime}\|+\|z^{\text{{T}}}\nabla_{y}g(x,\theta^{*}(w ^{\prime}))-z^{\text{{T}}}\nabla_{y}g(x^{\prime},\theta^{*}(w^{\prime}))\|$ $\displaystyle+\|z^{\text{{T}}}\nabla_{y}g(x^{\prime},\theta^{*}(w^{\prime}))-( z^{\prime})^{\text{{T}}}\nabla_{y}g(x^{\prime},\theta^{*}(w^{\prime}))\|+\|y-y ^{\prime}\|/\gamma_{1}$ $\displaystyle\leq$ $\displaystyle L_{f}\|x-x^{\prime}\|+rL_{g_{2}}\|x-x^{\prime}\|+L_{g}\|z-z^{ \prime}\|+\|y-y^{\prime}\|/\gamma_{1},$ $$

where the last inequality follows from the $L_{g_{2}}$-Lipschitz continuity of $\nabla_{y}g$ on $X\times Y$, $\max_{x\in X,y\in Y}\|\nabla_{y}g(x,y)\|\leq L_{g}$, $(x^{\prime},\theta^{*}(w^{\prime}))\in X\times Y$ and $z\in Z$.
Combining this inequality with (32) and (33), we have

$$ $\|\theta^{*}(w)-\theta^{*}(w^{\prime})\|\leq\left(\gamma_{1}L_{f}+\gamma_{1}rL _{g_{2}}\right)\|x-x^{\prime}\|+\|y-y^{\prime}\|+\gamma_{1}L_{g}\|z-z^{\prime}\|.$ $$

Further, exploiting the monotonicity of $\mathcal{N}_{\mathbb{R}^{p}_{+}}$, we obtain from (30) and (31) that

$$ $\left\langle g(x,y)+(z-\lambda^{*}(w))/\gamma_{2}-g(x^{\prime},y^{\prime})-(z^ {\prime}-\lambda^{*}(w^{\prime}))/\gamma_{2},\lambda^{*}(w)-\lambda^{*}(w^{ \prime})\right\rangle\geq 0.$ (34) $$

Then, invoking Assumption 4.3, we have that

$$ $\|\lambda^{*}(w)-\lambda^{*}(w^{\prime})\|\leq\gamma_{2}L_{g}\|(x,y)-(x^{ \prime},y^{\prime})\|+\|z-z^{\prime}\|.$ $$

∎

###### Lemma A.2 .

Under Assumptions 4.2 and 4.3, assume that $X$, $Y$ are compact sets, and $\gamma_{1}>0$, $\gamma_{2}>0$. Then, there exists $L_{\mathcal{G}}>0$ such that for any points $({x},{y},z),(x^{\prime},y^{\prime},z^{\prime})\in X\times Y\times Z$,

$$ $\|\nabla\mathcal{G}_{\gamma}(x,y,z)-\nabla\mathcal{G}_{\gamma}(x^{\prime},y^{ \prime},z^{\prime})\|\leq L_{\mathcal{G}}\|(x,y,z)-(x^{\prime},y^{\prime},z^{ \prime})\|,$ $$

and

$$ $-\mathcal{G}_{\gamma}(x^{\prime},y^{\prime},z^{\prime})\leq-\mathcal{G}_{ \gamma}(x,y,z)-\langle\nabla\mathcal{G}_{\gamma}(x,y,z),(x^{\prime},y^{\prime} ,z^{\prime})-(x,y,z)\rangle+\frac{L_{\mathcal{G}}}{2}\|(x,y,z)-(x^{\prime},y^{ \prime},z^{\prime})\|^{2}.$ $$

###### Proof.

For conciseness, we denote $(x,y,z),(x^{\prime},y^{\prime},z^{\prime})\in X\times Y\times Z$ by $w$ and $w^{\prime}$, respectively.
Recalling from Lemma A.1, we have

$$ $\|\theta^{*}(w)-\theta^{*}(w^{\prime})\|\leq L_{\theta}\|w-w^{\prime}\|,\qquad \|\lambda^{*}(w)-\lambda^{*}(w^{\prime})\|\leq L_{\lambda}\|w-w^{\prime}\|.$ (35) $$

As specified in (5), the norm of $\lambda^{*}(w)$ is bounded above by

$$ $\|\lambda^{*}(w)\|=\|\mathrm{Proj}_{\mathbb{R}^{p}_{+}}\left(z+\gamma_{2}g(x,y )\right)\|\leq\|z+\gamma_{2}g(x,y)\|\leq r+\gamma_{2}M_{g},\qquad\forall\,w\in X \times Y\times Z,$ $$

where $M_{g}:=\max_{x\in X,y\in Y}\|g(x,y)\|$.
Drawing upon Lemma 2.2, Assumptions 4.2 and 4.3, and (35), for any $w,w^{\prime}\in X\times Y\times Z$, we have

$$ $\displaystyle\|\nabla_{x}\mathcal{G}_{\gamma}(w)-\nabla_{x}\mathcal{G}_{\gamma }(w^{\prime})\|$ $\displaystyle\leq$ $\displaystyle\|\nabla_{x}f(x,y)-\nabla_{x}f(x^{\prime},y^{\prime})\|+\|\lambda ^{*}(w)^{\text{{T}}}\nabla_{x}g(x,y)-\lambda^{*}(w^{\prime})^{\text{{T}}} \nabla_{x}g(x^{\prime},y^{\prime})\|$ $\displaystyle+\|\nabla_{x}f(x,\theta^{*}(w))-\nabla_{x}f(x^{\prime},\theta^{*} (w^{\prime}))\|+\|z^{\text{{T}}}\nabla_{x}g(x,\theta^{*}(w))-(z^{\prime})^{ \text{{T}}}\nabla_{x}g(x^{\prime},\theta^{*}(w^{\prime}))\|$ $\displaystyle\leq$ $\displaystyle L_{f}\|(x,y)-(x^{\prime},y^{\prime})\|+L_{f}\|(x,\theta^{*}(w))- (x^{\prime},\theta^{*}(w^{\prime}))\|$ $\displaystyle+\|\lambda^{*}(w)^{\text{{T}}}\nabla_{x}g(x,y)-\lambda^{*}(w)^{ \text{{T}}}\nabla_{x}g(x^{\prime},y^{\prime})\|+\|\lambda^{*}(w)^{\text{{T}}} \nabla_{x}g(x^{\prime},y^{\prime})-\lambda^{*}(w^{\prime})^{\text{{T}}}\nabla_ {x}g(x^{\prime},y^{\prime})\|$ $\displaystyle+\|z^{\text{{T}}}\nabla_{x}g(x,\theta^{*}(w))-z^{\text{{T}}} \nabla_{x}g(x^{\prime},\theta^{*}(w^{\prime}))\|+\|z^{\text{{T}}}\nabla_{x}g(x ^{\prime},\theta^{*}(w^{\prime}))-(z^{\prime})^{\text{{T}}}\nabla_{x}g(x^{ \prime},\theta^{*}(w^{\prime}))\|$ $\displaystyle\leq$ $\displaystyle L_{f}\|(x,y)-(x^{\prime},y^{\prime})\|+L_{f}\|x-x^{\prime}\|+L_{ f}L_{\theta}\|w-w^{\prime}\|$ $\displaystyle+(r+\gamma_{2}M_{g})L_{g_{1}}\|(x,y)-(x^{\prime},y^{\prime})\|+L_ {g}L_{\lambda}\|w-w^{\prime}\|$ $\displaystyle+rL_{g_{1}}\|x-x^{\prime}\|+rL_{g_{1}}L_{\theta}\|w-w^{\prime}\|+ L_{g}\|z-z^{\prime}\|,$ $$

where the last inequality follows from $\theta^{*}(w^{\prime})\in Y$, $z\in Z$ and $\max_{x\in X,y\in Y}\|\nabla_{x}g(x,y)\|\leq L_{g}$. Similarly, for gradients with respect to $y$ and $z$, for any $w,w^{\prime}\in X\times Y\times Z$, we have

$$ $\displaystyle\|\nabla_{y}\mathcal{G}_{\gamma}(w)-\nabla_{y}\mathcal{G}_{\gamma }(w^{\prime})\|$ $\displaystyle\leq$ $\displaystyle\|\nabla_{y}f(x,y)-\nabla_{y}f(x^{\prime},y^{\prime})\|+\|\lambda ^{*}(w)^{\text{{T}}}\nabla_{y}g(x,y)-\lambda^{*}(w)^{\text{{T}}}\nabla_{y}g(x^ {\prime},y^{\prime})\|$ $\displaystyle+\|\lambda^{*}(w)^{\text{{T}}}\nabla_{y}g(x^{\prime},y^{\prime})- \lambda^{*}(w^{\prime})^{\text{{T}}}\nabla_{y}g(x^{\prime},y^{\prime})\|+\frac {1}{\gamma_{1}}\|y-y^{\prime}\|+\frac{1}{\gamma_{1}}\|\theta^{*}(w)-\theta^{*} (w^{\prime})\|$ $\displaystyle\leq$ $\displaystyle L_{f}\|(x,y)-(x^{\prime},y^{\prime})\|+(r+\gamma_{2}M_{g})L_{g_{ 2}}\|(x,y)-(x^{\prime},y^{\prime})\|$ $\displaystyle+L_{g}L_{\lambda}\|w-w^{\prime}\|+\frac{1}{\gamma_{1}}\|y-y^{ \prime}\|+\frac{1}{\gamma_{1}}L_{\theta}\|w-w^{\prime}\|,$ $$

where the last inequality follows from the fact that $\max_{x\in X,y\in Y}\|\nabla_{y}g(x,y)\|\leq L_{g}$, and

$$ $\displaystyle\|\nabla_{z}\mathcal{G}_{\gamma}(w)-\nabla_{z}\mathcal{G}_{\gamma }(w^{\prime})\|$ $\displaystyle\leq$ $\displaystyle\frac{1}{\gamma_{2}}\|z-z^{\prime}\|+\frac{1}{\gamma_{2}}\| \lambda^{*}(w)-\lambda^{*}(w^{\prime})\|+\|g(x,\theta^{*}(w))-g(x^{\prime}, \theta^{*}(w^{\prime}))\|$ $\displaystyle\leq$ $\displaystyle\frac{1}{\gamma_{2}}\|z-z^{\prime}\|+\frac{L_{\lambda}}{\gamma_{2 }}\|w-w^{\prime}\|+L_{g}\|x-x^{\prime}\|+L_{g}L_{\theta}\|w-w^{\prime}\|.$ $$

The above inequalities together yields the existence of constant $L_{\mathcal{G}}>0$ such that

$$ $\|\nabla\mathcal{G}_{\gamma}(x,y,z)-\nabla\mathcal{G}_{\gamma}(x^{\prime},y^{ \prime},z^{\prime})\|\leq L_{\mathcal{G}}\|(x,y,z)-(x^{\prime},y^{\prime},z^{ \prime})\|.$ $$

Then the conclusion follows from Lemma 5.7 of [].
∎

The update of $\theta^{k}$ constitutes a single gradient descent step to the minimization problem defined in (5). The progress of this update is quantified in the subsequent lemma.

###### Lemma A.3 .

Under Assumptions 4.2 and 4.3, let $\gamma_{1}>0$, $\gamma_{2}>0$ and $\eta_{k}\in(0,1/(L_{f}+rL_{g_{2}}+1/\gamma_{1}))$, then the sequence of $(x^{k},y^{k},z^{k},\theta^{k},\lambda^{k})$ generated by Algorithm 1 satisfies

$$ $\|\theta^{k+1}-\theta^{*}(x^{k},y^{k},z^{k})\|^{2}\leq\,(1-\eta_{k}/\gamma_{1} )\|\theta^{k}-\theta^{*}(x^{k},y^{k},z^{k})\|^{2},$ (36) $$

and

$$ $\lambda^{k+1}=\lambda^{*}(x^{k},y^{k},z^{k}).$ $$

###### Proof.

Consider $(x^{k},y^{k},z^{k})\in X\times Y\times Z$. For brevity, we denote $\theta^{*}(x^{k},y^{k},z^{k})$ and $\lambda^{*}(x^{k},y^{k},z^{k})$ by $\theta^{*}$ and $\lambda^{*}$, respectively. Under Assumptions 4.2 and 4.3, and given $\gamma_{1}>0$, the function

$$ $f(x^{k},\theta)+(z^{k})^{\text{{T}}}g(x^{k},\theta)+\frac{1}{2\gamma_{1}}\| \theta-y^{k}\|^{2},$ $$

is $\frac{1}{\gamma_{1}}$-strongly convex and $(L_{f}+rL_{g_{2}}+1/\gamma_{1})$-smooth with respect to $\theta$. invoking Theorem 10.29 of [], we obtain

$$ $\|\theta^{k+1}-\theta^{*}\|^{2}\leq\left(1-\eta_{k}/\gamma_{1}\right)\|\theta^ {k}-\theta^{*}\|^{2}.$ $$

Additionally, according to (5), it follows that $\lambda^{k+1}=\lambda^{*}$.
∎

The update of variables $(x,y,z)$ in (14) can be viewed as an inexact alternating proximal gradient step from $(x^{k},y^{k},z^{k})$ on solving $\min_{(x,y,z)\in X\times Y\times Z}\phi_{c_{k}}(x,y,z)$.
In the lemma below, we demonstrate that the function $\phi_{c_{k}}(x,y,z)$ exhibits a decreasing property with bounded errors at each iteration.

###### Lemma A.4 .

Under Assumptions 4.2 and 4.3, assume $X$, $Y$ are compact sets, and let $\gamma_{1}>0$, $\gamma_{2}>0$. Then the sequence of $(x^{k},y^{k},z^{k},\theta^{k},\lambda^{k})$ generated by Algorithm 1 satisfies

$$ $\displaystyle\phi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})\leq$ $\displaystyle\phi_{c_{k}}(x^{k},y^{k},z^{k})-\left(\frac{1}{2\alpha_{k}}-\frac {L_{\phi_{k}}}{2}\right)\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|^{2}$ (37) $\displaystyle+\frac{\alpha_{k}}{2}\left(L_{f}+u_{z}L_{g_{1}}+\frac{1}{\gamma_{ 1}}+L_{g}\right)\left\|\theta^{*}(x^{k},y^{k},z^{k})-\theta^{k+1}\right\|^{2}.$ $$

where $L_{\phi_{k}}:=L_{F}/c_{k}+L_{\mathcal{G}}$.

###### Proof.

For clarity, we denote $(x^{k},y^{k},z^{k}),(x^{k+1},y^{k+1},z^{k+1})$ as $w^{k}$ and $w^{k+1}$, respectively.
Under the Assumptions 4.2 and 4.3, where $\nabla F$ and $\nabla f$ are $L_{F}$- and $L_{f}$-Lipschitz continuous on $X\times Y$, and applying Lemma 5.7 of [] and Lemma A.2, we obtain

$$ $\displaystyle\phi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})\leq$ $\displaystyle\phi_{c_{k}}(x^{k},y^{k},z^{k})+\langle\nabla\phi_{c_{k}}(x^{k},y ^{k},z^{k}),(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\rangle$ (38) $\displaystyle+\frac{L_{\phi_{k}}}{2}\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z ^{k})\|^{2},$ $$

with $L_{\phi_{k}}:=L_{F}/c_{k}+L_{\mathcal{G}}$. Based on the update rule of variables $(x,y,z)$ in (14), the convexity of $X\times Y\times Z$ and the property of the projection operator $\mathrm{Proj}_{X\times Y\times Z}$, we have

$$ $\big{\langle}(x^{k},y^{k},z^{k})-\alpha_{k}(d_{x}^{k},d_{y}^{k},d_{z}^{k})-(x^ {k+1},y^{k+1},z^{k+1}),(x^{k},y^{k},z^{k})-(x^{k+1},y^{k+1},z^{k+1})\big{ \rangle}\leq 0,$ $$

yielding

$$ $\big{\langle}(d_{x}^{k},d_{y}^{k},d_{z}^{k}),(x^{k+1},y^{k+1},z^{k+1})-(x^{k}, y^{k},z^{k})\big{\rangle}\leq-\frac{1}{\alpha_{k}}\|(x^{k+1},y^{k+1},z^{k+1})- (x^{k},y^{k},z^{k})\|^{2}.$ $$

Integrating this inequality with (38), we infer that

$$ $\displaystyle\phi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})\leq$ $\displaystyle\phi_{c_{k}}(x^{k},y^{k},z^{k})-\left(\frac{1}{\alpha_{k}}-\frac{ L_{\phi_{k}}}{2}\right)\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|^{2}$ (39) $\displaystyle+\left\langle\nabla\phi_{c_{k}}(x^{k},y^{k},z^{k})-(d_{x}^{k},d_{ y}^{k},d_{z}^{k}),(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\right\rangle.$ $$

Furthermore, considering $w^{k}\in X\times Y\times Z$ and with $\theta^{*}(w^{k}),\theta^{k}\in Y$ for all $k$, and utilizing the formula of $\nabla\mathcal{G}_{\gamma}({x},{y},z)$ derived in Lemma 2.2, along with the definitions of $d_{x}^{k}$, $d_{y}^{k}$ and $d_{z}^{k}$ in (13) and $\lambda^{k+1}=\lambda^{*}(w^{k})$ from Lemma A.3, we can obtain from Assumptions 4.2 and 4.3 that

$$ $\displaystyle\left\|\nabla\phi_{c_{k}}(x^{k},y^{k},z^{k})-(d_{x}^{k},d_{y}^{k} ,d_{z}^{k})\right\|$ (40) $\displaystyle\leq$ $\displaystyle\left\|\nabla_{x}\mathcal{G}_{\gamma}(x^{k},y^{k},z^{k})-\nabla_{ x}f(x^{k},y^{k})-(\lambda^{k+1})^{\text{{T}}}\nabla_{x}g(x^{k},y^{k})+\nabla_{ x}f(x^{k},\theta^{k+1})+(z^{k})^{\text{{T}}}\nabla_{x}g(x^{k},\theta^{k+1})\right\|$ $\displaystyle+\left\|\nabla_{y}\mathcal{G}_{\gamma}(x^{k},y^{k},z^{k})-\nabla_ {y}f(x^{k},y^{k})-(\lambda^{k+1})^{\text{{T}}}\nabla_{y}g(x^{k},y^{k})+(y^{k}- \theta^{k+1})/\gamma_{1}\right\|$ $\displaystyle+\left\|\nabla_{z}\mathcal{G}_{\gamma}(x^{k},y^{k},z^{k})+(z^{k}- \lambda^{k+1})/\gamma_{2}+g(x^{k},\theta^{k+1})\right\|$ $\displaystyle\leq$ $\displaystyle\left\|\lambda^{*}(w^{k})^{\text{{T}}}\nabla_{x}g(x^{k},y^{k})-( \lambda^{k+1})^{\text{{T}}}\nabla_{x}g(x^{k},y^{k})\right\|+\left\|\nabla_{x}f (x^{k},\theta^{*}(w^{k}))-\nabla_{x}f(x^{k},\theta^{k+1})\right\|$ $\displaystyle+\left\|(z^{k})^{\text{{T}}}\nabla_{x}g(x^{k},\theta^{*}(w^{k}))- (z^{k})^{\text{{T}}}\nabla_{x}g(x^{k},\theta^{k+1})\right\|$ $\displaystyle+\left\|\lambda^{*}(w^{k})^{\text{{T}}}\nabla_{y}g(x^{k},y^{k})-( \lambda^{k+1})^{\text{{T}}}\nabla_{y}g(x^{k},y^{k})\right\|+\frac{1}{\gamma_{1 }}\left\|\theta^{*}(w^{k})-\theta^{k+1}\right\|$ $\displaystyle+\frac{1}{\gamma_{2}}\left\|\lambda^{*}(w^{k})-\lambda^{k+1} \right\|+\left\|g(x^{k},\theta^{*}(w^{k}))-g(x^{k},\theta^{k+1})\right\|$ $\displaystyle\leq$ $\displaystyle L_{f}\left\|\theta^{*}(w^{k})-\theta^{k+1}\right\|+rL_{g_{1}} \left\|\theta^{*}(w^{k})-\theta^{k+1}\right\|+\frac{1}{\gamma_{1}}\left\| \theta^{*}(w^{k})-\theta^{k+1}\right\|+L_{g}\left\|\theta^{*}(w^{k})-\theta^{k +1}\right\|$ $\displaystyle=$ $\displaystyle\left(L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+L_{g}\right)\left\| \theta^{*}(w^{k})-\theta^{k+1}\right\|.$ $$

This yields that

$$ $\displaystyle\left\langle\nabla\phi_{c_{k}}(x^{k},y^{k},z^{k})-(d_{x}^{k},d_{y }^{k},d_{z}^{k}),(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\right\rangle$ $\displaystyle\leq$ $\displaystyle\frac{\alpha_{k}}{2}\left(L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+L _{g}\right)\left\|\theta^{*}(w^{k})-\theta^{k+1}\right\|^{2}+\frac{1}{2\alpha_ {k}}\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|^{2},$ $$

Combining this with (39) leads to the following inequality

$$ $\displaystyle\phi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})\leq$ $\displaystyle\phi_{c_{k}}(x^{k},y^{k},z^{k})-\left(\frac{1}{2\alpha_{k}}-\frac {L_{\phi_{k}}}{2}\right)\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|^{2}$ (41) $\displaystyle+\frac{\alpha_{k}}{2}\left(L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+ L_{g}\right)\left\|\theta^{*}(w^{k})-\theta^{k+1}\right\|^{2}.$ $$

∎

###### Lemma A.1 .

###### Proof.

###### Lemma A.2 .

###### Proof.

###### Lemma A.3 .

###### Proof.

###### Lemma A.4 .

###### Proof.

#### A.5.2 Sufficient descent property of V k subscript 𝑉 𝑘 V_{k} italic_V start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT

Utilizing the auxiliary lemmas established previously, we now proceed to demonstrate the sufficient decreasing property of $V_{k}$.

###### Lemma A.5 .

Under Assumptions 4.1, 4.2 and 4.3, suppose $X$, $Y$ are compact sets, $\gamma_{1}>0$, $\gamma_{2}>0$ and $\eta_{k}\in(\underline{\eta},1/(L_{f}+rL_{g_{2}}+1/\gamma_{1}))$ with $\underline{\eta}>0$. Then there exists $c_{\alpha}>0$ such that when $0<\alpha_{k}\leq{c_{\alpha}}$, the sequence of $(x^{k},y^{k},z^{k},\theta^{k},\lambda^{k})$ generated by Algorithm 1 satisfies

$$ $\displaystyle V_{k+1}-V_{k}\leq$ $\displaystyle-\frac{1}{4\alpha_{k}}\left\|w^{k+1}-w^{k}\right\|^{2}-\frac{\eta _{k}C_{\theta}}{2\gamma_{1}}\|\theta^{k}-\theta^{*}(w^{k})\|^{2}.$ (42) $$

###### Proof.

For clarity and conciseness, we represent $(x^{k},y^{k},z^{k}),(x^{k+1},y^{k+1},z^{k+1})$ as $w^{k}$ and $w^{k+1}$, respectively.
Recall (37) from Lemma A.4 that

$$ $\displaystyle\phi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})\leq$ $\displaystyle\phi_{c_{k}}(x^{k},y^{k},z^{k})-\left(\frac{1}{2\alpha_{k}}-\frac {L_{\phi_{k}}}{2}\right)\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|^{2}$ (43) $\displaystyle+\frac{\alpha_{k}}{2}\left(L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+ L_{g}\right)\left\|\theta^{*}(x^{k},y^{k},z^{k})-\theta^{k+1}\right\|^{2}.$ $$

Given that $c_{k+1}\geq c_{k}$,
it follows that $(F(x^{k+1},y^{k+1})-\underline{F})/c_{k+1}\leq(F(x^{k+1},y^{k+1})-\underline{F
})/c_{k}$.
Combining this with (43) leads to

$$ $\displaystyle V_{k+1}-V_{k}=$ $\displaystyle\phi_{c_{k+1}}(w^{k+1})-\phi_{c_{k}}(w^{k})+C_{\theta}\left\| \theta^{k+1}-\theta^{*}(w^{k+1})\right\|^{2}-C_{\theta}\left\|\theta^{k}- \theta^{*}(w^{k})\right\|^{2}$ (44) $\displaystyle\leq$ $\displaystyle\phi_{c_{k}}(w^{k+1})-\phi_{c_{k}}(w^{k})+C_{\theta}\left\|\theta ^{k+1}-\theta^{*}(w^{k+1})\right\|^{2}-C_{\theta}\left\|\theta^{k}-\theta^{*}( w^{k})\right\|^{2}$ $\displaystyle\leq$ $\displaystyle-\left(\frac{1}{2\alpha_{k}}-\frac{L_{\phi_{k}}}{2}\right)\|w^{k+ 1}-w^{k}\|^{2}+C_{\theta}\left\|\theta^{k+1}-\theta^{*}(w^{k+1})\right\|^{2}-C _{\theta}\left\|\theta^{k}-\theta^{*}(w^{k})\right\|^{2}$ $\displaystyle+\frac{\alpha_{k}}{2}\left(L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+ L_{g}\right)\left\|\theta^{*}(w^{k})-\theta^{k+1}\right\|^{2}$ $\displaystyle=$ $\displaystyle-\left(\frac{1}{2\alpha_{k}}-\frac{L_{\phi_{k}}}{2}\right)\|w^{k+ 1}-w^{k}\|^{2}+C_{\theta}\left\|\theta^{k+1}-\theta^{*}(w^{k+1})\right\|^{2}-C _{\theta}\left\|\theta^{k}-\theta^{*}(w^{k})\right\|^{2}$ $\displaystyle+\frac{\alpha_{k}}{2}C_{\theta}\left\|\theta^{*}(w^{k})-\theta^{k +1}\right\|^{2}$ $$

where the last equation follows from defining $C_{\theta}:=L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+L_{g}$.
Further, we can demonstrate that

$$ $\displaystyle\left\|\theta^{k+1}-\theta^{*}(w^{k+1})\right\|^{2}-\left\|\theta ^{k}-\theta^{*}(w^{k})\right\|^{2}+\frac{\alpha_{k}}{2}\left\|\theta^{*}(w^{k} )-\theta^{k+1}\right\|^{2}$ $\displaystyle\leq$ $\displaystyle(1+\epsilon_{k}+\alpha_{k}/2)\left\|\theta^{k+1}-\theta^{*}(w^{k} )\right\|^{2}-\left\|\theta^{k}-\theta^{*}(w^{k})\right\|^{2}+(1+\frac{1}{ \epsilon_{k}})\|\theta^{*}(w^{k+1})-\theta^{*}(w^{k})\|^{2}$ $\displaystyle\leq$ $\displaystyle(1+\epsilon_{k}+\alpha_{k}/2)(1-\eta_{k}/\gamma_{1})\|\theta^{k}- \theta^{*}(w^{k})\|^{2}-\left\|\theta^{k}-\theta^{*}(w^{k})\right\|^{2}+(1+ \frac{1}{\epsilon_{k}})L_{\theta}^{2}\left\|w^{k+1}-w^{k}\right\|^{2},$ $$

for any $\epsilon_{k}>0$, where the second inequality is a consequence of Lemmas A.1 and A.3.
By setting $\epsilon_{k}=\frac{1}{4}\eta_{k}/\gamma_{1}$ in the above inequality, we obtain that when $\alpha_{k}\leq\frac{1}{2}\eta_{k}/\gamma_{1}$, the following inequalities hold

$$ $(1+\epsilon_{k}+\alpha_{k}/2)(1-\eta_{k}/\gamma_{1})\leq(1+\frac{1}{2}\eta_{k} /\gamma_{1})(1-\eta_{k}/\gamma_{1})\leq 1-\frac{\eta_{k}}{2\gamma_{1}}.$ $$

Consequently, we establish the inequality

$$ $\displaystyle\left\|\theta^{k+1}-\theta^{*}(w^{k+1})\right\|^{2}-\left\|\theta ^{k}-\theta^{*}(w^{k})\right\|^{2}+\frac{\alpha_{k}}{2}\left\|\theta^{*}(w^{k} )-\theta^{k+1}\right\|^{2}$ (45) $\displaystyle\leq$ $\displaystyle-\frac{\eta_{k}}{2\gamma_{1}}\|\theta^{k}-\theta^{*}(w^{k})\|^{2} +\left(1+\frac{4\gamma_{1}}{\eta_{k}}\right)L_{\theta}^{2}\left\|w^{k+1}-w^{k} \right\|^{2}.$ $$

Combining (44) and (45) implies

$$ $\displaystyle V_{k+1}-V_{k}\leq$ $\displaystyle-\left[\frac{1}{2\alpha_{k}}-\frac{L_{\phi_{k}}}{2}-\left(1+\frac {4\gamma_{1}}{\eta_{k}}\right)L_{\theta}^{2}C_{\theta}\right]\left\|w^{k+1}-w^ {k}\right\|^{2}-\frac{\eta_{k}C_{\theta}}{2\gamma_{1}}\|\theta^{k}-\theta^{*}( w^{k})\|^{2}.$ (46) $$

When $c_{k+1}\geq c_{k}$, $\eta_{k}\geq\underline{\eta}>0$, $\alpha_{k}\leq\frac{1}{2}\underline{\eta}/\gamma_{1}$, it holds that for any $k$, $\alpha_{k}\leq\frac{1}{2}\eta_{k}/\gamma_{1}$,

$$ $\displaystyle\frac{L_{\phi_{k}}}{2}+\left(1+\frac{4\gamma_{1}}{\eta_{k}}\right )L_{\theta}^{2}C_{\theta}\leq\frac{L_{\phi_{0}}}{2}+\left(1+\frac{4\gamma_{1}} {\underline{\eta}}\right)L_{\theta}^{2}C_{\theta}=:C_{\alpha}.$ $$

If $c_{\alpha}>0$ satisfies

$$ $c_{\alpha}\leq\min\left\{\frac{1}{2}\underline{\eta}/\gamma_{1},\frac{1}{4C_{ \alpha}}\right\},$ $$

then, for $0<\alpha_{k}\leq c_{\alpha}$, it is guaranteed that

$$ $\frac{1}{2\alpha_{k}}-\frac{L_{\phi_{k}}}{2}-\left(1+\frac{4\gamma_{1}}{\eta_{ k}}\right)L_{\theta}^{2}C_{\theta}\geq\frac{1}{4\alpha_{k}}.$ $$

Therefore, the conclusion follows directly from (46).
∎

###### Lemma A.5 .

###### Proof.

#### A.5.3 Proofs of Theorem 4.4 and Theorem 4.5

Indeed, Theorem 4.4 is a special case of Theorem 4.5 with $\rho=0$. Consequently, we present the proof of Theorem 4.5.

###### Theorem 4.5 .

Under Assumptions 4.1, 4.2 and 4.3, assume that $X$, $Y$ are compact sets, $\gamma_{1}>0$, $\gamma_{2}>0$, $c_{k}=c(k+1)^{\rho}$ with $c>0$, $\rho\in[0,1/2)$ and $\eta_{k}\in(\underline{\eta},1/(L_{f}+rL_{g_{2}}+1/\gamma_{1}))$ with $\underline{\eta}>0$. Then there exists $c_{\alpha}>0$ such that when $\alpha_{k}\in(\underline{\alpha},{c_{\alpha}})$ with $\underline{\alpha}>0$, the sequence of $(x^{k},y^{k},z^{k},\theta^{k},\lambda^{k})$ generated by Algorithm 1 satisfies

$$ $\min_{0\leq k\leq K}\left\|\theta^{k}-\theta^{*}(x^{k},y^{k},z^{k})\right\|=O \left(\frac{1}{K^{1/2}}\right),$ $$

and

$$ $\min_{0\leq k\leq K}R_{k}(x^{k+1},y^{k+1},z^{k+1})=O\left(\frac{1}{K^{(1-2\rho )/2}}\right).$ $$

Furthermore, if $\rho>0$ and $\psi_{c_{k}}(x^{k},y^{k},z^{k})$ is uniformly bounded above, then the sequence of $(x^{k},y^{k},z^{k})$ satisfies

$$ $0\leq\mathcal{G}_{\gamma}(x^{K},y^{K},z^{K})=O\left(\frac{1}{K^{p}}\right).$ $$

###### Proof.

Lemma A.5 establishes the existence of $c_{\alpha}>0$ such that (42) holds when $\alpha_{k}\leq c_{\alpha}$.
Summing (42) over $k=0,1,\ldots,K-1$, we obtain

$$ $\displaystyle\sum_{k=0}^{K-1}\Bigg{(}\frac{1}{4\alpha_{k}}\|(x^{k+1},y^{k+1},z ^{k+1})-(x^{k},y^{k},z^{k})\|^{2}+\frac{\underline{\eta}C_{\theta}}{2\gamma_{1 }}\|\theta^{k}-\theta^{*}(x^{k},y^{k},z^{k})\|^{2}\Bigg{)}$ (47) $\displaystyle\leq$ $\displaystyle V_{0}-V_{K}\leq V_{0},$ $$

where the last inequality follows from the nonnegativity of $V_{K}$. Consequently, it holds that

$$ $\sum_{k=0}^{\infty}\|\theta^{k}-\theta^{*}(x^{k},y^{k},z^{k})\|^{2}<\infty,$ $$

and

$$ $\min_{0\leq k\leq K}\left\|\theta^{k}-\theta^{*}(x^{k},y^{k},z^{k})\right\|=O \left(\frac{1}{K^{1/2}}\right).$ $$

According to the update rule for variables $(x,y,z)$ in (14), we derive

$$ $\displaystyle 0$ $\displaystyle\in c_{k}(d_{x}^{k},d_{y}^{k},d_{z}^{k})+\mathcal{N}_{X\times Y \times Z}(x^{k+1},y^{k+1},z^{k+1})+\frac{c_{k}}{\alpha_{k}}\left((x^{k+1},y^{k +1},z^{k+1})-(x^{k},y^{k},z^{k})\right).$ (48) $$

Following the definitions of $d_{x}^{k}$, $d_{y}^{k}$ and $d_{z}^{k}$ in (13), we obtain

$$ $e^{k}\in\nabla\psi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})+\mathcal{N}_{X\times Y \times Z}(x^{k+1},y^{k+1},z^{k+1}),$ $$

where

$$ $e^{k}:=\nabla\psi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})-c_{k}(d_{x}^{k},d_{y}^{k},d _{z}^{k})-\frac{c_{k}}{\alpha_{k}}\left((x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k} ,z^{k})\right).$ (49) $$

Next, we estimate $\|e^{k}\|$. We have

$$ $\displaystyle\|e^{k}\|\leq$ $\displaystyle\|\nabla\psi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})-\nabla\psi_{c_{k}}( x^{k},y^{k},z^{k})\|+\|\nabla\psi_{c_{k}}(x^{k},y^{k},z^{k})-c_{k}(d_{x}^{k},d _{y}^{k},d_{z}^{k})\|$ $\displaystyle+\frac{c_{k}}{\alpha_{k}}\left\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k}, y^{k},z^{k})\right\|.$ $$

For the first term in the right hand side of the above inequality, by using Assumptions 4.1, 4.2 and 4.3, along with Lemma A.1, we establish the existence of $L_{\psi}:=L_{F}+L_{\mathcal{G}}>0$ such that

$$ $\|\nabla\psi_{c_{k}}(x^{k+1},y^{k+1},z^{k+1})-\nabla\psi_{c_{k}}(x^{k},y^{k},z ^{k})\|\leq c_{k}L_{\psi}\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|.$ $$

Using (40) and Lemma A.3, we have

$$ $\displaystyle\|\nabla\psi_{c_{k}}(x^{k},y^{k},z^{k})-c_{k}(d_{x}^{k},d_{y}^{k} ,d_{z}^{k})\|=$ $\displaystyle c_{k}\left\|\nabla\phi_{c_{k}}(x^{k},y^{k},z^{k})-(d_{x}^{k},d_{ y}^{k},d_{z}^{k})\right\|$ (50) $\displaystyle\leq$ $\displaystyle c_{k}C_{\theta}\left\|\theta^{*}(x^{k},y^{k},z^{k})-\theta^{k+1}\right\|$ $\displaystyle\leq$ $\displaystyle c_{k}C_{\theta}\left\|\theta^{*}(x^{k},y^{k},z^{k})-\theta^{k} \right\|.$ $$

with $C_{\theta}=L_{f}+rL_{g_{1}}+\frac{1}{\gamma_{1}}+L_{g}$.
Consequently, we have

$$ $\displaystyle\|e^{k}\|\leq$ $\displaystyle c_{k}L_{\psi}\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\|+ \frac{c_{k}}{\alpha_{k}}\left\|(x^{k+1},y^{k+1},z^{k+1})-(x^{k},y^{k},z^{k})\right\|$ $\displaystyle+c_{k}C_{\theta}\left\|\theta^{*}(x^{k},y^{k},z^{k})-\theta^{k} \right\|.$ $$

Using this bound on $\|e^{k}\|$, we have that

$$ $\displaystyle R_{k}(x^{k+1},y^{k+1},z^{k+1})\leq$ $\displaystyle c_{k}\left(L_{\psi}+1/\alpha_{k}\right)\|(x^{k+1},y^{k+1},z^{k+1 })-(x^{k},y^{k},z^{k})\|$ $\displaystyle+c_{k}C_{\theta}\left\|\theta^{*}(x^{k},y^{k},z^{k})-\theta^{k} \right\|.$ $$

Utilizing this inequality, and given that $\alpha_{k}\in(\underline{\alpha},{c_{\alpha}})$ for some positive constant $\underline{\alpha}$, we can establish the existence of a constant $C_{R}>0$ such that

$$ $\displaystyle\frac{1}{c_{k}^{2}}R_{k}(x^{k+1},y^{k+1},z^{k+1})^{2}$ (51) $\displaystyle\leq$ $\displaystyle C_{R}\Bigg{(}\frac{1}{4\alpha_{k}}\|(x^{k+1},y^{k+1},z^{k+1})-(x ^{k},y^{k},z^{k})\|^{2}+\frac{\underline{\eta}C_{\theta}}{2\gamma_{1}}\|\theta ^{k}-\theta^{*}(x^{k},y^{k},z^{k})\|^{2}\Bigg{)}.$ $$

This inequality, combined with (47), implies that

$$ $\sum_{k=0}^{\infty}\frac{1}{c_{k}^{2}}R_{k}(x^{k+1},y^{k+1},z^{k+1})^{2}<\infty.$ (52) $$

Because $2\rho<1$, it follows that

$$ $\sum_{k=0}^{K}\frac{1}{c_{k}^{2}}=\sum_{k=0}^{K}\left(\frac{1}{k+1}\right)^{2 \rho}\frac{1}{c^{2}}\geq\left(\int_{1}^{K+2}\frac{1}{x^{{}^{2\rho}}}d\,x\right )\frac{1}{c^{2}}=\left(\frac{(K+2)^{1-2\rho}-1}{1-2\rho}\right)\frac{1}{c^{2}},$ $$

leading us to conclude from (52) that

$$ $\min_{0\leq k\leq K}R_{k}(x^{k+1},y^{k+1},z^{k+1})=O\left(\frac{1}{K^{(1-2\rho )/2}}\right).$ $$

Given that $\psi_{c_{k}}(x^{k},y^{k},z^{k})\leq M$ and $F(x^{k},y^{k})\geq\underline{F}$ for all $k$, it follows that

$$ $c_{k}\mathcal{G}_{\gamma}(x^{k},y^{k},z^{k})\leq M-\underline{F},\quad\forall k.$ $$

From $c_{k}={c}(k+1)^{\rho}$ and $\rho>0$, we can obtain that

$$ $\mathcal{G}_{\gamma}(x^{K},y^{K},z^{K})=O\left(\frac{1}{K^{\rho}}\right).$ $$

∎

###### Theorem 4.5 .

###### Proof.

### A.6 Extension to bilevel optimization with minimax lower-level problem

In this part, we explore the extension of our proposed gradient-based, single-loop, Hessian-free algorithm, originally designed for bilevel optimization problems with constrained lower-level problems, to bilevel optimization problems with minimax lower-level problem,

$$ $\displaystyle\min_{x\in X,y\in Y,z\in Z}F(x,y,z)\quad\mathrm{s.t.}\quad(y,z) \in\mathcal{SP}(x),$ $$

where $\mathcal{SP}(x)$ denotes the set of saddle points for the convex-concave minimax problem,

$$ $\displaystyle\min_{y\in Y}\max_{z\in Z}f(x,y,z),$ $$

where $x\in\mathbb{R}^{n}$, $y\in\mathbb{R}^{m}$ and $z\in\mathbb{R}^{p}$, the sets $X$, $Y$ and $Z$ are closed convex sets in $\mathbb{R}^{n}$, $\mathbb{R}^{m}$ and $\mathbb{R}^{p}$, respectively. The UL objective $F:X\times Y\times Z\rightarrow\mathbb{R}$, and the LL objective $f:X\times Y\times Z\rightarrow\mathbb{R}$ are continuously differentiable with $f$ being convex in $y$ and concave in $z$.
Building upon the idea applied in the development of the regularized gap function (3), we introduce the doubly regularized gap function for lower-level minimax problems, defined as:

$$ $\displaystyle\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z):=\max_{\theta\in Y, \lambda\in Z}\Big{\{}f(x,y,\lambda)-\frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}-f(x ,\theta,z)-\frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\Big{\}}.$ $$

By employing proof techniques analogous to those used in Lemma 2.1 or Theorem 3.3 from [], we can derive similar results for the doubly regularized gap function $\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)$.

###### Lemma A.6 .

Assume that $f(x,y,z)$ is convex in $y$ on $Y$ for any given $x\in X,z\in Z$ and concave in $z$ on $Z$ for any given $x\in X,y\in Y$. Let $\gamma_{1},\gamma_{2}>0$, we have $\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)\geq 0$ for any $(x,y,z)\in X\times Y\times Z$, and

$$ $\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)\leq 0,$ $$

if and only if $(y,z)\in\mathcal{SP}(x)$.

Similarly, utilizing the proof methods in Lemma 2.2, we establish the differentiability of $\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)$ and derive the formula for its gradient.

###### Lemma A.7 .

Assume that $f(x,y,z)$ is convex in $y$ on $Y$ for any given $x\in X,z\in Z$ and concave in $z$ on $Z$ for any given $x\in X,y\in Y$ and is continuously differentiable on an open set containing $X\times Y\times Z$. Then $\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)$ is continuously differentiable on $X\times Y\times Z$, and for any $(x,y,z)\in X\times Y\times Z$,

$$ $\nabla\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)=\begin{pmatrix}\nabla_{x}f (x,y,\lambda^{*})\\ \nabla_{y}f(x,y,\lambda^{*})\\ -\left(z-\lambda^{*}\right)/\gamma_{2}\end{pmatrix}-\begin{pmatrix}\nabla_{x}f (x,\theta^{*},z)\\ \left(y-\theta^{*}\right)/\gamma_{1}\\ \nabla_{z}f(x,\theta^{*},z)\end{pmatrix},$ (53) $$

where $\theta^{*}$ and $\lambda^{*}$ denote $\theta^{*}(x,y,z)$ and $\lambda^{*}(x,y,z)$, respectively, defined as

$$ $\displaystyle\theta^{*}(x,y,z)$ $\displaystyle:=\underset{\theta\in Y}{\mathrm{argmin}}\left\{f(x,\theta,z)+ \frac{1}{2\gamma_{1}}\|\theta-y\|^{2}\right\},$ (54) $\displaystyle\lambda^{*}(x,y,z)$ $\displaystyle:=\underset{\lambda\in Z}{\mathrm{argmax}}\left\{f(x,y,\lambda)- \frac{1}{2\gamma_{2}}\|\lambda-z\|^{2}\right\}.$ $$

This newly introduced gap function enables the following equivalent single-level reformulation of the problem (16),

$$ $\min_{(x,y,z)\in X\times Y\times Z}\ F(x,y,z)\quad\mathrm{s.t.}\quad\mathcal{G }^{\mathrm{saddle}}_{\gamma}(x,y,z)\leq 0.$ $$

Using Lemma A.6 and the proof techniques in Theorem 2.3, we can establish the equivalence between the reformulation (19) and the bilevel optimization problem (16).

###### Theorem A.8 .

Assume that $f(x,y,z)$ is convex in $y$ on $Y$ for any given $x\in X,z\in Z$ and concave in $z$ on $Z$ for any given $x\in X,y\in Y$. Let $\gamma_{1},\gamma_{2}>0$, the reformulation (19) is equivalent to the bilevel optimization problem (16).

Utilizing this reformulation, analogous to BiC-GAFFA, we can propose a gradient-based, single-loop, Hessian-free algorithm for problem (16).
At each iteration, we employ a single projected gradient descent step to update $\theta^{k+1},\lambda^{k+1}$ to approximate $\theta^{*}(x^{k},y^{k},z^{k})$ and $\lambda^{*}(x^{k},y^{k},z^{k})$, as follows:

$$ $\displaystyle\theta^{k+1}$ $\displaystyle=\mathrm{Proj}_{Y}\left(\theta^{k}-\eta_{k}d_{\theta}^{k}\right),$ $\displaystyle\lambda^{k+1}$ $\displaystyle=\mathrm{Proj}_{Z}\left(\lambda^{k}-\eta_{k}d_{\lambda}^{k}\right),$ $$

where $\eta_{k}>0$ is the step size, and

$$ $\displaystyle d_{\theta}^{k}$ $\displaystyle:=\nabla_{y}f(x^{k},\theta^{k},z^{k})+\frac{1}{\gamma_{1}}(\theta ^{k}-y^{k}),$ (55) $\displaystyle d_{\lambda}^{k}$ $\displaystyle:=-\nabla_{z}f(x^{k},y^{k},\lambda^{k})+\frac{1}{\gamma_{2}}( \lambda^{k}-z^{k}).$ $$

By substituting $(\theta^{k+1},\lambda^{k+1})$ for $(\theta^{*},\lambda^{*})$ in (53), we can approximate the gradients of the function

$$ $\frac{1}{c}F(x,y,z)+\mathcal{G}^{\mathrm{saddle}}_{\gamma}(x,y,z)$ $$

to define the update directions:

$$ $\displaystyle d_{x}^{k}:=$ $\displaystyle\frac{1}{c_{k}}\nabla_{x}F(x^{k},y^{k},z^{k})+\nabla_{x}f(x^{k},y ^{k},\lambda^{k+1})-\nabla_{x}f(x^{k},\theta^{k+1},z^{k}),$ (56) $\displaystyle d_{y}^{k}:=$ $\displaystyle\frac{1}{c_{k}}\nabla_{y}F(x^{k},y^{k},z^{k})+\nabla_{y}f(x^{k},y ^{k},\lambda^{k+1})-(y^{k}-\theta^{k+1})/\gamma_{1},$ $\displaystyle d_{z}^{k}:=$ $\displaystyle\frac{1}{c_{k}}\nabla_{y}F(x^{k},y^{k},z^{k})-(z^{k}-\lambda^{k+1 })/\gamma_{2}-\nabla_{z}f(x^{k},\theta^{k+1},z^{k}).$ $$

Finally, we implement an update for the variables $(x,y,z)$ using a step size $\alpha_{k}>0$:

$$ $(x^{k+1},y^{k+1},z^{k+1})=\mathrm{Proj}_{X\times Y\times Z}\left((x^{k},y^{k}, z^{k})-\alpha_{k}(d_{x}^{k},d_{y}^{k},d_{z}^{k})\right).$ $$

The complete algorithm is presented in Algorithm 2.

Figure: Algorithm 2 single-loop Hessian-free algorithm for bilevel optimization problems with minimax lower-level problem

While the primary focus of this paper is the bilevel optimization with constrained lower-level problems, we defer the convergence analysis of this algorithm to future work.

###### Lemma A.6 .

###### Lemma A.7 .

###### Theorem A.8 .

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: In Proposition 3.1, it is assumed that the optimal solution $(x^*,y^*,z^*)$ to (6) exists with finite $z^*$. How can this be guaranteed? Additionally, how should $r$ be chosen in practice? Setting $r$ too large may hinder algorithm convergence, while choosing $r$ too small could exclude the optimal point $z^*$, which is unknown. 

Additionally, several concurrent related works are missing: [1]--[2]. 

[1] A Primal-Dual-Assisted Penalty Approach to Bilevel Optimization with Coupled Constraints. L Jiang, et. al. arXiv:2406.10148. 

[2] First-Order Methods for Linearly Constrained Bilevel Optimization. G Kornowski, et.al. 	arXiv:2406.12771.

In Table 3, the term "Required accuracy" could be replaced with a word like "Time" to clarify the values reported. As it stands, it is unclear whether the numbers in this table represent accuracy or computation time.
soundness: 3
strengths: 1. This work tackles an important bilevel problem which involves coupled constraints at the lower level.
2. The design of the doubly regularized gap function is novel and effective. 
3. The numerical studies are comprehensive and thorough.
summary: This paper studies a bilevel problem where lower-level problem is convex with coupled constraints. To avoid joint projection onto the coupled constraint set which requires coupling convexity and is costly, this paper introduces a doubly regularized gap function to convert lower-level problem into a smooth constraint. This constrained problem is then solved by penalty method. Convergence analysis is provided and numerical results validate the effectiveness of the proposed method.
weaknesses: This paper addresses a relaxation of the original bilevel problem through a penalty reformulation and truncation approach. For the penalty reformulation, it is shown that by selecting a sufficiently large penalty constant $c$, the penalty problem becomes an $\epsilon$-approximate solution to the original problem. At the same time, the truncation parameter $r$ also appears to need a large value to ensure equivalence with the original problem. But in the final convergence theorem, $r$ is chosen as a constant, independent of the target error $\epsilon$.

---

confidence: 5
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. The convergence results for extending the proposed algorithm to bilevel optimization problems with minimax lower-level problems are lacking, which leads to that Section 5 feels somewhat abrupt overall.

2. In the experimental settings, we observe that the step sizes $\alpha_k$ and $\eta_k$ are fixed or vary in a specific form. Are the step sizes related to the specific formulations of different problems?

Some more detailed questions:

1. Can the conclusion of Lemma 2.1 be stated as "$\mathcal{G}_\gamma(x,y,z)=0$ if and only if $y\in S(x)$ and $z\in \mathcal{M}(x, y)$"?

2. In line 845, the loss function in lower-level problem should be $\mathcal{L}_{tr}$.

3. In line 1007, can you provide the specific mathematical format of the functions $\mathcal{L}_{gen/det}$
soundness: 3
strengths: 1.The paper introduces a smooth constrained reformulation for bilevel optimization with a constrained lower-level problem. This reformulation avoids the use of any implicit functions associated with the lower-level problem, which are typically known to necessitate additional iterations or inexact solutions.

2. The authors present the first-order single-loop algorithm BiC-GAFFA, which exhibits overall complexity that is lower than that of many existing first-order algorithms requiring double-loop iterations. Notably, it does not rely on the assumption of full convexity of the lower-level constraints.

3. When updating the proximal variable $\lambda$ associated with the multipliers $z$, there is no need to compute the gradient of the gap function. Instead, the update can be effectively executed solely by applying the first-order optimality conditions of the maximization problem.

4. BiC-GAFFA possesses a broader range of potential applications, making it suitable for solving bilevel optimization problems where the lower-level problem is a minimax problem. Existing bilevel optimization algorithms are unable to achieve this simultaneously.
summary: The paper introduces a method for transforming constrained bilevel optimization problems with coupled lower-level constraints into a single-level problem with a smooth constraint, employing a doubly regularized gap function. Based on this reformulation, the authors propose a first-order single-loop penalty-based algorithm called BiC-GAFFA. Theoretically, they establish non-asymptotic convergence results under general assumptions. Additionally, they investigate the extension of the proposed algorithm to bilevel optimization problems featuring minimax lower-level problems, a topic that has received limited attention in the existing literature. Extensive experiments are conducted to validate the practicality, robustness and efficiency of BiC-GAFFA, encompassing multiple complex problem models within the learning problems.
weaknesses: 1. From Section 6.2 and A.2.2., we observe that when dealing with nonsmooth problems, such as sparse group lasso and SVM, it is necessary to reformulate them into a smooth form for effective resolution. It may impact the applicability of the algorithms.

2. By introducing the gap function, the authors have designed a single-loop algorithm. However, the dimensionality of the introduced variables remains the same as that of the original problem, resulting in an increased iteration scale. Compared to LV-HBA proposed in [1], there is no significant improvement in mathematical representation or convergence results.

3. The necessity to control multiple iteration step sizes within the algorithm significantly affects the convergence theory and empirical performance, making them highly dependent on the adjustment of these step size parameters.

[1] Wei Yao, Chengming Yu, Shangzhi Zeng, and Jin Zhang. Constrained bi-level optimization: Proximal Lagrangian value function approach and Hessian-free algorithm. In ICLR, 2024.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: See weaknesses.
soundness: 3
strengths: 1. The introduction of the doubly regularized gap function is a novel (as far as I am aware in the context of bilevel optimization) solution to deal with constraints at the lower level in bilevel optimization. This regularization allows the transformation of the bilevel problem into a smooth optimization problem, which simplifies the numerical implementation and analysis. 

2. The proposed BiC-GAFFA algorithm avoids costly second-order computations, which is usually the case for methods based on constrained optimization reformulation. The single-loop nature of the algorithm is also computationally efficient, as supported by the non-asymptotic convergence analysis. 

3. The extension to minimax inner level problems provides flexibility in dealing with a broader class of problems beyond the traditional settings.
summary: The paper introduces a novel algorithm for constrained bilevel optimization (BiO) using a regularized gap function to transform bilevel problems with lower-level (LL) constraints into a single-level optimization problem with a smooth inequality constraint. The authors propose a first-order, single-loop algorithm, BiC-GAFFA, that operates without requiring Hessian evaluations or projections onto the LL constraint set, offering both theoretical convergence results and empirical validation on synthetic and real-world tasks, such as hyperparameter optimization and generative adversarial networks (GANs). The paper also presents an extension of the method to minimax lower level problems, which broadens its applicability.
weaknesses: 1. Some theoretical assumptions are insufficiently justified (i.e. very weak statement), particularly in Remark 2.5, where it’s stated that the existence of a multiplier in Theorem 2.3 “can be guaranteed” under certain constraints qualification conditions like MFCQ. Is the assumption in fact guaranteed when the constraints qualification conditions are satisfied? Or does those conditions directly imply that the assumption is satisfied? Please clarify. 

2. The update steps in BiC-GAFFA are quite similar to existing constrained optimization reformulation methods, such as primal-dual [1] and dynamic barrier [2] methods. While these approaches have been applied successfully in other constrained settings, the paper does not explain why they are inadequate for the additional lower-level constraint here, suggesting a need for BiC-GAFFA. A discussion on this distinction would help clarify the necessity of more complex methods such as the proposed BiC-GAFFA. 

3. The algorithm requires additional sequences (e.g.,  $\theta_k$  and  $\lambda_k$), adding to the computational and memory load, which may challenge its feasibility for large-scale applications. 

4. The experimental results focus on toy numerical problems and small-scale real-world problems. While these results are encouraging, a broader evaluation on large-scale or diverse domains would better demonstrate BiC-GAFFA’s practical applicability and scalability. 

Reference: 

[1] Sow et al. A Primal-Dual Approach to Bilevel Optimization with Multiple Inner Minima. 

[2] Ye et al. BOME! Bilevel Optimization Made Easy: A Simple First-Order Approach.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 2
questions: See summary
soundness: 3
strengths: See summary
summary: This work studies bilevel problems where the LL problem is constrained. 

The key contributions are:

1. The equivalence between main problem and the smoothed problem (6) under convexity.

2. The relation between (6) and the sequence of problems in Thm 3.3

3. A gradient based algorithm with nonasymptotic grad-norm rate for an instance of problem 3.3

Strength:

1.The paper is technically sound. 

2.The comparison against known results are clear, and hence the main contributions are clear.

3. I like the discussions on reformulations and the experiments.

Weakness:
I hope the authors could address my concerns below.

1. When and why would constraints in the LL be necessary? 

a. GAN would not require constraints.

b. SVM has uncontrained versions.

c. Hyperparameter tuning / data selection can be solved with more naive approach without bilevel problems.

2. There are more than one way to smoothen/regularize the lower level so that the limit is the original problem.

For example: one can add positive coefficients c1, c2, so that the bilevel is now F +  c1 f + c2 max(g, 0) . Then by sending c1, c2/c1 to infinity, one could also get an approximate solution. Why would one formulation better than another?


Minors: 
1. The notation N_Y(y) and other similar ones are not defined.
weaknesses: See summary

Meta-review: This paper studies bilevel optimization problems with coupled lower-level constraints beyond linear lower-level constraints. This is a challenging setting because the standard bilevel methods do not work. The novelty of this paper is that it introduces a new penalty function and the corresponding gradient update to solve the resultant problems.  This paper addresses an important problem and provide a comprehensive answer. Therefore, I will support for "accept."

Decision: Accept (Poster)