# An Introduction to Bi-level Optimization: Foundations and Applications in
Signal Processing and Machine Learning

Yihua Zhang1, Prashant Khanduri2, Ioannis Tsaknakis3, Yuguang Yao1,
Mingyi Hong3, Sijia Liu1,4
1CSE, Michigan State University, USA, 2CS, Wayne State University, USA,
3ECE, University of Minnesota, USA, 4MIT-IBM Watson AI Lab, USA

###### Abstract

Recently, bi-level optimization (BLO) has taken center stage in some very
exciting developments in the area of signal processing (SP) and machine
learning (ML). Roughly speaking, BLO is a classical optimization problem that
involves two levels of hierarchy (i.e., upper and lower levels), wherein
obtaining the solution to the upper-level problem requires solving the lower-
level one. BLO has become popular largely because it is powerful in modeling
problems in SP and ML, among others, that involve optimizing nested objective
functions. Prominent applications of BLO range from resource allocation for
wireless systems to adversarial machine learning. In this work, we focus on a
class of tractable BLO problems that often appear in SP and ML applications.
We provide an overview of some basic concepts of this class of BLO problems,
such as their optimality conditions, standard algorithms (including their
optimization principles and practical implementations), as well as how they
can be leveraged to obtain state-of-the-art results for a number of key SP and
ML applications. Further, we discuss some recent advances in BLO theory, its
implications for applications, and point out some limitations of the state-of-
the-art that require significant future research efforts. Overall, we hope
that this article can serve to accelerate the adoption of BLO as a generic
tool to model, analyze, and innovate on a wide array of emerging SP and ML
applications.

###### Index Terms:

Bi-level optimization (BLO), convex/non-convex optimization, signal processing
(SP), machine learning (ML), resource allocation, wireless communications,
adversarial robustness, generalization, datamodel efficiency,

##  I Introduction

Bi-level optimization (BLO) is a class of optimization problems involving two
nested levels (upper and lower-level), where the objective and variables of
the upper-level problem depend on the optimizer of the lower-level one. The
canonical formulation of BLO is given by:

$\displaystyle\overbrace{\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}\in\mathcal{U}}\leavevmode\nobreak\
\leavevmode\nobreak\ f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta})),}^{\text{Upper-
level optimization over $\bm{\theta}$}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname*{\text{subject to}}\leavevmode\nobreak\ \leavevmode\nobreak\
\underbrace{\displaystyle\bm{\phi}^{*}(\bm{\theta})\in\operatorname*{arg\,min}_{h(\bm{\theta},\bm{\phi})\leq
0}g(\bm{\theta},\bm{\phi})}_{\text{Lower-level optimization over
$\bm{\phi}$}},$ (BLO)

where we assume that $f$, $g$ and $h$ are bivariate smooth functions,
$\bm{\theta}\in\mathbb{R}^{m}$ denotes the upper-level variable subject to the
upper-level constraint set $\mathcal{U}$, $\bm{\phi}\in\mathbb{R}^{n}$ is the
lower-level variable subject to the constraint $h(\bm{\theta},\bm{\phi})\leq
0$ that couples both $\bm{\theta}$ and $\bm{\phi}$, and
$\bm{\phi}^{*}(\bm{\theta})$ is one lower-level optimal solution. It is
evident that the lower-level problem is an auxiliary problem since its
solution supports the upper-level problem in finding a better solution.

The study of BLO can be traced to that of Stackelberg games [1], where the
upper (resp. lower) problem optimizes the action taken by a leader (resp. the
follower). Early works in optimization utilized BLO to solve resource
allocation problems [2, 3, 4]; see [5] for a comprehensive survey of BLO
algorithms up until mid-2000, and some more recent surveys on discrete BLO
[6], BLO under uncertainty [7], and nonlinear and nonconvex aspects of BLO
[8]. In recent years, BLO has regained popularity because a subclass of BLO
has been used to formulate and solve various challenging problems in signal
processing (SP), machine learning (ML), and artificial intelligence (AI).
Notable applications in SP include resource management [9, 10, 11, 12], signal
demodulation [13, 14, 15], channel prediction [16, 17, 18], image
reconstruction [19] and image denoising [20]. In addition, BLO has also been
used to make ML models, especially deep neural networks (DNNs), robust [21,
22, 23, 24, 25, 26, 27], generalizable [28, 29, 30, 31, 32, 33, 34, 35, 36,
37, 38, 39], efficient [40, 41, 42, 43, 44], easier to train [45, 46, 47] and
scalable [48, 49, 50, 51, 52, 53, 54].

As can be easily imagined, the popularity of BLO in aforementioned
applications is largely attributed to its ability to handle (often implicit)
hierarchical structures. To better illustrate the challenges brought by the
hierarchical architecture, an example application of coreset selection for
model training [55, 40, 12] is discussed below.

Motivating application: Coreset selection for model training Many contemporary
SP and ML applications are facing significant challenges in data storage,
transportation, and computation because they have to deal with excessive
amounts of data. Consequently, the task of identifying the most informative
subset of data from a larger pool becomes crucial [55, 40, 12]. This leads to
the problem of coreset selection, which consists of two tasks: (T1) selecting
the most representative data samples to form the coreset, and (T2) validating
the performance of the selected coreset in model training. More specifically,
the problem can be formulated as follows:
$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\mathbf{w}\in\mathcal{U}}$
$\displaystyle\,\,\ell_{\mathrm{val}}({\bm{\theta}^{\ast}}(\mathbf{w}))$ (T1)
$\displaystyle\operatorname*{\text{subject to}}$
$\displaystyle\,\,\bm{\theta}^{\ast}(\mathbf{w})=\displaystyle\operatorname*{arg\,min}_{\bm{\theta}}\ell_{\mathrm{tr}}(\bm{\theta},\mathbf{w}),$
(T2) where $\mathbf{w}$ represents the weight vector for data selection, with
$w_{i}=0$ indicating that the $i^{\text{th}}$ data sample is not selected.
These weights are subject to the sparsity constraint $\mathcal{U}$, such as
$\|\mathbf{w}\|_{1}\leq k$, with $k$ being the selection budget. The model
parameters trained on the selected data points are denoted by $\bm{\theta}$.
The training loss for the model $\bm{\theta}$ with the data selection scheme
$\mathbf{w}$ is denoted by $\ell_{\mathrm{tr}}$, while the validation loss
measuring the performance of the learned model
$\bm{\theta}^{\ast}(\mathbf{w})$ over the coreset is denoted by
$\ell_{\mathrm{val}}$. The above BLO formulation is also related to the data
reweighting problem [56] and hyperparameter optimization [57].

Clearly, the coreset selection is a typical BLO problem, where the upper-level
and lower-level tasks are intertwined: without knowing the training result, it
is hard to gauge how representative the selected dataset is effective, while
without having the coreset, one cannot perform model training to evaluate the
utility of the learned model further. More importantly, these two tasks form a
hierarchy, with the model training problem being the main optimization
problem, while the data selection problem is an auxiliary problem that
supports the training.

Given the growing interest in BLO, this work presents an overview of a class
of tractable (BLO) problems that hold significant importance in SP and ML.
Roughly speaking, the class of BLO problems we consider have some desirable
properties (to be discussed shortly) that allow the development of efficient
and practical algorithms. We will discuss the basic concepts for this class of
BLO problems along with their optimality conditions, standard algorithms
(including their theoretical properties and practical implementations), as
well as how they can be used to obtain state-of-the-art results for a number
of key SP and ML applications.

In the existing literature, several recent surveys have been conducted on
general BLO problems [58, 7, 8]. However, these surveys primarily focus on the
mathematical foundations of BLO through a classical optimization lens. Other
works [59, 60, 61] aim to provide comprehensive reviews of BLO algorithms, but
they lack an in-depth discussion on the proper recognition and utilization of
BLO, specifically addressing ‘when’ and ‘how’ to apply it effectively. The
most relevant work to ours is [61], which examines complex learning and vision
problems from the BLO perspective. However, the theoretical component of BLO
is missing, and it overlooks a significant portion of emerging SP and ML
applications (e.g., those discussed in Sec. V-IX).

Different from the existing surveys on BLO [7, 8, 59, 60, 61] to provide a
broad overview of BLO in its most generic form, we focus on problems that are
relevant to SP and ML applications and strive to tightly integrate theory with
applications, especially those that have been largely left out of the existing
surveys. A few highlights of this article are listed below.

First, we distill the common structures and properties of BLO that emerge
across applications related to developing robust, parsimonious, and
generalizable data-driven models in SP and ML. Our goal is to provide insights
about when, where, and how BLO formulations and algorithms can be best used to
yield a significant performance boost, as compared with traditional, or
heuristic algorithms. In this process, we present some recent theoretical
results about BLOs and their associated algorithms to give a flavor of the
current advances in the research area, while discussing their practical and
scalable implementations.

Second, we dive deep to understand the performance of a selected subset of
state-of-the-art (SOTA) BLO algorithms on a number of representative
applications. Instead of relying on results reported in existing works, which
may not always be comprehensive, we designed an experiment plan and
implemented all benchmarking algorithms ourselves. The goal is not only to
showcase the effectiveness of the BLO-based algorithms but also to analyze the
pros (e.g., modeling flexibility and accuracy performance) and the cons (e.g.,
run-time efficiency) of different subclasses of BLO methods.

Overall, we hope that our balanced treatment of the subject will serve as the
cornerstone for the accelerated adoption of BLO in diverse application areas.
We hope that this tutorial can assist researchers with modeling, analyzing,
and innovating on the emerging SP and ML applications. And it can serve to
sparkle new theoretical and applied research to advance BLO. Fig. 1 provides
an overview of the topics to be covered in this article.

Notations: We use lower case letters (e.g., $a$), lower case bold face letters
(e.g., $\mathbf{a}$), and upper case bold face letters (e.g., $\mathbf{A}$) to
denote scalars, vectors, and matrices, respectively. For a vector
$\mathbf{a}$, we use $\|\mathbf{a}\|_{p}$ to denote its $\ell_{p}$ norm with
the typical choice $p\in\\{1,2,\infty\\}$. For a matrix $\mathbf{A}$, we use
the superscript $\top$ (or -1) to denote the transpose (or inverse) operation.
We use $\mathbf{I}$ to represent the identity matrix. For a function
$f(\mathbf{x},\mathbf{y})$ (with $\mathbf{x}\in\mathbb{R}^{m}$ and
$\mathbf{y}\in\mathbb{R}^{n}$), we use
$\nabla_{\mathbf{x}}f(\mathbf{x},\mathbf{y})\in\mathbb{R}^{m}$ (or
$\frac{\partial f}{\partial\mathbf{x}}$) and
$\nabla_{\mathbf{y}}f(\mathbf{x},\mathbf{y})\in\mathbb{R}^{n}$ (or
$\frac{\partial f}{\partial\mathbf{y}}$) to denote the partial derivatives of
$f$ w.r.t. the partial input argument $\mathbf{x}$ and $\mathbf{y}$,
respectively. In contrast, we use $\frac{df}{d\mathbf{x}}\in\mathbb{R}^{m}$ to
represent the full derivative of $f$ w.r.t. $\mathbf{x}$, namely,
$\frac{df}{d\mathbf{x}}=\nabla_{\mathbf{x}}f(\mathbf{x},\mathbf{y})+\frac{d\mathbf{y}}{d\mathbf{x}}^{\top}\nabla_{\mathbf{y}}f(\mathbf{x},\mathbf{y})$
following the chain rule, where
$\frac{d\mathbf{y}}{d\mathbf{x}}^{\top}\in\mathbb{R}^{m\times n}$ denotes the
Jacobian matrix of $\mathbf{y}$ w.r.t. $\mathbf{x}$. For ease of notation, the
transpose in $\frac{d\mathbf{y}}{d\mathbf{x}}^{\top}\in\mathbb{R}^{m\times n}$
will be omitted if its definition is clear from the context. We use
$\nabla_{\mathbf{x}}\nabla_{\mathbf{y}}f$ or
$\nabla_{\mathbf{x},\mathbf{y}}^{2}f\in\mathbb{R}^{m\times n}$ to denote the
second-order partial derivative of $f$.

Figure 1: An taxonomy of the solvers and application tasks for BLO.

##  II Warm-up: Introducing Basic Concepts of BLO

###  II-A A Class of Tractable BLOs

We start by discussing the challenges associated with the generic form of
(BLO). Even under the assumption that all involved functions are well-behaved,
such as the convexity or concavity of $f(\cdot,\cdot)$ and $g(\cdot,\cdot)$,
and the linearity of $h(\cdot,\cdot)$, solving the problem can still be highly
challenging (i.e., NP-hard). To see this, let us consider the following simple
example.

Example 1: Non-convex BLO Consider the following BLO, where
$f(\theta,{\phi})=-g(\theta,{\phi})=\theta^{2}-\theta\cdot{\phi}-{\phi}^{2}$:
$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\theta\in[-1,1]}\leavevmode\nobreak\
\theta^{2}-\theta\cdot{\phi}^{\ast}(\theta)-{\phi}^{\ast}(\theta)^{2};\leavevmode\nobreak\
\leavevmode\nobreak\ \operatorname*{\text{subject to}}\leavevmode\nobreak\
{\phi}^{\ast}(\theta)=\operatorname*{arg\,min}_{{\phi}\in[-1,1],\theta-{\phi}=0}-\left(\theta^{2}-\theta\cdot{\phi}-{\phi}^{2}\right).$
Notice that the objective $f(\theta,{\phi})$ (resp. $g(\theta,{\phi})$) is
strongly convex (resp. strongly concave) in $\theta$, strongly concave in
${\phi}$ (resp. strongly convex), and the above BLO problem is subject to
linear constraints in both the upper and lower level. In other words, both the
upper and lower-level problems are convex problems. Nonetheless, it can be
shown that solving the above BLO requires tackling a non-convex problem, which
in general is NP-hard. Indeed, it is not hard to see that
${\phi}^{\ast}(\theta)=\theta$. As a result, the outer function can be
expressed as $\ell(\theta)=-\theta^{2}$, which is a non-convex function.

We remark that the source of difficulty of the above problem is the coupling
constraint $\theta-\phi=0$. If this constraint is removed, the problem will
become a classical saddle-point problem, expressed below, whose global optimal
solution can be easily obtained:

$\displaystyle\operatorname*{\text{minimize}}\limits_{\theta\in[-1,1]}\operatorname*{\text{maximize}}_{\phi\in[-1,1]}\,\,\theta^{2}-\theta\cdot{\phi}-{\phi}^{2}.$
(1)

The aforementioned examples, along with numerous others in existing survey
papers like [61, 59, 60, 58], strongly motivate us to proceed with a focused
discussion for the subset of tractable (BLO) problems. This subset often
serves as a basis for developing practical BLO algorithms. By focusing on this
subset, we can address the specific needs of modeling SP and ML problems,
which frequently demand the development of efficient, and sometimes real-time,
algorithms.

To this end, we consider some special classes of (BLO), with the following
simplifications: 1) The lower-level constraint set, if present, is linear that
is only related to $\bm{\phi}$, that is,
$h(\bm{\theta},\bm{\phi})=\mathbf{A}\bm{\phi}-\mathbf{b}$ for some matrix
$\mathbf{A}$ and vector $\mathbf{b}$ of appropriate sizes; 2) The solution of
the lower-level problem is a singleton, and in most cases we assume an even
stronger condition that the objective function $g(\cdot,\cdot)$ is strongly
convex in the second argument. With these simplifications, it is possible to
show some nice properties, for example, the gradient of the upper-level
objective function may exist, making algorithm design and analysis tractable.
Note that there have been recent works that extend these conditions and are
still able to develop efficient algorithms. We will discuss these works in
Sec. IV.

In summary, depending on whether the lower-level problem has constraint or
not, we consider the following two classes of problems, referred to as the
lower-uncsontrained (LU) and lower-constrained (LC) BLO, respectively:

$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}\in\mathcal{U}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta})),\leavevmode\nobreak\
\leavevmode\nobreak\ \operatorname*{\text{subject to}}\leavevmode\nobreak\
\leavevmode\nobreak\
{\displaystyle\bm{\phi}^{*}(\bm{\theta})=\operatorname*{arg\,min}_{\bm{\phi}\in\mathbb{R}^{n}}g(\bm{\theta},\bm{\phi})}.$
(LU-BLO)
$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}\in\mathcal{U}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta})),\leavevmode\nobreak\
\leavevmode\nobreak\ \operatorname*{\text{subject to}}\leavevmode\nobreak\
\leavevmode\nobreak\
{\displaystyle\bm{\phi}^{*}(\bm{\theta})=\operatorname*{arg\,min}_{\bm{\phi}\in\mathcal{C}}g(\bm{\theta},\bm{\phi})}.$
(LC-BLO)

where the set
$\mathcal{C}:=\\{\bm{\phi}\mid\mathbf{A}\bm{\phi}-\mathbf{b}\leq\mathbf{0}\\}$.
As will be evident in Sec. III, the presence of lower-level constrains, even
in the form of linear and uncoupled constraints, can make (LU-BLO) much harder
to deal with than (LC-BLO).

###  II-B Connections of BLO with Game Theory

It is important to note that BLO problems have strong ties with Stackelberg or
leader-follower games [1], including Stackelberg congestion [62] and security
games [63]. These are sequential games involving two players, the leader and
the follower. The leader acts first, aiming to maximize its utility by
leveraging its knowledge of the follower’s anticipated response. The follower,
acting second, maximizes its utility based on the leader’s action. The
connection between BLO and Stackelberg games can be bifurcated as follows:
First, in certain Stackelberg games, the process of identifying a solution
(i.e., a Stackelberg equilibrium) can be framed as a BLO. Second, BLO allows a
(Stackelberg) game-theoretic interpretation, where the upper- and lower-level
problems correspond to the tasks of identifying the optimal actions for the
leader and the follower (i.e., the upper/lower-level variables), respectively.
A special case of Stackelberg game is min-max optimization (MMO) also referred
to as the saddle point problem. MMO follows a bi-level structure wherein the
lower-level objective $g$ in (BLO). is exactly opposite of the upper-level
objective function $f$ (i.e., $g=-f$) resulting in the following special case
of BLO:

$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}\in\mathcal{U}}\,\operatorname*{\text{maximize}}_{\bm{\phi}\in\mathcal{C}}&f(\bm{\theta},\bm{\phi}).\end{array}$
(MMO)

In fact, MMO is much simpler to deal with than BLO, and it has been heavily
studied in the SP and ML communities. We refer the interested reader to the
recent tutorial [64] for detailed discussion. In addition to MMO, it is also
worth noting that the algorithms reviewed in this paper are generally
applicable to Stackelberg games, provided that the game’s BLO formulation
adheres to the assumptions of the respective algorithms.

###  II-C Implicit Gradient (IG)

As alluded previously, one important reason to consider problems (LU-BLO) and
(LC-BLO) is that the objective functions of these problems are potentially
differentiable w.r.t. $\bm{\theta}$. Indeed, by applying the chain rule, and
suppose for now that the Jacobian matrix
$\frac{{d}\bm{\phi}^{*}(\bm{\theta})}{{d}\bm{\theta}}$ exists, we have:

$\displaystyle\frac{df(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))}{d\bm{\theta}}=\nabla_{\bm{\theta}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))+{\underbrace{\frac{{d}\bm{\phi}^{*}(\bm{\theta})^{\top}}{{d}\bm{\theta}}}_{\text{{{\text{IG}}}}}}\nabla_{\bm{\phi}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta})),$
(3)

where recall from our notation convention that
$\nabla_{\bm{\theta}}f(\bm{\theta},\bm{\phi})$ and
$\nabla_{\bm{\phi}}f(\bm{\theta},\bm{\phi})$ represent the partial derivatives
of $f$ w.r.t. the partial input arguments $\bm{\theta}$ and $\bm{\phi}$
respectively, and $\frac{df}{d\bm{\theta}}$ denotes the full derivative of $f$
w.r.t. $\bm{\theta}$. For ease of notation, the transpose operation ⊤ might be
omitted in the rest of the paper. We refer to the Jacobian matrix
$\frac{{d}\bm{\phi}^{*}(\bm{\theta})}{{d}\bm{\theta}}$ as IG following [65].
This term is introduced to characterize the gradient of the
$\operatorname*{arg\,min}$-based lower-level objective function w.r.t the
upper-level variable $\bm{\theta}$. However, IG does not always exist for
generic BLO problems. Even for (LU-BLO) and (LC-BLO), relatively strong
assumptions have to be imposed; For example, $g(\cdot,\cdot)$ needs to be
strongly convex in its second argument. Further, even IG exists, computing it
could be quite different for the two classes of problems (LU-BLO) and (LC-
BLO). For the former, we will show in Sec. III that the IG can be expressed in
closed-form using the Implicit Function Theorem [66] based on the first-order
stationary condition of the lower-level problem, i.e.,
$\nabla_{\bm{\phi}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))=\mathbf{0}$. And
this is the reason that it is referred to as “implicit” gradient. Yet, in (LC-
BLO), the stationary condition cannot be used since a stationary point might
violate the constraint $\bm{\phi}\in\mathcal{C}$. Therefore IG generally does
not admit any closed-form. Additionally, for a more restricted subset of
problems MMO, the influence of the IG-involved term
($\mathrm{IG}\cdot\nabla_{\bm{\phi}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))$)
in (3) can be neglected. To see this, assume that the inner problem in
unconstrained, i.e., $\mathcal{C}\equiv\mathbb{R}^{n}$, then
$\nabla_{\bm{\phi}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))=\mathbf{0}$ based
on the fact that
$\nabla_{\bm{\phi}}g(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))=\mathbf{0}$ and
$g=-f$ for MMO.

###  II-D BLO with Non-Singleton Lower-Level Solutions

As we have mentioned in Sec. II-A, throughout this article, we will mostly
focus on the case where the lower-level solution $\phi^{\ast}(\bm{\theta})$ is
unique, i.e., a singleton. Yet, if the lower-level problem involves a non-
singleton (NS) solution, the resulting BLO problem is typically cast as

$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta},\bm{\phi}^{\prime}\in\mathcal{S}(\bm{\theta})}f(\bm{\theta},\bm{\phi}^{\prime})\leavevmode\nobreak\
\leavevmode\nobreak\ \operatorname*{\text{subject to}}\leavevmode\nobreak\
\leavevmode\nobreak\
\mathcal{S}(\bm{\theta})\coloneqq\operatorname*{arg\,min}_{\bm{\phi}\in\mathcal{C}}g(\bm{\theta},\bm{\phi})$
(NS-BLO)

where $\mathcal{S}(\bm{\theta})$ denotes a solution set. The above formulation
also referred to as the optimistic BLO with non-singleton lower-level
solutions, has been discussed in the literature such as [6, 58]. Note that
problem (NS-BLO) presents significantly greater challenges from both practical
and theoretical perspectives compared to problem (BLO). This is because
optimization over $\bm{\phi}$ is coupled across both upper- and lower-level
objectives. While our work primarily focuses on BLO with a singleton lower-
level solution, we will also explore in Sec.,III the applicability of BLO
algorithms, derived from (BLO), to solve problem (NS-BLO).

Theory and Algorithms for Tractable BLOs

In the next two sections, we will delve into the essential optimization
principles employed in BLO, explore several popular classes of BLO algorithms,
and examine their theoretical properties.

##  III Algorithmic Foundations of BLO

This section presents an overview of three key optimization frameworks used to
solve the tractable BLO problems (LU-BLO) and (LC-BLO). The first two classes
both leverage (some approximated version of) the IG as defined in (3). The key
difference is how the approximation of the IG is conducted: One directly
assumes that there is some given procedure that can provide a high-quality
solution of the lower-level problem, while the other approximates the lower-
level solution by unrolling a given algorithm for a fixed number of steps. The
third class is referred to as the value function (VF)-based approach, which
reformulates BLO as a single-level regularized optimization problem. It is
worth mentioning that this approach offers flexibility in handling lower-level
constraints and solving non-singleton lower-level problems (NS-BLO).

###  III-A The Implicit Function (IF)-based Approach

####  III-A1 IF for lower-level unconstrained BLO

Let us examine the problem setup (LU-BLO) with a singleton lower-level
solution. For ease of theoretical analysis in Sec. IV, we further assume that
$g(\cdot)$ is strongly convex. In certain applications, one can explicitly add
a strongly convex regularization function, such as
$\gamma\times\|\bm{\phi}\|_{2}^{2}$ (with large enough $\gamma$), to satisfy
such an assumption.

The reason that we call this approach IF-based is that we will explicitly
utilize the Implicit Function Theorem [67] to calculate the IG as expressed in
(3). Recall from (LU-BLO) that $\bm{\phi}^{*}(\bm{\theta})$ is a lower-level
solution, thus it satisfies:

$\displaystyle\nabla_{\bm{\phi}}g(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))=\mathbf{0}.$
(4)

Following the implicit function theorem, we can take the first-order
derivative of (4) w.r.t. $\bm{\theta}$, yielding
$\frac{d}{d\bm{\theta}}[\nabla_{\bm{\phi}}g(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))]=\mathbf{0}$.
Further assume that the lower-level objective $g(\cdot)$ is second-order
differentiable, we can then obtain the IG in (3) in the following:

$\displaystyle{\frac{d\bm{\phi}^{*}(\bm{\theta})}{d\bm{\theta}}}=-\nabla^{2}_{\bm{\theta},\bm{\phi}}g(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))\nabla^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))^{-1}.$
(5)

As observed above, the computation of IG involves the mixed (second-order)
partial derivative $\nabla^{2}_{\bm{\theta},\bm{\phi}}g$ and the inverse of
the Hessian $\nabla^{2}_{\bm{\phi},\bm{\phi}}g$. Yet, computing these
quantities can be challenging in practice. Therefore, the class of IF approach
utilizes different kinds of approximation techniques to approximately compute
the IG as expressed (5). We summarize the IF approach in Algorithm 1.

Algorithm 1: IF-based approach for solving (LU-BLO) Given initialization
$\bm{\theta}_{0}$, learning rate $\alpha>0$, and iteration number $T$;
Iteration $t\geq 0$ yields: $\bullet$ Lower-level optimization: Given
$\bm{\theta}_{t}$, obtain an approximate solution of the lower-level problem,
denoted as $\widetilde{\bm{\phi}}(\bm{\theta}_{t})$; $\bullet$ Approximation:
Based on $\widetilde{\bm{\phi}}(\bm{\theta}_{t})$, compute approximated
versions of two second-order matrices in (5), denoted as
$\widetilde{\nabla}^{2}_{\bm{\theta},\bm{\phi}}g(\bm{\theta},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))$
and
$\widetilde{\nabla}^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))^{-1}$;
Compute an approximated IG following (3):
$\displaystyle{\tilde{\nabla}f(\bm{\theta}_{t})}\coloneqq\nabla_{\bm{\theta}}f(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))-\widetilde{\nabla}^{2}_{\bm{\theta}_{t},\bm{\phi}}g(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))\widetilde{\nabla}^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))^{-1}\nabla_{\bm{\phi}}f(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t})).$
(6) $\bullet$ Upper-level optimization: Utilize
$\tilde{\nabla}f(\bm{\theta}_{t})$ in (6) to update $\bm{\theta}_{t}$ through,
e.g., gradient descent (GD):
$\bm{\theta}_{t+1}\leftarrow\bm{\theta}_{t}-\alpha\tilde{\nabla}f(\bm{\theta}_{t})$.

####  III-A2 Practical considerations of IF

In Algorithm 1, the main computational overhead arises from the inverse-
Hessian gradient product $\mathbf{H}^{-1}\mathbf{g}$, where
$\mathbf{H}:=\widetilde{\nabla}^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))$
and
$\mathbf{v}:=\nabla_{\bm{\phi}}f(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))$.
Yet, in many contemporary applications, directly computing and storing the
Hessian is computationally prohibitive. To address the scalability challenge
in the IF method, we introduce four approaches to approximate the inverse-
Hessian gradient product (or the inverse-Hessian) $\mathbf{H}^{-1}\mathbf{g}$:
the conjugate gradient (CG) method [68, 57], the WoodFisher approximation
[69], the Neumann-series method to directly estimate the inverse-Hessian [66,
70], and a Hessian-free simplification [22, 44]. These methods offer different
trade-offs between computational costs, with the CG method being the most
computationally expensive and the Hessian-free simplification being the least
expensive.

First, the CG approach maps the product $\mathbf{H}^{-1}\mathbf{g}$ to the
solution of a quadratic program defined as
$\min_{\mathbf{x}}\mathbf{x}^{\top}\mathbf{H}\mathbf{x}/2-\mathbf{g}^{\top}\mathbf{x}$.
By utilizing the first-order gradient descent algorithm, we can numerically
approximate $\mathbf{H}^{-1}\mathbf{g}$. However, the convergence speed of the
CG method relies on the smallest eigenvalue of the positive definite matrix
$\mathbf{H}$. Therefore, if the lower-level problem is not well-conditioned,
the CG method can be slow. This approach has also been employed in the context
of model-agnostic meta-learning (MAML) [29] and adversarially robust training
[22].

Second, the WoodFisher approximation [69] expresses the Hessian as a
recurrence of a rank-one modified Hessian estimate and calls the Woodbury
matrix identity to compute the inverse of a rank-one modification to the given
Hessian matrix. The one-shot WoodFisher approximation is equivalent to the
quasi-Newton approximation [71],
$\mathbf{H}\approx\mathbf{v}\mathbf{v}^{T}+\gamma\mathbf{I}$, where $\gamma>0$
is the dampening term to render the invertibility of $\mathbf{H}$. We can then
readily obtain the inverse-Hessian vector product by the Woodbury matrix
identity
$\mathbf{H}^{-1}\mathbf{v}=\gamma^{-1}\mathbf{v}+\frac{\gamma^{-2}\mathbf{v}\mathbf{v}^{\top}}{1+\gamma^{-1}\mathbf{v}^{\top}\mathbf{v}}\mathbf{v}$.
Further, an iterative WoodFisher approximation for Hessian inverse proposed in
[69] enhances the estimation accuracy.

Third, one may utilize a Neumann-series approximation to estimate the inverse-
Hessian directly by the approximation111Assuming $\mathbf{H}$ is normalized to
ensure $\|\mathbf{H}\|\leq 1$.
$\mathbf{H}^{-1}\approx\sum_{i=0}^{K}[\mathbf{I}-\mathbf{H}]^{i}$. Note that
as $K\to\infty$ the approximation becomes more and more accurate [72]. This
technique is popular for approximating the inverse-Hessian in a stochastic
setting wherein the upper- and lower-level objectives are accessed via a
stochastic oracle [66, 70]. Here, we briefly describe the procedure to
approximate the inverse-Hessian stochastically using the Neumann-series
method. Let us choose $k$ uniformly randomly from the set
$\\{0,1,\ldots,K-1\\}$, access batch samples of $g(\bm{\theta},\bm{\phi})$
denoted by $\\{g(\bm{\theta},\bm{\phi};\zeta_{k})\\}_{k=1}^{k}$, and compute:

$\displaystyle\mathbf{H}^{-1}\approx\frac{k}{L_{g}}\prod_{i=1}^{k}\Big{(}I-{\nabla^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta},\bm{\phi};\zeta_{i})}/{L_{g}}\Big{)}$
(7)

where $L_{g}$ is the Lipschitz-smoothness constant of
$g(\bm{\theta},\bm{\phi};\zeta_{k})$. The above procedure requires the
computation of $k$ stochastic Hessians and their products. Importantly, this
estimator is a biased estimator of the inverse-Hessian with the bias
decreasing exponentially with $K$ [66, Lem. 3.2].

Finally, to ensure the local convexity, some quadratic regularization term is
usually added to the lower-level problem in BLO [22, 44]. This modifies (LU-
BLO) to

$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}\in\mathcal{U}}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta})),\leavevmode\nobreak\
\leavevmode\nobreak\ \operatorname*{\text{subject to}}\leavevmode\nobreak\
\leavevmode\nobreak\
\displaystyle\bm{\phi}^{*}(\bm{\theta})=\operatorname*{arg\,min}_{\bm{\phi}\in\mathbb{R}^{n}}g(\bm{\theta},\bm{\phi})+\frac{\lambda}{2}\|\bm{\phi}\|_{2}^{2},$
(8)

where recall that $\lambda>0$ is a regularization parameter. In this context,
the Hessian-free simplification is usually adopted, which assumes
$\widetilde{\nabla}^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))=\mathbf{0}$.
This assumption can be reasonable when the lower-level objective function $g$
involves deep model training. For instance, in the case of a neural network
with ReLU activation, the decision boundary is piece-wise linear in a tropical
hyper-surface, leading to an approximate Hessian of zeros [73]. This Hessian-
free simplification has been used for pruning deep neural networks [44]. Thus,
the Hessian matrix of the entire lower-level objective function in (8) will be
simplified to $\mathbf{H}\approx\lambda\mathbf{I}$.

####  III-A3 Extension to lower-level constrained BLO

Unlike the previous subsection, it turns out that when including constraints
to the lower-level problem, the IG no longer has the closed-form expression
because the stationary condition in (4) does not hold anymore. To see the
impact of having constraints (even linear ones) on the lower-level problem, we
present the following example, where the gradient $df/d\bm{\theta}$ is not
rigorously defined.

Example 2: (LC-BLO) can still be non-differentiable [74] Consider the
following special case of (LC-BLO), where the lower-level objective is
strongly convex in both scalar variable $\theta$ and $\phi$, the upper-level
is linear, and both levels are subject to linear constraints:
$\displaystyle\operatorname*{\text{minimize}}_{\theta\in[0,1]}\,\,\theta+\phi^{\ast}(\theta)\qquad\operatorname*{\text{subject
to}}\,\,\phi^{\ast}(\theta)\in\operatorname*{arg\,min}_{1/2\leq\phi\leq
1}(\theta-\phi)^{2}$ It follows that $\phi^{\ast}(\theta)=1/2$, for
$\theta\leq 1/2$, and $\phi^{\ast}(\theta)=\theta$, for $\theta>1/2$. We
notice that at the point $\theta=1/2$ the mapping $\phi^{\ast}(\theta)$ is
continuous, but not differentiable. As a result, the outer function
$\theta+\phi^{\ast}(\theta)$ is non-differentiable.

An immediate question is, can we still leverage the IF-based approach for this
subclass of problems? It turns out that if we make some additional assumptions
on the matrix $\mathbf{C}$ in the constraint set of (LC-BLO), one can still
apply the Implicit Function Theorem to the Karush-Kuhn-Tucker (KKT) condition
of the lower level problem to calculate the IG [22, 74]. It is also important
to note that IF-based approaches are typically not suitable for handling
general non-linear constraints in the lower-level problem. For problems with
complex constraints, value function or penalty-based approaches are often
employed [75, 76].

###  III-B The Gradient Unrolling (GU)-based Approach

GU-based approach is another class of popular algorithms for solving BLO
problems in practice [77, 78, 57]. Different from the IF-based framework, it
employs an unrolled lower-level optimizer as an intermediary step to connect
the lower-level solution with the upper-level optimization process. The
automatic differentiation (AD) technique is then used to compute gradients
w.r.t. the upper-level optimization variable $\bm{\theta}$. Consequently, the
computation of IG in GU is dependent on the choice of the lower-level
optimizer, and it no longer uses the implicit function-based expression
(4)-(5).

####  III-B1 GU-based approach for unconstrained BLO

In particular, the GU-based approach approximates
$\bm{\phi}^{\ast}(\bm{\theta})$ by running a given algorithm for a fixed
number of iterations, and then inserting the entire trajectory into the upper-
level objective; see Algorithm 2 for an illustration of the idea.

Algorithm 2: GU-based approach for solving (LU-BLO) Given initialization
$\bm{\phi}_{0}$ and $\bm{\theta}_{0}$, and iteration numbers $K$ and $T$; Let
$h(\cdot):\mathcal{U}\times\mathcal{C}\to\mathcal{C}$ denote one step of a
given algorithm, which takes both $\bm{\theta}$ and $\bm{\phi}$ as input, and
outputs an updated $\bm{\phi}$. At iteration $t\geq 0$, $\bullet$ Lower-level
optimization by $K$-step optimization:
$\displaystyle\begin{array}[]{ll}&\bm{\phi}_{k}=q(\bm{\theta}_{t},\bm{\phi}_{k-1}),\quad
k=1,\cdots K.\end{array}$ (10) Define
$\widetilde{\bm{\phi}}(\bm{\theta}_{t})\mathrel{\mathop{:}}=\bm{\phi}_{K}=q(\bm{\theta}_{t},q(\bm{\theta}_{t},\cdots,q(\bm{\theta}_{t},\bm{\phi}_{0})))$;
$\bullet$ Upper-level Optimization: Leverage AD to compute the approximated
gradient
$\widetilde{\nabla}f(\bm{\theta}_{t},\widetilde{\phi}(\bm{\theta}_{t}))\mathrel{\mathop{:}}=\frac{df(\bm{\theta}_{t},q(\bm{\theta}_{t},q(\bm{\theta}_{t},\cdots,q(\bm{\theta}_{t},\bm{\phi}_{0}))))}{d\bm{\theta}}$,
and use this gradient to update $\bm{\theta}_{t}$.

To see the difference between the AD and the IF-based approaches, let us
consider the simple case where $h(\cdot)$ is the gradient mapping,
$q(\bm{\theta}_{t},\bm{\phi}_{k-1})=\bm{\phi}_{k-1}-\beta\times\nabla_{\bm{\phi}}g(\bm{\theta}_{t},\bm{\phi}_{k-1})$
(for some constant stepsize $\beta>0$), and $K=1$ (i.e., a single-step
gradient descent step is performed for lower-level optimization). Further
assume that $\bm{\phi}_{0}$ is independent of $\bm{\theta}$, then the closed-
form expression of IG can be written as:

$\displaystyle\displaystyle\frac{d\widetilde{\bm{\phi}}(\bm{\theta}_{t})}{d\bm{\theta}}=\displaystyle\frac{d[\bm{\phi}_{0}-\beta\times\nabla_{\bm{\phi}}g(\bm{\theta}_{t},\bm{\phi}_{0})]}{d\bm{\theta}}=-\beta\nabla_{\bm{\theta},\bm{\phi}}^{2}g(\bm{\theta}_{t},\bm{\phi}_{0}).$
(11)

In some sense, the above computation is simpler than the computation of IG (5)
in the IF-based approach, since the Hessian inverse is no longer needed.
However, things can get much more complicated very quickly, as the total
number of inner iterations $K$ increases. Suppose that $K=2$, we obtain:

$\displaystyle\displaystyle\frac{d\widetilde{\bm{\phi}}(\bm{\theta}_{t})}{d\bm{\theta}}=\displaystyle\frac{d[\bm{\phi}_{1}-\beta\nabla_{\bm{\phi}}g(\bm{\theta}_{t},\bm{\phi}_{1})]}{d\bm{\theta}}=-\beta\left[\mathbf{I}+\beta\times\nabla^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta}_{t},\bm{\phi}_{1})\right]\nabla_{\bm{\theta},\bm{\phi}}^{2}g(\bm{\theta}_{t},\bm{\phi}_{0}).$
(12)

Clearly, the Hessian inverse is still not needed, but as the number of
unrolling steps increases, much higher computational and memory requirements
will be involved.

####  III-B2 Practical considerations

When the unrolling step $K$ becomes too large or the problem scale itself is
computationally expensive for GU, manual unrolling becomes necessary to save
memory costs and reduce computational overhead. Various GU approaches have
been proposed to achieve this goal efficiently. Notable techniques include
forward gradient unrolling (FGU) [47], backward gradient unrolling (BGU) [78,
57, 47, 46], and truncated gradient unrolling (TGU) [79, 57, 26]. In FGU,
unrolling is performed iteratively, and in the final step $K$, the Jacobi of
$\phi_{K}$ w.r.t. $\theta$ is given by:

$\displaystyle\underbrace{\frac{d\bm{\phi}_{K}}{d\bm{\theta}}}_{\mathbf{Z}_{K}}=\underbrace{\frac{\partial\bm{\phi}_{K}}{\partial\bm{\phi}_{K-1}}}_{\mathbf{A}_{K}}\underbrace{\frac{d\bm{\phi}_{K-1}}{d\bm{\theta}}}_{\mathbf{Z}_{K-1}}+\underbrace{\frac{\partial\bm{\phi}_{K}}{\partial\bm{\theta}}}_{\mathbf{B}_{K}}=\cdots=\sum_{i=1}^{K}\left(\prod_{t=i+1}^{K}\mathbf{A}_{t}\right)\mathbf{B}_{i}+\left(\prod_{t=1}^{K}\cdot\mathbf{A}_{t}\right)\mathbf{Z}_{0}.$
(13)

We also assume that $\bm{\phi}_{0}$ is independent of $\bm{\theta}$, which
implies $\frac{d\bm{\phi}_{0}}{d\bm{\theta}}=0$. Consequently, the expression
(13) can be rewritten as the following iterative form:

$\displaystyle\mathbf{Z}_{k}=\mathbf{A}_{k}\mathbf{Z}_{k-1}+\mathbf{B}_{k},\leavevmode\nobreak\
k=1,2,\cdots,K.$ (FGU)

Both $\mathbf{A}_{k}$ and $\mathbf{B}_{k}$ will be calculated along with the
$k$-th lower-level step $\bm{\phi}_{k}=q(\bm{\theta}_{t},\bm{\phi}_{k-1})$ and
will be discarded immediately after $\mathbf{Z}_{k}$ is obtained. Such an
iterative nature of (FGU) makes it particularly suitable for scenarios that
involve a large number of unrolling steps $K$, as the memory cost of
calculating $\mathbf{A}_{k}$ and $\mathbf{B}_{k}$ only involves any gradient
flow generated within the $k$th step. However, (FGU) requires keeping track of
the matrices $\mathbf{A}_{k}$, $\mathbf{B}_{k}$, and $\mathbf{Z}_{k-1}$.
Hence, it may not be suitable for problems with high-dimensional variables
$\bm{\theta}$ and $\bm{\phi}$.

To achieve more efficient computations when $\bm{\theta}$ and $\bm{\phi}$ are
of large scale, BGU is introduced, which eliminates the need for storing any
intermediate matrices [78, 57, 47, 46]. BGU explores the calculation of the IG
following (3):

$\displaystyle\frac{df(\bm{\theta},\bm{\phi}_{K})}{d\bm{\theta}}=\underbrace{\frac{\partial
f(\bm{\theta},\bm{\phi}_{K})}{\partial\bm{\theta}}}_{\mathbf{c}_{K}}+{\underbrace{\frac{d\bm{\phi}_{K}^{\top}}{d\bm{\theta}}}_{\mathbf{Z}_{K}^{\top}}}\underbrace{\frac{\partial
f(\bm{\theta},\bm{\phi}_{K})}{\partial\bm{\phi}_{K}}}_{\mathbf{d}_{K}}\overset{(\ref{eq:
FGU})}{=}\mathbf{c}_{K}+(\mathbf{Z}_{K-1}^{\top}\mathbf{A}_{K}^{\top}+\mathbf{B}_{K}^{\top})\mathbf{d}_{K}$
$\displaystyle=$
$\displaystyle\underbrace{(\mathbf{c}_{K}+\mathbf{B}_{K}^{\top}\mathbf{d}_{K})}_{\mathbf{c}_{K-1}}+\mathbf{Z}_{K-1}^{\top}\cdot\underbrace{\mathbf{A}_{K}^{\top}\mathbf{d}_{K}}_{\mathbf{d}_{K-1}}=\mathbf{c}_{K-1}+\mathbf{Z}_{K-1}^{\top}\mathbf{d}_{K-1}=\cdots=\mathbf{c}_{0}+\mathbf{Z}_{0}^{\top}\mathbf{d}_{0}=\mathbf{c}_{-1}.$
(14)

Instead of calculating IG explicitly like (FGU) does, (14) directly obtains
the gradient of the upper-level variable, which can be further simplified the
following recursive formulas:

$\displaystyle\begin{array}[]{l}\mathbf{c}_{k-1}=\mathbf{c}_{k}+\mathbf{B}_{k}^{\top}\mathbf{d}_{k},\leavevmode\nobreak\
k=0,1,\cdots,K,\text{ with }\mathbf{c}_{K}=\displaystyle\frac{\partial
f(\bm{\theta},\bm{\phi}_{K})}{\partial\bm{\theta}}\\\
\mathbf{d}_{k-1}=\mathbf{A}_{t}^{\top}\mathbf{d}_{k},\leavevmode\nobreak\
k=0,1,\cdots,K,\text{ with }\mathbf{d}_{K}=\displaystyle\frac{\partial
f(\bm{\theta},\bm{\phi}_{K})}{\partial\bm{\phi}_{K}}.\end{array}$ (BGU)

It can be observed that (BGU) only requires storing vectors ($\mathbf{c}_{k}$
and $\mathbf{d}_{k}$) throughout the recursion by utilizing the Jacobian-
vector product trick. As a result, BGU is particularly advantageous for
problems with large-scale variables compared to (FGU). Yet, due to the
recursive nature, BGU can be conducted only after all the $K$ lower-level
steps are finished. Thus, (BGU) needs to store all the unrolling steps
$\\{\bm{\phi}_{k}\in\mathbb{R}^{n}\\}_{k=1}^{K}$, compared with FGU.
Consequently, it may not be efficient to handle BLO as the number of unrolling
steps grows $K$.

It should also be noted that GU differs from IF as its computation relies on
the choice of the lower-level optimizer. For instance, employing sign-based
gradient descent (signGD) [80, 81] as the lower-level optimizer leads to a
computationally efficient GU variant referred to as signGD-based GU [82].
Specifically, the modified lower-level update rule (10) becomes:

$\displaystyle\widetilde{\bm{\phi}}(\bm{\theta}_{t})=\bm{\phi}_{K};\leavevmode\nobreak\
\leavevmode\nobreak\
\bm{\phi}_{k}=\bm{\phi}_{k-1}-\beta\mathrm{sign}\left(\nabla_{\bm{\phi}}g(\bm{\theta}_{t},\bm{\phi}_{k-1})\right),\leavevmode\nobreak\
k=1,2,\ldots,K,$ (signGD)

where $\mathrm{sign}(\cdot)$ denotes the element-wise sign operation, and
$\beta>0$ is a certain learning rate. Given the approximation
$\frac{d\,\mathrm{sign(\mathbf{x})}}{d\mathbf{x}}=\mathbf{0}$ (holding almost
everywhere), IG can be simplified to:

$\displaystyle\frac{d\widetilde{\bm{\phi}}(\bm{\theta}_{t})}{d\bm{\theta}}=\frac{d\bm{\phi}_{K}}{d\bm{\theta}}=\frac{d\bm{\phi}_{K-1}}{d\bm{\theta}}=\ldots=\frac{d\bm{\phi}_{0}}{d\bm{\theta}}.$
(17)

In the case that $\bm{\phi}_{0}$ is independent of $\bm{\theta}$, we can
achieve the IG-free variant of the GU approach.

###  III-C Value Function (VF)-based Approach

VF-based methods [83, 84, 85, 76, 75, 86] can also avoid the computation of
the inverse of the Hessian required in IF. The key technique is to reformulate
a standard BLO problem into a constrained single-level optimization problem.
This reformulation involves transforming the lower-level problem into an
upper-level inequality constraint. The resulting VF-based variants can then be
solved using algorithms for constrained optimization. Furthermore, in
comparison to IF and GU-based methods, the VF-based approach has broader
applicability in solving complex BLO problems. Not only can it handle lower-
level objectives with non-singleton solutions (including both convex and non-
convex objectives), but more importantly, it can accommodate lower-level
constraints as well. However, the VF-based approach has not been popular in
practical SP and ML applications yet, mainly because this approach has not
been able to deal with large-scale stochastic problems. This point will be
illustrated shortly in Sec. IV.

To understand VF-based approach, consider the following equivalent
reformulation of (LC-BLO):

$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta},\bm{\phi}\in\mathcal{C}}f(\bm{\theta},\bm{\phi}),\quad\operatorname*{\text{subject
to}}\leavevmode\nobreak\ g(\bm{\theta},\bm{\phi})\leq g^{\ast}(\bm{\theta}),$
(18)

where
$g^{\ast}(\bm{\theta})\coloneqq\min_{\bm{\phi}\in\mathcal{C}}g(\bm{\theta},\bm{\phi})$
is referred to as the VF (value function) of the lower-level problem. However,
solving (18) is highly non-trivial, partly because $g^{\ast}(\bm{\theta})$ is
not necessarily smooth, and can be non-convex. To address these challenges, a
relaxed version of problem (18) is typically considered by replacing
$g^{\ast}(\bm{\theta})$ with a smooth surrogate [83, 75]:

$\displaystyle g^{\ast}_{\mu}(\bm{\theta})$
$\displaystyle=\displaystyle\operatorname*{\text{minimize}}_{\bm{\phi}\in\mathcal{C}}g(\bm{\theta},\bm{\phi})+\frac{\mu_{1}}{2}\|\bm{\phi}\|_{2}^{2}+\mu_{2},$
(19)

where $\mu\coloneqq(\mu_{1},\mu_{2})$ is a pair of positive coefficients that
are introduced to guarantee the smoothness of $g^{\ast}_{\mu}(\bm{\theta})$
and to ensure the feasibility of the inequality constraint
$g(\bm{\theta},\bm{\phi})\leq g^{\ast}_{\mu}(\bm{\theta})$.

Given the relaxed VF formulation, one can adopt standard non-linear
optimization algorithms, such as a penalty-based algorithm to solve the
constrained optimization problem (18). For example, a log-barrier interior-
point method named BVFIM is leveraged in [83] to solve a sequence of smooth
approximated single-level problems of (18). Other related methods to solve
(18) include PDBO [75], BOME [86], and V-PBGD [76].

##  IV Theoretical Results of BLO

In this section, we examine the theoretical guarantees of various BLO
optimization methods. The section is divided into two parts. In the first
part, we discuss the convergence results of popular algorithms for solving the
(LU-BLO) problem, while the second part will focus on the more general
formulations such as (LC-BLO) and (NS-BLO). We list specific algorithms that
can handle both stochastic and deterministic BLO problems. Given the
differences in theoretical analysis between stochastic and deterministic
optimization, we also consider a generalized stochastic version of BLO, whose
upper- and lower-level objectives are

$\displaystyle\begin{array}[]{ll}\displaystyle
f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta}))\coloneqq\frac{1}{N}\sum_{i=1}^{N}f(\bm{\theta},\bm{\phi}^{*}(\bm{\theta});\xi_{i}),&\displaystyle
g(\bm{\theta},\bm{\phi})\coloneqq\frac{1}{N}\sum_{i=1}^{N}g(\bm{\theta},\bm{\phi};\zeta_{i}),\end{array}$
(21)

where $\xi_{i}\sim\mathcal{D}_{f}$ (resp. $\zeta_{i}\sim\mathcal{D}_{g}$)
represents the data sample of the upper-level (resp. lower-level) objective
from the distribution $\mathcal{D}_{f}$ (resp. $\mathcal{D}_{g}$), and $N$ is
the total number of data samples.

###  IV-A Convergence Measures of BLO

In what follows, we introduce the convergence measures utilized for evaluating
the performance of BLO algorithms. These measures serve to assess the quality
of solutions obtained by these algorithms. First, we focus on the concept of
an $\epsilon$-stationary point for (LU-BLO), which plays a crucial role in
characterizing the convergence properties of BLO when the upper-level problem
is unconstrained, i.e., $\mathcal{U}=\mathbb{R}^{m}$. In the deterministic
setting of (LU-BLO), a point $\bar{\bm{\theta}}\in\mathbb{R}^{m}$ is
considered an $\epsilon$-stationary point if it satisfies $\|\nabla
f(\bar{\bm{\theta}},\bm{\phi}^{\ast}(\bar{\bm{\theta}}))\|_{2}^{2}\leq\epsilon$.
In the stochastic setting (21), where the algorithm incorporates randomness,
the expectation is taken over the stochasticity of the algorithm. Thus, an
$\epsilon$-stationary point is defined as $\mathbb{E}\|\nabla
f(\bar{\bm{\theta}},\bm{\phi}^{\ast}(\bar{\bm{\theta}}))\|_{2}^{2}\leq\epsilon$.
It is important to note that when the upper-level problem in (LU-BLO) is
constrained, i.e., $\mathcal{U}\subset\mathbb{R}^{m}$ then the upper-level
objective $f(\bar{\bm{\theta}},\bm{\phi}^{\ast}(\bar{\bm{\theta}}))$ may not
be differentiable over $\bar{\bm{\theta}}$ in general. When solving (LC-BLO)
using the IF-based approach, if the IF is differentiable, similar measures of
stationarity as in standard optimization can be utilized. However, if the IF
is non-differentiable, alternative stationarity measures are commonly
employed. These include subgradient optimality [74], proximal gradient methods
[87, 76], and Moreau envelope techniques [74, 70, 88]. Furthermore, when
employing VF-based approaches to solve (LC-BLO), a widely used measure of
stationarity is to evaluate the convergence of the algorithms towards the KKT
points of the constrained reformulation of the BLO problem [89, 75].

In addition, the concept of oracle complexity is employed to quantify the
number of gradient evaluations needed to obtain an $\epsilon$-stationary
solution, as defined earlier. We denote $\mathcal{G}(f,\epsilon)$ (resp.
$\mathcal{G}(g,\epsilon)$) as the total number of (stochastic) gradients of
$f$ (resp. $g$) evaluated to achieve an $\epsilon$-stationary solution. This
measure provides insights into the computational requirements of BLO
algorithms and their scalability with respect to the problem size and desired
solution accuracy.

###  IV-B Convergence Guarantees for LU-BLO

As discussed in Sec. III, solving (LU-BLO) requires computing the Hessian or
its inverse. However, the stochastic formulation (21) leads to some
challenges, especially for convergence analysis. For example, in the case of
solving the stochastic BLO problem using Algorithm 1 (IF), the gradient
estimates would be replaced with appropriate stochastic gradient estimates for
both the upper and lower-level updates. However, obtaining an unbiased
estimator for the Hessian inverse term in IG (as defined in (5)) is
challenging. To overcome this challenge, a biased stochastic gradient
estimator based on Neumann-series approximation as discussed in Sec. III-A2
has been used in [66, 70]. We note that the bias of the estimator can be
easily controlled by choosing a larger batch to compute the Hessian of the
lower-level objective [66, Lem. 3.2]. Moreover, we point out that stochastic
CG (discussed in Sec. III-A2) can also be utilized to obtain the inverse-
Hessian gradient vector product to approximate the IG [90, 91].

In addition, a key design choice for BLO algorithms is whether the inner
problem is solved accurately or not. In a single-loop algorithm, one only
performs a fixed number of steps for the lower-level updates before every
upper-level updates [70, 92, 93], while in a double-loop, many lower-level
updates are carried out to obtain a very accurate approximation of
$\bm{\phi}^{\ast}(\bm{\theta})$ [66, 91, 90]. Typically, the former is simpler
to implement in practice, while the latter is easier to analyze since the
error caused by approximating $\bm{\phi}^{\ast}(\bm{\theta})$ can be well
managed. In addition, the stochastic descent direction to solve both (or
either) of upper and lower-level problems can be constructed using either
vanilla stochastic gradient descent (SGD) [94] or variance reduced (VR)
algorithms [95, 96]. Specifically, it is well known that VR-based algorithms
can lead to improved theoretical convergence of stochastic algorithms compared
to vanilla SGD-based algorithms to solve standard optimization problems. The
VR-based algorithms accomplish this improved convergence by computing
additional stochastic gradients on optimization variables computed in
consecutive iterations [95, 96]. Similar behavior is observed in solving BLO
algorithms using VR-based gradient constructions, as discussed next.

Algorithm 3: (S)GD and VR for solving stochastic (LU-BLO) and (LC-BLO) Given
initialization $\bm{\theta}_{0}$ and iteration number $T$; Iteration $t\geq 0$
yields: $\bullet$ Lower-level optimization: Given $\bm{\theta}_{t}$, call SGD
(or VR based on both $\bm{\theta}_{t}$ and $\bm{\theta}_{t-1}$) to obtain a
lower-level solution $\widetilde{\bm{\phi}}(\bm{\theta}_{t})$; $\bullet$
Approximation: Given $\widetilde{\bm{\phi}}(\bm{\theta}_{t})$, compute a
stochastic gradient estimate of (6) as follows: – Estimate stochastic versions
of
$\nabla_{\bm{\theta}}f(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))$,
$\nabla_{\bm{\phi}}f(\bm{\theta}_{t},\\!\widetilde{\bm{\phi}}(\bm{\theta}_{t}))$,
$\widetilde{\nabla}^{2}_{\bm{\theta},\bm{\phi}}g(\bm{\theta}_{t},\\!\widetilde{\bm{\phi}}(\bm{\theta}_{t}))$;
– Approximate Hessian inverse,
$\widetilde{\nabla}^{2}_{\bm{\phi},\bm{\phi}}g(\bm{\theta}_{t},\widetilde{\bm{\phi}}(\bm{\theta}_{t}))^{-1}$,
following (7) in Sec. III-A2 – Obtain stochastic estimate of (6) and construct
a descent direction ${\tilde{\nabla}f(\bm{\theta}_{t};\tilde{\xi})}$ for
$\bm{\theta}$ $\bullet$ Upper-level optimization: Call SGD (or VR) to update
$\bm{\theta}_{t}$ using ${\tilde{\nabla}f(\bm{\theta}_{t};\tilde{\xi})}$ (or
${\tilde{\nabla}f(\bm{\theta}_{t};\tilde{\xi})}$ and
${\tilde{\nabla}f(\bm{\theta}_{t-1};\tilde{\xi})}$).

In Algorithm 3, we provide a generic stochastic algorithm to solve BLO
problems using the IF-based approach. As pointed out earlier, for double-loop
algorithms $\widetilde{\bm{\phi}}(\bm{\theta}_{t})$ will approximate
${\bm{\phi}^{\ast}}(\bm{\theta}_{t})$ closely, while for a single-loop
algorithm $\widetilde{\bm{\phi}}(\bm{\theta}_{t})$ could be given by a crude
approximation of ${\bm{\phi}^{\ast}}(\bm{\theta}_{t})$. Table I provides a
summary of the oracle complexities of existing BLO algorithms for achieving an
$\epsilon$-stationary point in solving problem (LU-BLO). The theoretical
results are categorized based on three algorithmic families: stochastic BLO
(21), VR-based BLO, and deterministic BLO. The convergence performance is
evaluated in terms of the oracle complexities $\mathcal{G}(f,\epsilon)$ and
$\mathcal{G}(g,\epsilon)$, as introduced in Sec. IV-A. In Table I, we also
illustrate the BLO solver employed by different optimization principles (i.e.,
IF, GU, or VF) and the algorithmic design choices, such as double versus
single-loop.

Table I: Summary of convergence results of representative BLO algorithms for
solving (LU-BLO). The convergence performance is measured by the oracle
complexity $\mathcal{G}(f,\epsilon)$ and $\mathcal{G}(g,\epsilon)$ for upper-
and lower-level optimization, respectively. ‘Double’ or ‘Single’ refers to the
choice of the algorithm design, either double loop or single loop.

Stochastic BLO | Variance Reduced BLO | Deterministic BLO
---|---|---
Algorithm | Principle | Design | Oracle Complexity | Algorithm | Principle | Design | Oracle Complexity | Algorithm | Principle | Design | Oracle Complexity
$\mathcal{G}(f,\epsilon)$ | $\mathcal{G}(g,\epsilon)$ | $\mathcal{G}(f,\epsilon)$ | $\mathcal{G}(g,\epsilon)$ | $\mathcal{G}(f,\epsilon)$ | $\mathcal{G}(g,\epsilon)$
BSA [66] | IF | Double | $\mathcal{O}(1/\epsilon^{2})$ | $\mathcal{O}(1/\epsilon^{3})$ |  STABLE [88] | IF | Single | $\mathcal{O}(1/\epsilon^{2})$ | $\mathcal{O}(1/\epsilon^{2})$ |  BA [66] | IF | Double | $\mathcal{O}(1/\epsilon)$ | $\mathcal{O}(1/\epsilon^{5/4})$
TTSA [70] | IF | Single | $\mathcal{O}(1/\epsilon^{5/2})$ | $\mathcal{O}(1/\epsilon^{5/2})$ |  SUSTAIN [92] | IF | Single | $\mathcal{O}(1/\epsilon^{3/2})$ | $\mathcal{O}(1/\epsilon^{3/2})$ |  AID-BiO [91] | IF | Single | $\mathcal{O}(1/\epsilon)$ | $\mathcal{O}(1/\epsilon)$
stocBiO [91] | IF | Double | $\mathcal{O}(1/\epsilon^{2})$ | $\mathcal{O}(1/\epsilon^{2})$ |  VRBO [97] | IF | Double | $\mathcal{O}(1/\epsilon^{3/2})$ | $\mathcal{O}(1/\epsilon^{3/2})$ |  ITD-BiO [91] | IF | Double | $\mathcal{O}(1/\epsilon)$ | $\mathcal{O}(1/\epsilon)$
ALSET [93] | IF | Single | $\mathcal{O}(1/\epsilon^{2})$ | $\mathcal{O}(1/\epsilon^{2})$ |  SABA [98] | IF | Double | $\mathcal{O}(N^{2/3}/\epsilon)$ | $\mathcal{O}(N^{2/3}/\epsilon)$ |  MSTSA [99] | IF | Single | $\mathcal{O}(1/\epsilon)$ | $\mathcal{O}(1/\epsilon)$
F2SA [100] | VF | Single | $\mathcal{O}(1/\epsilon^{7/2})$ | $\mathcal{O}(1/\epsilon^{7/2})$ |  F3SA [100] | VF | Single | $\mathcal{O}(1/\epsilon^{5/2})$ | $\mathcal{O}(1/\epsilon^{5/2})$ |  K-RMD [57] | GU | Double | $\mathcal{O}(1/\epsilon^{2})$ | $\mathcal{O}(K/\epsilon^{2})$
AmIGO [90] | IF | Double | $\mathcal{O}(1/\epsilon^{2})$ | $\mathcal{O}(1/\epsilon^{2})$ |  SBFW [101] | IF | Single | $\mathcal{O}(1/\epsilon^{4})$ | $\mathcal{O}(1/\epsilon^{2})$ |  FGU/BGU [46] | GU | Double | N/A | N/A

##### Stochastic BLO

This set of algorithms updates the lower- and upper-level variables using SGD,
following the IF or GU optimization principle. The stochastic gradient
estimates are evaluated as discussed in Algorithm 3 while for the
deterministic setting, the gradients are approximated using the techniques
discussed in Sec. III. Bi-level Stochastic Approximation (BSA) [66] was the
first algorithm that offered finite-time convergence guarantees for solving
unconstrained stochastic BLO problems. It employed a double-loop algorithm
where the lower-level variable is iteratively estimated with multiple SGD
updates, resulting in a larger oracle complexity of
$\mathcal{O}(1/\epsilon^{3})$ for the inner-level optimization compared to
$\mathcal{O}(1/\epsilon^{2})$ for the upper-level optimization. Two-Timescale
Stochastic Approximation (TTSA) [70], a fully single-loop algorithm (with
projected SGD update for upper-level constrained optimization) improved the
lower-level oracle complexity to $\mathcal{O}(1/\epsilon^{5/2})$, however, at
the cost of worsening the upper-level oracle complexity to
$\mathcal{O}(1/\epsilon^{5/2})$. More recently, Stochastic Bi-level Optimizer
(stocBiO) [91], Stochastic Bi-level Algorithm (SOBA) [98], and Alternating
Stochastic Gradient Descent (ALSET) [93] algorithms are developed to achieve
$\mathcal{O}(1/\epsilon^{2})$ complexity for both the upper- and lower-level
optimization. Note that stocBiO [91] requires a batch size of
$\mathcal{O}(1/\epsilon)$ to achieve this complexity while ALSET [93] and SOBA
[98] rely on only $\mathcal{O}(1)$ batches to achieve the same performance. In
[100], the authors developed Fully First-order Stochastic Approximation
(F2SA), a VF-based algorithm to solve (LU-BLO). The algorithm achieved an
oracle complexity of $\mathcal{O}(1/\epsilon^{7/2})$ for both upper and lower-
level objectives while circumventing the need to compute Hessians (or Hessian
vector products) during the execution of the algorithm.

##### Variance-reduced (VR) stochastic BLO

Several VR algorithms have been proposed to improve the performance of vanilla
stochastic algorithms by computing additional stochastic gradients in each
iteration (see Algorithm 3). Examples including Single-Timescale Stochastic
Bilevel Optimization (STABLE) [88] and Momentum-Assisted Single-Timescale
Stochastic Approximation (MSTSA) [99] utilized momentum-based variance
reduction techniques [95, 102] for upper-level optimization, improving the
performance of TTSA [70] and BSA [66] and achieving a oracle complexity of
$\mathcal{O}(1/\epsilon^{2})$ for both upper- and lower-level objectives.
Single-timescale Double-momentum Stochastic Approximation (SUSTAIN) [92],
Momentum-Based Recursive Bilevel Optimizer (MRBO) [97], and Stochastic
Variance-Reduced Bilevel Method (SVRB) [103] further improved the performance
by applying variance reduction to both upper- and lower-level optimization,
achieving the oracle complexity of $\mathcal{O}(1/\epsilon^{3/2})$. Variance-
Reduced Bilevel Optimizer (VRBO) [97] employed Stochastic Path-Integrated
Differential Estimator (SPIDER) [96], a double-loop variance reduced gradient
estimator, to achieve the same complexity. More recently, Stochastic Average
Bilevel Algorithm (SABA) [98] was developed, applying SAGA (an incremental
gradient estimator) [104] achieved an oracle complexity of
$\mathcal{O}(N^{2/3}/\epsilon)$, where $N$ is the number of empirical data
points. In addition to F2SA, the authors of [100] also proposed Faster Fully
First-order Stochastic Approximation (F3SA) that improved on the oracle
complexity of F2SA. F3SA utilized momentum-based variance reduction to achieve
an oracle complexity of $\mathcal{O}(1/\epsilon^{5/2})$ for both upper and
lower-level objectives.

##### Deterministic BLO

Popular approaches for solving deterministic BLO problems include Iterative
Differentiation (ITD) and Approximate Implicit Differentiation (AID) methods.
ITD-based approaches, proposed in [47, 46], established asymptotic convergence
guarantees and were extended to TGU (truncated gradient unrolling) in [57],
achieving oracle complexities of $\mathcal{O}(1/\epsilon^{2})$. Improved
guarantees were later shown in [91], achieving oracle complexities of
$\mathcal{O}(1/\epsilon)$, comparable to solving a deterministic single-level
optimization problem. AID-based approaches for solving deterministic BLO
problems include Bi-level Approximation (BA) [66], Approximate Implicit
Differentiation Bi-level optimizer (AID-BiO) [91], and MSTSA [99]. BA, a
double-loop algorithm, was the first to establish convergence guarantees for
deterministic BLO using AID. BA achieved an oracle complexity of
$\mathcal{O}(1/\epsilon)$ for the upper-level optimization and
$\mathcal{O}(1/\epsilon^{5/4})$ for the lower-level optimization. The
performance of BA was improved in AID-BiO [91] and MSTSA [99], which achieved
an oracle complexity of $\mathcal{O}(1/\epsilon)$ for both upper- and lower-
level objectives. Please refer to Table I for a summary of the discussed
approaches.

###  IV-C Convergence Guarantees for (LC-BLO) and (NS-BLO)

Obtaining convergence guarantees for BLO algorithms become more challenging
when solving more complex problems, such as those involving lower-level
constraints in (LC-BLO) or non-singleton lower-level solutions in (NS-BLO). In
the previous section, the majority of the algorithms discussed employed an IF-
based approach to solve the BLO problem. However, in this section, only the
algorithms designed to solve (LC-BLO) utilize IF-based approaches. It is also
worth mentioning that IF-based approaches are not applicable for solving (NS-
BLO) (or (LC-BLO) with general constraints) due to the inapplicability of the
implicit function theorem in this context. Instead, standard approaches to
solve these more complex problems include interior-point methods [83], primal-
dual methods [75], dynamic barrier gradient descent [86], and penalty based
gradient descent [76]. As pointed out earlier, a major drawback of these
algorithms is that they are exclusively developed for deterministic problems
and lack efficient implementations for stochastic formulations (21).
Consequently, they are not well-suited for large-scale SP and ML applications
which often involve learning over large volumes of data.

In the following, we present a summary of recent theoretical advancements for
solving highly complex BLO problems such as (LC-BLO) and (NS-BLO) in Table II.
We list the oracle complexities of representative methods, design principles,
and the type of constraints present in the lower-level objective function.
Note that in Table II, we list the oracle complexity for only upper-level
objective with the notion of stationarity defined according to either squared
norm of the projected gradient [76, 89, 22, 74] or KKT conditions [86, 75].

####  IV-C1 Theoretical results for (LC-BLO)

BLO problems of the form (LC-BLO) involving linear constraints of the form
$\mathcal{C}\coloneqq\\{\bm{\phi}\mid
h(\bm{\theta},\bm{\phi})\leq\mathbf{0}\\}$ with
$h(\bm{\theta},\bm{\phi})\coloneqq\mathbf{A}\bm{\phi}-\mathbf{b}$, in the
lower-level have gained popularity in both theory and practice [22, 87, 74].
Under some regularity assumptions on upper- and lower-level objectives and the
constraint set of the lower-level problem, IF-based methods can be developed
for solving such problems. These algorithms follow the same structure as the
one presented in Algorithm 1, with the key difference that the construction of
the (stochastic) gradient estimate depends on the lower-level constraints.

Algorithm 4: SIGD, an IF-based approach for (LC-BLO) Given initialization
$\bm{\theta}_{0}$ and iteration number $T$; Iteration $t\geq 0$ yields:
$\bullet$ Call Algorithm 1 and use the following procedure to compute the IG
$\bullet$ Notations: $\overline{\mathbf{A}}$ be the matrix that contains the
rows of $\mathbf{A}$ that correspond to the active constraints of inequality
$\mathbf{A}\bm{\phi}-\mathbf{b}\leq\mathbf{0}$,
$\overline{\bm{\lambda}}^{\ast}(\bm{\theta})$ is the Lagrange multipliers
vector that corresponds to the active constraints at
${\bm{\phi}}^{\ast}(\bm{\theta})$ $\displaystyle\text{IG}:\leavevmode\nobreak\
\frac{d\bm{\phi}^{\ast}(\bm{\theta})^{\top}}{d\bm{\theta}}=$
$\displaystyle\big{[}\nabla_{\bm{\phi}\bm{\phi}}^{2}g(\bm{\theta},\bm{\phi}^{\ast}(\bm{\theta}))\big{]}^{-1}\cdot\big{[}-\nabla_{\bm{\theta}\bm{\phi}}^{2}g(\bm{\theta},\bm{\phi}^{\ast}(\bm{\theta}))-\overline{\mathbf{A}}^{T}\nabla\overline{\bm{\lambda}}^{\ast}(\bm{\theta})\big{]}$
(22) $\displaystyle\nabla\overline{\bm{\lambda}}^{\ast}(\bm{\theta})=-$
$\displaystyle\big{[}\overline{\mathbf{A}}\big{[}\nabla_{\bm{\phi}\bm{\phi}}^{2}g(\bm{\theta},\bm{\phi}^{\ast}(\bm{\theta}))\big{]}^{-1}\overline{\mathbf{A}}^{\top}\big{]}^{-1}\big{[}\overline{\mathbf{A}}\big{[}\nabla_{\bm{\phi}\bm{\phi}}^{2}g(\bm{\theta},\bm{\phi}^{\ast}(\bm{\theta}))\big{]}^{-1}\nabla_{\bm{\theta}\bm{\phi}}^{2}g(\bm{\theta},\bm{\phi}^{\ast}(\bm{\theta}))\big{]}.$

In [74], the Smoothed Implicit Gradient (SIGD) approach was developed to
handle (LC-BLO) with linear inequality constraints in the lower level. SIGD is
an implicit gradient descent algorithm that ensures the differentiability of
the implicit function through perturbation-based smoothing. The SIGD algorithm
and the expression for IG used in [74] is stated in (22). SIGD guarantees
asymptotic convergence to a stationary point. A similar IF-based approach was
also utilized in [22], termed Fast Bi-level Adversarial Training (Fast-BAT).
Fast-BAT aims to enhance the robustness of deep learning models against
adversarial attacks using BLO. And it achieves an oracle complexity of
$\mathcal{O}(1/\epsilon)$ under certain smoothness assumptions. In [87], the
stochastic BLO problem with linear equality constraints in both the upper- and
lower-level problems is considered. The authors proposed an IF-based approach
by constructing an approximate stochastic implicit gradient for linearly
constrained BLO. They also proposed the Alternating Implicit Projected SGD
(AiPOD) algorithm, an alternating projection method that achieves an oracle
complexity of $\mathcal{O}(1/\epsilon^{2})$ for both upper- and lower-level
objectives. The work [75] considered BLO with general constraints in both
upper- and lower-level objectives and non-singleton lower-level solutions. It
utilized the VF-based approach developed in Bi-level Value-Function-based
Interior-point Method (BVFIM) [83] for solving (NS-BLO) and proposed a Primal-
Dual Bi-level Optimizer (PDBO), a primal-dual algorithm for solving (18) when
the value function is approximated using (19). Under the assumptions of convex
and compact constraint sets and convex lower-level objectives, PDBO achieves
an oracle complexity of $\mathcal{O}(1/\epsilon^{3/2})$. Recently, the authors
in [76] proposed Penalty-Based Bilevel Gradient Descent (PBGD) for BLO with
general constraints and non-singleton lower-level solutions. The authors
established the equivalence of BLO and its penalty-based reformulations based
on VF and KKT conditions. PBGD achieved an oracle complexity of
$\mathcal{O}(1/\epsilon^{3/2})$ for solving constrained BLO with lower-level
objectives satisfying the PL inequality. Note that the algorithms PDBO [75]
and PBGD [76] can be utilized to solve both (LC-BLO) and (NS-BLO) problems.
Next, we discuss specific algorithms for solving (NS-BLO).

Table II: Convergence results of representative algorithms for solving (LC-
BLO) and (NS-BLO) problems. Similar to Table I, convergence is measured by
oracle complexity. Other algorithmic details, including optimization
principles (IF, GU, and VF), problem setups (LC-BLO and NS-BLO), objective
function types, and lower-level constraint types.

Algorithm | Principle | Objective functions | Constraints |  Oracle Complexity
---|---|---|---|---
Upper-Level | Lower-Level
BDA [45] | GU | Strongly-Convex | Convex | No | Asymptotic
BVFIM [83] | VF | Smooth | Non-Convex | No | Asymptotic
BOME! [86] | VF | Smooth | PL | No | $\mathcal{O}(1/\epsilon^{4})$
SIGD [74] | IF | Smooth | Strongly-Convex | Linear Inequality | Asymptotic
AiPOD [87] | IF | Smooth | Strongly-Convex | Linear Equality | $\mathcal{O}(1/\epsilon^{2})$
PDBO [75] | VF | Smooth | Convex | Non-Linear | $\mathcal{O}(1/\epsilon^{3/2})$
PBGD [76] | VF | Smooth | PL | Non-Linear | $\mathcal{O}(1/\epsilon^{3/2})$

####  IV-C2 Theoretical results for (NS-BLO)

An attempt to relax the lower-level singleton assumption for the lower-level
problem was made in [45] with the introduction of Bi-level Descent Aggregation
(BDA), a bi-level descent framework for solving (NS-BLO). BDA assumes
convexity of the lower-level objective and strong convexity of the upper-level
objective with respect to $\bm{\phi}$. The framework updates the lower-level
variable, $\bm{\phi}$, using a convex combination of the upper- and lower-
level partial gradients, and then updates the upper-level variable,
$\bm{\theta}$, using standard FGU/BGU techniques. The authors established the
asymptotic convergence of BDA in [45]. In [83], the authors relaxed the
convexity assumptions on the lower and upper-level objectives in (NS-BLO) and
proposed BVFIM, a VF-based approach to solving the problem. BVFIM solves a
sequence of penalty-based reformulations of the VF problem using the interior-
point method, with the asymptotic convergence. In [86], the authors introduced
Bi-level Optimization Made Easy! (BOME!), an alternative approach to directly
solve the VF problem using a dynamic barrier gradient descent algorithm. Under
the assumption of PL inequality for the lower-level objective, BOME! achieves
a finite-time sample complexity of $\mathcal{O}(1/\epsilon^{4})$ in the worst
case.

BLO-Enabled SP and ML Applications

Table III: An overview of emerging applications of BLO in SP and ML (
indicates applications we studied).

Representative Applications |  Application Areas | Problem Description |  Selected References
---|---|---|---
Wireless resource allocation  |  SP | To allocate wireless resources optimally and maximize their utilities | [9, 10, 11, 12]
Signal demodulation  | To accurately estimate transmitted symbols from received baseband signals | [13, 14]
Channel prediction | To predict the states of a communication channel by leveraging previous observations | [16, 17, 18]
Image reconstruction | To recover images from their sparse measurements | [20, 19]
Adversarial training  |  Robust ML | To train an ML model with adversarial robustness against adversarial attacks | [21, 22, 25]
Poisoning attack generation | To generate malicious data into the training set, creating vulnerabilities/backdoors in ML models. | [105, 26, 23]
Model pruning  |  Efficient ML | To find sparse sub-networks from a dense DNN w/o generalization loss | [106, 44]
Dataset condensation | To select a subset or distill a condensed version of the training set w/o generalization loss | [42, 40, 41, 60, 107]
Model-agnostic meta-learning  |  Generalized ML | To train an ML model that can quickly adapt to new tasks using limited data | [32, 28, 33, 29, 34]
Invariant risk minimization  | To train an ML model with invariant features against distribution shift | [35, 36, 37, 38, 39]
Neural architecture search |  Automated ML | To automatically optimize the architecture of DNNs for improved performance | [48, 52, 53, 49, 50, 54]
Hyper-parameter optimization | To optimize the hyper-parameters and model selection schemes in an ML pipeline | [47, 78, 46, 57, 45]

In the following sections, we will showcase how BLO can be leveraged to obtain
state-of-the-art results for a number of key SP and ML applications, such as
wireless resource allocation (Sec. V), signal demodulation (Sec. VI),
adversarial training for robustifying ML models (Sec. VII), weight pruning for
enhancing model efficiency (Sec. VIII), and invariant representation learning
for improving domain generalization (Sec. IX). Table III summarizes a number
of emerging BLO application areas, together with some representative
references.

##  V BLO for Wireless Resource Allocation

In this section, we explore the application of BLO in wireless communications,
specifically in the context of wireless optimal resource allocation (power
control) [9, 12]. The objective is to allocate power efficiently among
multiple transmitter-receiver pairs to maximize some system-level performance.
We consider a dynamic environment where wireless channel statistics change
episodically, reflecting the variations in user behavior. This scenario is
commonly observed in cellular networks, where user density and behavior vary
throughout the day, such as more users in downtown areas during peak hours
compared to night time.

To solve the resource allocation problem, we employ a neural network trained
to predict the optimal power allocation for users based on channel
information. However, neural networks often struggle when evaluated on data
that deviate from the training distribution, which is the case in our
scenario. To address this challenge and ensure the model’s performance across
episodes, we adopt a continual learning framework [108]. We maintain a memory
set, containing a representative subset of samples encountered so far, to
facilitate adaptation to new episodes while preserving satisfactory
performance on previous ones. The model is trained not only on the current
batch of data but also using the current memory set. In what follows, we
provide a detailed description of the problem and its formulation, following
the formulation and approach presented in [12].

Formulation. Consider a dynamic wireless environment comprising $T$ episodes,
where the channel state information (CSI) statistics remain stationary within
each episode. Consider a supervised learning setting, where the $i$th data
pair $(\mathbf{h}^{(i)},\mathbf{p}^{(i)})$ consists of the CSI vector
$\mathbf{h}^{(i)}$ (the feature vector) capturing the channel characteristics,
and the corresponding optimal power allocation $\mathbf{p}^{(i)}$ across the
users (the label). We train a neural network
$\pi(\bm{\theta};\mathbf{h}^{(i)})$ on these data pairs, where $\bm{\theta}$
represents the model parameters and $\mathbf{h}^{(i)}$ serves as the network
input, with the output being the power allocation prediction. Assume that the
data arrives sequentially in multiple batches. Let $D_{t}$ denote the batch we
receive at time $t$. Assume that there is a fixed-size memory set $M_{t}$
available, which stores representative historical data to be combined with
$D_{t}$ for training, and it is updated when a new batch arrives. The
performance of a power allocation scheme $\mathbf{p}$ (for a given CSI
$\mathbf{h}$) is measured by the weighted sum-rate loss function
$R(\mathbf{p};\mathbf{h})$ [12, eq. (1)].

At each time $t$, our problem involves two tasks. The first task is to train
the neural network on a (weighted) set of training samples, aiming to find the
optimal model parameter $\bm{\theta}$. The second task is to select the most
representative subset from the available training data, which includes the
samples in memory $M_{t}$ and the current data batch $D_{t}$ (denoted as
$M_{t}\cup D_{t}$). This selected subset will be then used for training, as
well as to form the new memory to be used in the next time $t+1$. To achieve
this, we introduce the variable $\bm{\lambda}$, which represents the weights
associated with each sample. Higher weights are assigned to samples that are
more representative or challenging, as determined by the system performance
metric $R(\mathbf{p};\mathbf{h})$. These weighted samples are then selected to
form the updated memory set, and they contribute more during the model
training. The idea is that by focusing on training the model on these
challenging samples, we can expect better performance on the remaining easier
samples. It is important to note that the model training process relies on the
weights obtained from the second task. Consequently, the problem described
above can be naturally formulated as the following BLO problem:

$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}&\displaystyle\sum_{i\in
M_{t}\cup
D_{t}}\lambda^{(i)}(\bm{\theta})\ell(\bm{\theta};\mathbf{h}^{(i)},\mathbf{p}^{(i)})\\\
\operatorname*{\text{subject
to}}&\bm{\lambda}(\bm{\theta})=\displaystyle\operatorname*{arg\,min}_{\mathbf{1}^{T}\bm{\lambda}=1}\sum_{i\in
M_{t}\cup
D_{t}}\lambda^{(i)}R(\pi(\bm{\theta};\mathbf{h}^{(i)});\mathbf{h}^{(i)}),,\end{array}$
(25)

where
$\ell(\bm{\theta};\mathbf{h}^{(i)},\mathbf{p}^{(i)})=\|\mathbf{p}^{(i)}-\pi(\bm{\theta};\mathbf{h}^{(i)})\|_{2}^{2}$
is the MSE loss over the $i$th sample, and $M_{t}\cup D_{t}$ are the available
training samples at time $t$. At the upper-level, we conduct supervised
training by assigning weights to the per-sample losses based on the weights
obtained from the lower-level task, which is similar to (T2). On the other
hand, at the lower-level, we tackle the task of assigning weights to the
samples based on their achieved rates, where higher weights are assigned to
samples with lower rates. This is similar to (T1).

Methods. Problem (25) is a constrained BLO problem with linear equality
constraints w.r.t. the lower-level variable $\bm{\lambda}$. Based on the
problem structure, and the optimization principles and algorithms introduced
in Sec. III and Sec. IV-C, we utilize the IF-based SIGD (Smoothed Implicit
Gradient) method [74] to solve (25). However, it should be noted that the SIGD
method assumes a strongly convex lower-level problem (see Table II). To ensure
this property, we introduce a regularization term
$(\gamma/2)\|\bm{\lambda}\|_{2}^{2}$, where $\gamma$ is the regularization
parameter as described in Sec. III-A3. As a classical baseline approach, we
also consider transfer learning (TL). In TL, when a new data batch $D_{t}$
arrives, the current model trained on data up to time $t-1$ is fine-tuned
using only $D_{t}$. This approach is motivated by the expectation that
previous knowledge can be transferred to the new environment, enabling quick
adaptation of the model. However, updating the model may result in a loss of
prior knowledge, leading to performance degradation on the prior episodes.

|  |
---|---|---
(a) | (b) | (c)
Figure 2:  The average sum-rate achieved on the combined test set (across all
four episodes) is plotted as a function of the total no. of samples used for
model training, which arrive sequentially in batches. Three experiments are
conducted with different types of channel statistics across the four episodes.
(a) The channel sequence is Rayleigh-Rician-Geometry10-Geometry50. (b) The
channel sequence is Rician-Geometry10-Rayleigh-Geometry50. (c) The channel
sequence is Geometry10-Geometry20-Geometry50-Rayleigh. The numbers after the
‘Geometry’ channel indicate the spatial arrangement of the nodes, such as a
$20m\times 20m$ area for ‘Geometry20’.

Experiment Results. We consider an experiment setting with $T=4$ episodes and
$10$ users, and with three different types of communication channels: Rayleigh
fading, Rician fading, and Geometry channels; see [12, pg. 15] for more
details. The neural network trained for power allocation consists of three
hidden layers with sizes $200,80,80$, respectively. Fig. 2 illustrates the
performance of power allocation, measured by the sum rate; in the horizontal
axis we have the total number of samples used for model training, as these
arrive sequentially in batches. Here the power allocation schemes are obtained
using the BLO-based SIGD method and the baseline approach (TL), respectively.
As we can see, the SIGD method exhibits smoother adaptation to each episode
(note that the gray line marks the boundary between episodes). It also
experiences less deterioration in performance compared to the baseline
approach. These results demonstrate the advantage of using BLO for power
allocation.

##  VI BLO for Wireless Signal Demodulation

In this section, we explore the application of BLO in wireless signal
demodulation by associating it with another BLO application, model-agnostic
meta-learning (MAML) [28]. Thus, we begin by introducing the fundamental
concepts of MAML within the framework of BLO, and then establish the
connection between MAML and signal demodulation.

###  VI-A Fundamentals of BLO in MAML

MAML, as an optimization-based meta-learning approach, has gained significant
popularity in various fields, especially in scenarios with limited resources
[28, 82, 109, 110, 29, 111, 112]. Specifically, MAML learns a meta-
initialization of optimizee variables (e.g., model weights $\bm{\theta}$) to
enable fast adaptation to new tasks when fine-tuning the model from the
learned initialization with only a few new data points [28, 113]. With $N$
learning tasks $\\{\mathcal{T}_{i}\\}_{i=1}^{N}$, (a) a fine-tuning set
$\mathcal{D}^{\texttt{tr}}_{i}$ is used in $\mathcal{T}_{i}$ for the task-
specific lower-level optimization over the task-agnostic model initialization
$\bm{\theta}$, and (b) a validation set $\mathcal{D}^{\texttt{val}}_{i}$ is
used in the upper-level optimization for evaluating the fine-tuned model
$\bm{\theta}_{i}^{*}$ from $\bm{\theta}$. Thus, MAML can be formulated as the
following BLO problem:

$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}&\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}_{(\mathcal{D}^{\texttt{tr}}_{i},\mathcal{D}^{\texttt{val}}_{i})\sim\mathcal{T}_{i}}\left[\ell_{i}(\bm{\theta}_{i}^{*}(\bm{\theta});\mathcal{D}^{\texttt{val}}_{i})\right]\\\
\operatorname*{\text{subject
to}}&\bm{\theta}_{i}^{*}(\bm{\theta})=\displaystyle\operatorname*{arg\,min}_{\bm{\theta}_{i}\in\mathbb{R}^{n}}\ell_{i}(\bm{\theta}_{i};\mathcal{D}^{\texttt{tr}}_{i},\bm{\theta}),\end{array}$
(MAML-BLO)

where $\bm{\theta}_{i}^{*}(\bm{\theta})$ signifies the fine-tuned model
weights using the initialization $\bm{\theta}$ under the task
$\mathcal{T}_{i}$, and $\ell_{i}$ denotes the model training (or validation)
loss over $\mathcal{D}^{\texttt{tr}}_{i}$ (or
$\mathcal{D}^{\texttt{val}}_{i}$) with initialization $\bm{\theta}$ (or fine-
tuned model $\bm{\theta}_{i}^{*}(\bm{\theta})$).

The MAML-BLO problem falls into the category of unconstrained BLO. Thereby,
existing works such as [28, 82, 114, 29, 115] commonly employ the IF or GU-
based approaches to solve it. The vanilla MAML algorithm [28] utilizes GU-
based BLO solver, which carries out the upper- and lower-level updates using
the following steps:

$\displaystyle\begin{array}[]{l}\text{Lower:}\leavevmode\nobreak\
\bm{\theta}_{i}^{*}(\bm{\theta})=\bm{\theta}_{i}^{(M)},\leavevmode\nobreak\
\bm{\theta}_{i}^{(m)}\leftarrow\bm{\theta}_{i}^{(m-1)}-\beta\nabla_{\bm{\theta}_{i}}\ell_{i}(\bm{\theta}_{i}^{(m-1)}),\text{
$m=1,\ldots,M$, given $\bm{\theta}_{i}^{(0)}=\bm{\theta}$}\\\
\text{Upper:}\leavevmode\nobreak\
\bm{\theta}\leftarrow\bm{\theta}-\alpha\displaystyle\frac{1}{N}\sum_{i=1}^{N}\nabla_{\bm{\theta}}\ell_{i}(\bm{\theta}_{i}^{*}(\bm{\theta})),\\\
\end{array}$ (MAML)

where $\alpha$, $\beta>0$ represent the learning rates for the SGD updates in
the upper and lower level respectively. As mentioned in Sec. III-B1, the
choice of the lower-level optimizer will greatly influence the GU-based BLO
solvers.

As described in Sec. III-B2, the usage of the sign-based SGD in the lower
level can lead the vanilla MAML to a first-order BLO solver, known as Sign-
MAML [82]. In contrast to GU-based MAML methods discussed above, the implicit
MAML (iMAML) [29] utilizes IF-based approach, where the CG (conjugate
gradient) method is used to compute the inverse-Hessian gradient product.
Compared to the vanilla MAML, iMAML shares the same lower-level updating rule,
while adopting the IF-based upper-level iteration, similar to (6). Different
from IF and GU, the first-order (FO) MAML (FO-MAML) operates by alternating
between SGD-based lower-level optimization and SGD-based upper-level
optimization, without explicitly considering the implicit gradient. We refer
to this optimization procedure as AO (Alternating Optimization).

Table IV: Performance comparison of various MAML methods on the commonly used
datasets for few-shot learning tasks Omniglog [116] and Mini-Imagenet [117].
The best performance in each setting is marked in bold. The standard
deviations are reported based on 5 random trials. Rows marked in gray indicate
BLO-enabled algorithms.

Method |  BLO (Solver) |  Test Acc. (%) |  Time (min) |  Test Acc. (%) |  Time (min) |  Test Acc. (%) |  Time (min) |  Test Acc. (%) |  Time (min)
---|---|---|---|---|---|---|---|---|---
Omniglot 20-way-1-shot | Omniglot 20-way-5-shot | Mini-ImageNet 5-way-1-shot | Mini-ImageNet 5-way-5-shot
FO-MAML [28] | AO | 90.62$\pm 0.29$ | 2.71 | 96.44$\pm 0.23$ | 2.99 | 46.39$\pm 0.44$ | 4.32 | 54.45$\pm 0.29$ | 4.93
Sign-MAML [82] | GU | 91.75$\pm 0.26$ | 2.85 | 97.79$\pm 0.18$ | 2.91 | 47.73$\pm 0.55$ | 4.51 | 55.12$\pm 0.33$ | 4.72
vanilla MAML [28] | GU | 95.65$\pm 0.25$ | 4.71 | 98.42$\pm 0.23$ | 4.94 | 48.77$\pm 0.65$ | 15.5 | 55.72$\pm 0.36$ | 15.9
iMAML [29] | IF | 95.99$\pm 0.19$ | 3.64 | 98.63$\pm 0.14$ | 3.85 | 49.31$\pm 0.41$ | 11.6 | 54.71$\pm 0.27$ | 13.3

We next demonstrate the effectiveness of BLO in MAML by applying it to
benchmark few-shot learning tasks on the Omniglot [116] and Mini-ImageNet
[117] datasets, where the generalization of the learned meta-initialization is
evaluated on the new tasks, each with only a few examples. We follow the
standard experimental setting [28, 29, 118, 82], considering 20-way-1-shot and
20-way-5-shot learning on Omniglot, and 5-way-1-shot and 5-way-5-shot learning
on Mini-ImageNet. Here $P$-way-$Q$-shot refers to training a model using a
small set of data points sampled from $P$ classes, with each class containing
$Q$ examples. Tab. IV provides an overview of the accuracy and runtime
efficiency of different methods for solving MAML-BLO, including AO-based
first-order MAML (FO-MAML) [28], GU-based Sign-MAML [82], GU-based vanilla
MAML [28], and IF-based iMAML [29]. As we can see, the vanilla MAML and iMAML
generally achieve higher testing accuracy than other MAML variants, but they
require more computation time. This is expected for their fewer implementation
assumptions, resulting in more precise meta initialization states. In the FO
optimization category, we observe that Sign-MAML outperforms FO-MAML,
demonstrating the advantage of using GU to solve BLO problems.

###  VI-B BLO for Wireless Signal Demodulation

Background and connection to MAML. We next examine the application of BLO in
the context of wireless signal demodulation through the lens of MAML. Wireless
signal demodulation is a critical process in the physical layer of a wireless
communication system, which seeks to recover the transmitted symbols
$\mathbf{s}$ from the received signals $\mathbf{y}$. Our investigation aligns
with prior research [14] that focuses on a scenario where wireless devices
frequently transmit short packets consisting of a few pilot symbols through a
varying channel. Notably, when working with any new wireless transmitter
device in the field, only a limited number of pilot symbols are available to
optimize the demodulator. At the same time, the historical data pairs for
previous devices and channel conditions can hardly transfer to new ones.
Consequently, we view this demodulation modeling as a few-shot meta-learning
problem, aiming to derive a meta initialization from historical data that is
able to quickly adapt to future devices. In this context, similar to MAML, a
meta-demodulator characterized by a learnable parameter $\bm{\theta}$ is
trained to quickly adjust to new devices with only a handful of new pilot
symbols.

Formulation. Specifically, $N$ supervised datasets are collected to train the
meta-demodulator, each associated with a specific device, which can be treated
as $N$ tasks. The $i$-th dataset with $K$ samples is given by
$\mathcal{D}_{i}=\\{(\mathbf{s}_{i}^{(k)},\mathbf{y}_{i}^{(k)})\\}_{k=1}^{K}$,
where $\mathbf{y}_{i}^{(k)}$ represents the $k$-th received signal and
$\mathbf{s}^{(k)}_{i}$ is its corresponding ground-truth transmitted symbol.
The demodulation task for each single device can be formulated as a
classification problem, as each symbol $s_{i}^{(k)}$ can only be one of the
several binary encodings, e.g., ranging from $0000$ to $1111$ following the
16-Quadrature Amplitude Modulation (16-QAM) [14]. Thus, the cross-entropy loss
between the predicted transmitted symbol
$\hat{\mathbf{s}}(\mathbf{y},\bm{\theta})$ and the true symbol $\mathbf{s}$ is
used to train the demodulator model:
$\ell_{i}(\bm{\theta};D_{i})=\mathbb{E}_{(\mathbf{s}_{i}^{(k)},\mathbf{y}_{i}^{(k)})\sim\mathcal{D}_{i}}\mathcal{L}_{\texttt{CE}}(\hat{\mathbf{s}}(\mathbf{y}_{i}^{(k)},\bm{\theta}),\mathbf{s}_{i}^{(k)})$,
where we use a multi-layer neural network as the demodulation model
$\bm{\theta}$ following [14], in order to predict the transmitted symbol
$\hat{\mathbf{s}}$ using the received signal $\mathbf{y}$. Given $N$ devices
(datasets) within the meta-training dataset, each contains $K$ data samples,
which are divided into $K^{\texttt{tr}}$ for fine-tuning $\bm{\theta}$ and
$K^{\texttt{val}}$ for validating the performance of the fine-tuned model. In
line with the notions used in (MAML-BLO), the demodulation of the $i$-th
device can be regarded as learning task $\mathcal{T}_{i}$, which consists of a
fine-tuning dataset $\mathcal{D}_{i}^{\texttt{tr}}$ and a validation dataset
$\mathcal{D}_{i}^{\texttt{val}}$. To this end, we can apply the previously
introduced MAML methods [28, 82, 29] to address the problem of wireless signal
demodulation.

Experiment Results. During the meta-training phase, we consider $N=1000$
different devices, each of which has $K^{\texttt{tr}}\in\\{1,5,10,20\\}$
training samples designated as the fine-tuning set. For the meta-testing
phase, we use another set of $100$ devices, each of which has
$K^{\texttt{tr}}$ pairs for the few-shot demodulator fine-tuning and an
additional $10000$ pilot symbol pairs for symbol classification accuracy
evaluation. Tab. V shows the average classification accuracy and running time
for different MAML methods, providing insights into how the choice of
$K^{\texttt{tr}}$ affects the performance of these methods. As we can see,
MAML and iMAML achieved the highest symbol classification accuracy, although
they required more computational time, consistent with the findings in Tab.
IV. Sign-MAML outperformed FO-MAML in the 1-shot scenario, benefiting from the
effectiveness of the GU solver. However, in scenarios with more shots, FO-MAML
can achieve the performance on par with MAML.

Table V: Performance comparison of different MAML methods for few-shot
demodulation of 16-QAM modulated wireless signals. The best performance in
each setting is indicated in bold. Standard deviations are reported based on 5
random trials. Rows marked in gray indicate BLO-enabled algorithms.

Method |  BLO (Solver) |  Test Acc. (%) |  Time (min) |  Test Acc. (%) |  Time (min) |  Test Acc. (%) |  Time (min) |  Test Acc. (%) |  Time (min)
---|---|---|---|---|---|---|---|---|---
16-way-1-shot | 16-way-5-shot | 16-way-10-shot | 16-way-20-shot
FO-MAML [28] | AO | 94.46$\pm 0.35$ | 8.59 | 97.93$\pm 0.11$ | 8.73 | 99.13$\pm 0.13$ | 8.87 | 99.17$\pm 0.04$ | 9.01
Sign-MAML [82] | GU | 96.91$\pm 0.10$ | 8.10 | 97.35$\pm 0.28$ | 8.14 | 97.42$\pm 0.30$ | 8.23 | 97.46$\pm 0.24$ | 8.32
vanilla MAML [28] | GU | 98.84$\pm 0.10$ | 13.04 | 99.49$\pm 0.03$ | 13.27 | 99.63$\pm 0.01$ | 13.29 | 99.66$\pm 0.02$ | 13.35
iMAML [29] | IF | 97.66$\pm 0.07$ | 14.14 | 98.58$\pm 0.04$ | 14.52 | 99.35$\pm 0.04$ | 14.81 | 99.43$\pm 0.02$ | 14.88

##  VII BLO for Adversarially Robust Training

Adversarial attacks demonstrate the vulnerability of ML models to small
perturbations in their inputs, posing a challenge to their reliability [119,
120, 121]. As a result, developing robust models that withstand such attacks
has become a critical focus in the field of robust ML. In what follows, we
delve into the BLO application to improve the adversarial robustness of DNNs.

The lack of adversarial robustness in ML models has prompted extensive
research on adversarial defense mechanisms [122, 123, 124, 125, 126]. While
most of existing defenses rely on min-max optimization (MMO) to minimize
worst-case training loss by incorporating a synthesized adversary, this
approach requires completely opposing objectives for the defender and
attacker. This limits its applicability to scenarios where differing
objectives are desired. Recent works [22, 25] have demonstrated the use of BLO
with customizable attack objectives to improve the efficiency and robustness
of robust model training across a wide range of adversarial attack strengths.

Formulation. We formulate the BLO-based robust training for defending against
adversarial attacks. Using variables $\bm{\theta}$ for model parameters and
$\bm{\delta}$ for input perturbation, and loss functions $\ell_{\mathrm{tr}}$
for training and $\ell_{\mathrm{atk}}$ for attacks, we define the task of
robust training (RT) as a BLO problem. The upper-level problem involves
training the model $\bm{\theta}$, while the lower-level problem optimizes
$\bm{\delta}$ for adversarial attack generation to produce worst-case input in
model training. This yields:

$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}&\mathbb{E}_{(\mathbf{x},y)\in\mathcal{D}}[\ell_{\mathrm{tr}}(\bm{\theta},\mathbf{x}+\bm{\delta}^{*}(\bm{\theta};\mathbf{x},y),y)]\\\
\operatorname*{\text{subject
to}}&\bm{\delta}^{*}(\bm{\theta};\mathbf{x},y)=\displaystyle\operatorname*{arg\,min}_{\bm{\delta}\in\mathcal{C}}\ell_{\mathrm{atk}}(\bm{\delta},\bm{\theta};\mathbf{x},y),\end{array}$
(RT-BLO)

where $(\mathbf{x},y)$ is a data pair with feature $\mathbf{x}$ and label $y$
drawn from the training dataset $\mathcal{D}$, $(\mathbf{x}+\bm{\delta})$ is
an adversarial example w.r.t. $\mathbf{x}$, and $\bm{\delta}\in\mathcal{C}$
denotes a perturbation constraint, e.g.,
$\mathcal{C}=\\{\bm{\delta}\,|\,\|\bm{\delta}\|_{\infty}\leq\epsilon,\mathbf{x}+\bm{\delta}\in[\mathbf{0},\mathbf{1}]\\}$
for $\epsilon$-tolerated $\ell_{\infty}$-norm constrained attack with
normalized input in $[\mathbf{0},\mathbf{1}]$.

If we choose $\ell_{\rm atk}=-\ell_{\rm tr}$, then problem (RT-BLO) reduces to
the MMO-based adversarial training [122]. However, the flexibility of
independently choosing the lower-level attack objective allows for a broader
range of robust training scenarios. In particular, it enables the development
of a fast robust training variant called Fast Bi-level Adversarial Training
(Fast-BAT) [22]. Fast-BAT formulation, given below, specifies the lower-level
attack generation problem of (RT-BLO) as a constrained convex quadratic
program:

$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}&\mathbb{E}_{(\mathbf{x},y)\in\mathcal{D}}[\ell_{\mathrm{CE}}(\bm{\theta},\mathbf{x}+\bm{\delta}^{*}(\bm{\theta};\mathbf{x},y),y)]\\\
\operatorname*{\text{subject
to}}&\bm{\delta}^{*}(\bm{\theta};\mathbf{x},y)=\displaystyle\operatorname*{arg\,min}_{\bm{\delta}\in\mathcal{C}}\leavevmode\nobreak\
\underbrace{-(\bm{\delta}-\bm{\delta}_{0})^{\top}\nabla_{\bm{\delta}=\bm{\delta}_{0}}\ell_{\mathrm{CE}}(\bm{\theta},\mathbf{x}+\bm{\delta},y)+(\gamma/2)\|\bm{\delta}-\bm{\delta}_{0}\|_{2}^{2}}_{\text{$\mathrel{\mathop{:}}=\ell_{\mathrm{atk}}(\bm{\theta},\bm{\delta};\mathbf{x},y)$}},\end{array}$
(Fast-BAT)

where $\ell_{\mathrm{CE}}(\bm{\theta},\mathbf{x},y)$ is the CE (cross-entropy)
loss for training model weights $\bm{\theta}$ evaluated at the data point
$(\mathbf{x},y)$, the lower-level attack objective is given by a first-order
Taylor expansion of $-\ell_{\mathrm{CE}}$ (at the linearization point
$\bm{\delta}_{0}$) plus a quadratic residual with the regularization parameter
$\gamma>0$. Since the lower-level problem becomes a convex quadratic program,
it leads to the closed-form projected gradient descent (PGD) solution
$\bm{\delta}^{*}(\bm{\theta})=\mathcal{P}_{\mathcal{C}}\left(\bm{\delta}_{0}-(1/\gamma)\nabla_{\bm{\delta}}\ell_{\mathrm{CE}}(\bm{\theta},\bm{\delta})\right)\left.\right|_{\bm{\delta}=\bm{\delta}_{0}}$.

Methods. As problem (Fast-BAT) falls into the category of (LC-BLO), we can
solve it using optimization methods introduced in Sec. IV-C. Specifically, we
consider the KKT-oriented IF approach (Sec.III-A3) as our BLO solver (see
Algorithm 4). In our experiments, we refer to this method as Fast-BAT-IF.
Furthermore, we compare Fast-BAT-IF, with non-BLO representative robust
training baselines, such as Fast-AT [125], Fast-AT-GA [127], and BackSmooth
[128].

Table VI:  Performance comparison of different robust training methods using
PreActResNet-18 on CIFAR-10 and Tiny-ImageNet datasets. The training phase
includes adversarial perturbations with two budgets: $\epsilon=8/255$ and
$16/255$ over 20 epochs. Results are presented as mean $\pm$ standard
deviation over 10 random trials. Rows marked in gray indicate BLO-enabled
algorithms.

Method |  | BLO
---
(Solver)
|  SA (%)
---
($\epsilon=8/255$)
|  RA-PGD (%)
---
($\epsilon=8/255$)
|  RA-AA (%)
---
($\epsilon=8/255$)
|  SA (%)
---
($\epsilon=16/255$)
|  RA-PGD (%)
---
($\epsilon=16/255$)
|  RA-AA (%)
---
($\epsilon=16/255$)
| Time
---
(s/epoch)
CIFAR-10, PreActResNet-18
Fast-AT [125] | N/A | 82.39$\pm 0.14$ | 45.49$\pm 0.21$ | 41.87$\pm 0.15$ | 44.15$\pm 7.27$ | 21.83$\pm 1.32$ | 12.49$\pm 0.33$ | 23.1
Fast-AT-GA [127] | 79.71$\pm 0.24$ | 47.27$\pm 0.22$ | 43.24$\pm 0.27$ | 58.29$\pm 1.32$ | 26.01$\pm 0.16$ | 17.97$\pm 0.33$ | 75.3
BackSmooth [129] | 79.31$\pm 0.17$ | 48.06$\pm 0.07$ | 44.55$\pm 0.17$ | 64.88$\pm 1.75$ | 24.18$\pm 1.37$ | 15.47$\pm 0.92$ | 60.6
Fast-BAT-IF [22] | IF | 79.97$\pm 0.12$ | 48.83$\pm 0.17$ | 45.19$\pm 0.12$ | 68.16$\pm 0.25$ | 27.69$\pm 0.16$ | 18.79$\pm 0.24$ | 61.4
Tiny-ImageNet, PreActResNet-18
Fast-AT [125] | N/A | 41.37$\pm 3.08$ | 17.05$\pm 3.25$ | 12.31$\pm 2.73$ | 31.38$\pm 0.19$ | 5.42$\pm 2.17$ | 3.13$\pm 0.24$ | 284.6
Fast-AT-GA [127] | 45.52$\pm 0.24$ | 20.39$\pm 0.19$ | 16.25$\pm 0.17$ | 29.17$\pm 0.32$ | 6.79$\pm 0.27$ | 4.27$\pm 0.15$ | 592.7
BackSmooth [129] | 45.32$\pm 0.24$ | 20.49$\pm 0.16$ | 17.13$\pm 0.14$ | 31.17$\pm 0.25$ | 7.32$\pm 0.21$ | 4.71$\pm 0.16$ | 511.3
Fast-BAT-IF [22] | IF | 45.80$\pm 0.22$ | 21.97$\pm 0.21$ | 17.64$\pm 0.15$ | 33.78$\pm 0.23$ | 8.83$\pm 0.22$ | 5.52$\pm 0.14$ | 572.4

Experiment Results. In Tab. VI, we empirically show the performance of
different robust training methods to robusify PreActResNet-18 [130] on the
CIFAR-10 [131] and Tiny-ImageNet [132] datasets. The evaluation metrics
include ① the test-time robust accuracy (RA) of the learned model against
50-step PGD attacks [122] with 10 restarts (RA-PGD) using the perturbation
budgets $\epsilon=8/255$ and $16/255$; ② RA against AutoAttack (RA-AA) [133]
in the setup similar to RA-PGD; ③ the standard accuracy (SA) of the learned
model on natural examples; ④ the time consumption required for robust
training. As observed, Fast-BAT-IF exhibits higher robustness compared to non-
BLO baselines, highlighting the effectiveness of BLO in robust training.

##  VIII BLO for Model Pruning

While over-parameterized structures are key to the improved generalization of
DNNs, they create new problems - the millions or even billions of parameters
not only increase computational costs during inference but also pose serious
deployment challenges on resource-limited devices. Thus, the problem of model
pruning arises, aiming to reduce the sizes of an ML model by identifying and
removing redundant model weights. In this section, we investigate the
application of BLO in the context of model pruning [134, 135, 44, 136].

Formulation. The study of model pruning through the lens of BLO was first
explored in [44]. Specifically, there exist two main tasks in model pruning:
pruning and re-training. Pruning involves determining the sparse pattern of
model weights, while re-training focuses on recovering model accuracy using
the remaining non-zero weights [135, 136]. To facilitate these tasks, one can
introduce the binary pruning mask variable $\mathbf{m}\in\\{0,1\\}^{m}$ and
the model weight variable $\bm{\phi}\in\mathbb{R}^{m}$, where $m$ represents
the total number of model parameters. Accordingly, the pruned model is given
by $(\mathbf{m}\odot\bm{\phi})$, where $\odot$ denotes element-wise
multiplication. To achieve a pruning ratio of $p\%$, we impose a sparsity
constraint on $\mathbf{m}$, where $\mathbf{m}\in\Omega$ and
$\Omega=\\{\mathbf{m},|,\mathbf{m}\in{0,1}^{n},\mathbf{1}^{T}\mathbf{m}\leq
k\\}$, with $k=(1-p\%)n$. Our goal is to prune the original dense model to the
targeted pruning ratio of $p\%$ and obtain the optimal sparse model
$(\mathbf{m}\odot\bm{\phi})$. To achieve this, we interpret the pruning task
(①) and the model retraining task (②) as two optimization levels, leading to
the formulation of bi-level pruning (BiP):

$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\mathbf{m}\in\Omega}\underbrace{\ell_{\mathrm{tr}}(\mathbf{m}\odot{\bm{\phi}}^{*}(\mathbf{m}))}_{\text{
\char 172: Pruning task}};\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname*{\text{subject to}}\leavevmode\nobreak\ \leavevmode\nobreak\
\underbrace{\displaystyle{\bm{\phi}}^{*}(\mathbf{m})=\operatorname*{arg\,min}_{\bm{\phi}\in\mathbb{R}^{n}}\tilde{\ell}_{\mathrm{tr}}(\mathbf{m}\odot\bm{\phi})+\frac{\gamma}{2}\|\bm{\phi}\|_{2}^{2}}_{\text{\char
173: Sparsity-fixed model retraining}},$ (BiP)

where $\ell_{\mathrm{tr}}$ and $\tilde{\ell}_{\mathrm{tr}}$ denote the
training losses under different data batches, $\mathbf{m}$ and $\bm{\phi}$ are
the upper-level and lower-level optimization variables respectively, and
${\bm{\phi}}^{*}(\mathbf{m})$ signifies the retrained model weights given the
pruning mask $\mathbf{m}$. In (BiP), the lower-level training objective was
regularized using a strongly convex regularizer
$\frac{\gamma}{2}\|\bm{\phi}\|_{2}^{2}$ like (Fast-BAT).

Methods. Since BiP is an unconstrained BLO problem, it can be solved using BLO
algorithms, e.g., IF and GU, introduced in Secs. III-A1 and III-B1. Moreover,
since the Hessian of the lower-level objective function w.r.t. model
parameters is of high dimension, we impose the Hessian-free assumption
$\nabla_{\bm{\phi},\bm{\phi}}\ell_{\mathrm{tr}}=\mathbf{0}$ to make the BLO
implementation computationally feasible. Following (5), one can then obtain
the closed form of IG [44]:
$\frac{d\bm{\phi}^{*}(\mathbf{m})}{d\mathbf{m}}=-\frac{1}{\gamma}\nabla^{2}_{\mathbf{m},\bm{\phi}}\ell_{\mathrm{tr}}(\mathbf{m}\odot\bm{\phi}^{*})$,
where $\bm{\phi}^{*}$ signifies a lower-level solution. Furthermore, the bi-
linearity of the pruning mask $\mathbf{m}$ and the model weights $\bm{\phi}$
allows us to further simplify the IG to:

$\displaystyle\displaystyle\frac{d\bm{\phi}^{*}(\mathbf{m})}{d\mathbf{m}}=-\frac{1}{\gamma}\mathrm{diag}(\nabla_{\mathbf{z}}\ell_{\mathrm{tr}}(\mathbf{z})\left|{}_{\mathbf{z}={\mathbf{m}}\odot{\bm{\phi}}^{*}}\right.),$
(34)

where the Hessian-free assumption is adopted and $\mathrm{diag}(\mathbf{a})$
denotes the diagonal matrix with $\mathbf{a}$ being the main diagonal vector.
Detailed proof can be found in [44, Sec. 3].

Experiments Results. To implement BiP, we adopt two BLO methods: IF (Sec.
III-A1) and GU (Sec. III-B1). We term the resulting BLO-inspired model pruning
approaches as BiP-IF and BiP-GU. For comparison, we also consider two
commonly-used non-BLO based pruning methods, the state-of-the-art iterative
magnitude pruning (IMP) [136] and the most efficient one-shot magnitude
pruning (OMP) [137]. We remark that the notable lottery ticket hypothesis
(LTH) [136] stated that IMP is able to identify a trainable sparse subnetwork
(known as a ‘winning ticket’), with test accuracy surprisingly on par or even
better than the original model.

(a) CIFAR-10 (b) CIFAR-100 (c) Tiny-ImageNet (d) Efficiency comparison

Figure 3: Experiment results of model pruning on different datasets under
ResNet-18. (a-c) Pruning trajectory is given by test accuracy (%) vs. sparsity
(%) under different datasets. (d) The entire time consumption vs. the pruning
ratio.

Fig. 3 illustrates the pruning accuracy and the run-time efficiency of BLO-
based pruning methods vs. non-BLO approaches across diverse image
classification datasets (including CIFAR-10, CIFAR-100, and Tiny-ImageNet)
under ResNet-18. As we can see, BiP-IF yields the best performance in all the
dataset settings evidenced in Fig. 3(a)-(c). This also the pruning recipe used
in [44]. Thanks to the closed-form expression of IG in (34), the computation
of BiP-IF is also more efficient than BiP-GU, as shown in Fig. 3(d). In
addition, BiP-GU can also provide competitive pruning accuracy to IMP and
takes less computation time than IMP. Furthermore, we observe that OMP yields
the least computation time but the worst pruning accuracy. This is not
surprising since OMP adopts a non-iterative pruning scheme to find the model’s
sparse pattern.

##  IX BLO for Invariant Representation Learning

In this section, we explore the application of BLO in improving the
generalization capabilities of ML models. Specifically, we investigate the use
of BLO for acquiring training environment-agnostic data representations
through invariant risk minimization (IRM) [35]. The goal is to enhance the
domain generalization of ML models.

Formulation. IRM [35] is proposed to acquire invariant data representations
and to enforce invariant predictions against distribution shifts. Different
from the conventional environment risk minimization (ERM)-based training, IRM
yields the BLO-like formulation: The upper-level optimization task of IRM is
to train a network backbone to capture environment-agnostic data
representations, and the lower-level optimization task is to find an invariant
prediction head (one top of the learned representation network) to produce a
global optima to all the training environments. Formally, IRM can be cast as

$\displaystyle\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}\sum_{i=1}^{E}\ell_{i}(\bm{\phi}^{*}(\bm{\theta})\circ\bm{\theta})\leavevmode\nobreak\
\leavevmode\nobreak\ \operatorname*{\text{subject to}}\leavevmode\nobreak\
\leavevmode\nobreak\
\bm{\phi}^{*}(\bm{\theta})\in\operatorname*{arg\,min}_{\bm{\phi}}\ell_{i}(\bm{\phi}\circ\bm{\theta}),$
(IRM)

where $\bm{\phi}\circ\bm{\theta}$ denotes the representation-acquisition model
$\bm{\theta}$ and the predictor $\bm{\phi}$, $\ell_{i}$ is the training loss
associated with the $i$th training environment, and $E$ is the total number of
training environments. The rationale behind (IRM) is that given the invariant
representation extractor $\bm{\theta}$, there exists an invariant predictor
$\bm{\phi}^{*}(\bm{\theta})$ which is optimal across all the training
environments.

Methods. Solving problem (IRM) is highly non-trivial since the lower-level
solution $\bm{\phi}(\bm{\theta})$ should be universal across all $E$ training
environments. To circumvent this difficulty, IRM is typically relaxed to a
single-level optimization problem, known as IRMv1 [35]:

$\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}\sum_{i=1}^{E}[\ell_{i}(\bm{\phi}^{\ast}(\bm{\theta})\circ\bm{\theta})+\gamma\|\nabla_{w|w=1.0}\ell_{i}(w\circ\bm{\theta})\|_{2}^{2}]$
(35)

where $\gamma>0$ is a regularization parameter and
$\nabla_{w|w=1.0}\ell_{i}(w\odot\bm{\theta})$ denotes the gradient of
$\ell_{i}$ w.r.t. w, computed at $w=1.0$. However, the above is restricted to
linear invariant prediction and penalizes the deviation of individual
environment losses from stationarity to approach the lower-level optimality in
IRM. Despite the practical simplicity of IRMv1, it may fail to achieve the
desired invariance [138, 139].

Beyond IRMv1, a consensus-constrained BLO method is developed in [140] to
solve problem (IRM). The key idea is to introduce $E$ auxiliary predictors
$\\{\bm{\phi}_{i}\\}$ and explicitly enforce prediction invariance by infusing
a consensus prediction constraint
$\mathcal{C}=\\{\\{\bm{\phi}_{i}\\}|\bm{\phi}_{1}=\ldots=\bm{\phi}_{E}\\}$ to
the lower-level problem of (IRM), and promote the per-environment stationarity
in its upper-level problem. This modifies (IRM) to an ordinary BLO problem:

$\displaystyle\begin{array}[]{ll}\displaystyle\operatorname*{\text{minimize}}_{\bm{\theta}}&\sum_{i=1}^{E}[\ell_{i}(\bm{\phi}_{i}^{*}(\bm{\theta})\circ\bm{\theta})+\gamma\|\nabla_{\bm{\phi}_{i}}\ell_{i}(\bm{\phi}_{i}^{*}(\bm{\theta})\circ\bm{\theta})\|_{2}^{2}]\\\
\operatorname*{\text{subject
to}}&\\{\bm{\phi}_{i}^{*}(\bm{\theta})\\}_{i=1}^{E}=\displaystyle\operatorname*{arg\,min}_{\\{\bm{\phi}_{i}\\}\in\mathcal{C}}\sum_{i=1}^{E}\ell_{i}(\bm{\phi}_{i}\circ\bm{\theta}),\leavevmode\nobreak\
\forall i\in[E],\end{array}$ (IRM-BLO)

where $\gamma>0$ is a regularization parameter, and $[E]$ denotes the integer
set $\\{1,2,\ldots,E\\}$. The advantage of converting (IRM) into the
consensus-constrained (IRM-BLO) is that projection onto the consensus
constraint yields a closed-form solution, i.e.,
$\mathcal{P}_{\mathcal{C}}(\mathbf{a})=\operatorname*{arg\,min}_{\\{\bm{\phi}_{i}\\}\in\mathcal{C}}\sum_{i=1}^{E}\|\bm{\phi}_{i}-\mathbf{a}_{i}\|_{2}^{2}=\frac{1}{E}\sum_{i}\mathbf{a}_{i}$,
where $\mathcal{P}_{\mathcal{C}}(\mathbf{a})$ denotes the projection operation
to project the point $\mathbf{a}$ onto the constraint $\mathcal{C}$. It has
been shown in [140] that problem (IRM-BLO) can be effectively solved using the
GU approach, which approximates each individual lower-level solution using
$K$-step GD unrolling together with the consensus projection. Thus, the lower-
level solution becomes

$\displaystyle\forall i,\leavevmode\nobreak\ \leavevmode\nobreak\
\bm{\phi}_{i}^{*}(\bm{\theta})\approx\frac{1}{E}\sum_{i=1}^{E}\bm{\phi}_{i}^{(K)},\leavevmode\nobreak\
\leavevmode\nobreak\
\bm{\phi}_{i}^{(k)}=\bm{\phi}_{i}^{(k-1)}-\beta\nabla_{\bm{\phi}_{i}}\ell_{i}(\bm{\phi}_{i}^{(k-1)}\circ\bm{\theta}),\text{
for }k\in[K],$ (38)

where $\beta>0$ is the lower-level learning rate. Based on the above,
automatic differentiation (AD) can be called to compute the implicit gradient
from $\bm{\phi}_{i}^{(K)}(\bm{\theta})$ to the variable $\bm{\theta}$ in the
GU process.

Table VII: Performance of different IRM training methods. The best performance
per evaluation metric is highlighted in bold and the performance of the BLO-
enabled method is marked in gray.

Method |  BLO (Solver) | Colored-MNIST | Colored-FashionMNIST
---|---|---|---
Avg. Acc. | Acc. Gap | Avg. Acc. | Acc. Gap
ERM | N/A | 49.19$\pm 1.89$ | 90.72$\pm 2.08$ | 49.77$\pm 1.71$ | 88.62$\pm 2.49$
IRM-v1 [35] | N/A | 68.33$\pm 0.31$ | 2.04$\pm 0.05$ | 68.76$\pm 0.31$ | 1.45$\pm 0.09$
IRM-Game [36] | N/A | 67.73$\pm 0.24$ | 1.67$\pm 0.14$ | 67.49$\pm 0.32$ | 1.82$\pm 0.13$
IRM-BLO [140] | GU | 69.47$\pm 0.24$ | 1.04$\pm 0.07$ | 69.43$\pm 0.21$ | 1.14$\pm 0.11$

Experiment Results. We evaluate the performance of IRM-BLO with two commonly-
used image classification datasets, Colored-MNIST [35] and Colored-
FashionMNIST [36], where spurious correlation is imposed in the datasets which
makes the conventional ERM training ineffective. To capture both the accuracy
and the variance of invariant predictions across multiple testing
environments, the average accuracy and the accuracy gap (the difference
between best and worst-case accuracy) are measured for IRM methods. Tab. VII
presents the resulting performance of IRM-BLO and compares it with that of ERM
and two IRM baselines, IRMv1 [35]and IRM-Game [36]. Note that all IRM variants
outperform ERM, which justifies the importance of IRM training to improve
model generalization across diverse environments. Within the IRM training
family, IRM-BLO outperforms others by achieving the highest average accuracy
and the smallest accuracy gap across both datasets. This superior performance
underscores the value of BLO compared to the sub-optimal design IRM-v1.

##  X Discussion

BLO is a challenging but rapidly developing subject. Despite recent progresses
discussed in this article, significant work is yet to be done to address
various challenges ranging from developing scalable algorithms, to extending
its applicability to a wider range of problems. Below we highlight several
worthwhile future directions.

$\bullet$ BLO algorithms. First, the development and analysis of BLO
algorithms for more general BLO problems, including those with complex lower-
level constraints (e.g., non-linear constraints), require additional
exploration. The current focus has predominantly been on problems with linear
constraints, and extending the BLO framework to handle non-linear constraints
is an important and challenging task; some recent work towards this direction
in the context of the min-max problems can be found in [141]. Additionally,
exploring scenarios with coupled constraints between lower and upper levels,
such as resource sharing among adversarial and normal agents, presents a
complex area that is under-explored. Second, BLO formulations with non-
singleton lower-level solutions, such as (NS-BLO), lack theoretically
grounded, scalable, and easy-to-implement algorithms. This aspect of BLO has
received less attention from the community, making it an open topic for future
investigation. Developing efficient algorithms that can handle non-singleton
lower-level solutions and provide convergence guarantees is a key research
direction. Third, beyond the scope of bi-level optimization, exploring
problems involving more than two levels, such as dataset pruning for transfer
learning, represents an exciting frontier. These multi-level problems
introduce additional complexity and challenges, requiring the design of novel
algorithms to tackle the inherent hierarchical structure effectively. Finally,
large-scale and distributed availability of data demands the development of
decentralized and federated algorithms to solve these complex BLO problems
which also presents a compelling research direction for future exploration.

$\bullet$ BLO theories. First, while significant progress has been made in
establishing theoretical guarantees for solving the basic BLO problem (LU-
BLO), more attention is needed on exploring practical settings including
(coupled) lower-level constraints, non-convex lower-level problems, and/or
black-box settings where one may not have access to upper/lower-level
parameters. These scenarios present unique challenges and complexities, making
it difficult to analyze the problem and derive theoretical guarantees.
Investigating the convergence properties and establishing theoretical
foundations for solving BLO problems under these practical settings is an
important avenue for future research. Second, when datasets become massive, it
is crucial to develop theoretically-grounded BLO algorithms that can adhere to
practical requirements. Therefore, developing theoretical frameworks and
analyzing the convergence properties of algorithms for solving large-scale BLO
problems under realistic assumptions is also an important research topic.
Third, with the discovery of the phenomenon of double descent [142],
theoretical analysis of standard ML algorithms on overparameterized neural
networks has received significant attention from the research community.
Theoretical investigation of BLO algorithms for such overparameterized
problems is certainly an interesting research direction. Finally, developing a
thorough theoretical understanding of distributed/federated algorithms for
solving complex BLO problems will certainly be at the forefront of future BLO
research.

$\bullet$ BLO applications. First, in the context of Mixture-of-Experts (MoE)
training, there is a complex interplay between the training of the gating
network that selects experts and the training of the expert-oriented pathways
used for final predictions. Exploring BLO techniques to effectively optimize
the coupling between these two training processes in MoE can lead to improved
performance and better utilization of emerging ML models like MoE. Second,
prompt learning, a key technique used in today’s foundation models, involves a
crucial coupling between prompt pattern learning and label/feature mapping
optimization. Leveraging BLO methods to model and optimize the interactions
between prompt pattern learning and label/feature mapping can enhance the
learning process and enable accurate and robust prompt generations. Third, BLO
is highly applicable in (inverse) reinforcement learning. For instance, the
Actor/Critic Algorithm can be formulated as a BLO problem, with separate
agents evaluating and optimizing the policy. In inverse reinforcement
learning, the tasks involve inferring the agents’ reward function and finding
the optimal policy based on it. Applying BLO frameworks to these scenarios
offers potential for novel insights and improved efficiency.

In summary, the interplay between the theoretical underpinnings of BLO and its
practical applications promises a fertile ground for future exploration and
innovation, pushing the boundaries of optimization theory and applications in
SP and ML.

## References

  * [1] Heinrich Von Stackelberg,  The theory of the market economy,  Oxford University Press, 1952.
  * [2] Jerome Bracken and James T. McGill,  “Mathematical programs with optimization problems in the constraints,”  Operations Research, vol. 21, no. 1, pp. 37–44, 1973.
  * [3] Jerome Bracken and James T McGill,  “Defense applications of mathematical programs with optimization problems in the constraints,”  Operations Research, vol. 22, no. 5, pp. 1086–1096, 1974.
  * [4] Jerome Bracken and JT McGill,  “Production and marketing decisions with multiple objectives in a competitive environment,”  Journal of Optimization Theory and Applications, vol. 24, no. 3, pp. 449–458, 1978.
  * [5] Benoît Colson, Patrice Marcotte, and Gilles Savard,  “An overview of bilevel optimization,”  Annals of operations research, vol. 153, no. 1, pp. 235–256, 2007\.
  * [6] Ankur Sinha, Pekka Malo, and Kalyanmoy Deb,  “A review on bilevel optimization: from classical to evolutionary approaches and applications,”  IEEE Transactions on Evolutionary Computation, vol. 22, no. 2, pp. 276–295, 2017.
  * [7] Yasmine Beck, Ivana Ljubić, and Martin Schmidt,  “A survey on bilevel optimization under uncertainty,”  European Journal of Operational Research, 2023.
  * [8] Yasmine Beck, Martin Schmidt, Johannes Thürauf, and Daniel Bienstock,  “On a computationally ill-behaved bilevel problem with a continuous and nonconvex lower level,”  arXiv preprint arXiv:2202.01033, 2022.
  * [9] Haoran Sun, Xiangyi Chen, Qingjiang Shi, Mingyi Hong, Xiao Fu, and Nicholas D Sidiropoulos,  “Learning to optimize: Training deep neural networks for interference management,”  IEEE Transactions on Signal Processing, vol. 66, no. 20, pp. 5438–5453, 2018.
  * [10] Yifan Chen, Zhiyong Li, Bo Yang, Ke Nai, and Keqin Li,  “A stackelberg game approach to multiple resources allocation and pricing in mobile edge computing,”  Future Generation Computer Systems, vol. 108, pp. 273–287, 2020\.
  * [11] Yulan Gao, Chao Yong, Zehui Xiong, Dusit Niyato, Yue Xiao, and Jun Zhao,  “A stackelberg game approach to resource allocation for intelligent reflecting surface aided communications,”  arXiv preprint arXiv:2003.06640, 2020.
  * [12] Haoran Sun, Wenqiang Pu, Xiao Fu, Tsung-Hui Chang, and Mingyi Hong,  “Learning to continuously optimize wireless resource in a dynamic environment: A bilevel optimization perspective,”  IEEE Transactions on Signal Processing, vol. 70, pp. 1900–1917, 2022.
  * [13] Osvaldo Simeone,  “A very brief introduction to machine learning with applications to communication systems,”  IEEE Transactions on Cognitive Communications and Networking, vol. 4, no. 4, pp. 648–664, 2018.
  * [14] Sangwoo Park, Hyeryung Jang, Osvaldo Simeone, and Joonhyuk Kang,  “Learning to demodulate from few pilots via offline and online meta-learning,”  IEEE Transactions on Signal Processing, vol. 69, pp. 226–239, 2020\.
  * [15] Sangwoo Park and Osvaldo Simeone,  “Predicting flat-fading channels via meta-learned closed-form linear filters and equilibrium propagation,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 8817–8821.
  * [16] Ali Abdi and Mostafa Kaveh,  “A space-time correlation model for multielement antenna systems in mobile fading channels,”  IEEE Journal on Selected Areas in communications, vol. 20, no. 3, pp. 550–560, 2002.
  * [17] Osvaldo Simeone and Umberto Spagnolini,  “Lower bound on training-based channel estimation error for frequency-selective block-fading rayleigh mimo channels,”  IEEE Transactions on Signal Processing, vol. 52, no. 11, pp. 3265–3277, 2004.
  * [18] Marcello Cicerone, Osvaldo Simeone, and Umberto Spagnolini,  “Channel estimation for mimo-ofdm systems by modal analysis/filtering,”  IEEE transactions on communications, vol. 54, no. 11, pp. 2062–2074, 2006.
  * [19] Caroline Crockett, Jeffrey A Fessler, et al.,  “Bilevel methods for image reconstruction,”  Foundations and Trends® in Signal Processing, vol. 15, no. 2-3, pp. 121–289, 2022.
  * [20] Caroline Crockett and Jeffrey A Fessler,  “Motivating bilevel approaches to filter learning: A case study,”  in 2021 IEEE International Conference on Image Processing (ICIP). IEEE, 2021, pp. 2803–2807.
  * [21] Simiao Zuo, Chen Liang, Haoming Jiang, Xiaodong Liu, Pengcheng He, Jianfeng Gao, Weizhu Chen, and Tuo Zhao,  “Adversarial regularization as stackelberg game: An unrolled optimization approach,”  arXiv preprint arXiv:2104.04886, 2021.
  * [22] Yihua Zhang, Guanhua Zhang, Prashant Khanduri, Mingyi Hong, Shiyu Chang, and Sijia Liu,  “Revisiting and advancing fast adversarial training through the lens of bi-level optimization,”  in International Conference on Machine Learning, 2022, pp. 26693–26712.
  * [23] Peng Yang, Yingjie Lao, and Ping Li,  “Robust watermarking for deep neural networks via bi-level optimization,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 14841–14850.
  * [24] Daniel Zügner, Amir Akbarnejad, and Stephan Günnemann,  “Adversarial attacks on neural networks for graph data,”  in Proceedings of the 24th ACM SIGKDD international conference on knowledge discovery & data mining, 2018, pp. 2847–2856.
  * [25] Alexander Robey, Fabian Latorre, George J Pappas, Hamed Hassani, and Volkan Cevher,  “Adversarial training should be cast as a non-zero-sum game,”  arXiv preprint arXiv:2306.11035, 2023.
  * [26] W Ronny Huang, Jonas Geiping, Liam Fowl, Gavin Taylor, and Tom Goldstein,  “Metapoison: Practical general-purpose clean-label data poisoning,”  Advances in Neural Information Processing Systems, vol. 33, pp. 12080–12091, 2020.
  * [27] Luis Muñoz-González, Battista Biggio, Ambra Demontis, Andrea Paudice, Vasin Wongrassamee, Emil C Lupu, and Fabio Roli,  “Towards poisoning of deep learning algorithms with back-gradient optimization,”  in Proceedings of the 10th ACM workshop on artificial intelligence and security, 2017, pp. 27–38.
  * [28] Chelsea Finn, Pieter Abbeel, and Sergey Levine,  “Model-agnostic meta-learning for fast adaptation of deep networks,”  in International conference on machine learning. PMLR, 2017, pp. 1126–1135.
  * [29] Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine,  “Meta-learning with implicit gradients,”  Advances in neural information processing systems, vol. 32, 2019\.
  * [30] Marcin Andrychowicz, Misha Denil, et al.,  “Learning to learn by gradient descent by gradient descent,”  in NeurIPS, 2016.
  * [31] Tianlong Chen, Xiaohan Chen, Wuyang Chen, Zhangyang Wang, Howard Heaton, Jialin Liu, and Wotao Yin,  “Learning to optimize: A primer and a benchmark,”  The Journal of Machine Learning Research, vol. 23, no. 1, pp. 8562–8620, 2022.
  * [32] Sachin Ravi and Hugo Larochelle,  “Optimization as a model for few-shot learning,”  in International conference on learning representations, 2016.
  * [33] Alex Nichol, Joshua Achiam, and John Schulman,  “On first-order meta-learning algorithms,”  arXiv preprint arXiv:1803.02999, 2018.
  * [34] Harkirat Singh Behl, Atılım Güneş Baydin, and Philip HS Torr,  “Alpha maml: Adaptive model-agnostic meta-learning,”  arXiv preprint arXiv:1905.07435, 2019.
  * [35] Martin Arjovsky, Léon Bottou, Ishaan Gulrajani, and David Lopez-Paz,  “Invariant risk minimization,”  arXiv preprint arXiv:1907.02893, 2019.
  * [36] Kartik Ahuja, Karthikeyan Shanmugam, Kush Varshney, and Amit Dhurandhar,  “Invariant risk minimization games,”  in International Conference on Machine Learning. PMLR, 2020, pp. 145–155.
  * [37] Elan Rosenfeld, Pradeep Ravikumar, and Andrej Risteski,  “The risks of invariant risk minimization,”  arXiv preprint arXiv:2010.05761, 2020.
  * [38] Yong Lin, Hanze Dong, Hao Wang, and Tong Zhang,  “Bayesian invariant risk minimization,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2022, pp. 16021–16030.
  * [39] Xiao Zhou, Yong Lin, Weizhong Zhang, and Tong Zhang,  “Sparse invariant risk minimization,”  in International Conference on Machine Learning. PMLR, 2022, pp. 27222–27244.
  * [40] Zalán Borsos, Mojmir Mutny, and Andreas Krause,  “Coresets via bilevel optimization for continual learning and streaming,”  Advances in Neural Information Processing Systems, vol. 33, pp. 14879–14890, 2020.