# A Primal-Dual-Assisted Penalty Approach to Bilevel Optimization with Coupled
Constraints

Liuyuan Jiang†, Quan Xiao†, Victor M. Tenorio⋆, Fernando Real-Rojas⋆
Antonio Marques⋆, Tianyi Chen†
†Rensselaer Polytechnic Institute, Troy, NY, United States
⋆King Juan Carlos University, Madrid, Spain
{jiangl7, xiaoq5<EMAIL_ADDRESS>
{victor.tenorio<EMAIL_ADDRESS><EMAIL_ADDRESS>The work was supported by NSF project 2134168, NSF CAREER project 2047177,
Cisco Research Award, Amazon Research Award, and the IBM-Rensselaer Future of
Computing Research Collaboration.

###### Abstract

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

###  1 Introduction

Bilevel optimization (BLO) approaches are pertinent in various machine
learning problems, including hyperparameter optimization [48, 24], meta-
learning [22], and reinforcement learning [65, 62]. Moreover, the ability to
handle BLO with constraints is particularly important, as these constraints
appear in applications such as pricing [16], transportation [49, 1], and
kernelized SVM [30]. Although there is extensive research on BLO problems
without constraints or with uncoupled constraints [31, 12, 14, 38, 41, 61],
solutions for BLO problems with coupled constraints (CCs) remain limited; see
details in Table 1. However, it is of particular interest to investigate BLO
with lower-level CCs. Taking infrastructure planning in a transportation
network as an example, the lower-level problem seeks to optimize the total
utility given a network configuration (determined by an upper-level variable),
and the upper-level problem is to optimize the network configuration itself.

Motivated by this, we consider the coupled-constrained BLO problem in the
following form

$\displaystyle\min_{x\in\mathcal{X}}$ $\displaystyle~{}f(x,y^{*}_{g}(x))$ (1a)
s.t. $\displaystyle\quad
y^{*}_{g}(x):=\arg\min_{y\in\mathcal{Y}(x)}g(x,y)\quad\text{with
}\quad\mathcal{Y}(x):=\\{y\in\mathcal{Y}:g^{c}(x,y)\leq 0\\}$ (1b)

where $f:\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{y}}\rightarrow\mathbb{R}$ and
$g:\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{y}}\rightarrow\mathbb{R}$ are
respectively the upper-level and lower-level objectives,
$g^{c}(x,y):\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{y}}\rightarrow\mathbb{R}^{d_{c}}$
defines the lower-level CCs, and
$\mathcal{X}\subseteq\mathbb{R}^{d_{x}},\mathcal{Y}\subseteq\mathbb{R}^{d_{y}}$
are the domain of $x$ and $y$ that are easy to project, such as the Euclidean
ball.

The difficulty of solving (1) lies in the coupling between the upper and
lower-level problems. To develop efficient methods, previous literature
started from the unconstrained BLO problem and solved it by implicit gradient
descent (IGD) method [27, 31, 34, 12, 13, 39, 60, 42, 64] and penalty-based
methods [44, 43, 60, 40, 41], to name a few. To solve BLO with CCs, AiPOD [68]
and GAM [69] have investigated the IGD method under different constraint
settings. However, AiPOD only considered equality constraints, and GAM lacked
finite-time convergence guarantees. Leveraging a penalty reformulation, [72]
developed a Hessian-free method with finite-time convergence. However, a key
algorithm step in [72] is a joint projection of the current iterate $(x,y)$
onto the constraint set. This projection, which needs to be carried out at
every iteration, can be challenging when the number of variables
($d_{x},d_{y}$) or constraints ($d_{c}$) is large or when $g^{c}(x,y)$ is not
_jointly_ convex.

To this end, this paper aims to address the following question

Can we develop an efficient algorithm that bypasses joint projections on
$g^{c}(x,y)$ and quickly solves the BLO problem with coupled inequality
constraints in (1)?

We address this question affirmatively, focusing on the setting where the
lower-level objective, $g(x,y)$, is strongly convex in $y$, and the
constraints $g^{c}(x,y)$ are convex in $y$. To avoid implementing a joint
projection, we put forth a novel single-level primal-dual-assisted penalty
reformulation that decouples $x$ and $y$. Specifically, with
$\mu\in\mathbb{R}^{d_{c}}$ denoting the Lagrange multiplier of (1b), we
propose solving

$\displaystyle\min_{x\in\mathcal{X}}~{}~{}F_{\gamma}(x)$
$\displaystyle:=\max_{\mu\in\mathbb{R}^{dx}_{+}}\min_{y\in\mathcal{Y}}f(x,y)+\underbrace{\gamma(g(x,y)-v(x))}_{\text{penalty
term}}+\underbrace{\langle\mu,g^{c}(x,y)\rangle}_{\text{Lagrangian term}}$
(2a) where $\displaystyle\quad v(x):=\min_{y\in\mathcal{Y}(x)}g(x,y)$ (2b)

where the penalty term with the constant $\gamma$ controls the distance
between $y$ and $y^{*}_{g}(x)$ by penalizing the distance between the lower-
level objective $g(x,y)$ and its value function $v(x)$, and the Lagrangian
term penalizes the constraint violation of $g^{c}(x,y)$.

However, recognizing the max-min subproblems involved in (2a), it becomes
computationally costly to evaluate the penalty function $F_{\gamma}(x)$ and
its gradient. To this end, we pose the following question

Can we develop efficient algorithms to solve the max-min subproblem and
evaluate $\nabla F_{\gamma}(x)$?

We answer this question by proving that this reformulation exhibits several
favorable properties, including smoothness. These properties are critical for
designing gradient-based algorithms and characterizing their performance.
However, the presence of the CCs renders the calculation of the gradient
$\nabla F_{\gamma}(x)$ more challenging than for its unconstrained counterpart
[61]. Building upon this, we design a primal-dual gradient method with
rigorous convergence guarantees for the BLO with general inequality CCs, and
provide an improved result for the case of $g^{c}$ being affine in $y$.

####  1.1 Main contributions

In a nutshell, our main contributions are outlined below.

  * C1)

In Section 2, leveraging the Lagrangian duality theorem, we introduce the
function $F_{\gamma}(x)$ in (2a) as a penalty-based reformulation of (1),
establish the continuity and smoothness of $F_{\gamma}(x)$, and develop a
novel way to compute its gradient.

  * C2)

In Section 3, we develop BLOCC, a fully first-order algorithm to tackle BLO
problems with CCs. With $\epsilon$ being the target error in terms of the
generalized gradient norm square of $F_{\gamma}(x)$, we establish that the
iteration complexity under the generic constraint $g^{c}(x,y)$ in (1b) is
$\tilde{\cal O}(\epsilon^{-2.5})$, which reduces to $\tilde{\cal
O}(\epsilon^{-1.5})$ when the constraint $g^{c}(x,y)$ is affine in $y$.

  * C3)

In Section 4, we apply the proposed BLOCC algorithm to two real-world
applications: SVM model training and planning of a transportation network. By
comparison with LV-HBA [72] and GAM [69], we demonstrate the algorithm’s
effectiveness.

| LL constraint | First Order | Complexity
---|---|---|---
BLOCC |  $g^{c}(x,y)\leq 0$ convex in $y$ and LICQ holds; | ✓ |  $\tilde{\cal O}(\epsilon^{-2.5})$;
Special case: $g^{c}(x,y)$ affine in $y$ | $\tilde{\cal O}(\epsilon^{-1.5})$
LV-HBA |  $g^{c}(x,y)\leq 0$ convex in $x\times y$ | ✓ | ${\cal O}(\epsilon^{-3})$
GAM |  $g^{c}(x,y)\leq 0$ convex in $x\times y$ and LICQ holds | ✘ | ✘
BVFSM |  $g^{c}(x,y)\leq 0$ satisfying other requirements | ✓ | ✘
AiPOD | $g^{c}(x,y)=Ay-b(x)=0$ | ✘ | ${\cal O}(\epsilon^{-1.5})$
Table 1: Comparison of our work with LV-HBA [72], BVFSM [45], AiPOD [68], and
Gradient Approximation method (GAM) [69]. LL convergence is on metric the
squared distance of $y_{t}$ to its optimal solution, and UL convergence is on
squared (generalized) gradient norm.

####  1.2 Related works

BLO has a long history, dating back to the seminal work of [8]. It has
inspired a rich body of literature, e.g., [73, 66, 15, 63]. The recent focus
on BLO is centered on developing efficient gradient-based approaches with
provable finite-time guarantees.

Methods for BLO without constraints. A cluster of BLO gradient-based
approaches gravitates around the implicit gradient descent (IGD) method [54],
where the key idea is to approximate the hypergradient by the implicit
function theorem. The finite-time convergence of IGD was first established in
[27] for unconstrained strongly-convex lower-level problems. Subsequent works
[31, 34, 12, 13, 39, 60, 42, 64, 70, 35] improved the convergence rates under
various settings. Another cluster of works is based on iterative
differentiation (ITD) [48, 23, 52, 59], which estimates the hypergradient by
differentiating the entire iterative algorithm used to solve the lower-level
problem with respect to the upper-level variables. The finite-time guarantee
was first established in [28, 46, 33, 6]. Penalty-based methods have also
emerged as a promising approach for BLO. Dated back to [74], this line of
works [44, 50, 43, 40, 61, 25, 47] reformulated the original BLO as the
single-level problems with various penalty terms and leveraged first-order
methods to solve them.

Figure 1: Calculation of $\nabla v(x)$. The blue line is $v(x)$, the red
dashed line is calculated by the formulation given in [61, 40], while the
yellow line is derived by our BLOCC. It can be seen that $\nabla v(x)$ without
the Lagrange multiplier is in the opposite direction of the true gradient.

BLO with constraints. While substantial progress has been made for
unconstrained BLO, the analysis for constrained BLO is more limited. Upper-
level constraints of the form $x\in\mathcal{X}$ were considered in [31, 12,
14]. For the lower-level uncoupled constraint, SIGD [38] considered the
uncoupled constraint $Ay\leq b$ and achieved asymptotic convergence, [41, 61]
employed penalty reformulation and considered both upper and lower uncoupled
constraints. However, the literature on BLO with CCs is scarce. BVFSM [45]
conducted a penalty-based method to avoid the calculation of the Hessian, as
IGD methods do. However, only asymptotic convergence was achieved. GAM [69]
investigated the IGD method under inequality constraints while failing to
provide finite-time convergence results as well. AiPOD [68] also applied IGD
and successfully achieved finite-time convergence, but it only considered
equality constraints. LV-HBA [72] considered inequality constraints and
constructed a penalty-based reformulation. However, it employed a joint
projection of $(x,y)$ onto $\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq
0\\}$ which is computationally inefficient when there are many constraints or
when $g^{c}(x,y)$ is not jointly convex. After our initial submission, we
found a concurrent work [71] posted on ArXiv, which used Lagrange duality
theory differently from ours, applying it to construct a new smoothed penalty
term. However, it does not quantify the relationship between the relaxation of
this penalty and the relaxation of lower-level optimality, and it does not
guarantee lower-level feasibility. We summarized prior works on BLO with
lower-level CCs in Table 1.

###  2 Primal-dual Penalty-based Reformulation

In this section, our goal is to construct a primal-dual-assisted penalty
reformulation for our BLOCC problem. The technical challenge comes from
finding a suitable penalty function for BLO with CCs.

####  2.1 The challenges in BLO with coupled constraints

Here, we will elaborate on the two technical challenges of BLO with CCs.

The first challenge associated with the presence of CCs is that the closed-
form expression of the gradient $\nabla v(x)$ without CCs, which is
$\nabla_{x}g(x,y_{g}^{*}(x))$ provided in [61, 41, 40, 43], is not applicable.
For example, when $g(x,y)=(y-2x)^{2}$ and $g^{c}(x,y)=3x-y$, the optimal
lower-level solution is $y_{g}^{*}(x)=3x$ and thus, $v(x)=x^{2}$ with $\nabla
v(x)=2x$. However, $\nabla_{x}g(x,y_{g}^{*}(x))=-4x\neq\nabla v(x)$.

The closed form gradient for BLO with CCs should be (8), which considers a
Lagrange term that is going to be illustrated later in this paper. In Figure
1, we present the gradient $\nabla v(x)$ without the Lagrange multiplier in
[61, 41, 40, 43] and ours with the Lagrange term. In this example, the
gradient without the Lagrange term leads to the opposite direction to the true
gradient.

The second challenge associated with the presence of CCs is the difficulty of
performing the joint projection. If we directly extend the penalty
reformulation in [61, 40], we could treat the coupling constraint set
$\mathcal{Y}(x)$ as a joint constraint $\\{(x,y):g^{c}(x,y)\leq 0\\}$, and
employ a joint projection of $(x,y)$ to ensure the feasibility. However, this
can be computationally inefficient when either the number of variables or
constraints is high. Moreover, when $g^{c}(x,y)$ is complex, the projection
may not have a closed form and may not even be well-defined if $g^{c}(x,y)$ is
not jointly convex.

####  2.2 The Lagrangian duality-based penalty reformulation

Before proceeding, we summarize the assumptions as follows.

######  Assumption 1 (Lipschitz Continuity).

Assume that $f$, $\nabla f$, $\nabla g$, $g^{c}$ and $\nabla g^{c}$, and are
respectively $l_{f,0}$, $l_{f,1}$, $l_{g,1}$, $l_{g^{c},0}$, and
$l_{g^{c},1}$-Lipschitz jointly over $(x,y)$, and $g$ is
$l_{g_{x},0}$-Lipschitz in $x$.

######  Assumption 2 (Convexity in $y$).

For any given $x\in\mathcal{X}$, $g(x,y)$ and $g^{c}(x,y)$ are
$\alpha_{g}$-strongly convex and convex in $y$, respectively.

######  Assumption 3 (Domain Feasibility).

Domain $\mathcal{X}\subset\mathbb{R}^{d_{x}}$ and
$\mathcal{Y}\subset\mathbb{R}^{d_{y}}$ are non-empty, closed, and convex. For
any given $x\in\mathcal{X}$,
$\mathcal{Y}(x):=\\{y\in\mathcal{Y}:g^{c}(x,y)\leq 0\\}$ is non-empty.

######  Assumption 4 (Constraint Qualification).

For any given $x\in\mathcal{X}$, $g^{c}$ satisfies the Linear Independence
Constraint Qualification (LICQ) condition in $y$.

The Lipschitz continuity in Assumption 1 and the strong convexity of $g$ over
$y$ in Assumption 2 are conventional [27, 31, 40, 68, 35, 13]. Moreover, we
only require $g^{c}(x,y)$ to be convex in $y$ rather than: i) jointly convex
in $(x,y)$ as in [72], or ii) linear as in [38, 68]. Assumption 3 pertains to
the convexity and closeness of the domain, which is also conventional, and
Assumption 4 is a standard constraint qualification condition.

To build a penalty reformulation for BLO with lower-level CCs, the first
challenge is to find a penalty that regulates $\|y-y_{g}^{*}(x)\|$. In the
following lemma, we show that $g(x,y)-v(x)$ is a good choice.

######  Lemma 1.

Suppose that Assumptions 2-4 hold and $v(x)$ is defined as in (2b). Then, it
holds that

  1. c1)

$g(x,y)-v(x)\geq\frac{\alpha_{g}}{2}\|y-y^{*}_{g}(x)\|^{2}$, for all
$y\in\mathcal{Y}(x)$; and

  2. c2)

$g(x,y)=v(x)$ if and only if $y=y^{*}_{g}(x)$, for all $y\in\mathcal{Y}(x)$.

The technical challenge of proving this lemma lies in showing the quadratic
growth property as in c1) of Lemma 1. For BLO without the CCs, this naturally
holds as the lower-level objective $g(x,y)$ is strongly convex in $y$. For BLO
with CCs, one needs to use the Lagrangian duality theorem. The full proof of
Lemma 1 is given in Appendix B.1, where the key challenge is to show the
invariance of the modulus of strong convexity in the Lagrangian reformulated
lower-level objective.

With $\gamma$ denoting a penalty constant, we consider the following penalty
reformulation:

$\displaystyle\min_{(x,y)\in\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq
0\\}}f(x,y)+\gamma(g(x,y)-v(x))$ (3)

where the penalty term confines the squared Euclidean distance from $y$ to
$y_{g}^{*}(x)$. Furthermore, to avoid projecting $(x,y)$ onto the
$\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq 0\\}$, we propose the primal-
dual-assisted penalty reformulation, which was defined in (2). The approximate
equivalence of the reformulation (2) to the original BLO problem (1) is
established in the following theorem.

######  Theorem 1 (Equivalence).

Suppose that $f$ is Lipschitz in $y$ and $\nabla f$ is $l_{f,1}$-Lipschitz in
$(x,y)$ in Assumption 1 hold, and Assumptions 2-4 hold. Then, solving the
$\epsilon$-approximation problem of (1):

$\displaystyle\min_{(x,y)\in\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq
0\\}}f(x,y)~{}~{}~{}\text{ s.t. }~{}~{}\|y-y^{*}_{g}(x)\|^{2}\leq\epsilon,$
(4)

is equivalent to solve the primal-dual penalty reformulation in (2) with
$\gamma={\cal O}(\epsilon^{-0.5})$ and $\gamma>\frac{l_{f,1}}{\alpha_{g}}$.

The detailed proof is provided in Appendix B.2. By setting $\gamma={\cal
O}(\epsilon^{-0.5})$, we effectively state the equivalence between the penalty
reformulation in (3) and the approximated original problem (4). We can then
decouple the joint minimization on $(x,y)$ to a min-min problem on $x$ and
$g^{c}$ constrained $y$. For the inner minimization problem in $y$, we choose
$\gamma$ in a way that $\gamma\alpha_{g}-l_{f,1}>0$ holds. This ensures the
objective in (2a) being strongly convex in $y$, as $l_{f,1}$-smoothness
ensures a lower bound for negative curvature of $f(x,y)$. Furthermore, the
convexity of $g^{c}(x,y)$ in $y$ validates the strong duality theorem in [58,
32], thereby enabling the max-min primal-dual reformulation in (2).

####  2.3 Smoothness of the penalty reformulation

To evaluate $F_{\gamma}(x)$ defined in (2a), we can find the solution to the
inner max-min problem:

$\displaystyle(\mu_{F}^{*}(x),y_{F}^{*}(x)):=\arg\max_{\mu\in\mathbb{R}^{dx}_{+}}\min_{y\in\mathcal{Y}}\underbrace{f(x,y)+\gamma(g(x,y)-v(x))+\langle\mu,g^{c}(x,y)\rangle}_{=:L_{F}(\mu,y;x)}.$
(5)

The uniqueness of $y_{F}^{*}(x)$ and $\mu_{F}^{*}(x)$ is guaranteed under
Assumptions 2 and 4 (see Lemma 5 in Appendix A). Therefore, $F_{\gamma}(x)$ in
(2a) can be evaluated using the unique optimal solutions by

$\displaystyle F_{\gamma}(x)=L_{F}(\mu_{F}^{*}(x),y_{F}^{*}(x);x).$ (6)

Similarly, we can evaluate $v(x)=L_{g}(\mu_{g}^{*}(x),y_{g}^{*}(x);x)$, where

$\displaystyle(\mu_{g}^{*}(x),y_{g}^{*}(x)):=\arg\max_{\mu\in\mathbb{R}^{d_{x}}_{+}}\min_{y\in\mathcal{Y}}\underbrace{g(x,y)+\langle\mu,g^{c}(x,y)\rangle}_{=:L_{g}(\mu,y;x)}.$
(7)

The penalty reformulation $F_{\gamma}(x)$ can hardly be convex, as $-v(x)$ may
be concave, even when $g(x,y)=g(y)$ and $g^{c}(x,y)=A^{\top}y-x$; see Lemma
4.24 in [58]. Instead, in the subsequent lemma, we not only show that $v(x)$
is differentiable, but also provide a closed-form expression of $\nabla v(x)$.

######  Lemma 2 (Danskin-like theorem for $v(x)$).

Suppose that Assumptions 1–4 hold, and let $B_{g}<\infty$ be a constant such
that $\|\mu_{g}^{*}(x)\|<B_{g}$ for all $x\in\mathcal{X}$. Then, it holds that

1\. $y_{g}^{*}(x)$ and $\mu_{g}^{*}(x)$ defined in (7) are $L_{g}$-Lipschitz
for some finite constant $L_{g}\geq 0$.

2\. $v(x)$ defined in (2b) is $l_{v,1}$-smooth where
$l_{v,1}\leq(l_{g,1}+B_{g}l_{g^{c},1})(1+L_{g})+l_{g^{c},0}L_{g}$ and

$\displaystyle\nabla
v(x)=\nabla_{x}g(x,y_{g}^{*}(x))+\langle\mu_{g}^{*}(x),\nabla_{x}g^{c}(x,y_{g}^{*}(x))\rangle.$
(8)

The assumption of the existence of the upper bound for the Lagrange multiplier
is mild; see e.g., [32, Theorem 1.6] and [72]. Finding $\nabla v(x)$ is
crucial for the design of a gradient-based method to solve
$\min_{x\in\mathcal{X}}F_{\gamma}(x)$. Leveraging Lemma 2, the gradient
$\nabla F_{\gamma}(x)$ can be obtained by the next lemma.

######  Lemma 3 (Danskin-like theorem for $F_{\gamma}(x)$).

Suppose that the conditions in Lemma 2 hold. Moreover, assume that
$\gamma>\frac{l_{f,1}}{\alpha_{g}}$, and there exist $B_{F}\\!<\\!\infty$ such
that $\|\mu_{F}^{*}(x)\|\\!<\\!B_{F},\,\forall x\in\mathcal{X}$. Then, it
holds that

1\. $y_{F}^{*}(x)$ and $\mu_{F}^{*}(x)$ defined in (5) are $L_{F}$-Lipschitz
for some constant $L_{F}\geq 0$.

2\. $F_{\gamma}(x)$ is $l_{F,1}$-smooth with $l_{F,1}\leq(l_{f,1}+\gamma
l_{g,1}+B_{F}l_{g^{c},1})(1+L_{F})+\gamma l_{v,1}+l_{f^{c},0}L_{F}$, and

$\displaystyle\nabla
F_{\gamma}(x)=\nabla_{x}f(x,y_{F}^{*}(x))+\gamma\left(\nabla_{x}g(x,y_{F}^{*}(x))-\nabla
v(x)\right)+\langle\mu_{F}^{*}(x),\nabla_{x}g^{c}(x,y^{*}_{F}(x))\rangle.$ (9)

The proof of Lemmas 2 and 3 is provided in Appendix C. Similar to [61], the
Danskin-like theorems in Lemmas 2 and 3 rely on the Lipschitzness of the
solutions in (5) and (7). Different from BLO without CCs [61], the results
here also hinge on the Lagrange multipliers $\mu_{g}^{*}(x)$ and
$\mu_{F}^{*}(x)$.

###  3 Main Results

We will first introduce an algorithm tailored for BiLevel Optimization
problems with inequality Coupled Constraints, and present its convergence
analysis. In Section 3.2, we will propose a primal-dual solver for the inner
max-min problems and characterize the overall convergence.

####  3.1 BLOCC algorithm

As $F_{\gamma}(x)$ features differentiability and smoothness, we can apply a
projected gradient descent (PGD)-based method to solve
$\min_{x\in\mathcal{X}}F_{\gamma}(x)$. At iteration $t$, update

$x_{t+1}=\operatorname{Proj}_{\mathcal{X}}(x_{t}-\eta g_{F,t})$ (10)

where $\eta$ is the step size and $g_{F,t}$ is an estimate of $\nabla
F_{\gamma}(x_{t})$.

In the previous section, we have obtained the closed-form expressions of
$\nabla v(x)$ in (8) and $\nabla F_{\gamma}(x)$ in (9). Evaluating the closed-
form expression requires finding $(y_{g}^{*}(x_{t}),\mu_{g}^{*}(x_{t}))$ and
$(y_{F}^{*}(x_{t}),\mu_{F}^{*}(x_{t}))$, the solutions to the max-min problem
$L_{g}(\mu,y;x)$ in (7) and $L_{F}(\mu,y;x)$ in (5) respectively.

Given $x_{t}\in\mathcal{X}$, we can use some minmax optimization solvers with
$T_{g}$ iterations on (7) and $T_{F}$ iterations on (5) to find an
$\epsilon_{g}$-solution $(y_{g,t}^{T_{g}},\mu_{g,t}^{T_{g}})$ and an
$\epsilon_{F}$-solution $(y_{F,t}^{T_{F}},\mu_{F,t}^{T_{F}})$ satisfying

$\displaystyle\|(y_{g,t}^{T_{g}},\mu_{g,t}^{T_{g}})-(y_{g}^{*}(x_{t}),\mu_{g}^{*}(x_{t}))\|^{2}={\cal
O}(\epsilon_{g});~{}~{}~{}\|(y_{F,t}^{T_{F}},\mu_{F,t}^{T_{F}})-(y_{F}^{*}(x_{t}),\mu_{F}^{*}(x_{t}))\|^{2}={\cal
O}(\epsilon_{F})$ (11)

for target estimation accuracy $\epsilon_{g},\epsilon_{F}>0$. Such an
effective minmax optimization solver will be introduced in the following
Section 3.2. In this way, $\nabla v(x_{t})$ in (9) can be estimated as

$\displaystyle
g_{v,t}=\nabla_{x}g(x_{t},y_{g,t}^{T_{g}})+\langle\mu_{g,t}^{T_{g}},\nabla_{x}g^{c}(x_{t},y_{g,t}^{T_{g}})\rangle.$
(12)

Leveraging $g_{v,t}$, the gradient $\nabla F_{\gamma}(x)$ can be estimated via
(8) as

$\displaystyle g_{F,t}$
$\displaystyle=\nabla_{x}f(x_{t},y_{F,t}^{T_{F}})+\gamma\left(\nabla_{x}g(x_{t},y_{F,t}^{T_{F}})-g_{v,t}\right)+\langle\mu_{F}^{T_{F}},\nabla_{x}g^{c}(x,y^{T_{F}}_{F,t})\rangle.$
(13)

We summarize the methods for finding $\min F_{\gamma}(x)$ as in Algorithm 1.
In the following, we present the convergence result of it allowing estimation
error $\epsilon_{g},\epsilon_{F}>0$ from the MaxMin Solver.

######  Theorem 2.

Suppose that the assumptions in Lemma 3 hold. Run Algorithm 1 with $T_{F}$ and
$T_{g}$ sufficiently large such that $(y_{g,t}^{T_{g}},\mu_{g,t}^{T_{g}})$ and
$(y_{F,t}^{T_{F}},\mu_{F,t}^{T_{F}})$ are respectively ${\cal
O}(\epsilon_{g})$ and ${\cal O}(\epsilon_{F})$-optimal in squared distance as
in (11). Set $\eta\leq\frac{1}{l_{F,1}}$ with some $l_{F,1}$ defined in Lemma
3. It then holds that

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\|G_{\eta}(x_{t})\|^{2}:=\frac{1}{T\eta^{2}}\sum_{t=0}^{T-1}\|(x_{t+1}-x_{t})\|^{2}=$
$\displaystyle{\cal O}(\gamma
T^{-1}+\gamma^{2}\epsilon_{F}+\gamma^{2}\epsilon_{g}).$ (14)

Using the projected gradient $G_{\eta}=\eta^{-1}(x_{t+1}-x_{t})$ as the
convergence metric for constrained optimization problems is standard [26]. The
proof of this theorem is available in Appendix D.1. The term ${\cal
O}(\gamma^{2}\epsilon_{F}+\gamma^{2}\epsilon_{g})$ arises from estimation
errors using the MaxMin Solver. In the next section, we develop a particular
MaxMin Solver for Algorithm 1.

####  3.2 MaxMin Solver for the BLO with inequality CCs

Algorithm 1 Our main algorithm BLOCC

1:inputs: initial points $x_{0}$, $y_{g,0}$, $\mu_{g,0}$, $y_{F,0}$,
$\mu_{F,0}$; step size $\eta$, $\eta_{g}$, $\eta_{F}$; counters $T_{g}$,
$T_{F}$.

2:for $t=0,1,\ldots,T$ do

3: $(y_{g,t}^{T_{g}},\mu_{g,t}^{T_{g}})=\textsf{MaxMin}(T_{g},T_{y}^{g})$.

4: $(y_{F,t}^{T_{F}},\mu_{F,t}^{T_{F}})=\textsf{MaxMin}(T_{F},T_{y}^{F})$.

5: update $x_{t+1}$ via (10) with $g_{F,t}$ in (13).

6:end for

7:outputs: $(x_{T},y_{F,T}^{T_{F}})$

Algorithm 2 Subroutine on MaxMin$(T,T_{y})$

1:inputs: initial values $\mu_{-1}=\mu_{0}\in\mathbb{R}^{d_{c}}_{+}$,
$y_{0}\in\mathcal{Y}$, step sizes $\eta_{1},\eta_{2}>0$, counters $T,T_{y}$.
Blue part is the version with acceleration; Red part is without acceleration.

2:for $t=0,\ldots,T-1$ do

3: update $\mu_{t+\frac{1}{2}}$ via (17) or $\mu_{t+\frac{1}{2}}=\mu_{t}$

4: for $t_{y}=0,\ldots,T_{y}-1$ do

5: update $y_{t,t_{y}+1}$ via (18) $\triangleright$ set $y_{t,0}=y_{t}$

6: end for

7: update $\mu_{t+1}$ via (19) $\triangleright$ set $y_{t+1}=y_{t,T_{y}}$

8:end for

9:outputs: $(y_{T},\mu_{T})$

In this section, we specify the MaxMin solver in Algorithm 1. By viewing
$L_{g}(\mu,y;x)$ in (7) and $L_{F}(\mu,y;x)$ in (5) as $L(\mu,y)$ for fixed
given $x\in\mathcal{X}$, we consider the following max-min problem

$\displaystyle\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}\min_{y\in\mathcal{Y}}L(\mu,y)$
(15)

which is concave (linear) in $\mu$ and strongly convex in $y$.

We can evaluate the dual function of $L(\mu,y)$ defined below by finding
$y_{\mu}^{*}(\mu)$.

$\displaystyle
D(\mu):=\min_{y\in\mathcal{Y}}L(\mu,y)=L(\mu,y_{\mu}^{*}(\mu))\quad\text{where}\quad
y_{\mu}^{*}(\mu):=\arg\min_{y\in\mathcal{Y}}L(\mu,y).$ (16)

According to Danskin’s theorem, we have $\nabla
D(\mu)=\nabla_{\mu}L(\mu,y^{*}_{\mu}(\mu))$. Taking either $L_{g}(\mu,y;x)$ or
$L_{F}(\mu,y;x)$ as $L(\mu,y)$ for given $x\in\mathcal{X}$, $D(\mu)$ defined
as (16) exhibits favorable properties namely smoothness and concavity, with
details illustrated in Lemma 10 in Appendix D.2. We can, therefore, apply
accelerated gradient methods designed for smooth and convex functions such as
[51, 4].

At each iteration $t$, we first perform a momentum-based update step to update
$\mu_{t+\frac{1}{2}}$ as

$\displaystyle\mu_{t+\frac{1}{2}}=\mu_{t}+\frac{t-1}{t+2}(\mu_{t}-\mu_{t-1}).$
(17)

To evaluate $\nabla
D(\mu_{t+\frac{1}{2}})=\nabla_{\mu}L(\mu_{t+\frac{1}{2}},y^{*}_{\mu}(\mu_{t+\frac{1}{2}})$,
with an arbitrary small target accuracy $\epsilon>0$, we can run $T_{y}={\cal
O}(\ln(\epsilon^{-1}))$ PGD steps on $L(\mu_{t+\frac{1}{2}},y)$ in $y$ via

$\displaystyle
y_{t,t_{y}+1}=\operatorname{Proj}_{\mathcal{Y}}\left(y_{t,t_{y}}-\eta_{1}\nabla_{y}L(\mu_{t+\frac{1}{2}},y_{t,t_{y}})\right).$
(18)

Defining the output after $T_{y}$ iterations as $y_{t+1}$, since strongly
convexity of $L(\mu,\cdot)$ ensures that PGD converges linearly [9, Theorem
3.10], it implies that $\|y_{t+1}-y_{\mu}^{*}(\mu_{t+\frac{1}{2}})\|={\cal
O}(\epsilon)$. We can conduct

$\displaystyle\mu_{t+1}=\operatorname{Proj}_{\mathbb{R}^{d_{c}}_{+}}(\mu_{t+\frac{1}{2}}+\eta_{2}\nabla_{\mu}L(\mu_{t+\frac{1}{2}},y_{t+1})).$
(19)

We summarize this oracle in Algorithm 2 based on an accelerated method [51].
When we skip the momentum update, i.e. setting $\mu_{t+\frac{1}{2}}=\mu_{t}$,
it is a simple PGD method on $D(\mu)$.

However, for a convex function $-D(\mu)$, the standard results only provide
convergence of the function value, i.e.
$\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D(\mu)-D(\mu_{t})\stackrel{{\scriptstyle
t\rightarrow\infty}}{{\longrightarrow}}0$. To establish the convergence of
$\|\mu_{g,t}^{T_{g}}-\mu_{g}^{*}(x_{t})\|$ and
$\|\mu_{F,t}^{T_{F}}-\mu_{F}^{*}(x_{t})\|$, we define the dual functions
associated with inner problems (7) and (5) as

$\displaystyle D_{g}(\mu)=\min_{y\in\mathcal{Y}}L_{g}(\mu,y;x)\quad\text{ and
}\quad D_{F}(\mu)=\min_{y\in\mathcal{Y}}L_{F}(\mu,y;x)$ (20)

and make the following curvature assumption near the optimum.

######  Assumption 5.

There exist $\delta_{g},\delta_{F}>0$ and $C_{\delta_{g}},C_{\delta_{F}}>0$
such that $\displaystyle\langle-\nabla D_{g}(\mu)+\nabla
D_{g}(\mu_{g}^{*}(x)),\mu-\mu_{g}^{*}(x)\rangle\geq$ $\displaystyle
C_{\delta_{g}}\|\mu-\mu_{g}^{*}(x)\|^{2},\quad\forall\mu\in\mathcal{B}(\mu_{g}^{*}(x);\delta_{g}),$
(21a) $\displaystyle\langle-\nabla D_{F}(\mu)+\nabla
D_{F}(\mu_{F}^{*}(x)),\mu-\mu_{F}^{*}(x)\rangle\geq$ $\displaystyle
C_{\delta_{F}}\|\mu-\mu_{F}^{*}(x)\|^{2},\quad\forall\mu\in\mathcal{B}(\mu_{F}^{*}(x);\delta_{F}).$
(21b)

It is worth noting that $\langle-\nabla D_{g}(\mu)+\nabla
D_{g}(\mu_{g}^{*}(x)),\mu-\mu_{g}^{*}(x)\rangle\geq 0$ holds for all $\mu$ due
to concavity and in a neighborhood of the optimal, the equality only happens
at the optimal due to the uniqueness of $\mu_{g}^{*}(x)$. The same argument
applies to $D_{F}$ due to the concavity of the dual functions. Therefore,
Assumption 5 essentially asserts a positive lower bound on the curvature of
the left-hand side term, which is mild as it only applies to the neighborhood
of the optima $\mu_{g}^{*}(x)$ and $\mu_{F}^{*}(x)$. It is also weaker than
the local strong concavity or global restricted secant inequality (RSI)
conditions.

By choosing Algorithm 2 with acceleration as the MaxMin solver, we provide the
convergence analysis of Algorithm 1 next, the proof of which can be found in
Appendix D.2.

######  Theorem 3.

Suppose that Assumptions 1–5 and the conditions in Theorem 2 hold. Let
$\gamma>\frac{l_{f,1}}{\alpha_{g}}$,
$\epsilon_{g}\leq\frac{C_{\delta_{g}}}{2}\delta_{g}$ and
$\epsilon_{F}\leq\frac{C_{\delta_{F}}}{2}\delta_{F}$. If we choose the inner
loop $T_{g}={\cal O}(\epsilon_{g}^{-0.5}),T_{F}={\cal
O}(\epsilon_{F}^{-0.5})$, and $T_{y}^{g}={\cal
O}(\log(\epsilon_{g}^{-1})),T_{y}^{F}={\cal O}(\log(\epsilon_{F}^{-1}))$, then
the iterates generated by Algorithm 1 satisfy

$\displaystyle\|\mu_{t,g}^{T_{g}}-\mu^{*}_{g}(x_{t})\|^{2}={\cal
O}(\epsilon_{g})\quad\text{and}\quad\|y_{t,g}^{T_{g}}-y^{*}_{g}(x_{t})\|^{2}={\cal
O}(\epsilon_{g})$
$\displaystyle\|\mu_{t,F}^{T_{F}}-\mu^{*}_{F}(x_{t})\|^{2}={\cal
O}(\epsilon_{F})\quad\text{and}\quad\|y_{t,F}^{T_{F}}-y^{*}_{F}(x_{t})\|^{2}={\cal
O}(\epsilon_{F}).$

For solving an $\epsilon$-approximation problem of BLO defined in (4), we
solve $\min_{x\in\mathcal{X}}F_{\gamma}(x)$ with $\gamma={\cal
O}(\epsilon^{-0.5})$ according to Theorem 1. In this way, to achieve
$\frac{1}{T}\sum_{t=0}^{T-1}\|G_{\eta}(x_{t})\|^{2}\leq\epsilon$ where
$\frac{1}{T}\sum_{t=0}^{T-1}\|G_{\eta}(x_{t})\|^{2}={\cal O}(\gamma
T^{-1}+\gamma^{2}\epsilon_{F}+\gamma^{2}\epsilon_{g})$ by (14), we need
$\epsilon_{g},\epsilon_{F}={\cal O}(\epsilon^{2})$, i.e. complexity
$\tilde{\cal O}(\epsilon^{-1})$ for the MaxMin solver in Algorithm 2, and
$T={\cal O}(\epsilon^{-1.5})$ for the number of iteration in BLOCC (Algorithm
1). Therefore, the overall complexity is $\tilde{\cal O}(\epsilon^{-2.5})$,
where $\tilde{\cal O}$ omits the $\log$ terms.

####  3.3 Special case of the MaxMin Solver: $g^{c}(x,y)$ being affine in $y$
and $\mathcal{Y}=\mathbb{R}^{d_{y}}$

In this section, we investigate a special case of BLO with CCs where
Assumption 5 automatically holds. Specifically, we focus on the case where
$g^{c}$ is affine in $y$ and $\mathcal{Y}=\mathbb{R}^{d_{y}}$, i.e.

$\displaystyle g^{c}(x,y)=g^{c}_{1}(x)^{\top}y-g^{c}_{2}(x).$ (22)

In this case, fixing $x$, taking $L_{g}(\mu,y;x)$ in (7) and $L_{F}(\mu,y;x)$
as $L(\mu,y)$, (16) gives

$\displaystyle D_{g}(\mu)=$
$\displaystyle\min_{y\in\mathbb{R}^{d_{y}}}-\langle
g^{c}_{2}(x),\mu\rangle+\langle
y,g^{c}_{1}(x)\mu\rangle+g(x,y)\quad\text{and}$ (23a) $\displaystyle
D_{F}(\mu)=$ $\displaystyle\min_{y\in\mathbb{R}^{d_{y}}}-\langle
g^{c}_{2}(x),\mu\rangle+\langle
y,g^{c}_{1}(x)\mu\rangle+f(x,y)+\gamma(g(x,y)-v(x)).$ (23b)

When $g^{c}_{1}(x)$ is of full column rank, both $D_{g}(\mu)$ and $D_{F}(\mu)$
are strongly concave according to Lemma 13 in Appendix so that Assumption 5
holds globally. Moreover, when applying PGD on $-D_{g}(\mu)$ and
$-D_{F}(\mu)$, strongly convexity guarantees linear convergence (Theorem 3.10
in [9]). Therefore, even without acceleration, Algorithm 2 with $T_{y}$
sufficiently large performs PGD on $-D(\mu)$ and it converges linearly up to
inner loop accuracy. Moreover, PGD on $L(\mu,y)$ in $y$ also converges
linearly. This motivates us to implement a single-loop version ($T_{y}=1$) of
Algorithm 2.

When both $y$ and $\mu$ are unconstrained, the analysis has been established
in [19]. However, as $\mu$ is constrained to $\mathbb{R}_{+}^{d_{c}}$ in the
Lagrangian formulation for inequality constraints, the convergence analysis in
[19] is not applicable, and the extension is nontrivial due to the non-
differentiability of the projection. We address this technical challenge by
treating $\mathbb{R}^{d_{c}}_{+}$ as an inequality constraint and
reformulating it as an unconstrained problem using Lagrange duality theory.
This approach demonstrates that a single-loop $T_{y}=1$ update in Algorithm 2,
without acceleration, achieves linear convergence for the max-min problem
(15). The detailed proof is provided in Appendix D.3.

Thus, by selecting the single-loop version ($T_{y}=1$) of Algorithm 2 without
acceleration as the MaxMin solver in Algorithm 1, we establish the following
theorem.

######  Theorem 4 (Inner linear convergence).

Consider BLO with $\mathcal{Y}=\mathbb{R}^{d_{y}}$ and $g^{c}(x,y)$ defined in
(22). Suppose Assumptions 1–4 and the conditions in Theorem 2 hold. Suppose
for any $x\in\mathcal{X}$, there exist constants $s_{\min}$ and $s_{\max}$
such that
$0<s_{\min}\leq\sigma_{\min}(g^{c}_{1}(x))\leq\sigma_{\max}(g^{c}_{1}(x))\leq
s_{\max}<\infty$. Let $\gamma>\frac{l_{f,1}}{\alpha_{g}}$. If we choose the
inner loop $T_{g}={\cal O}(\log(\epsilon_{g}^{-1})),T_{F}={\cal
O}(\log(\epsilon_{F}^{-1}))$, and $T_{y}^{g}=T_{y}^{F}=1$ with proper constant
stepsizes, then the iterates generated by Algorithm 1 satisfy

$\displaystyle\|\mu_{t,g}^{T_{g}}-\mu^{*}_{g}(x_{t})\|^{2}={\cal
O}(\epsilon_{g})\quad\text{and}\quad\|y_{t,g}^{T_{g}}-y^{*}_{g}(x_{t})\|^{2}={\cal
O}(\epsilon_{g})$
$\displaystyle\|\mu_{t,F}^{T_{F}}-\mu^{*}_{F}(x_{t})\|^{2}={\cal
O}(\epsilon_{F})\quad\text{and}\quad\|y_{t,F}^{T_{F}}-y^{*}_{F}(x_{t})\|^{2}={\cal
O}(\epsilon_{F}).$

The proof is available in Appendix D.3. Similar to the previous analysis, we
choose $\gamma={\cal O}(\epsilon^{-0.5})$ to solve the equivalent (2a) to
$\epsilon$-approximation problem of BLO (4). To achieve
$\frac{1}{T}\sum_{t=0}^{T-1}\|G_{\eta}(x_{t})\|^{2}\leq\epsilon$, the number
of iteration in BLOCC $T={\cal O}(\epsilon^{-1.5})$ by (14). As the inner
MaxMin solver has linear convergence, the overall complexity is $\tilde{\cal
O}(\epsilon^{-1.5})$ where $\tilde{\cal O}$ omits the $\log$ terms. We
summarized the overall complexity of our BLOCC algorithm in different settings
in Table 1.

###  4 Numerical Experiments

This section reports the results of numerical experiments for three different
problems: a toy example used to validate our method, an SVM training
application, and a network design problem in both synthetic and real-world
transportation scenarios. In the last two experiments, we compare the proposed
algorithm with two baselines, LV-HBA [72] and GAM [69]. The code is available
at https://github.com/Liuyuan999/Penalty_Based_Lagrangian_Bilevel.

####  4.1 Toy example

Figure 2: 3-D plot of the upper-level objective $f(x,y)$ of the toy example,
with the line $f(x,y)|_{y=x}$ shown in dashed red and the convergence points
marked as red dots.

Consider the BLO problem with an inequality coupled constraint
$\mathcal{Y}(x):=\\{y\in\mathcal{Y}:y-x\leq 0\\}$, given by

$\displaystyle\min_{x\in[0,3]}\;$
$\displaystyle\;f(x,y^{*}_{g}(x))=\frac{e^{-y^{*}_{g}(x)+2}}{2+\cos(6x)}+\frac{1}{2}\log((4x-2)^{2}+1)$
with
$\displaystyle\;y^{*}_{g}(x)\in\text{arg}\min_{y\in\mathcal{Y}(x)}\;g(x,y)=(y-2x)^{2}.$
(24)

Problem (4.1) satisfies all assumptions for Theorem 2 and Theorem 4. The
lower-level problem is strongly convex and $y^{*}_{g}(x)=x$. Therefore, the
BLO problem with inequality constraint in (4.1) reduces to
$\min_{x\in[0,3]}f(x,y)|_{y=x}$. In Figure 2, we plot the dashed line as the
intersected line of the surface $f(x,y)$ and the plane $f(x,y_{g}^{*}(x))$,
and the red points as the converged points by running BLOCC with $\gamma=5$
and with 200 different initialization values. It can be seen that BLOCC
consistently finds the local minima, verifying the effectiveness.

Methods | diabetes | fourclass
---|---|---
BLOCC | $\bf 0.767\pm 0.039$ | $\bf 0.761\pm 0.014$
(Ours) | $(\bf 1.729\pm 0.529)$ | $(\bf 1.922\pm 0.108)$
LV-HBA | $0.765\pm 0.039$ | $0.748\pm 0.060$
$(2.899\pm 1.378)$ | $(2.404\pm 0.795)$
GAM | $0.721\pm 0.047$ | $0.715\pm 0.056$
$(8.752\pm 4.736)$ | $(13.481\pm 0.970)$

Table 2: Numerical results on the training outcome of our BLOCC with LV-HBA
[72] and GAM [69]. The first row represents accuracy $\pm$ standard deviation,
and the second row between brackets represents the running time until the
upper-level objective’s update is smaller than 1e-5.

####  4.2 Hyperparameter optimization for SVM

We test the performance of our algorithm BLOCC with $\gamma=10$ when training
a linear SVM model on the diabetes [20] and fourclass datasets [29]. The model
is trained via the BLO formulation (1) with inequality CCs; see the details in
Appendix E. We compared performance with two baselines, LV-HBA [72] and GAM
[69].

Table 2 shows that the model trained by our BLOCC algorithm outperforms that
of GAM [69] significantly and is of a similar level as that of LV-HBA [72].
Looking into the test accuracy in Figure 3, our algorithm achieves more than
0.76 accuracy in the first 2 iterations, which is significantly better than
the other ones. For the upper-level objective, in the first few iterations,
the loss decreases significantly under all algorithms, in which our BLOCC
achieves the lowest results. Moreover, in the right plot of Figure 3, the
lower-level optimum for LV-HBA was not attained until the very end. This
indicates that the decrease of upper loss between 0-40 iterations in the
middle figure may be due to the suboptimality of lower-level variables. For
our BLOCC and GAM, the lower-level minimum is attained as lower-level
objective $g\geq 0$ in this case.

Figure 3: Test accuracy (left), upper loss $f(x,y)$ (middle), and lower loss
$g(x,y)$ (right) for the SVM on the diabetes dataset. The experiments are
executed for 50 different random train-validation-test splits, with the bold
line representing the mean, and the shaded regions being the standard
deviation.

####  4.3 Transportation network design problem

| NN = 3, NV = 48 | NN = 9, NV = 2,262 | NN = 26, NV = 49,216
---|---|---|---
| NC = 24, NZ = 126 | NC = 678, NZ = 6,222 | NC = 13,336, NZ = 129,256
Methods | Runtime (s) | UL utility | Runtime (s) | UL utility | Runtime (s) | UL utility
LV-HBA | $3.51e2$ | 1.53 | / | / | / | /
BLOCC (Ours)-$\gamma=2$ | $2.02e1$ | 1.07 | $7.27e2$ | 8.09 | $6.82e4$ | 98.40
BLOCC (Ours)-$\gamma=3$ | $2.00e1$ | 1.69 | $8.50e2$ | 10.37 | $6.42e4$ | 111.39
BLOCC (Ours)-$\gamma=4$ | $2.01e1$ | 1.71 | $8.68e2$ | 11.04 | $6.70e4$ | 138.78
Table 3: Results of the transportation experiment, both in terms of running
time (Runtime) and convergenced upper-level objective value (UL utility,
larger there better), with stepsize $\eta=1.6e{-4}$. We use ”/” for algorithms
that cannot converge within 24 hours of execution. NN (Number of Nodes) is the
number of stations in the network. Analogously, NV (Number of Variables), NC
(Number of Constraints), and NNZ (Number of non-zero elements) are the number
of optimization variables, constraints, and non-zero elements of the
constraints matrix, respectively.

BLO is particularly relevant in transportation, where network planning must
consider different time horizons and actors. Those problems have a large
number of CCs, challenging the use of traditional BLO techniques. As a result,
our final experiment considers a network design problem, where we act as an
operator (whose utility is modeled as the upper-level objective) that
considers the passengers’ behavior, modeled in the lower level. We considered
three networks: two synthetic networks of 3 and 9 nodes, respectively, and a
real-world network of 26 nodes in the city of Seville, Spain. Further details
about the formulation and the experiment can be found in Appendix F.

In this experiment, we only compare our BLOCC with LV-HBA, as GAM involves
calculating the inverse of the Hessian, which is not tractable in this
setting. Table 3 summarizes the experiment results. We can see that LV-HBA
failed to work efficiently, especially for large networks. This is mainly
because the increased constraints render the projection step impracticable.
Our BLOCC, in contrast, is much faster, and it successfully converges even
with large real-world networks.

###  5 Conclusions and Future Work

This paper proposed a novel primal-dual-assisted penalty reformulation for BLO
problems with coupled lower-level constraints, and developed a new first-order
method BLOCC to solve the resultant problem. The non-asymptotic convergence
rate of our algorithm is $\tilde{\cal O}(\epsilon^{-2.5})$, tightening to
$\tilde{\cal O}(\epsilon^{-1.5})$ when the lower-level constraints are affine
in $y$ without any other constraints. Our method achieves the best-known
convergence rate and is projection-free, making it more favorable for large-
scale, high-dimensional, constrained BLO problems. Experiments on SVM model
training and transportation network planning showcased the effectiveness of
our algorithm.

### References

  * [1] Seyed Mehdi Alizadeh, Patrice Marcotte, and Gilles Savard.  Two-stage stochastic bilevel programming over a transportation network.  Transportation Research Part B: Methodological, 58:92–105, 2013.
  * [2] Aram V Arutyunov, Evgeniy R Avakov, and Alexey F Izmailov.  Directional regularity and metric regularity.  SIAM Journal on Optimization, 18(3):810–833, 2007.
  * [3] Dominique Azé and Jean-Paul Penot.  Uniformly convex and uniformly smooth convex functions.  In Annales de la Faculté des sciences de Toulouse: Mathématiques, 1995.
  * [4] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences, 2(1):183–202, 2009.
  * [5] Moshe E Ben-Akiva and Steven R Lerman.  Discrete choice analysis: theory and application to travel demand, volume 9.  MIT press, 1985.
  * [6] J. Bolte, E. Pauwels, and S. Vaiter.  Automatic differentiation of nonsmooth iterative algorithms.  In Advances in Neural Information Processing Systems, 2022.
  * [7] J Frédéric Bonnans and Alexander Shapiro.  Perturbation analysis of optimization problems.  Springer Science & Business Media, 2013.
  * [8] Jerome Bracken and James T McGill.  Mathematical programs with optimization problems in the constraints.  Operations Research, 21(1):37–44, 1973.
  * [9] Sébastien Bubeck.  Convex Optimization: Algorithms and Complexity.  Foundations and Trends® in Machine Learning, 2015.
  * [10] Luis Cadarso and Angel Marin.  Combining robustness and recovery in rapid transit network design.  Transportmetrica A: Transport Science, 12:1–26, 11 2015.
  * [11] Ennio Cascetta.  Transportation systems analysis: models and applications, volume 29.  Springer Science & Business Media, 2009.
  * [12] Tianyi Chen, Yuejiao Sun, Quan Xiao, and Wotao Yin.  A single-timescale method for stochastic bilevel optimization.  In International Conference on Artificial Intelligence and Statistics, virtual, 2022.
  * [13] Tianyi Chen, Yuejiao Sun, and Wotao Yin.  Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems.  In Advances in Neural Information Processing Systems, Virtual, 2021.
  * [14] Ziyi Chen, Bhavya Kailkhura, and Yi Zhou.  A fast and convergent proximal algorithm for regularized nonconvex and nonsmooth bi-level optimization.  arXiv preprint arXiv:2203.16615, 2022.
  * [15] Benoît Colson, Patrice Marcotte, and Gilles Savard.  An overview of bilevel optimization.  Annals of operations research, 153(1):235–256, 2007.
  * [16] Jean-Philippe Côté, Patrice Marcotte, and Gilles Savard.  A bilevel modelling approach to pricing and fare optimisation in the airline industry.  Journal of Revenue and Pricing Management, 2:23–36, 2003.
  * [17] Olivier Devolder, François Glineur, and Yurii Nesterov.  First-order methods of smooth convex optimization with inexact oracle.  Mathematical Programming, 146:37–75, 2014.
  * [18] Asen L Dontchev and R Tyrrell Rockafellar.  Implicit functions and solution mappings, volume 543.  Springer, 2009.
  * [19] Simon S Du and Wei Hu.  Linear convergence of the primal-dual gradient method for convex-concave saddle point problems without strong convexity.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 196–205, 2019.
  * [20] Dheeru Dua and Casey Graff.  Uci machine learning repository, 2017.  Accessed: 2024-05-21.
  * [21] Laureano Escudero and Susana Muñoz.  An approach for solving a modification of the extended rapid transit network design problem.  Top, 17:320–334, 12 2009.
  * [22] Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In International Conference on Machine Learning, pages 1126–1135, Sydney, Australia, 2017.
  * [23] Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil.  Forward and reverse gradient-based hyperparameter optimization.  In International Conference on Machine Learning, pages 1165–1173, 2017.
  * [24] Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil.  Bilevel programming for hyperparameter optimization and meta-learning.  In International Conference on Machine Learning, Stockholm, Sweden, 2018.
  * [25] Lucy L Gao, Jane Ye, Haian Yin, Shangzhi Zeng, and Jin Zhang.  Value function based difference-of-convex algorithm for bilevel hyperparameter selection problems.  In International Conference on Machine Learning, pages 7164–7182, Baltimore, MD, 2022.
  * [26] Saeed Ghadimi, Guanghui Lan, and Hongchao Zhang.  Mini-batch stochastic approximation methods for nonconvex stochastic composite optimization.  Mathematical Programming, 155(1):267–305, 2016.
  * [27] Saeed Ghadimi and Mengdi Wang.  Approximation methods for bilevel programming.  arXiv preprint arXiv:1802.02246, 2018.
  * [28] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo.  On the iteration complexity of hypergradient computation.  In International Conference on Machine Learning, pages 3748–3758, virtual, 2020.
  * [29] Tin Kam Ho and Eugene M Kleinberg.  Building projectable classifiers of arbitrary complexity.  In Proceedings of 13th International Conference on Pattern Recognition, volume 2, pages 880–885, 1996.
  * [30] Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola.  Kernel methods in machine learning.  2008\.
  * [31] Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang.  A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic.  SIAM Journal on Optimization, 33(1), 2023.
  * [32] Kazufumi Ito and Karl Kunisch.  Lagrange Multiplier Approach to Variational Problems and Applications, volume 15 of Advances in Design and Control.  Not specified in the provided information, 2008.  Includes bibliographical references and index.
  * [33] Kaiyi Ji, Mingrui Liu, Yingbin Liang, and Lei Ying.  Will bilevel optimizers benefit from loops.  arXiv preprint arXiv:2205.14224, 2022.
  * [34] Kaiyi Ji, Junjie Yang, and Yingbin Liang.  Provably faster algorithms for bilevel optimization and applications to meta-learning.  In Advances in Neural Information Processing Systems, 2020.
  * [35] Kaiyi Ji, Junjie Yang, and Yingbin Liang.  Bilevel optimization: Convergence analysis and enhanced design.  In International Conference on Machine Learning, virtual, 2021.
  * [36] Sham Kakade, Shai Shalev-Shwartz, Ambuj Tewari, et al.  On the duality of strong convexity and strong smoothness: Learning applications and matrix regularization.  Unpublished Manuscript, http://ttic. uchicago. edu/shai/papers/KakadeShalevTewari09. pdf, 2(1):35, 2009.
  * [37] Hamed Karimi, Julie Nutini, and Mark Schmidt.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16, pages 795–811, 2016.
  * [38] Prashant Khanduri, Ioannis Tsaknakis, Yihua Zhang, Jia Liu, Sijia Liu, Jiawei Zhang, and Mingyi Hong.  Linearly constrained bilevel optimization: A smoothed implicit gradient approach.  In International Conference on Machine Learning, pages 16291–16325, 2023.
  * [39] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang.  A near-optimal algorithm for stochastic bilevel optimization via double-momentum.  In Advances in Neural Information Processing Systems, Virtual, 2021.
  * [40] Jeongyeol Kwon, Dohyun Kwon, Stephen Wright, and Robert D Nowak.  A fully first-order method for stochastic bilevel optimization.  In International Conference on Machine Learning, pages 18083–18113, 2023.
  * [41] Jeongyeol Kwon, Dohyun Kwon, Steve Wright, and Robert Nowak.  On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation.  In International Conference on Learning Representations, Vienna, Austria, 2024.
  * [42] Junyi Li, Bin Gu, and Heng Huang.  A fully single loop algorithm for bilevel optimization without hessian inverse.  In Association for the Advancement of Artificial Intelligence, pages 7426–7434, virtual, 2022.
  * [43] Bo Liu, Mao Ye, Stephen Wright, Peter Stone, and Qiang Liu.  Bome! bilevel optimization made easy: A simple first-order approach.  Advances in neural information processing systems, 35:17248–17262, 2022.
  * [44] Risheng Liu, Xuan Liu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang.  A value-function-based interior-point method for non-convex bi-level optimization.  In International conference on machine learning, pages 6882–6892, 2021.
  * [45] Risheng Liu, Xuan Liu, Shangzhi Zeng, Jin Zhang, and Yixuan Zhang.  Value-function-based sequential minimization for bi-level optimization.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023.
  * [46] Risheng Liu, Yaohua Liu, Shangzhi Zeng, and Jin Zhang.  Towards gradient-based bilevel optimization with non-convex followers and beyond.  In Advances in Neural Information Processing Systems, volume 34, pages 8662–8675, 2021.
  * [47] Songtao Lu.  Slm: A smoothed first-order lagrangian method for structured constrained nonconvex optimization.  Advances in Neural Information Processing Systems, 36, 2023.
  * [48] Dougal Maclaurin, David Duvenaud, and Ryan Adams.  Gradient-based hyperparameter optimization through reversible learning.  In International Conference on Machine Learning, pages 2113–2122, Lille, France, 2015.
  * [49] Patrice Marcotte.  Network design problem with congestion effects: A case of bilevel programming.  Mathematical programming, 34(2):142–162, 1986.
  * [50] Akshay Mehra and Jihun Hamm.  Penalty method for inversion-free deep bilevel optimization.  In Asian conference on machine learning, pages 347–362, 2021.
  * [51] Yu Nesterov.  Smooth minimization of non-smooth functions.  Mathematical programming, 103:127–152, 2005.
  * [52] Alex Nichol, Joshua Achiam, and John Schulman.  On first-order meta-learning algorithms.  arXiv preprint arXiv:1803.02999, 2018.
  * [53] Norbert Oppenheim.  Equilibrium trip distribution/assignment with variable destination costs.  Transportation Research Part B: Methodological, 27(3):207–217, 1993.
  * [54] Fabian Pedregosa.  Hyperparameter optimization with approximate gradient.  In International conference on machine learning, pages 737–746, 2016.
  * [55] Fernando Real-Rojas.  Diseño de redes de transporte: integrando la competencia modal con técnicas de optimización convexa, BSc Thesis, King Juan Carlos University. June 2024.
  * [56] Fernando Real-Rojas, Victor M Tenorio, and Antonio G Marques.  A sparse nonlinear approach for designing hub-and-spoke air transportation networks.  In Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, 2024.
  * [57] R. Tyrrell Rockafellar.  Convex Analysis.  Princeton University Press, Princeton, NJ, USA, 1970s.
  * [58] Andrzej Ruszczyński.  Nonlinear Optimization.  Princeton University Press, 2006.
  * [59] Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots.  Truncated back-propagation for bilevel optimization.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1723–1732, 2019.
  * [60] Han Shen and Tianyi Chen.  A single-timescale analysis for stochastic approximation with multiple coupled sequences.  Advances in Neural Information Processing Systems, 35:17415–17429, 2022.
  * [61] Han Shen, Quan Xiao, and Tianyi Chen.  On penalty-based bilevel gradient descent method.  In International Conference on Machine Learning, Honolulu, HI, 2023.
  * [62] Han Shen, Zhuoran Yang, and Tianyi Chen.  Principled penalty-based methods for bilevel reinforcement learning and RLHF.  In International Conference on Machine Learning, Vienna, Austria, 2024.
  * [63] Ankur Sinha, Pekka Malo, and Kalyanmoy Deb.  A review on bilevel optimization: from classical to evolutionary approaches and applications.  IEEE Transactions on Evolutionary Computation, 22(2):276–295, 2017.
  * [64] Daouda Sow, Kaiyi Ji, and Yingbin Liang.  On the convergence theory for hessian-free bilevel algorithms.  volume 35, pages 4136–4149, 2022.
  * [65] Bradly Stadie, Lunjun Zhang, and Jimmy Ba.  Learning intrinsic rewards as a bi-level optimization problem.  In Conference on Uncertainty in Artificial Intelligence, virtual, 2020.
  * [66] Luis N Vicente and Paul H Calamai.  Bilevel and multilevel programming: A bibliography review.  Journal of Global optimization, 5(3):291–306, 1994.
  * [67] Gerd Wachsmuth.  On licq and the uniqueness of lagrange multipliers.  Operations Research Letters, 41(1):78–80, 2013.
  * [68] Quan Xiao, Han Shen, Wotao Yin, and Tianyi Chen.  Alternating implicit projected sgd and its efficient variants for equality-constrained bilevel optimization.  In International Conference on Artificial Intelligence and Statistics, 2023.
  * [69] Siyuan Xu and Minghui Zhu.  Efficient gradient approximation method for constrained bilevel optimization.  In Proceedings of the AAAI Conference on Artificial Intelligence, 2023.
  * [70] Haikuo Yang, Luo Luo, Chris Junchi Li, Michael Jordan, and Maryam Fazel.  Accelerating inexact hypergradient descent for bilevel optimization.  In OPT 2023: Optimization for Machine Learning, 2023.
  * [71] Wei Yao, Haian Yin, Shangzhi Zeng, and Jin Zhang.  Overcoming lower-level constraints in bilevel optimization: A novel approach with regularized gap functions.  arXiv preprint arXiv:2406.01992, 2024.
  * [72] Wei Yao, Chengming Yu, Shangzhi Zeng, and Jin Zhang.  Constrained bi-level optimization: Proximal lagrangian value function approach and hessian-free algorithm.  arXiv preprint arXiv:2401.16164, 2024.
  * [73] Jane J Ye and Daoli Zhu.  Optimality conditions for bilevel programming problems.  Optimization, 33(1):9–27, 1995.
  * [74] Jane J Ye, Daoli Zhu, and Qiji Jim Zhu.  Exact penalization and necessary optimality conditions for generalized bilevel programming problems.  SIAM Journal on optimization, 7(2):481–507, 1997.

Supplementary Material for “A Primal-Dual-Assisted Penalty Approach to Bilevel
Optimization with Coupled Constraints”

##

###  Appendix A Preliminaries

This section will provide some preliminaries for our subsequent analysis.

######  Definition 5.

For a convex function $h:\mathbb{R}^{d_{q}}\rightarrow\mathbb{R}$ whose domain
is $\mathcal{Q}\subseteq\mathbb{R}^{d_{q}}$, the Legendre conjugate of
$h^{*}:\mathcal{Q}^{*}\rightarrow\mathbb{R}$ is defined as:

$\displaystyle h^{*}(q):=\sup_{q^{\prime}\in\mathcal{Q}}\\{\langle
q^{\prime},q\rangle-h(q^{\prime})\\}=-\inf_{q^{\prime}\in\mathcal{Q}}\\{-\langle
q^{\prime},q\rangle+h(q^{\prime})\\},$ $\displaystyle\forall
q\in\mathcal{Q}^{*}:=\\{q\in\mathbb{R}^{d_{q}}:\sup_{q^{\prime}\in\mathcal{Q}}\\{\langle
q^{\prime},q\rangle-h(q^{\prime})\\}<\infty\\}.$

######  Remark 1.

When $h$ is strongly convex in $\mathbb{R}^{d_{q}}$, it is lower bounded and
therefore $\mathcal{Q}^{*}=\mathbb{R}^{d_{y}}$.

######  Lemma 4.

Suppose $h:\mathbb{R}^{d_{y}}\rightarrow\mathbb{R}$ is $l_{h,1}$-smooth and
$\alpha_{h}$-strongly convex and its domain
$\mathcal{Q}\subseteq\mathbb{R}^{d_{q}}$ is convex, closed and non-empty.

  1. 1.

If $\mathcal{Q}=\mathbb{R}^{d_{q}}$, the gradient mappings $\nabla h$ and
$\nabla h^{*}$ are inverse of each other ([57]); and
$h^{*}:\mathbb{R}^{d_{q}}\rightarrow\mathbb{R}$ is
$\frac{1}{\alpha_{h}}$-smooth and $\frac{1}{l_{h,1}}$-strongly convex
(Proposition 2.6 [3]).

  2. 2.

If $\mathcal{Q}\subset\mathbb{R}^{d_{q}}$, $h^{*}$ is
$\frac{1}{\alpha_{h}}$-smooth ([36]) and and convex (Theorem 4.43 [32]).

######  Lemma 5.

Suppose $\mathcal{Q}\subseteq\mathbb{R}^{d_{q}}$ is convex, closed and non-
empty, $h:\mathbb{R}^{d_{q}}\rightarrow\mathbb{R}$ is strongly convex on
$\mathcal{Q}$, $h^{c}:\mathbb{R}^{d_{q}}\rightarrow\mathbb{R}^{d_{c}}$ is
convex on $\mathcal{Q}$ and $d_{c}$ is finite, and
$\\{q\in\mathcal{Q}:h^{c}(q)\leq 0\\}$ is non-empty.

  1. 1.

The problem $\min_{q\in\\{q\in\mathcal{Q}:h^{c}(q)\leq 0\\}}h(q)$ has a unique
feasible solution.

  2. 2.

When linear independence constraint qualification (LICQ) condition holds for
$h^{c}(q)$, the Lagrange multiplier for the problem
$\min_{q\in\\{q\in\mathcal{Q}:h^{c}(q)\leq 0\\}}h(q)$, i.e. solution to the
problem
$\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}\min_{q\in\mathcal{Q}}h(q)+\langle\mu,h^{c}(q)\rangle$
is unique [67].

######  Lemma 6 (Lemma 3.1 in [9]).

Suppose $\mathcal{Q}\subseteq\mathbb{R}^{d_{q}}$ is convex, closed, and
nonempty. For any $q_{1}\in\mathbb{R}^{d_{q}}$ and any $q_{2}\in\mathcal{Q}$,
it follows that

$\displaystyle\langle\operatorname{Proj}_{\mathcal{Q}}(q_{1})-q_{2},\operatorname{Proj}_{\mathcal{Q}}(q_{1})-q_{1}\rangle\leq
0.$ (25)

In this way, take $q_{1}=q_{3}-\eta g$ for any $q_{3}\in\mathcal{Q}$, and
denote $q_{3}^{\eta,g}=\operatorname{Proj}_{\mathcal{Q}}(q_{3}-\eta g)$ as a
projected gradient update in direction $g$ with step size $\eta$, we have,

$\displaystyle\langle g,q_{3}^{\eta,g}-q_{2}\rangle\leq-\frac{1}{\eta}\langle
q_{3}^{\eta,g}-q_{2},q_{3}^{\eta,g}-q_{3}\rangle,\quad\forall
q_{2},q_{3}\in\mathcal{Q}.$ (26)

######  Lemma 7 (Theorem 3.10 [9]).

Suppose a differentiable function $h$ is $l_{h,1}$-smooth and
$\alpha_{h_{2}}$-strongly convex. Consider the constrained problem
$\min_{q\in\mathcal{Q}}h(q)$ where $\mathcal{Q}$ is non-empty, closed and
convex. Projected Gradient Descent with $\eta\leq\frac{1}{l_{h,1}}$ converges
linearly to the unique $q^{*}=\arg\min_{q\in\mathcal{Q}}h(q)$:

$\displaystyle\|\operatorname{Proj}_{\mathcal{Q}}(q-\eta\nabla
h(q))-q^{*}\|\leq(1-\alpha\eta)^{1/2}\|q-q^{*}\|\leq(1-\alpha\eta/2)\|q-q^{*}\|,\quad\forall
q\in\mathcal{Q}.$ (27)

###  Appendix B Analysis of the Penalty-Based Lagrangian Reformulation

####  B.1 Proof of Lemma 1

According to Lemma 5, for any fixed $x$, there exists a unique
$\mu_{g}^{*}(x)$. Therefore, according to Lagrange duality theorem, for any
fixed $x$, the primal problem

$\displaystyle\min_{y\in\mathcal{Y}}g(x,y)\quad\text{s.t.}\quad g^{c}(x,y)\leq
0~{}~{}~{}\Leftrightarrow~{}~{}~{}~{}\min_{y\in\mathcal{Y}}g(x,y)+\langle\mu_{g}^{*}(x),g^{c}(x,y)\rangle.$

As $g(x,y)$ is $\alpha_{g}$-strongly convex in $y$ and $g^{c}(x,y)$ is convex
in $y$, we know $g(x,y)+\langle\mu_{g}^{*}(x),g^{c}(x,y)\rangle$ is
$\alpha_{g}$-strongly convex in $y$, where the modulus $\alpha_{g}$ is
independent of $x$. Therefore, according to Theorem 2 in [37], the quadratic
growth in statement 1 can be concluded.

As $g(x,y)$ is strongly convex in $y$, and $\mathcal{Y}(x)$ is a non-empty,
closed and convex set under assumption 3, there exists a unique solution
$y^{*}_{g}(x)$ such that $g(x,y_{g}^{*}(x))=v(x)$ by Lemma 5. In this way, if
$y\neq y^{*}_{g}(x)$ and $y\in\mathcal{Y}(x)$, we have $g(x,y)>v(x)$. This
completes the proof of statement 2.

####  B.2 Proof of Theorem 1

We know from Lemma 1 that
$g(x,y)-v(x)\geq\frac{\alpha_{g}}{2}\|y-y^{*}_{g}(x)\|^{2}$ and $g(x,y)=v(x)$
if and only if $y=y^{*}_{g}(x)$. This is squared-distance bound follows
Definition 1 in [61]. Under Lipschitzness of $f(x,y)$ with respect to $y$,
finding the solutions to the $\epsilon$-approximate problem in (4) is
equivalent to finding the solutions to its penalty reformulation

$\displaystyle\min_{(x,y)\in\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq
0\\}}f(x,y)+\gamma(g(x,y)-v(x))$ (28)

with $\gamma={\cal O}(\epsilon^{-0.5})$ following Theorems 1 and 2 in [61].

Moreover, jointly finding solution for $(x,y)$ in (28) is in equivalence to
finding solutions in

$\displaystyle\min_{x\in\mathcal{X}}\min_{y\in\mathcal{Y}(x)}f(x,y)+\gamma(g(x,y)-v(x)).$
(29)

The proof of this equivalence is as follows. Suppose
$(x_{0},y_{0})\in\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq 0\\}$ being a
solution to (28). Suppose for any $x\in\mathcal{X}$,
$y_{F}^{*}(x)\in\arg\min_{y\in\mathcal{Y}(x)}f(x,y)+\gamma(g(x,y)-v(x))$. We
know that for any $x\in\mathcal{X}$, $y\in\mathcal{Y}(x)$, it follows that

$\displaystyle f(x_{0},y_{0})+\gamma(g(x_{0},y_{0})-v(x_{0}))\leq$
$\displaystyle f(x,y_{F}^{*}(x))+\gamma(g(x,y_{F}^{*}(x))-v(x))$
$\displaystyle\leq$ $\displaystyle f(x,y)+\gamma(g(x,y)-v(x)).$

This means any solution to (28) is a solution to (29). On the other hand,
suppose $x_{0}\in\mathcal{X}$, $y_{F}^{*}(x_{0})\in\mathcal{Y}(x_{0})$ is a
solution to (29). We know that for any
$(x,y)\in\\{\mathcal{X}\times\mathcal{Y}:g^{c}(x,y)\leq 0\\}$,

$\displaystyle
f(x_{0},y_{F}^{*}(x_{0}))+\gamma(g(x_{0},y_{F}^{*}(x_{0}))-v(x_{0}))\leq$
$\displaystyle f(x,y_{F}^{*}(x))+\gamma(g(x,y_{F}^{*}(x))-v(x))$
$\displaystyle\leq f(x,y)+\gamma(g(x,y)-v(x)).$

This means any solution to (29) is a solution to (28).

Besides, we know $f(x,y)$ is $l_{f,1}$-smooth, $g(x,y)$ is
$\alpha_{g}$-strongly convex in $y$. By the definitions of strongly convexity
and smoothness, we know for fixed $x$, for any $y_{1},y_{2}\in\mathcal{Y}$,

$\displaystyle
f(x,y_{1})+\gamma(g(x,y_{1})-v(x))-f(x,y_{2})+\gamma(g(x,y_{2})-v(x))$
$\displaystyle=$ $\displaystyle
f(x,y_{1})-f(x,y_{2})+\gamma(g(x,y_{1})-g(x,y_{2}))$ $\displaystyle\geq$
$\displaystyle\langle\nabla_{y}f(x,y_{2}),y_{1}-y_{2}\rangle-\frac{l_{f,1}}{2}\|y_{1}-y_{2}\|^{2}+\gamma\langle\nabla_{y}g(x,y_{2}),y_{1}-y_{2}\rangle+\gamma\frac{\alpha_{g}}{2}\|y_{1}-y_{2}\|^{2}$
$\displaystyle=$
$\displaystyle\langle\nabla_{y}f(x,y_{2})+\gamma\nabla_{y}g(x,y_{2}),y_{1}-y_{2}\rangle+\frac{\gamma\alpha_{g}-l_{f,1}}{2}\|y_{1}-y_{2}\|^{2}.$
(30)

This proves that $f(x,y)+\gamma(g(x,y)-v(x))$ is
$(\gamma\alpha_{g}-l_{f,1})$-strongly convex in $y$. Moreover, according to
Assumption 2, the constraint $g^{c}(x,y)$ is convex in $y$, and
$\min_{y\in\mathcal{Y}(x)}f(x,y)+\gamma(g(x,y)-v(x))$ is equivalent to its
equivalent Lagrangian Dual Form [58]

$\displaystyle\max_{\mu\in\mathbb{R}_{+}^{d_{c}}}\min_{y\in\mathcal{Y}}f(x,y)+\gamma(g(x,y)-v(x))+\langle\mu,g^{c}(x,y)\rangle.$
(31)

Therefore, (28) can be recovered to (2a) and this completes the proof.

###  Appendix C Analysis of the Differentiability of Value Functions

######  Lemma 8 (Theorem 2.16 in [32]).

Suppose $h(x,y)$ is strongly convex in $y\in\mathcal{Y}$ and is Lipschitz with
respect to $x\in\mathcal{X}$, $h^{c}(x,y)$ is convex in $y$ and is Lipschitz
with respect to $x$, and both $\mathcal{Y}$ and
$\\{y\in\mathcal{Y}:h^{c}(x,y)\leq 0\\}$ are non-empty, closed, and convex.
For the problem $\min_{y\in\\{y\in\mathcal{Y}:h^{c}(x,y)\leq 0\\}}h(x,y)$, the
unique solution $y_{h}^{*}(x)$ and unique Lagrange multiplier
$\mu_{h}^{*}(x)$, defined as

$\displaystyle(y_{h}^{*}(x),\mu_{h}^{*}(x)):=\arg\max_{\mu\in\mathbb{R}^{d_{x}}_{+}}\min_{y\in\mathcal{Y}}h(x,y)+\langle\mu,h^{c}(x,y)\rangle,$
(32)

is Lipschitz in $x$. In other words, there exist $L_{h}\geq 0$ that, for all
$x_{1},x_{2}\in\mathcal{X}$,

$\displaystyle\|(y_{h}^{*}(x_{1});\mu_{h}^{*}(x_{1}))-(y_{h}^{*}(x_{2});\mu_{h}^{*}(x_{2}))\|\leq
L_{h}\|x_{1}-x_{2}\|.$

Before proving Lemmas 2 and 3, we would like to introduce a more general form.

######  Lemma 9.

Suppose $\mathcal{Y}$ and $\\{y\in\mathcal{Y}:h^{c}(x,y)\leq 0\\}$ are both
non-empty, closed and convex, $h(x,y)$ is jointly smooth in $(x,y)$ and is
strongly convex in $y$, $h^{c}(x,y)$ is convex in $y$, and both $h(x,y)$ and
$h^{c}(x,y)$ are Lipschitz with respect to $x$. Then we have

$\displaystyle v_{h}(x)=\min_{y\in\mathcal{Y}}h(x,y)\quad\text{s.t.}\quad
h^{c}(x,y)\leq 0$

is differentiable with its gradient as

$\displaystyle\nabla
v_{h}(x)=\nabla_{x}h(x,y^{*}_{h}(x))+\langle\mu^{*}_{h}(x),h^{c}(x,y^{*}_{h}(x))\rangle,$
(33)

where $(y^{*}_{h}(x),\mu^{*}_{h}(x))$ defined in (32) are unique.

###### Proof.

We prove this using Theorem 4.24 in [7].

i) As $h(x,y)$ being strongly convex in $y$, Condition 1 in Theorem 4.24 in
[7] is satisfied and the solution sets are of singleton value
$(y_{h}^{*}(x),\mu_{h}^{*}(x))$ according to Lemma 5.

ii) Moreover, the smoothness of $h(x,y)$ guarantees Robinson’s constraint
qualification [2], which implies the directional regularity condition
(Definition 4.8 in [7]) for any direction $d$ (Theorem 4.9. (ii) in [7]). This
guarantees Condition 2 in Theorem 4.24 in [7] can be satisfied for all
directions $d$.

iii) Additionally, under the Lipschitzness of $h(x,y)$ and $h^{c}(x,y)$ with
respect to $x$, $y^{*}_{h}(x),\mu^{*}_{h}(x)$ are Lipschitz according to Lemma
8. This implies condition 3 in Theorem 4.24 in [7] holds.

In this way, all conditions in Theorem 4.24 in [7] hold and it gives the
gradient as in (33) for unique $(y_{h}^{*}(x),\mu_{h}^{*}(x))$. This completes
the proof. ∎

####  C.1 Proof of Lemma 2

###### Proof.

The problem $\min_{y\in\mathcal{Y}}g(x,y)\text{ s.t. }g^{c}(x,y)\leq 0$ fits
in the setting of Lemma 9 by taking $h(x,y)=g(x,y)$ and
$h^{c}(x,y)=g^{c}(x,y)$. Therefore the derivative (8) can be obtained
accordingly. Moreover, for any $x_{1},x_{2}\in\mathcal{X}$,

$\displaystyle\|\nabla v(x_{1})-\nabla v(x_{2})\|$ $\displaystyle=$
$\displaystyle\|\nabla_{x}g(x_{1},y_{g}^{*}(x_{1}))+\langle\mu_{g}^{*}(x_{1}),\nabla_{x}g^{c}(x_{1},y_{g}^{*}(x_{1}))\rangle-\nabla_{x}g(x_{2},y_{g}^{*}(x_{2}))$
$\displaystyle-\langle\mu_{g}^{*}(x_{2}),\nabla_{x}g^{c}(x_{2},y_{g}^{*}(x_{2}))\rangle\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}$
$\displaystyle\|\nabla_{x}g(x_{1},y_{g}^{*}(x_{1}))-\nabla_{x}g(x_{2},y_{g}^{*}(x_{2}))\|$
$\displaystyle+\|\langle\mu_{g}^{*}(x_{1}),\nabla_{x}g^{c}(x_{1},y_{g}^{*}(x_{1}))\rangle-\langle\mu_{g}^{*}(x_{1}),\nabla_{x}g^{c}(x_{2},y_{g}^{*}(x_{2}))\rangle\|$
$\displaystyle+\|\langle\mu_{g}^{*}(x_{1}),\nabla_{x}g^{c}(x_{2},y_{g}^{*}(x_{2}))\rangle-\langle\mu_{g}^{*}(x_{2}),\nabla_{x}g^{c}(x_{2},y_{g}^{*}(x_{2}))\rangle\|$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}$
$\displaystyle(l_{g,1}+B_{g}l_{g^{c},1})(\|x_{1}-x_{2}\|+\|y_{g}^{*}(x_{1})-y_{g}^{*}(x_{2})\|)+l_{g^{c},0}\|\mu_{g}^{*}(x_{1})-\mu_{g}^{*}(x_{2})\|$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}$
$\displaystyle((l_{g,1}+B_{g}l_{g^{c},1})(1+L_{g})+l_{g^{c},0}L_{g})\|x_{1}-x_{2}\|,$

where $(a)$ follows triangle inequality; $(b)$ leverage on the Lipschitzness
of $\nabla g$, $g^{c}$ and $\nabla g^{c}$ in $x$, and the upper bound for
$\|\mu_{g}^{*}(x)\|$; and $(c)$ uses the Lipschitzness of $y_{g}^{*}(x)$ and
$\mu_{g}^{*}(x)$ from Lemma 8. As the bound is loose due to the use of
triangle inequality, we can conclude that $v(x)$ is $l_{v,1}$-smooth where
$l_{v,1}\leq((1+B_{g})(1+L_{g})l_{g^{c},1}+l_{g^{c},0}L_{g})$. ∎

####  C.2 Proof of Lemma 3

By assumption, $f(x,y)$ is $l_{f,1}$-smooth and $g(x,y)$ is
$\alpha_{g}$-strongly convex in $y$. We know $f(x,y)+\gamma(g(x,y)-v(x))$ is
$(\gamma\alpha_{g}-l_{f,1})$-strongly convex when
$\gamma>\frac{l_{f,1}}{\alpha_{g}}$ as discussed in (30). Moreover, constraint
$g^{c}(x,y)$ is convex in $y$ by Assumption 2. In this way, the problem

$\displaystyle\min_{y\in\mathcal{Y}}~{}f(x,y)+\gamma(g(x,y)-v(x))~{}~{}~{}~{}~{}~{}\text{s.t.}$
$\displaystyle~{}g^{c}(x,y)\leq 0$

features strong convexity according to Chapter 4 in [58] and it equals to

$\displaystyle
F_{\gamma}(x)=\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}\min_{y\in\mathcal{Y}}$
$\displaystyle~{}f(x,y)+\gamma(g(x,y)-v(x))+\langle\mu,g^{c}(x,y)\rangle.$

Considering the smoothness of $v(x)$ as presented in Lemma 2, all assumptions
in Lemma 9 are satisfied. Therefore the derivative (9) can be obtained. In
addition, for any $x_{1},x_{2}\in\mathcal{X}$,

$\displaystyle\|\nabla F(x_{1})-\nabla F(x_{2})\|$ $\displaystyle=$
$\displaystyle\|\nabla_{x}f(x_{1},y_{F}^{*}(x_{1}))+\gamma\left(\nabla_{x}g(x_{1},y_{F}^{*}(x_{1}))-\nabla
v(x_{1})\right)+\langle\mu_{F}^{*}(x_{1}),\nabla_{x}g^{c}(x_{1},y^{*}_{F}(x_{1}))\rangle$
$\displaystyle-\nabla_{x}f(x_{2},y_{F}^{*}(x_{2}))-\gamma\left(\nabla_{x}g(x_{2},y_{F}^{*}(x_{2}))-\nabla
v(x_{2})\right)-\langle\mu_{F}^{*}(x_{2}),\nabla_{x}g^{c}(x_{2},y^{*}_{F}(x_{2}))\rangle\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}$
$\displaystyle\|\nabla_{x}f(x_{1},y_{F}^{*}(x_{1}))-\nabla_{x}f(x_{2},y_{F}^{*}(x_{2}))\|+\gamma\|\nabla_{x}g(x_{1},y_{F}^{*}(x_{1}))-\nabla_{x}g(x_{2},y_{F}^{*}(x_{2}))\|$
$\displaystyle+\gamma\|\nabla v(x_{1})-\nabla
v(x_{2})\|+\|\langle\mu_{F}^{*}(x_{1}),\nabla_{x}g^{c}(x_{1},y_{F}^{*}(x_{1}))\rangle-\langle\mu_{F}^{*}(x_{1}),\nabla_{x}g^{c}(x_{2},y_{F}^{*}(x_{2}))\rangle\|$
$\displaystyle+\|\langle\mu_{F}^{*}(x_{1}),\nabla_{x}g^{c}(x_{2},y_{F}^{*}(x_{2}))\rangle-\langle\mu_{F}^{*}(x_{2}),\nabla_{x}g^{c}(x_{2},y_{F}^{*}(x_{2}))\rangle\|$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}$
$\displaystyle(l_{f,1}+\gamma
l_{g,1}+B_{F}l_{g^{c},1})(\|x_{1}-x_{2}\|+\|y_{F}^{*}(x_{1})-y_{F}^{*}(x_{2})\|)+\gamma
l_{v,1}\|x_{1}-x_{2}\|$
$\displaystyle+l_{g^{c},0}\|\mu_{F}^{*}(x_{1})-\mu_{F}^{*}(x_{2})\|$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}$
$\displaystyle((l_{f,1}+\gamma l_{g,1}+B_{F}l_{g^{c},1})(1+L_{F})+\gamma
l_{v,1}+l_{f^{c},0}L_{F})\|x_{1}-x_{2}\|,$

where $(a)$ follows triangle inequality; $(b)$ leverage on the Lipschitzness
of $\nabla f$, $\nabla g$, $g^{c}$ and $\nabla g^{c}$ in $x$, and the upper
bound for $\|\mu_{F}^{*}(x)\|$; and $(c)$ uses the Lipschitzness of
$y_{F}^{*}(x)$ and $\mu_{F}^{*}(x)$ from Lemma 8. As the bound is loose due to
the use of triangle equality, we can conclude that $F(x)$ is $l_{F,1}$-smooth
where $l_{F,1}\leq(l_{f,1}+\gamma l_{g,1}+B_{F}l_{g^{c},1})(1+L_{F})+\gamma
l_{v,1}+l_{f^{c},0}L_{F}$.

###  Appendix D Convergence Analysis of the Main Result

####  D.1 Proof of Theorem 2

Define the bias term $b(x_{t})$ of the gradient $\nabla F_{\gamma}(x_{t})$ as

$\displaystyle b(x_{t}):=$ $\displaystyle\nabla F_{\gamma}(x_{t})-g_{F,t}$
$\displaystyle=$
$\displaystyle\bigg{(}\nabla_{x}f(x_{t},y_{F}^{*}(x_{t}))+\gamma\big{(}\nabla_{x}g(x,y_{F}^{*}(x_{t}))-\big{(}\nabla_{x}g(x_{t},y_{g}^{*}(x_{t}))+\left\langle\mu_{g}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{g}^{*}(x_{t}))\right\rangle\big{)}\big{)}$
$\displaystyle~{}~{}~{}~{}+\left\langle\mu_{F}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{F}^{*}(x_{t}))\right\rangle\bigg{)}$
$\displaystyle-\biggl{(}\nabla_{x}f(x_{t},y_{F,t}^{T_{F}})+\gamma\left(\nabla_{x}g(x_{t},y_{F,t}^{T_{F}})-\left(\nabla_{x}g(x_{t},y_{g,t}^{T_{g}})+\left\langle\mu_{g,t}^{T_{g}},\nabla_{x}g^{c}(x_{t},y_{g,t}^{T_{g}})\right\rangle\right)\right)$
$\displaystyle~{}~{}~{}~{}+\langle\mu_{F}^{T_{F}},\nabla_{x}g^{c}(x_{t},y_{F,t}^{T_{F}})\rangle\biggr{)}.$

In this way, we have

$\displaystyle\|b(x_{t})\|\stackrel{{\scriptstyle(a)}}{{\leq}}$
$\displaystyle\|\nabla_{x}f(x_{t},y_{F,t}^{T_{F}})-\nabla_{x}f(x_{t},y_{F}^{*}(x_{t}))\|$
$\displaystyle+\gamma\biggl{(}\|\nabla_{x}g(x_{t},y_{F,t}^{T_{F}})-\nabla_{x}g(x_{t},y_{F}^{*}(x_{t}))\|+\|\nabla_{x}g(x_{t},y_{g,t}^{T_{g}})-\nabla_{x}g(x_{t},y_{g}^{*}(x_{t}))\|$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left\|\left\langle\mu_{g}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{g}^{*}(x_{t}))\right\rangle-\left\langle\mu_{g}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{g,t}^{T_{g}})\right\rangle\right\|$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left\|\left\langle\mu_{g}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{g,t}^{T_{g}})\right\rangle-\left\langle\mu_{g,t}^{T_{g}},\nabla_{x}g^{c}(x_{t},y_{g,t}^{T_{g}})\right\rangle\right\|\biggr{)}$
$\displaystyle+\|\langle\mu_{F}^{T_{F}},\nabla_{x}g^{c}(x_{t},y_{F,t}^{T_{F}})\rangle-\langle\mu_{F}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{F,t}^{T_{F}})\rangle\|$
$\displaystyle+\|\langle\mu_{F}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{F,t}^{T_{F}})\rangle-\langle\mu_{F}^{*}(x_{t}),\nabla_{x}g^{c}(x_{t},y_{F}^{*}(x_{t}))\rangle\|$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}$ $\displaystyle
l_{f,1}\|y_{F,t}^{T_{F}}-y_{F}^{*}(x_{t})\|+\gamma\bigg{(}l_{g,1}\|y_{F,t}^{T_{F}}-y_{F}^{*}(x_{t})\|+l_{g,1}\|y_{g,t}^{T_{g}}-y_{g}^{*}(x_{t})\|$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+l_{g^{c},0}\|\mu_{g,t}^{T_{g}}-\mu_{g}^{*}(x_{t})\|+B_{g}l_{g^{c},1}\|y_{g,t}^{T_{g}}-y_{g}^{*}(x_{t})\|\bigg{)}$
$\displaystyle+l_{g^{c},0}\|\mu_{F,t}^{T_{F}}-\mu_{F}^{*}(x_{t})\|+B_{F}l_{g^{c},1}\|y_{F,t}^{T_{F}}-y_{F}^{*}(x_{t})\|$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}$ $\displaystyle(l_{f,1}+\gamma
l_{g,1}+B_{F}l_{g^{c},0})\|y_{F,t}^{T_{F}}-y_{F}^{*}(x_{t})\|+l_{g^{c},0}\|\mu_{F,t}^{T_{F}}-\mu_{F}^{*}(x_{t})\|$
$\displaystyle+\gamma\left((l_{g,1}+B_{g}l_{g^{c},1})\|y_{g,t}^{T_{g}}-y_{g}^{*}(x_{t})\|+l_{g^{c},0}\|\mu_{g,t}^{T_{g}}-\mu_{g}^{*}(x_{t})\|\right),$

where $(a)$ uses triangle inequality, $(b)$ relies on the Lipschitzness of
$\nabla f$, $\nabla g$, $g^{c}$, and $\nabla g^{c}$ in $x$, the upper bounds
for $\|\mu_{F}^{*}(x)\|$ and $\|\mu_{g}^{*}(x)\|$, and Cauchy-Schwartz
inequality, and $(c)$ is by rearrangement.

Furthermore, according to Young’s inequality, it follows that

$\displaystyle\|b(x_{t})\|^{2}\leq$ $\displaystyle 2\left((l_{f,1}+\gamma
l_{g,1}+B_{F}l_{g^{c},0})\|y_{F,t}^{T_{F}}-y_{F,t}^{*}\|+l_{g^{c},0}\|\mu_{F,t}^{T_{F}}-\mu_{F,t}^{*}\|\right)^{2}$
$\displaystyle+2\gamma^{2}\left((l_{g,1}+B_{g}l_{g^{c},1})\|y_{g,t}^{T_{g}}-y_{g}^{*}(x_{t})\|+l_{g^{c},0}\|\mu_{g,t}^{T_{g}}-\mu_{g}^{*}(x_{t})\|\right)^{2}$
$\displaystyle=$ $\displaystyle{\cal
O}(\gamma^{2}\epsilon_{F}+\gamma^{2}\epsilon_{g}).$

According to Lemma 3, $F_{\gamma}(x)$ is $l_{F,1}$-smooth in $\mathcal{X}$. In
this way, by the smoothness, we have

$\displaystyle F(x_{t+1})\leq$ $\displaystyle F(x_{t})+\langle\nabla
F(x_{t}),x_{t+1}-x_{t}\rangle+\frac{l_{F,1}}{2}\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle\leq$ $\displaystyle F(x_{t})+\langle
g_{F,t},x_{t+1}-x_{t}\rangle+\frac{1}{2\eta}\|x_{t+1}-x_{t}\|^{2}+\langle
b(x_{t}),x_{t+1}-x_{t}\rangle,$ (34)

where the second inequality is by $\eta\leq\frac{1}{l_{F,1}}$ and $\nabla
F(x_{t})=g_{F_{t}}+b(x_{t})$.

The projection guarantees that $x_{t+1}$ and $x_{t}$ are in $\mathcal{X}$.
Following Lemma 6, we know that

$\displaystyle\langle
g_{F,t},x_{t+1}-x_{t}\rangle\leq-\frac{1}{\eta}\|x_{t+1}-x_{t}\|^{2}.$

Plugging this back to (34), it follows

$\displaystyle F(x_{t+1})\leq$ $\displaystyle
F(x_{t})-\frac{1}{2\eta}\|x_{t+1}-x_{t}\|^{2}+\langle
b(x_{t}),x_{t+1}-x_{t}\rangle$ $\displaystyle\leq$ $\displaystyle
F(x_{t})-\frac{1}{2\eta}\|x_{t+1}-x_{t}\|^{2}+\eta\|b(x_{t})\|^{2}+\frac{1}{4\eta}\|x_{t+1}-x_{t}\|^{2}$
$\displaystyle=$ $\displaystyle
F(x_{t})-\frac{1}{4\eta}\|x_{t+1}-x_{t}\|^{2}+\eta\|b(x_{t})\|^{2},$

where the second inequality is from Young’s inequality. Telescoping therefore
gives

$\displaystyle\frac{1}{T}\sum_{t=0}^{T-1}\|G_{\eta}(x_{t})\|^{2}\leq$
$\displaystyle\frac{4}{\eta
T}(F(x_{0})-F(x_{T}))+\frac{4}{T}\sum_{t=0}^{T-1}\|b(x_{t})\|^{2}$
$\displaystyle=$ $\displaystyle{\cal
O}(\eta^{-1}T^{-1})+O(\gamma^{2}\epsilon_{F}+\gamma^{2}\epsilon_{g})$
$\displaystyle=$ $\displaystyle{\cal O}(\gamma
T^{-1}+\gamma^{2}\epsilon_{F}+\gamma^{2}\epsilon_{g})$

where last equality comes from $\eta={\cal O}(\gamma^{-1})$ as
$\eta\leq\frac{1}{l_{F,1}}$ and $l_{F,1}\leq(l_{f,1}+\gamma
l_{g,1}+B_{F}l_{g^{c},1})(1+L_{F})+\gamma l_{v,1}+l_{f^{c},0}L_{F}={\cal
O}(\gamma)$. This completes the proof.

####  D.2 Proof of Theorem 3

To restate, we are viewing $L_{g}(\mu,y;x)$ in (7) and $L_{F}(\mu,y;x)$ in (5)
respectively as $L(\mu,y)$ and considering the following max-min problem

$\displaystyle\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}\min_{y\in\mathcal{Y}}L(\mu,y).$

In (16), we defined the minimization part as

$\displaystyle
D(\mu):=\min_{y\in\mathcal{Y}}L(\mu,y)=L(\mu,y_{\mu}^{*}(\mu))\quad\text{where}\quad
y_{\mu}^{*}(\mu):=\arg\min_{y\in\mathcal{Y}}L(\mu,y).$

To evaluate $D(\mu)$ for $L_{g}(\mu,y;x)$ in (7) and $L_{F}(\mu,y;x)$ in (5)
as $L(\mu,y)$, we define the following mappings for fixed $x\in\mathcal{X}$:

$\displaystyle y^{*}_{\mu,g}(\mu):=$
$\displaystyle\arg\min_{y\in\mathcal{Y}}L_{g}(\mu,y;x),$ (35) $\displaystyle
y^{*}_{\mu,F}(\mu):=$ $\displaystyle\arg\min_{y\in\mathcal{Y}}L_{F}(\mu,y;x).$
(36)

In this way, for $L_{g}(\mu,y;x)$ in (7) and $L_{F}(\mu,y;x)$ in (5)
respectively as $L(\mu,y)$, $D(\mu)$ defined in (16) equals to $D_{g}(\mu)$
and $D_{F}(\mu)$ respectively, where

$\displaystyle D_{g}(\mu)=$
$\displaystyle\min_{y\in\mathcal{Y}}L_{g}(\mu,y;x)=L_{g}(\mu,y^{*}_{\mu,g}(\mu);x),$
(37) $\displaystyle D_{F}(\mu)=$
$\displaystyle\min_{y\in\mathcal{Y}}L_{F}(\mu,y;x)=L_{F}(\mu,y^{*}_{\mu,F}(\mu);x).$
(38)

In the following lemma, we show that $D(\mu)$ exhibits concavity and
smoothness, which are favorable properties for conducting gradient-based
algorithm.

######  Lemma 10 (Smoothness and Concavity of $D(\mu)$).

Suppose all the assumptions in Theorem 3 hold. For fixed $x\in\mathcal{X}$,
the following holds

  1. 1.

$y^{*}_{\mu,g}(\mu)$ in (35) and $y^{*}_{\mu,F}(\mu)$ in (36) are respectively
$\frac{1}{\alpha_{g}}$ and $\frac{1}{\gamma\alpha_{g}-l_{f,1}}$-Lipschitz to
$\mu$.

  2. 2.

$D_{g}(\mu)$ in (37) is concave and $\frac{l_{g^{c},0}}{\alpha_{g}}$-smooth,
and $D_{F}(\mu)$ in (38) is concave and
$\frac{l_{g^{c},0}}{\gamma\alpha_{g}-l_{f,1}}$-smooth.

###### Proof.

To restate,

$\displaystyle L_{g}(\mu,y;x)=$ $\displaystyle
g(x,y)+\langle\mu,g^{c}(x,y)\rangle,$ $\displaystyle L_{F}(\mu,y;x)=$
$\displaystyle f(x,y)+\gamma(g(x,y)-v(x))+\langle\mu,g^{c}(x,y)\rangle.$

Under Assumption 1 and 2, for fixed $x$ and given $\mu$, $L_{g}(\mu,y;x)$ is
$\alpha_{g}$-strongly convex and $(l_{g,1}+\|\mu\|l_{g^{c},1})$-smooth in $y$,
and $L_{F}(\mu,y;x)$ is $(\gamma\alpha_{g}-l_{f,1})$-strongly convex and
$(l_{f,1}+\gamma l_{g,1}+\|\mu\|l_{g^{c},1})$-smooth in $y$. Therefore, we
know $y^{*}_{\mu,g}(\mu)$ in (35) and $y^{*}_{\mu,F}(\mu)$ in (36) are
respectively $\frac{1}{\alpha_{g}}$ and
$\frac{1}{\gamma\alpha_{g}-l_{f,1}}$-Lipschitz to $\mu$ by directly quoting
Theorem F.10 in [18] or Theorem 4.47 in [32]. This proves the first part of
the Lemma.

For the second part, the concavity of $D_{g}(\mu)$ and $D_{F}(\mu)$ can be
directly obtained by Lemma 2.58 in [58] as $L_{g}(\mu,y;x)$ and
$L_{g}(\mu,y;x)$ are both convex in $y$.

Moreover, following Theorem 4.24 in [7], we have

$\displaystyle\nabla D_{g}(\mu)=$
$\displaystyle\nabla_{\mu}L_{g}(\mu,y^{*}_{\mu,g}(\mu))=g^{c}(x,y^{*}_{\mu,g}(\mu)),$
$\displaystyle\nabla D_{F}(\mu)=$
$\displaystyle\nabla_{\mu}L_{F}(\mu,y^{*}_{\mu,F}(\mu))=g^{c}(x,y^{*}_{\mu,F}(\mu)).$

As $g^{c}(x,y)$ is $l_{g^{c},0}$-Lipschitz by Assumption 1, for any
$\mu_{1},\mu_{2}\in\mathbb{R}^{d_{c}}_{+}$:

$\displaystyle\|\nabla D_{g}(\mu_{1})-\nabla D_{g}(\mu_{2})\|=$
$\displaystyle\|g^{c}(x,y^{*}_{\mu,g}(\mu_{1}))-g^{c}(x,y^{*}_{\mu,g}(\mu_{2}))\|$
$\displaystyle\leq$ $\displaystyle
l_{g^{c},0}\|y^{*}_{\mu,g}(\mu_{1})-y^{*}_{\mu,g}(\mu_{2})\|\leq\frac{l_{g^{c},0}}{\alpha_{g}}\|\mu_{1}-\mu_{2}\|,$
(39)

and similarly,

$\displaystyle\|\nabla D_{F}(\mu_{1})-\nabla D_{F}(\mu_{2})\|=$
$\displaystyle\|g^{c}(x,y^{*}_{\mu,F}(\mu_{1}))-g^{c}(x,y^{*}_{\mu,F}(\mu_{2}))\|$
$\displaystyle\leq$ $\displaystyle
l_{g^{c},0}\|y^{*}_{\mu,F}(\mu_{1})-y^{*}_{\mu,F}(\mu_{2})\|\leq\frac{l_{g^{c},0}}{\gamma\alpha_{g}-l_{f,1}}\|\mu_{1}-\mu_{2}\|.$
(40)

We can conclude that $D_{g}(\mu)$ and $D_{F}(\mu)$ are respectively
$\frac{l_{g^{c},0}}{\alpha_{g}}$ and
$\frac{l_{g^{c},0}}{\gamma\alpha_{g}-l_{f,1}}$-smooth. ∎

In Algorithm 2, we are implementing an accelerated projected gradient descent
on $-D(\mu)$ where the gradient bias is controlled by $T_{y}$. Following [17],
the following lemma presents the convergence analysis of the accelerated
method on smooth and convex functions.

######  Lemma 11 (Section 5 and 6 in [17]).

Suppose $D(\mu)$ is concave and $l_{D,1}$-smooth. Consider the constrained
problem $\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D(\mu)$. At iteration
$t=0,\dots,T-1$, perform accelerated projected gradient update with step size
$\eta\leq\frac{1}{l_{D,1}}$ and initial value $\mu_{0}=\mu_{-1}$:

$\displaystyle\mu_{t+\frac{1}{2}}=$
$\displaystyle\mu_{t}+\frac{t-1}{t+2}(\mu_{t}-\mu_{t-1})$ (41a)
$\displaystyle\mu_{t+1}=$
$\displaystyle\operatorname{Proj}_{\mathbb{R}^{d_{c}}_{+}}(\mu_{t+\frac{1}{2}}+\eta
g_{t+\frac{1}{2}})$ (41b)

where $g_{t+\frac{1}{2}}$ is an $\epsilon^{1.5}$-approximate to $\nabla
D(\mu_{t+\frac{1}{2}})$ satisfying $\|g_{t+\frac{1}{2}}\ -\nabla
D(\mu_{t+\frac{1}{2}})\|={\cal O}(\epsilon^{1.5})$, for a given accuracy
$\epsilon>0$. Denote $D^{*}=\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D(\mu)$,
performing $T={\cal O}(\epsilon^{-0.5})$ iterations leads to

$\displaystyle D^{*}-D(\mu_{T})={\cal O}(\epsilon).$

In this way, we are ready to proceed to the proof of Theorem 3.

###### proof of Theorem 3.

Algorithm 2 solves (7) and (5) by taking $L_{g}(\mu,y;x)$ and $L_{F}(\mu,y;x)$
respectively as $L(\mu,y)$ and run respectively iterations $T$ equals $T_{g}$
and $T_{F}$

We begin our proof with analysis on (7). The accelerated version of Algorithm
2 solves this by taking $L_{g}(\mu,y;x)$ in (7) with fixed $x\in\mathcal{X}$
as $L(\mu,y)$.

Fixing $\mu_{t+\frac{1}{2}}$, steps 4-6 are $T_{y}$-step projected gradient
descent in $y$. As $L_{g}(\mu,y;x)$ is $(l_{g,1}+l_{g^{c},1})$-smooth and
$\alpha_{g}$-strongly convex in $y$, taking the inner loop step size
$\eta_{1}$ as $\eta_{g,1}\leq\frac{1}{l_{g,1}+l_{g^{c},1}}$ ensures linear
convergence according to Lemma 7. Choosing $T_{y}={\cal
O}(\ln((\epsilon^{1.5})^{-1}))={\cal O}(\ln(\epsilon_{g}^{-1}))$ leads to

$\displaystyle\|y_{t+1}-y^{*}_{\mu,g}(\mu_{t+\frac{1}{2}})\|={\cal
O}(\epsilon_{g}^{1.5})$

for target accuracy $\epsilon_{g}>0$ where $y^{*}_{\mu,g}(\mu)$ is defined in
(35).

In step 7, we update $\mu$ with a projected gradient descent step which takes
$\nabla_{\mu}L(\mu_{t+\frac{1}{2}},y_{t+1})=g^{c}(x,y_{t+1})$ as an estimate
of $\nabla D_{g}(\mu_{t+\frac{1}{2}})=g^{c}(x,y^{*}_{\mu,g})$. The estimation
bias is bounded by

$\displaystyle\|\nabla_{\mu}L(\mu_{t+\frac{1}{2}},y_{t+1})-\nabla
D_{g}(\mu_{t+\frac{1}{2}})\|=$
$\displaystyle\|g^{c}(x,y_{t+1})-g^{c}(x,y^{*}_{\mu,g}(\mu_{t+\frac{1}{2}}))\|$
$\displaystyle\leq$ $\displaystyle
l_{g^{c},0}\|y_{t+1}-y^{*}_{\mu,g}(\mu_{t+\frac{1}{2}})\|={\cal
O}(\epsilon_{g}^{1.5}).$

By Lemma 11, we can conclude the complexity is $\tilde{\cal
O}(\epsilon_{g}^{-0.5})$ for conducting the accelerated version of Algorithm 2
on $L_{g}(\mu,y;x)$ in (7) as $L(\mu,y)$ to achieve

$\displaystyle D_{g}(\mu^{*}_{g}(x))-D_{g}(\mu_{T_{g}})={\cal
O}(\epsilon_{g}).$ (42)

Similarly, the complexity of the accelerated version of Algorithm 2 for (5) is
$\tilde{\cal O}(\epsilon_{F}^{-0.5})$ , i.e.,

$\displaystyle D_{F}(\mu^{*}_{F}(x))-D_{F}(\mu_{T_{F}})={\cal
O}(\epsilon_{F}).$ (43)

In the following, we are going to show that (42) and (43) and respectively
sufficient to bound $\|\mu_{T_{g}}-\mu^{*}_{g}(x)\|$ and
$\|\mu_{T_{F}}-\mu^{*}_{F}(x)\|$ considering Assumption 5 is satisfied.

As $D_{g}(\mu)$ and $D_{F}(\mu)$ are both concave in $\mu$ and
$\mu\in\mathbb{R}^{d_{c}}_{+}$ is equivalent to $\mu\geq 0$, the problems

$\displaystyle\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D_{g}(\mu)\quad\text{and}\quad\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D_{F}(\mu)$

are respectively equivalent to the unconstrained problems

$\displaystyle\max_{\mu\in\mathbb{R}^{d_{c}}}\tilde{D}_{g}(\mu):=D_{g}(\mu)+\lambda_{g}^{\top}\mu\quad\text{and}\quad\max_{\mu\in\mathbb{R}^{d_{c}}}\tilde{D}_{F}(\mu):=D_{F}(\mu)+\lambda_{F}^{\top}\mu$

with the Lagrange multipliers $\lambda_{g},\lambda_{F}$ being non-negative and
finite in all dimension, i.e. $0\leq\lambda_{g}<\infty$,
$0\leq\lambda_{F}<\infty$ according to Lagrange duality Theorem. Moreover, it
is well known (Chapter 4 in [58]) that the Largrangian terms respectively
equal zero when the problems attain the optimals. i.e.

$\displaystyle\lambda_{g}^{\top}\mu^{*}_{g}(x)=0\quad\text{and}\quad\lambda_{F}^{\top}\mu^{*}_{F}(x)=0.$
(44)

Moreoever, the first-order stationary condition requires
$\nabla\tilde{D}_{g}(\mu_{g}^{*}(x))=\nabla
D_{g}(\mu_{g}^{*}(x))+\lambda_{g}=0$ and
$\nabla\tilde{D}_{F}(\mu_{F}^{*}(x))=\nabla
D_{F}(\mu_{F}^{*}(x))+\lambda_{F}=0$ and therefore

$\displaystyle\nabla
D_{g}(\mu_{g}^{*}(x))=-\lambda_{g}\quad\text{and}\quad\nabla
D_{F}(\mu_{F}^{*}(x))=-\lambda_{F}.$ (45)

In this way, for all
$\mu\in\mathcal{B}(\mu_{g}^{*}(x);\delta_{g})\cap\mathbb{R}^{d_{c}}_{+}$.

$\displaystyle D_{g}(\mu_{g}^{*}(x))-D_{g}(\mu)=$
$\displaystyle\int_{\tau=0}^{1}\langle\nabla
D_{g}(\mu+\tau(\mu_{g}^{*}(x)-\mu)),\mu_{g}^{*}(x)-\mu\rangle d\tau$
$\displaystyle=$ $\displaystyle\int_{\tau=0}^{1}\frac{1}{\tau}\langle\nabla
D_{g}(\mu_{g}^{*}(x))-D_{g}(\mu+\tau(\mu_{g}^{*}(x)-\mu)),\tau(\mu-\mu_{g}^{*}(x))\rangle
d\tau$ $\displaystyle-\langle\nabla
D_{g}(\mu_{g}^{*}(x)),\mu-\mu_{g}^{*}(x)\rangle$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}$
$\displaystyle\int_{0}^{1}C_{\delta_{g}}\|\mu-\mu_{g}^{*}(x)\|^{2}\tau
d\tau-\langle\nabla D_{g}(\mu_{g}^{*}(x)),\mu-\mu_{g}^{*}(x)\rangle$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}$
$\displaystyle\frac{C_{\delta_{g}}}{2}\|\mu-\mu_{g}^{*}(x)\|^{2}+\langle\lambda_{g},\mu-\mu_{g}^{*}(x)\rangle$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\geq}}$
$\displaystyle\frac{C_{\delta_{g}}}{2}\|\mu-\mu_{g}^{*}(x)\|^{2},$

where $(a)$ uses (21a) in Assumption 5 and the fact that the
$\mu,\mu_{g}^{*}(x)\in\mathcal{B}(\mu_{g}^{*}(x);\delta_{g})\cap\mathbb{R}^{d_{c}}_{+}$
implies
$\mu+\tau(\mu_{g}^{*}(x)-\mu)\in\mathcal{B}(\mu_{g}^{*}(x);\delta_{g})\cap\mathbb{R}^{d_{c}}_{+}$;
$(b)$ solves the integral and uses $\lambda_{g}=-\nabla D_{g}(\mu_{g}^{*}(x))$
in (45); and $(c)$ follows from the fact that
$\langle\lambda,\mu_{g}^{*}(x)\rangle=0$ in (44) and $\mu,\lambda_{g}\geq 0$.

Analogously, for all
$\mu\in\mathcal{B}(\mu_{F}^{*}(x);\delta_{F})\cap\mathbb{R}^{d_{c}}_{+}$, it
follows that

$\displaystyle D_{F}(\mu_{F}^{*}(x))-D_{F}(\mu)\geq$
$\displaystyle\frac{C_{\delta_{F}}}{2}\|\mu-\mu_{F}^{*}(x)\|^{2}.$

In this way, for arbitrary $\epsilon_{g}<\frac{C_{\delta_{g}}}{2}\delta_{g}$,
the complexity of Algorithm 2 to solve (7) is $\tilde{\cal
O}(\epsilon_{F}^{-0.5})$, i.e.,

$\displaystyle\|\mu_{T_{g}}-\mu_{g}^{*}(x)\|^{2}=$ $\displaystyle~{}{\cal
O}(\epsilon_{g}),$ $\displaystyle\text{and
}\|y_{T_{g}}-y_{g}^{*}(x)\|^{2}\leq$
$\displaystyle\|y_{T_{g}}-y_{g}^{*}(\mu_{T_{g}};x)\|^{2}+\|\mu_{T_{g}}-\mu_{g}^{*}(x)\|^{2}$
$\displaystyle\leq$
$\displaystyle(1/\alpha_{g}+1)\|\mu_{T_{g}}-\mu_{g}^{*}(x)\|^{2}={\cal
O}(\epsilon_{g}).$

At each iteration $t$ in Algorithm 1, $x=x_{t}$, and the output
$(y_{T_{g}},\mu_{T_{g}})$ is chosen as $(y_{g,t}^{T_{g}},\mu_{g,t}^{T_{g}})$.

Similarly, to solve (5), for arbitrary
$\epsilon_{g}<\frac{C_{\delta_{g}}}{2}\delta_{g}$, applying Algorithm 2 with
complexity $\tilde{\cal O}(\epsilon_{F}^{-0.5})$ to achieve (43) can achieve

$\displaystyle\|\mu_{T_{F}}-\mu_{F}^{*}(x)\|^{2}=$ $\displaystyle~{}{\cal
O}(\epsilon_{F}),$ $\displaystyle\text{and
}\|y_{T_{F}}-y_{F}^{*}(x)\|^{2}\leq$
$\displaystyle\|y_{T_{F}}-y_{F}^{*}(\mu_{T_{F}};x)\|^{2}+\|\mu_{T_{F}}-\mu_{F}^{*}(x)\|^{2}$
$\displaystyle\leq$
$\displaystyle(1/\alpha_{F}+1)\|\mu_{T_{F}}-\mu_{F}^{*}(x)\|^{2}={\cal
O}(\epsilon_{F}).$

At each iteration $t$ in Algorithm 1, $x=x_{t}$, and the output
$(y_{T_{F}},\mu_{T_{F}})$ is chosen as $(y_{F,t}^{T_{F}},\mu_{F,t}^{T_{F}})$.

This completes the proof. ∎

####  D.3 Proof of Theorem 4

In this section, we consider

$\displaystyle g^{c}(x,y)=g^{c}_{1}(x)^{\top}y-g^{c}_{2}(x)$ (46)

being affine in $y$, and $\mathcal{Y}=\mathbb{R}^{d_{y}}$.

Therefore, for a fixed $x$, taking either $L_{g}(\mu,y;x)$ in (7) or
$L_{F}(\mu,y;x)$ in (5) as $L(\mu,y)$ fits into a special case of strongly-
convex-concave saddle point problems in the following form:

$\displaystyle\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}\min_{y\in\mathbb{R}^{d_{y}}}L(\mu,y)=-h_{1}(\mu)+y^{\top}A\mu+h_{2}(y)$
(47)

where $h_{1}(\mu)$ is smooth and linear (concave) in $\mu$ and $h_{2}(y)$ is
smooth and strongly convex in $y$. Specifically, for $L_{g}(\mu,y;x)$ as
$L(\mu,y)$, $h_{1}(\mu)=\langle g^{c}_{2}(x),\mu\rangle$ is $0$-smooth and
linear (concave), $A=g^{c}_{1}(x)$, and $h_{2}(y)=g(x,y)$ is $l_{g,1}$-smooth
and $\alpha_{g}$-strongly convex.

In this context, performing PGD on $L(\mu,y)$ on $y$ is equivalent to a
gradient descent step since $\mathcal{Y}=\mathbb{R}^{d_{y}}$. The following
lemma summarizes the error of $\|y_{t}-y_{\mu}^{*}(\mu_{t})\|$ where
$y_{\mu}^{*}(\mu)$ is defined in (16) and the update of $\mu_{t+1}$ and
$y_{t+1}$ is the non-accelerated version of Algorithm 2 with $T_{y}=1$.

######  Lemma 12 (Update of $\|y_{t}-\nabla h_{2}^{*}(-A\mu_{t})\|$).

Consider the problem (47) where $A$ is of full column rank and $h_{2}(y)$ is
$\alpha_{h_{2}}$-strongly convex and $l_{f,1}$-smooth. ${y}^{*}_{\mu}(\mu)$
defined in (16) satisfies

$\displaystyle{y}^{*}_{\mu}(\mu)=\nabla h_{2}^{*}(-A\mu)$ (48)

where $h_{2}^{*}(y)$ is the conjugate function of $h_{2}(y)$ by Definition 5.
At iteration $t$, $\mu_{t},y_{t}$ are known, and conduct
$y_{t+1}=y_{t}-\eta_{1}\nabla_{y}L(\mu_{t},y_{t})$, a gradient descent step
for $L(\mu_{t},y)$ in $y$. This gives

$\displaystyle\|y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t})\|\leq(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|$ (49)

when $\eta_{1}\leq\frac{1}{l_{h_{2},1}}$. Additionally, given another
$\mu_{t+1}$, we know

$\displaystyle\|y_{t+1}-\nabla h_{2}^{*}(-A\mu_{t+1})\|\leq$
$\displaystyle(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\|\mu_{t+1}-\mu_{t}\|.$
(50)

###### Proof.

Recall the definition ${y}^{*}_{\mu}(\mu)=\arg\min_{y}L(\mu,y)$ following
(16). The first-order stationary optimality condition requires that for any
given $\mu$, it holds

$\displaystyle\nabla_{y}L(\mu,{y}^{*}_{\mu})=A\mu+\nabla{h}_{2}({y}^{*}_{\mu})=0\quad\Leftrightarrow\quad\nabla{h}_{2}({y}^{*}_{\mu})=-A\mu.$

As the mapping $\nabla h_{2}$ and $\nabla h_{2}^{*}$ are the inverse of each
other according to Lemma 4, for any $\mu$:

$\displaystyle{y}^{*}_{\mu}=\nabla h_{2}^{*}(-A\mu).$

At iteration $t$, conducting a gradient descent step on $L(\mu_{t},y)$ gives
$y_{t+1}=y_{t}-\eta_{1}\nabla_{y}L(\mu_{t},y_{t})$. As $L(\mu_{t},y)$ is
$\alpha_{h_{2}}$-strongly convex and $l_{h_{2},1}$-smooth in $y$, following
Lemma 7, take $\eta_{1}\leq\frac{1}{l_{h_{2},1}}$, we have

$\displaystyle\|y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t})\|\leq(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|.$

Following triangle inequality, we also have

$\displaystyle\|y_{t+1}-\nabla h_{2}^{*}(-A\mu_{t+1})\|$ $\displaystyle\leq$
$\displaystyle\|y_{t+1}-\nabla h_{2}^{*}(-A\mu_{t})\|+\|\nabla
h_{2}^{*}(-A\mu_{t})-\nabla h_{2}^{*}(-A\mu_{t+1})\|$ $\displaystyle\leq$
$\displaystyle(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\|\mu_{t+1}-\mu_{t}\|$
(51)

where the second term comes from the smoothness of the conjugate function (see
Lemma 4). ∎

In (50), the update behavior of $\|y_{t}-\nabla h_{2}^{*}(-A\mu_{t})\|$
depends on $\|\mu_{t+1}-\mu_{t}\|$. Therefore, we are interested in the the
update behavior of $\|\mu_{t+1}-\mu_{t}\|$ where
$\mu_{t+1}=\operatorname{Proj}_{\mathbb{R}^{d_{c}}}(\mu_{t}+\eta\nabla_{\mu}L(\mu_{t},y_{t+1}))$
as in the non-accelerated version of Algorithm 2 with $T_{y}=1$. Before
proceeding, we would like to look into the properties of $D(\mu)$ defined in
(16) under the setting (47):

$\displaystyle D(\mu)=\min_{y}-h_{1}(\mu)+\langle y,A\mu\rangle+h_{2}(y)$ (52)

where $h_{1}(\mu)$ is smooth and linear (concave) in $\mu$, $h_{2}(y)$ is
smooth and strongly convex in $y$, and $A$ is of full rank in column. The next
lemma shows that $D(\mu)$ features strong concavity and smoothness.

######  Lemma 13 (Smoothness and strongly concavity of $D(\mu)$).

Suppose $h_{1}$ is concave and $l_{h_{1},1}$-smooth, $h_{2}$ is
$\alpha_{h_{2}}$-strongly convex and $l_{h_{2},1}$-smooth, and $A$ is full
column rank. Then $D(\mu)$ in (52) satisfies

$\displaystyle D(\mu)=-h_{1}(\mu)-h_{2}^{*}(-A\mu),$

and is $\frac{\sigma^{2}_{\min}(A)}{l_{h_{2},1}}$-strongly concave and
$(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})$-smooth with
respect to $\mu$.

###### Proof.

Following Definition 5, we have

$\displaystyle D(\mu)=-h_{1}(\mu)-h_{2}^{*}(-A\mu)$

where $h_{2}^{*}(y)$ is $\frac{1}{l_{h_{2},1}}$-strongly convex and
$\frac{1}{\alpha_{h_{2}}}$-smooth according to Lemma 4.

For all $\mu_{1},\mu_{2}$,

$\displaystyle-D(\mu_{1})-(-D(\mu_{2}))=$ $\displaystyle
h^{*}_{2}(-A\mu_{1})-h^{*}_{2}(-A\mu_{2})+h_{1}(\mu_{1})-h_{1}(\mu_{2})$
$\displaystyle\geq$ $\displaystyle\langle\frac{\partial
h_{2}^{*}(-A\mu_{2})}{\partial-A\mu_{2}},-A\mu_{1}+A\mu_{2}\rangle+\frac{1/{l_{h_{2},1}}}{2}\|A\mu_{1}-A\mu_{2}\|^{2}$
$\displaystyle+\langle\nabla h_{1}(\mu_{2}),\mu_{1}-\mu_{2}\rangle\rangle$
$\displaystyle\geq$ $\displaystyle\langle\nabla
D(\mu_{2}),\mu_{1}-\mu_{2}\rangle+\frac{\sigma_{\min}^{2}(A)/l_{h_{2},1}}{2}\|\mu_{1}-\mu_{2}\|^{2}.$

where the first inequality follows the strong convexity of $h_{2}^{*}(y)$ and
the fact that $-h_{1}(\mu)$ is convex as $h_{1}(y)$ is concave. and the second
inequality follows the chain rule to formulate $\nabla D(\mu_{2})$. Therefore,
$-D(\mu)$ is $\frac{\sigma^{2}_{\min}(A)}{l_{h_{2},1}}$-strongly convex, and
$D(\mu)$ is $\frac{\sigma^{2}_{\min}(A)}{l_{h_{2},1}}$-strongly concave.

Moreover $D(\mu)$ is
$(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})$-smooth as

$\displaystyle D(\mu_{1})-D(\mu_{2})=$
$\displaystyle-h^{*}_{2}(-A\mu_{1})-(-h^{*}_{2}(-A\mu_{2}))-h_{1}(\mu_{1})+h_{1}(\mu_{2})$
$\displaystyle\leq$ $\displaystyle\langle\frac{\partial-
h_{2}^{*}(-A\mu_{2})}{\partial-A\mu_{2}},-A\mu_{1}-(-A\mu_{2})\rangle+\frac{1/{\alpha_{h_{2}}}}{2}\|-A\mu_{1}-(-A\mu_{2})\|^{2}$
$\displaystyle+\langle-\nabla
h_{1}(\mu_{2}),\mu_{1}-\mu_{2}\rangle\rangle+\frac{l_{h_{1},1}}{2}\|\mu_{1}-\mu_{2}\|^{2}$
$\displaystyle\leq$ $\displaystyle\langle\nabla
D(\mu_{2}),\mu_{1}-\mu_{2}\rangle+\frac{l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}}{2}\|\mu_{1}-\mu_{2}\|^{2}.$

The first inequality holds as $h_{2}^{*}(y)$ and $h_{1}(\mu)$ are smooth. The
second follows the chain rule.

Note $\sigma_{\max}(A)\geq\sigma_{\min}(A)>0$ as $A$ is full column rank. This
completes the proof. ∎

Knowing $D(\mu)$ has such favorable properties, we next analyze the update of
$\|\mu_{t+1}-\mu_{t}\|$, where $\\{\mu_{t}\\}$ is the sequence generated in
the non-accelerated version of Algorithm 2 with $T_{y}=1$.

######  Lemma 14 (Update of $\|\mu_{t+1}-\mu_{t}\|$).

Consider the problem in (47) where $h_{1}$ is concave and
$l_{h_{1},1}$-smooth, $h_{2}$ is $\alpha_{h_{2}}$-strongly convex and
$l_{h_{2},1}$-smooth, and $A$ is full column rank. Running the non-accelerated
version of Algorithm 2 with $T_{y}=1$ and $\eta_{1}\leq{l_{h_{2},1}}^{-1}$
gives

$\displaystyle\frac{1}{\eta_{2}}\|\mu_{t+1}-\mu_{t}\|\leq$
$\displaystyle\left(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\right)\|\mu_{t}-\mu^{*}\|$
$\displaystyle+\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\|\lambda\|$ (53)

where the constant $\lambda$ satisfies $0\leq\lambda<\infty$ and
$\mu^{*}=\arg\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D(\mu)$ with $D(\mu)$ defined
in (52).

###### Proof.

According to Lemma 13,

$\displaystyle
D(\mu)=\min_{y\in\mathbb{R}^{d_{y}}}-h_{1}(\mu)+y^{\top}A\mu+h_{2}(y)=-h_{1}(\mu)-h_{2}^{*}(-A\mu)$

is $\frac{\sigma^{2}_{\min}(A)}{l_{h_{2},1}}$-strongly concave and
$(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})$-smooth with
respect to $\mu$. Moreover, the problem
$\max_{\mu\in\mathbb{R}^{d_{c}}_{+}}D(\mu)$ is equivalent to the unconstrained
problem with the Lagrange multiplier

$\displaystyle\max_{\mu\in\mathbb{R}^{d_{c}}}\tilde{D}(\mu):=D(\mu)+\lambda^{\top}\mu$

where unique $\lambda$ is non-negative and finite in all dimension, i.e.
$0\leq\lambda<\infty$, as $D(\mu)$ is strongly convex and
$\mu\in\mathbb{R}^{d_{c}}_{+}$ is equivalent to $\mu\geq 0$ satisfying the
LICQ condition (Lemma 5). In this way,

$\displaystyle\nabla\tilde{D}(\mu)=\nabla D(\mu)+\lambda=-\nabla
h_{1}(\mu)+A^{\top}\nabla h_{2}^{*}(-A\mu)+\lambda.$ (54)

We can see that $\tilde{D}(\mu)$ is smooth and strongly concave with the same
modulus as $D(\mu)$. The first-order stationary condition requires

$\displaystyle\nabla\tilde{D}(\mu^{*})=-\nabla h_{1}(\mu^{*})+A^{\top}\nabla
h_{2}^{*}(-A\mu^{*})+\lambda=0.$ (55)

In this way,

$\displaystyle\frac{1}{\eta_{2}}\|\mu_{t+1}-\mu_{t}\|=\frac{1}{\eta_{2}}\|\operatorname{Proj}_{\mathbb{R}^{d_{c}}_{+}}\left(\mu_{t}+\eta_{2}(-\nabla
h_{1}(\mu_{t})+A^{\top}y_{t+1})\right)-\mu_{t}\|$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}$ $\displaystyle\|-\nabla
h_{1}(\mu_{t})+A^{\top}y_{t+1}\|$ $\displaystyle=$ $\displaystyle\|-\nabla
h_{1}(\mu_{t})+A^{\top}\nabla
h_{2}^{*}(-A\mu_{t})+\lambda+A^{\top}y_{t+1}-A^{\top}\nabla
h_{2}^{*}(-A\mu_{t})-\lambda\|$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}$
$\displaystyle\|\nabla\tilde{D}(\mu_{t})\|+\sigma_{\max}(A)\|y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\|\lambda\|$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}$
$\displaystyle\|\nabla\tilde{D}(\mu_{t})-\nabla\tilde{D}(\mu^{*})\|+\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\|\lambda\|$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}$
$\displaystyle\left(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\right)\|\mu_{t}-\mu^{*}\|+\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\|\lambda\|$

Inequality $(a)$ comes from the non-expansiveness (1-Lipschitzness) of the
projection operation, $(b)$ follows triangle inequality and uses (54), $(c)$
uses (55) and (49) in Lemma 12, and $(d)$ comes from the smoothness of
$\tilde{D}(\mu)$, which is of the same modulus as $D(\mu)$. This completes the
proof. ∎

In (53), the update behavior of $\|\mu_{t+1}-\mu_{t}\|$ depends on
$\|\mu_{t}-\mu^{*}\|$. We further look into the update of
$\|\mu_{t}-\mu^{*}\|$ and summarize in the following lemma the bound of the
update of $\|\mu_{t}-\mu^{*}\|$.

######  Lemma 15 (Update of $\|\mu_{t}-\mu^{*}\|$).

Consider the problem in (47) where $h_{1}$ is concave and
$l_{h_{1},1}$-smooth, $h_{2}$ is $\alpha_{h_{2}}$-strongly convex and
$l_{h_{2},1}$-smooth, and $A$ is full column rank. Conduct the non-accelerated
version of Algorithm 2 with $T_{y}=1$, gives

$\displaystyle\|{\mu}_{t+1}-\mu^{*}\|\leq$
$\displaystyle\left(1-\eta_{2}\frac{\sigma^{2}_{\min}(A)}{2l_{h_{2},1}}\right)\|\mu_{t}-\mu^{*}\|+\eta_{2}\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|.$ (56)

when $\eta_{1}\leq{l_{h_{2},1}}^{-1}$ and
$\eta_{2}\leq\left(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\right)^{-1}$.
Here $\mu^{*}=\arg\min_{\mu\in\mathbb{R}^{d_{c}}_{+}}D(\mu)$ where $D(\mu)$ is
defined in (52).

###### Proof.

Define an auxiliary update as

$\displaystyle\tilde{\mu}_{t+1}:=\operatorname{Proj}_{\mathbb{R}^{d_{c}}_{+}}\left(\mu_{t}+\eta_{2}\nabla
D(\mu_{t})\right)=\operatorname{Proj}_{\mathbb{R}^{d_{c}}_{+}}\left(\mu_{t}+\eta_{2}(-\nabla
h_{1}(\mu_{t})+A^{\top}\nabla h_{2}^{*}(-A\mu_{t}))\right).$ (57)

This is a projected gradient descent on strongly convex $-D(\mu)$. As
$\mathbb{R}^{d_{c}}_{+}$ is closed and convex, following Lemma 7, for
$\eta_{2}\leq\left(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\right)^{-1}$,
where $\left(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\right)$
is the modulus for smoothness of $D(\mu)$ by Lemma 13, we have

$\displaystyle\|\tilde{\mu}_{t+1}-\mu^{*}\|\leq$
$\displaystyle\left(1-\eta_{2}\frac{\sigma^{2}_{\min}(A)}{2l_{h_{2},1}}\right)\|\mu_{t}-\mu^{*}\|.$

As the real update is
$\mu_{t+1}=\operatorname{Proj}_{\mathbb{R}^{d_{c}}_{+}}\left((\mu_{t}+\eta_{2}(-\nabla
h_{1}(\mu_{t})+A^{\top}y_{t})\right)$, by the non-expansiveness
(1-Lipschitzness) of projection operation, we have

$\displaystyle\|\tilde{\mu}_{t+1}-\mu_{t+1}\|\leq$
$\displaystyle\|\eta_{2}A^{\top}(y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t}))\|\leq\eta_{2}\sigma_{\max}(A)\|y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t})\|$

By triangle inequality and (49), we have

$\displaystyle\|{\mu}_{t+1}-\mu^{*}\|\leq\left(1-\eta_{2}\frac{\sigma^{2}_{\min}(A)}{2l_{h_{2},1}}\right)\|\mu_{t}-\mu^{*}\|+\eta_{2}\sigma_{\max}(A)\|y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t})\|$ $\displaystyle\leq$
$\displaystyle\left(1-\eta_{2}\frac{\sigma^{2}_{\min}(A)}{2l_{h_{2},1}}\right)\|\mu_{t}-\mu^{*}\|+\eta_{2}\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|.$ (58)

This completes the proof. ∎

We are ready to proceed with the convergence analysis for the single-loop
algorithm (Algorithm 2) without acceleration and $T_{y}=1$, on the problems
(47), which is a general form to (7) and (5). In this way, Theorem 4 follows
directly from the following theorem.

######  Theorem 6.

Suppose $L(\mu,y)$ is in the form of (47) where $A$ is full column rank,
$h_{1}$ is concave and $l_{h_{1},1}$-smooth, $h_{2}$ is
$\alpha_{h_{2}}$-strongly convex and $l_{h_{2},1}$-smooth satisfying
$l_{h_{1},1}={\cal O}(1)$, $l_{h_{2},1},l_{\alpha_{2}}\geq{\cal O}(1)$, and
$\frac{l_{h_{2},1}}{\alpha_{h_{2}}}={\cal O}(1)$. Conduct the non-accelerated
version of Algorithm 2 with $T_{y}=1$. For arbitrary small positive
$\epsilon\leq\left(\frac{4l_{h_{2},1}\sigma_{\max}(A)}{\alpha_{h_{2}}\sigma_{\min}^{2}(A)}(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})\right)^{-1}$,
when $\eta_{1}={\cal O}(\frac{1}{l_{h_{2},1}})\leq\frac{1}{l_{h_{2},1}}$ and
$\eta_{2}={\cal
O}(\epsilon)\leq\frac{1}{l_{h_{1},1}+\sigma_{\max}^{2}(A)/\alpha_{h_{2}}}$,
the algorithm yields output $(\mu_{T},y_{T})$ such that

$\displaystyle\|\mu_{T}-\mu^{*}\|^{2}<\epsilon,\quad\text{and}\quad\|y_{T}-y^{*}\|^{2}<\epsilon$

with complexity $T={\cal O}(\log(\epsilon^{-1}))$. Here,
$(\mu^{*},y^{*})=\arg\max_{\mu\in\mathbb{R}^{d_{c}}}\min_{y\in\mathbb{R}^{d_{y}}}L(\mu,y)$.

###### Proof.

For some positive constant $\rho>0$, denote

$\displaystyle P_{t}:=\rho\|\mu_{t}-\mu^{*}\|+\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|.$ (59)

Plugging (50) in Lemma 12 , (53) in Lemma 14, and (58) in Lemma 15 to (59), we
know

$\displaystyle P_{t+1}=\rho\|\mu_{t+1}-\mu^{*}\|+\|y_{t+1}-\nabla
h_{2}^{*}(-A\mu_{t+1})\|$ $\displaystyle\leq$
$\displaystyle\rho\left(\left(1-\eta_{2}\frac{\sigma^{2}_{\min}(A)}{2l_{h_{2},1}}\right)\|\mu_{t}-\mu^{*}\|+\eta_{2}\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|\right)$
$\displaystyle+(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}$
$\displaystyle\times\left((l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})\|\mu_{t}-\mu^{*}\|+\sigma_{\max}(A)(1-\eta_{1}\alpha_{h_{2}}/2)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\|\lambda\|\right)$ $\displaystyle=$
$\displaystyle\left(1-\eta_{2}\frac{\sigma_{\min}^{2}(A)}{2l_{h_{2},1}}+\frac{1}{\rho}\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})\right)\rho\|\mu_{t}-\mu^{*}\|$
$\displaystyle+(1-\eta_{1}\alpha_{h_{2}}/2)\left(1+\rho\eta_{2}\sigma_{\max}(A)+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}\right)\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}\|\lambda\|.$

To construct
$P_{t+1}\leq(1-c)P_{t}+\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}\|\lambda\|$
for some constant $0<c<1$, it is sufficient to find
$\eta_{1}\leq\frac{1}{l_{h_{2},1}},\eta_{2}\leq\frac{1}{(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})}$,
and $\rho>0$ such that

$\displaystyle\begin{cases}0<\left(1-\eta_{2}\frac{\sigma_{\min}^{2}(A)}{2l_{h_{2},1}}+\frac{1}{\rho}\frac{\sigma_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})\right)\leq
1-\eta_{2}\frac{\sigma_{\min}^{2}(A)}{4l_{h_{2},1}}<1\\\
0<(1-\eta_{1}\alpha_{h_{2}}/2)\left(1+\rho\eta_{2}\sigma_{\max}(A)+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\eta_{2}\right)\leq(1-\eta_{1}\alpha_{h_{2}}/2)(1+\eta_{1}\alpha_{h_{2}}/2)<1\end{cases}$

This can be obtained when

$\displaystyle\begin{cases}\rho\geq\frac{4l_{h_{2},1}\sigma_{\max}(A)}{\alpha_{h_{2}}\sigma_{\min}^{2}(A)}(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})\\\
\eta_{2}\leq\frac{\eta_{1}\alpha_{h_{2}}}{2\left(\rho\sigma_{\max}(A)+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}}\right)}\end{cases}$
(60)

Conditions in (60) can be satisfied when $\epsilon>0$ is sufficiently small
such that
$\rho=\epsilon^{-1}\geq\frac{4l_{h_{2},1}\sigma_{\max}(A)}{\alpha_{h_{2}}\sigma_{\min}^{2}(A)}(l_{h_{1},1}+\frac{\sigma^{2}_{\max}(A)}{\alpha_{h_{2}}})$,
$\eta_{1}={\cal O}(\frac{1}{l_{h_{2},1}})$ and $\eta_{2}={\cal
O}(\frac{\alpha_{h_{2}}}{l_{h_{2},1}}\rho^{-1})={\cal O}(\epsilon^{-1})$. In
this way,

$\displaystyle P_{t+1}\leq$ $\displaystyle(1-c)P_{t}+{\cal
O}(\alpha_{h_{2}}^{-1}\epsilon)$

where $c>0$ is of the order ${\cal O}(\epsilon)$. Iteration gives

$\displaystyle P_{t}\leq$ $\displaystyle(1-c)^{t}P_{0}+{\cal
O}(\alpha_{h_{2}}^{-1}).$ (61)

Notice ${\cal O}(\alpha_{h_{2}}^{-1})<{\cal O}(1)$ and $P_{0}={\cal
O}(\epsilon^{-1})$ as $\rho=\epsilon^{-1}$. In this way, there exist
$T_{1}={\cal O}(\log(\epsilon^{-1}))$ such that for all $t>T_{1}$,
$(1-c)^{t}P_{0}={\cal O}(1)$ and ${\cal O}(\alpha_{h_{2}}^{-1})\leq{\cal
O}(1)$. Accordingly, we can achieve

$\displaystyle P_{t}={\cal O}(1),\quad\forall t>T_{1}.$

Moreover, as $P_{t}=\epsilon^{-1}\|\mu_{t}-\mu^{*}\|+\|y_{t}-\nabla
h^{*}_{2}(-A\mu_{t})\|$,

$\displaystyle\|\mu_{t}-\mu^{*}\|\leq\epsilon P_{t}={\cal
O}(\epsilon),\quad\forall t>T_{1}.$ (62)

Furthermore, choose $\eta_{1}={\cal O}(\frac{1}{l_{h_{2},1}})$ satisfying
$\eta_{1}\leq\frac{1}{l_{h_{2},1}}$, for $t>T_{1}$,

$\displaystyle\|y_{t}-\nabla h_{2}^{*}(-A\mu_{t})\|\leq$
$\displaystyle(1-\eta_{1}\alpha_{h_{2}}/2)^{t-T_{1}}\|y_{T_{1}}-\nabla
h_{2}^{*}(-A\mu_{T_{1}})\|+{\cal O}(\epsilon).$

This is an iteration outcome using (50) in Lemma 12, (62), and the fact that
$\eta_{1}\alpha_{h_{2}}/2={\cal O}(\frac{\alpha_{h_{2}}}{l_{h_{2},1}})={\cal
O}(1)$. In this way, for another $T_{2}={\cal O}(\log(\epsilon^{-1}))$ steps,
we have

$\displaystyle\|y_{t}-\nabla h_{2}^{*}(-A\mu_{t})\|={\cal O}(\epsilon),$ and
$\displaystyle\|y_{t}-y^{*}\|\leq\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\|\nabla h_{2}^{*}(-A\mu_{t})-\nabla
h_{2}^{*}(-A\mu^{*})\|$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\leq\|y_{t}-\nabla
h_{2}^{*}(-A\mu_{t})\|+\frac{\sigma_{\max}^{2}(A)}{\alpha_{h_{2}}}\|\mu_{t}-\mu^{*}\|$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}={\cal
O}\left(\epsilon+\alpha_{h_{2}}^{-1}\epsilon\right)={\cal
O}\left(\epsilon\right),\quad\forall t>T_{1}+T_{2}.$

We can see that the algorithm converges linearly with complexity ${\cal
O}(T_{1}+T_{2})={\cal O}(\log(\epsilon^{-1}))$. In this way, obtaining

$\displaystyle\|y_{T}-y^{*}\|^{2}={\cal
O}(\epsilon)\quad\text{and}\quad\|\mu_{T}-\mu^{*}\|^{2}={\cal O}(\epsilon),$

requires complexity $T={\cal O}(\log((\sqrt{\epsilon})^{-1}))={\cal
O}(\log(\epsilon^{-1}))$. This completes the proof. ∎

###  Appendix E Applications to Hyperparameter Optimization for SVM

In this section, we provide additional details about the SVM model training
experiment for the linear SVM model, including the problem formulation and
analysis of the results.

####  E.1 Problem formulation

SVMs train a machine learning model by finding the optimal hyperplane that
separates data points of different classes with the maximum margin $w$.
Misclassification is not tolerated for hard-margin SVMs. In contrast, some
samples are allowed to be misclassified for soft-margin SVMs. Specifically,
one first introduces variables $\xi_{i}$, which measure the violation
associated with the classification of sample $i$, and then augments the
original SVM objective with the norm of $\xi$, which is a vector collecting
the values of $\xi_{i}$ for all the samples in the training set.

BLO can be applied to the hyperparameter selection task of SVM. For example,
it can be used to choose the value of the regularization parameters during
soft-margin linear SVM training. Let us consider a classification problem and
define $\mathcal{D}_{\rm tr}:=\\{(z_{{\rm tr},i},l_{{\rm
tr},i})\\}_{i=1}^{|\mathcal{D}_{\rm tr}|}$ as the training set, with $z_{{\rm
tr},i}$ being the input feature vector for sample $i$ and $l_{{\rm tr},i}$
being its associated binary label. Similarly, let $\mathcal{D}_{\rm
val}:=\\{(z_{{\rm val},},l_{{\rm val},})\\}_{i=1}^{|\mathcal{D}_{\rm val}|}$
be the validation set. For a linear classification problem, the parameters of
the SVM are $w$, a vector of coefficients with the same size as $z_{{\rm
tr},i}$, and the intercept $b$.

In short, we are interested in the following constrained BLO problem

$\displaystyle\min_{c}~{}~{}~{}\mathcal{L}_{\mathcal{D}_{\rm
val}}(w^{*},b^{*})=$ $\displaystyle\sum_{(z_{{\rm val},},l_{{\rm
val},})\in\mathcal{D}_{\rm val}}\exp\left({1-l_{{\rm val},}\left(z_{{\rm
val},}^{\top}w^{*}+b^{*}\right)}\right)+\frac{1}{2}\|c\|^{2}$ (63a) with
$\displaystyle w^{*},b^{*},\xi^{*}=\arg\min_{w,b,\xi}\frac{1}{2}\|w\|^{2}$
(63b) $\displaystyle\qquad\qquad\qquad\text{s.t. }~{}l_{{\rm tr},i}(z_{{\rm
tr},i}^{\top}w+b)\geq 1-\xi_{i}\quad\forall\;i\in\\{1,\ldots,|\mathcal{D}_{\rm
tr}|\\}$ (63c) $\displaystyle\qquad\qquad\qquad\quad\;\;\;\xi_{i}\leq
c_{i}\qquad\qquad\qquad\qquad\;\forall\;i\in\\{1,\ldots,|\mathcal{D}_{\rm
tr}|\\}.$ (63d)

The first term of upper-level objective (63a) is a validation loss, evaluated
on the validation set $\mathcal{D}_{\rm val}$, the second is a regulation term
on the upper-level variable $c$; the lower-level problem is designed to train
the parameters of the hyperplane on the training set $\mathcal{D}_{\rm tr}$,
with the soft margin violation $\xi_{i}$ being upper bounded by the
hyperparameter $c_{i}$ [cf. (63d), which are the CCs]. The lower-level
objective focuses on maximizing the margin by minimizing $\|w\|^{2}$ while
allowing for violations $\xi$ to the separating hyperplane, which are
regulated by the hyperparameter $c$. The BLO formulation aims to adjust the
hyperparameter $c$ using the validation loss (upper-level objective), ensuring
that the model parameters (lower-level variables) are optimal for the training
dataset.

####  E.2 Experiment details

In this section, we present detailed experimental results for the training of
the SVM model in (63) using our BLOCC algorithm. We compare our algorithm to
two baselines, LV-HBA [72] and GAM [69], both designed for BLO problems with
inequality CCs. The GAM algorithm often encounters issues with matrix
inversion, as the matrix required for the GAM algorithm to solve the problem
in (63) can be singular, preventing it from finding a solution. Consequently,
only BLOCC and LV-HBA solve the problem in (63), while GAM uses a different
formulation as introduced in [69]. We compare the algorithms based on
validation loss (the first term in the upper-level objective) and
classification accuracy. These metrics are computed for both validation and
test datasets. We evaluate the algorithms on two datasets: diabetes [20] and
fourclass [29].

The detailed results are shown in Figure 4 and 5. Ten realizations of the
problem are run, each with different training and validation sets. The plots
show both the mean (line) and standard deviation (shaded region). The results
reveal that our algorithm converges faster in terms of accuracy and loss,
achieving a lower loss value than the alternatives for both datasets, in
validation and test sets. Furthermore, in terms of accuracy, we observe that
for the diabetes dataset, our algorithm reaches a higher accuracy value faster
than the alternatives in both validation and test sets.

(a) Test accuracy.

(b) Validation loss.

(c) Lower level objective.

Figure 4: Plots for diabetes dataset

(a) Test accuracy.

(b) Validation loss.

(c) Lower level objective.

Figure 5: Plots for fourclass dataset

###  Appendix F Applications to Transportation Network Planning

This section applies the proposed BLOCC algorithm to a transportation network
design problem, comparing it with the baselines from [69] and [72].

####  F.1 Problem formulation

In transportation network planning, the planning operator is to construct a
new transportation network (upper level) connecting a collection of stations
$\mathcal{S}$, and the passengers will decide whether to use the new network
(lower level), considering the options given by the new network and existing
alternatives (constraints). The goal is to design a network that maximizes the
operator’s benefit, knowing the passengers will make rational choices
depending on the given design.

The operator considers building a new network on a collection of links
$\mathcal{A}\subseteq\mathcal{S}\times\mathcal{S}$. For any link
$(i,j)\in\mathcal{A}$ connecting stations $i\in\mathcal{S}$ to
$j\in\mathcal{S}$, the operator needs to design the capacity $x_{ij}$. If the
link’s capacity is set to zero, then the link is not constructed. The larger
the capacity, the larger the number of travelers, thus the larger the revenue
while the higher the construction cost $c_{ij}$.

The passengers are in demand to travel in the network
$\mathcal{K}\subseteq\mathcal{S}\times\mathcal{S}$. For every origin-
destination pair $(o,d)$, the traffic demand $w^{od}$ and the traffic time on
the existing route $t_{\rm ext}^{od}$ is known. We assume that there is only
one existing network. The proportion of passengers choosing the new network
$y^{od}$ and $y^{od}_{ij}$ the fraction of passengers using link $(i,j)$ to
travel from $o$ to $d$ will be determined following passengers’ rational
choice, modeled by logit choice model [5, 11] that will be explained later. If
$y_{ij}^{od}=0$, then the link does not belong to the route that passengers
follow to travel from $o$ to $d$.

In a) of Figure 6, the station $\mathcal{S}$ are presented as the dots, and
links $\mathcal{A}$ are as the dashed lines (input topology); b) of Figure 6
shows a heatmap for the demand matrix for each $(o,d)\in\mathcal{K}$ (input to
the design); c) is an example of the constructed network connecting some of
the stations (output of the design); and d) illustrates the number of
passengers using the constructed network (design output).

Figure 6: Example of a transportation network design. (a) represents the set
of stations $\mathcal{S}=\\{1,2,3,4,5,6,7,8,9\\}$ and the set of links
$\mathcal{A}\subseteq\mathcal{S}\times\mathcal{S}$, where the number of links
is $|\mathcal{A}|=30$ (15 segments with two orientations each) and the value
on each edge represents the travel time. (b) represents the demand
$\\{w^{od}\\}$ between all $(o,d)$ pairs, where $|\mathcal{K}|=9\times 8=72$
values are provided (those in the off-diagonal elements of the heatmap). (c)
represents the constructed network, where, to facilitate visualization, we
have assumed that the capacity is symmetric and used the width of the edge to
represent the capacity of every link. (d) represents
$\\{w^{od}y^{od}\\}_{\\{\forall(o,d)\in\mathcal{K}\\}}$, the number of
passengers served by the constructed network.

To summarize, in this experiment, we assume that

  1. A1)

Only one existing alternative network exists;

  2. A2)

Passengers decide whether to use our network rationally, which is modeled by
the logit choice model considering a utility depends on travel attributes
(travel time) [5, 11]; and,

  3. A3)

The demand per market, the trip prices, and the travel times per link are
known to operators.

We summarize the optimization variables as follows

  * •

$x_{ij}\in\mathbb{R}_{+}$, the capacity constructed for the link
$(i,j)\in\mathcal{A}$. The number of capacity variables is $|\mathcal{A}|$.
The link is not constructed if $x_{ij}=0$.

  * •

$y^{od}\in[0,1]$, the fraction (proportion) of passengers from market
$(o,d)\in\mathcal{K}$ choosing the new network for their travel. The number of
flow variables is $|\mathcal{K}|$. Since we consider only one competitor,
$1-y^{od}$ represents the fraction of incumbent network passengers. If
$y^{od}=0$, then the passengers of market $(o,d)$ do not use the new network.

  * •

$y_{ij}^{od}\in[0,1]$, the fraction of passengers from market
$(o,d)\in\mathcal{K}$ that, when choosing the new network to travel from $o$
to $d$, use the link $(i,j)\in\mathcal{A}$ in their route to the destination.
The number of link flow variables is $|\mathcal{A}||\mathcal{K}|$. With this
definition, it holds that $y_{ij}^{od}\leq y^{od}$; and, if $y_{ij}^{od}=0$,
then link $(i,j)$ does not belong to the route when traveling from $o$ to $d$.

To simplify the notation and to make it consistent with the notations used in
the paper, we denote

  * •

$x=\\{x_{ij}\\}_{\forall(i,j)\in\mathcal{A}}$ represents the upper-level
variables to be optimized.

  * •

$\mathcal{X}=\mathbb{R}^{|\mathcal{A}|}_{+}$ represents the domain of $x$.

  * •

$y=\\{y^{od},\\{y_{ij}^{od}\\}_{\forall(i,j)\in\mathcal{A}}\\}_{\forall(o,d)\in\mathcal{K}}$
represents the lower-level variables to be optimized

  * •

$\mathcal{Y}=[0,1]^{|\mathcal{K}|}\times[0,1]^{|\mathcal{A}||\mathcal{K}|}$
represents the domain of $y$.

Besides the optimization variables, our objective and constraints include the
following parameters

  * •

$w^{od}$, the total estimated demand (number of passengers) for the market
$(o,d)\in\mathcal{K}$.

  * •

$m^{od}$, the revenue obtained by the operator from a passenger in the market
$(o,d)\in\mathcal{K}$.

  * •

$c_{ij}$, the construction cost per passenger associated with link
$(i,j)\in\mathcal{A}$.

  * •

$t_{ij}$, the travel time for link $(i,j)\in\mathcal{A}$.

  * •

$t_{\rm ext}^{od}$, travel time on the alternative network for passengers in
the market $(o,d)\in\mathcal{K}$.

  * •

$\omega_{t}<0$, the coefficient associated with the travel time in passengers’
utility function.

In transportation network design, a bilevel formulation is essential due to
the interaction between two players with different levels of influence: the
operator, who constructs the network, and the passengers, who choose their
routes based on the network’s characteristics. The operator’s goal is to
maximize their benefit by minimizing construction costs and maximizing
attracted demand, with link capacity as the optimization variable. Conversely,
passengers aim to maximize their trip utility, determining the proportion of
demand using each link. This dual optimization requires that passenger choices
comply with link capacity constraints set by the operator, coupling the
variables at both levels. This necessitates a bilevel optimization approach,
specifically using BLOCC, to address the interdependent decisions and
constraints effectively.

Now, we are ready to introduce the objective formulations of our BLO problem.
For the upper level, the network operator aims to maximize profits and
minimize costs, and therefore, its interest is

$\min_{x\in\mathcal{X}}f(x,y^{*}_{g}(x)):=-\Bigg{(}\sum_{\forall(o,d)\in\mathcal{K}}m^{od}y^{\rm
od*}(x)-\sum_{\forall(i,j)\in\mathcal{A}}c_{ij}x_{ij}\Bigg{)},$ (64)

where $y^{\rm od*}(x)$ are optimal lower-level passenger flows associated with
the network design $x$.

For the lower level, we model the passenger’s behavior by finding the flow
variables that maximize utility and minimizes flow entropy cost.

$\displaystyle\min_{y\in\mathcal{Y}}g(x,y):=$
$\displaystyle-\Bigg{(}\underbrace{\sum_{(o,d)\in\mathcal{K}}\sum_{(i,j)\in\mathcal{A}}w^{od}\omega_{t}t_{ij}y_{ij}^{od}+\sum_{(o,d)\in\mathcal{K}}w^{od}\omega_{t}t_{\rm
ext}^{od}(1-y^{od})}_{\rm passengers~{}utility}$ (65)
$\displaystyle\qquad+\underbrace{\sum_{(o,d)\in\mathcal{K}}w^{od}y^{od}(\log(y^{od})-1)+\sum_{(o,d)\in\mathcal{K}}w^{od}(1-y^{od})(\log(1-y^{od})-1)}_{\rm
flow~{}entropy~{}cost}\Bigg{)}.$

The passengers’ utility considers the time cost of choosing the new network
and the existing network. The users set the flow variables $y$ so that the
transportation network yielding a higher utility is preferred. Here, the
probability of chosing new network is modeled by a logistic (softmax) model.
However, setting the objective as a simple linear utility maximization would
lead to an all-or-nothing policy, which is not the behavior observed in
practice. Hence, the second term is brought to consider. The approach
considered here is to formulate an objective given by the Legendre transform
(cf. Definition 5) of the softmax sharing (see [53, 56, 55] for additional
details). Intuitively, this means that rather than imposing the softmax
sharing a fortiori, we formulate a convex problem whose KKT conditions lead to
the softmax sharing. Using the fact that the Legendre transform of an
exponential $e^{y}$ is the negative entropy function $y(\log(y)-1)$, the
lower-level objective is traditional as in [56].

Having introduced the optimization variables, parameters, and objective
functions, we next formulate our BLO problem, where we also incorporate the
network constraints:

$\displaystyle\qquad\min_{x\in\mathcal{X}}\;\;-\sum_{\forall(o,d)\in\mathcal{K}}m^{od}y^{\rm
od*}+\sum_{\forall(i,j)\in\mathcal{A}}c_{ij}x_{ij}$ (66a)
$\displaystyle\text{s.t.}~{}(y^{\rm od*},y_{ij}^{\rm
od*})=\arg\min_{y\in\mathcal{Y}}{-\sum_{(o,d)\in\mathcal{K}}\sum_{(i,j)\in\mathcal{A}}w^{od}\omega_{t}t_{ij}y_{ij}^{od}-\sum_{(o,d)\in\mathcal{K}}w^{od}\omega_{t}t_{\rm
ext}^{od}(1-y^{od})}$ (66b)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+{\sum_{(o,d)\in\mathcal{K}}w^{od}y^{od}(\log(y^{od})-1)+\sum_{(o,d)\in\mathcal{K}}w^{od}(1-y^{od})(\log(1-y^{od})-1)}$
$\displaystyle\qquad\qquad\text{s.t.}\;\;\sum_{\forall
j|(i,j)\in\mathcal{A}}y_{ij}^{od}-\sum_{\forall
j|(j,i)\in\mathcal{A}}y_{ji}^{od}=\begin{cases}y^{od}\;\quad\text{if}\;i=o\\\
-y^{od}\;\ \text{if}\;i=d\\\ 0\;\;\qquad\text{otherwise}\end{cases}\forall
i,(o,d)\in\mathcal{S}\times\mathcal{K}$ (66c)
$\displaystyle\qquad\qquad\;\;\;\;\;\;\sum_{\forall(o,d)\in\mathcal{K}}w^{od}y_{ij}^{od}\leq
x_{ij}\qquad\qquad\;\;\forall(i,j)\in\mathcal{A}$ (66d)

where (66c) are the flow-conservation constraints, (66d) are the capacity
constraints that involve both upper and lower-level variables, coupling the
optimization and motivating the use of BLOCC. Note that, for networks with $n$
stations (nodes): i) the number of CCs is approximately $n^{2}$; and ii) each
of the constraints involves approximately $n^{2}$ variables. Hence, even for a
moderate-size network (say 30-50 nodes), we may have thousands of CCs
involving millions of variables.

##### Experiment roadmap.

To provide numerical results illustrating the behavior of our algorithm, we
solve the optimization in (66) for three scenarios:

  1. S1)

The design of a 3-node simple synthetic network;

  2. S2)

The design of a 9-node synthetic network from the prior transportation
literature; and

  3. S3)

The design of a (real-world) subway network for the city of Seville, Spain,
with 24 nodes

where S1) involves 6 CCs and 48 variables, and S3) involves around 100 CCs and
50,000 variables.

For the 3-node network scenario, we will conduct a comparative analysis
against other algorithms to evaluate the efficacy of our approach. In the
other two scenarios, the baseline algorithms cannot find a solution; hence,
for the 9-node and Seville networks, we will focus on providing insights into
the performance and behavior of our algorithm under varying parameters,
shedding light on the versatility and adaptability of our approach to real-
world transportation networks.

Before delving into the presentation and analysis of the results, two
additional remarks are in order:

  * R1)

While one of the goals of these experiments was to compare our BLOCC algorithm
against LV-HBA [72] and GAM [69], for the scenario at hand, the GAM algorithm
cannot be implemented, since the inverse of a matrix at each iteration for the
problem in (66) is not tractable. In this way, we only conducted the
experiments using our BLOCC and LV-HBA.

  * R2)

The BLOCC algorithm produces two lower level variables: $y_{F,T}^{T_{F}}$ and
$y_{g,T}^{T_{g}}$. While Theorems 1 and 2 guarantee that
$(x_{T},y_{F,T}^{T_{F}})$ is an $\epsilon$-approximate solution, the pair
$(x_{T},y_{g,T}^{T_{g}})$ is strictly feasible, meaning that
$y_{g,T}^{T_{g}}=y_{g}^{*}(x_{T})$. Since strict feasibility is important in
the transportation network, this section will report the results using both
$y_{F,T}^{T_{F}}$ and $y_{g,T}^{T_{g}}$.

Figure 7:  Negative upper-level objective $-f(x_{t},y_{t})$ for a 3-node
network design problem; Solid lines show mean value of
$-f(x_{t},y_{g,t}^{T_{g}})$ with the shaded region as a standard deviation;
Dashed lines show the mean value of $-f(x_{t},y_{F,t}^{T_{F}})$ with the
shaded region as a standard deviation; Three different $\gamma$ values (red,
purple, blue) and fixed stepsize $\eta=1.6e-4$; The orange color represents
the result of the LV-HBA algorithm.

Figure 8:  Optimality gap of the lower-level problem for a 3-node network
design problem $g(x_{t},y_{t})-g(x_{t},y^{*}_{t})$. Solid lines represent the
mean value of the 10 realizations of the upper-level variables, dashed lines
represent $g(x_{t},y_{F,t}^{T_{F}})-g(x_{t},y^{*}_{t})$, and the shaded region
is the standard deviation. Three different $\gamma$ values are represented in
our algorithm, and fixed stepsize $\eta=1.6e-4$.

####  F.2 Numerical results for the 3-node network

| Market $(o,d)$
---|---
Parameter | (1,2) | (1,3) | (2,1) | (2,3) | (3,1) | (3,2)
demand $w^{od}$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$
revenue $m^{od}$ | $2$ | $6$ | $2$ | $1$ | $6$ | $1$
travel time incumbent $t_{\rm ext}^{od}$ | $3$ | $3$ | $3$ | $3$ | $3$ | $3$