# First-Order Methods for
Linearly Constrained Bilevel Optimization††thanks: GK, SP, KW, and ZZ
contributed equally; authors listed in alphabetical order.

Guy Kornowski Weizmann Institute of Science<EMAIL_ADDRESS>Swati Padmanabhan Massachusetts Institute of Technology<EMAIL_ADDRESS>Kai Wang
Georgia Institute of Technology<EMAIL_ADDRESS>Zhe Zhang Purdue
University<EMAIL_ADDRESS>Suvrit Sra Massachusetts Institute of
Technology<EMAIL_ADDRESS>
###### Abstract

Algorithms for bilevel optimization often encounter Hessian computations,
which are prohibitive in high dimensions. While recent works offer first-order
methods for unconstrained bilevel problems, the constrained setting remains
relatively underexplored. We present first-order linearly constrained
optimization methods with finite-time hypergradient stationarity guarantees.
For linear equality constraints, we attain $\epsilon$-stationarity in
$\widetilde{O}(\epsilon^{-2})$ gradient oracle calls, which is nearly-optimal.
For linear inequality constraints, we attain $(\delta,\epsilon)$-Goldstein
stationarity in $\widetilde{O}(d{\delta^{-1}\epsilon^{-3}})$ gradient oracle
calls, where $d$ is the upper-level dimension. Finally, we obtain for the
linear inequality setting dimension-free rates of
$\widetilde{O}({\delta^{-1}\epsilon^{-4}})$ oracle complexity under the
additional assumption of oracle access to the optimal dual variable. Along the
way, we develop new nonsmooth nonconvex optimization methods with inexact
oracles. We verify these guarantees with preliminary numerical experiments.

##  1 Introduction

Bilevel optimization [1, 2, 3, 4], an important problem in optimization, is
defined as follows:

$\displaystyle\mbox{minimize}_{x\in\mathcal{X}}~{}F(x)\coloneqq
f(x,y^{*}(x))\text{ subject to }~{}y^{*}(x)\in\arg\min\nolimits_{y\in
S(x)}g(x,y).$ (1.1)

Here, the value of the upper-level problem at any point $x$ depends on the
solution of the lower-level problem. This framework has recently found
numerous applications in meta-learning [5, 6, 7, 8], hyperparameter
optimization [9, 10, 11], and reinforcement learning [12, 13, 14, 15]. Its
growing importance has spurred increasing efforts towards designing
computationally efficient algorithms for it.

As demonstrated by [16], a key computational step in algorithms for bilevel
optimization is estimating $\tfrac{dy^{*}(x)}{dx}$, the gradient of the lower-
level solution. This gradient estimation problem has been extensively studied
in differentiable optimization [17, 18] by applying the implicit function
theorem to the KKT system of the given problem [19, 20, 21, 22, 23, 24].
However, this technique typically entails computing (or estimating) second-
order derivatives, which can be prohibitive in high dimensions [25, 26, 27].

Recently, [28] made a big leap forward towards addressing this computational
bottleneck. Restricting themselves to the class of unconstrained bilevel
optimization, they proposed a first-order method with finite-time stationarity
guarantees. While a remarkable breakthrough, [28] does not directly extend to
the important setting of constrained bilevel optimization. This motivates the
question:

> _Can we develop a first-order algorithm for constrained bilevel
> optimization?_

Besides being natural from the viewpoint of complexity theory, this question
is well-grounded in applications such as mechanism design [29, 30], resource
allocation [31, 32, 33, 34], and decision-making under uncertainty [35, 36].
Our primary contribution is an affirmative answer to the highlighted question
for bilevel programs with linear constraints, an important problem class often
arising in adversarial training, decentralized meta learning, and sensor
networks (see [37] and the discussion therein). While there have been some
other recent works [37, 38, 39] on this problem, our work is first-order (as
opposed to [37]) and offers, in our view, a stronger guarantee on stationarity
(compared to [38, 39])— cf. Section 1.2. We hope our work paves the path to
progress for first-order methods for general constrained bilevel programs. We
now summarize our contributions.

###  1.1 Our contributions

We provide first-order algorithms (with associated finite-time convergence
guarantees) for linearly constrained bilevel programs (Eq. 1.1). By “first-
order”, we mean that we use only zeroth and first-order oracle access to $f$
and $g$. Our problem assumptions for each of the following contributions are
detailed in LABEL:{sec:differentiable-optimization}.

(1)

Linear equality constraints. As our first contribution, we design first-order
algorithms for solving Eq. 1.1 where the lower-level constraint set
$S(x):=\left\\{y:Ax-By-b=0\right\\}$ comprises linear equality constraints,
and $\mathcal{X}$ a convex compact set. With appropriate regularity
assumptions on $f$ and $g$, we show in this case smoothness in $x$ of the
hyperobjective $F$. Inspired by ideas from [40], we use implicit
differentiation of the KKT matrix of a slightly perturbed version of the
lower-level problem to design a first-order approximation to $\nabla F$.
Constructing our first-order approximation entails solving a strongly convex
optimization problem on affine constraints, which can be done efficiently.
With this inexact gradient oracle in hand, we then run projected gradient
descent, which converges in $\widetilde{O}(\epsilon^{-2})$ iterations for
smooth functions.

######  Theorem 1.1 (Informal; cf. Theorem 3.1).

Consider Eq. 1.1 with linear equality constraints $S(x)\coloneqq\left\\{y:Ax-
By-b=0\right\\}$ and $\mathcal{X}$ a convex compact set. Under mild regularity
assumptions on $f$ and $g$ (Assumptions 2.2 and 2.3), there exists an
algorithm, which in $\widetilde{O}(\epsilon^{-2})$ zeroth and first-order
oracle calls to $f$ and $g$, converges to an $\epsilon$-stationary point of
the hyperobjective $F$.

For linear equality constrained bilevel optimization, this is the first first-
order result attaining $\epsilon$-stationarity of $F$ with assumptions solely
on the constituent functions $f$ and $g$ and none on $F$ — cf. Section 1.2 for
a discussion of the results of [37] for this setting.

(2)

Linear inequality constraints. Next, we provide first-order algorithms for
solving Eq. 1.1 where the lower-level constraint set $S(x):=\left\\{y:Ax-
By-b\leq 0\right\\}$ comprises linear inequality constraints, and the upper-
level variable is unconstrained.

Our measure of convergence of algorithms in this case is that of
$(\delta,\epsilon)$-stationarity [41]: for a Lipschitz function, we say that a
point $x$ is $(\delta,\epsilon)$-stationary if within a $\delta$-ball around
$x$ there exists a convex combination of subgradients of the function with
norm at most $\epsilon$ (cf. Definition 2.1).

To motivate this notion of convergence, we note that the hyperobjective $F$
(in Eq. 1.1) as a function of $x$ could be nonsmooth and nonconvex (and is
Lipschitz, as we later prove). Minimizing such a function in general is well-
known to be intractable [42], necessitating local notions of stationarity.
Indeed, not only is it impossible to attain $\epsilon$-stationarity in finite
time [43], even getting near an approximate stationary point of an arbitrary
Lipschitz function is impossible unless the number of queries has an
exponential dependence on the dimension [44]. Consequently, for this function
class, $(\delta,\epsilon)$-stationarity has recently emerged, from the work of
[43], to be a natural and algorithmically tractable notion of stationarity. We
give the following guarantee under regularity assumptions on $f$ and $g$.

######  Theorem 1.2 (Informal; LABEL:{{thm:izo_complete_guarantees}}).

Consider Eq. 1.1 with linear inequality constraints
$S(x)\coloneqq\left\\{y:Ax-By-b\leq 0\right\\}$. Under mild regularity
assumptions on $f$ and $g$ (Assumption 2.2) and the lower-level primal
solution $y^{*}$ (Assumption 2.4), there exists an algorithm, which in
$\widetilde{O}(d\delta^{-1}\epsilon^{-3})$ oracle calls to $f$ and $g$,
converges to a $(\delta,\epsilon)$-stationary point of $F$, where $d$ is the
upper-level variable dimension.

To the best of our knowledge, this is the first result to offer a first-order
finite-time stationarity guarantee on the hyperobjective for linear inequality
constrained bilevel optimization (cf. Section 1.2 for a discussion of related
work of [37, 38, 39]). We obtain our guarantee in Theorem 1.2 by first
invoking a result by [45] to obtain inexact hyperobjective values of $F$ using
only $\widetilde{O}(1)$ oracle calls to $f$ and $g$. We also show
(LABEL:{lem:LipscConstrBilevel}) that this hyperobjective $F$ is Lipschitz. We
then employ our inexact zeroth-order oracle for $F$ in an algorithm (Algorithm
2) designed to minimize Lipschitz nonsmooth nonconvex functions (in
particular, $F$) with the following convergence guarantee.

######  Theorem 1.3 (Informal; LABEL:{thm:_Lipschitz-min-with-inexact-zero-
oracle}).

Given $L$-Lipschitz $F:\mathbb{R}^{d}\to\mathbb{R}$ and that
$|\widetilde{F}(\cdot)-F(\cdot)|\leq\epsilon$. Then there exists an algorithm,
which in $\widetilde{O}(d\delta^{-1}\epsilon^{-3})$ calls to
$\widetilde{F}(\cdot)$ outputs a point $x^{\mathrm{out}}$ such that
$\mathbb{E}[\mathrm{dist}(0,{\partial}_{\delta}F(x^{\mathrm{out}}))]\leq
2\epsilon$.

While such algorithms using _exact_ zeroth-order access already exist [46],
extending them to the inexact gradient setting is non-trivial; we leverage
recent ideas connecting online learning to nonsmooth nonconvex optimization by
[47] (cf. LABEL:{sec:nonsmooth}).

(3)

Linear inequality under assumptions on dual variable access. For the
inequality setting (i.e., Eq. 1.1 with the lower-level constraint set
$S(x):=\left\\{y:Ax-By-b\leq 0\right\\}$), we obtain dimension-free rates
under an additional assumption (Assumption 2.5) on oracle access to the
optimal dual variable $\lambda^{*}$ of the lower-level problem. We are not
aware of a method to obtain this dual variable in a first-order fashion
(though in practice, highly accurate approximations to $\lambda^{*}$ are
readily available), hence the need for imposing this assumption. We believe
that removing this assumption and obtaining dimension-free first-order rates
in this setting would be an important direction for future work. Our guarantee
for this setting is summarized below.

######  Theorem 1.4 (Informal; Theorem 4.5 combined with Theorem 5.3).

Consider Eq. 1.1 with linear inequality constraints
$S(x)\coloneqq\left\\{y:Ax-By-b\leq 0\right\\}$ and unconstrained upper-level
variable. Under mild regularity assumptions on $f$ and $g$ (Assumption 2.2),
on $y^{*}$ (Assumption 2.4), and assuming oracle access to the optimal dual
variable $\lambda^{*}$ (Assumption 2.5), there exists an algorithm, which in
$\widetilde{O}(\delta^{-1}\epsilon^{-4})$ oracle calls to $f$ and $g$
converges to a $(\delta,\epsilon)$-stationary point for $F$.

We obtain this result by first reformulating Eq. 1.1 via the penalty method
and constructing an inexact gradient oracle for the hyperobjective $F$ (cf.
Section 5). We then employ this inexact gradient oracle within an algorithm
(LABEL:{alg:_OIGRM}) designed to minimize Lipschitz nonsmooth nonconvex
functions (in particular, $F$), with the following convergence guarantee.

######  Theorem 1.5 (Informal; LABEL:{thm:Lipschitz-min-with-inexact-grad-
oracle}).

Given Lipschitz $F:\mathbb{R}^{d}\to\mathbb{R}$ and
$\|\widetilde{\nabla}F(\cdot)-\nabla F(\cdot)\|\leq\epsilon$, there exists an
algorithm that, in $T=O(\delta^{-1}\epsilon^{-3})$ calls to
$\widetilde{\nabla}F$, outputs a $(\delta,2\epsilon)$-stationary point of $F$.

Our Algorithm 3 is essentially a “first-order” version of Algorithm 2. As was
the case for Algorithm 2, despite the existence of algorithms with these
guarantees with access to exact gradients [48], their extensions to the
_inexact_ gradient setting are not trivial and also make use of the new
framework of [47]. With the ubiquity of nonsmooth nonconvex optimization
problems associated with training modern neural networks, we believe our
analysis for this general task can be of independent interest to the broader
optimization community. Lastly, we also design a more implementation-friendly
variant of Algorithm 3 (with slightly worse theoretical guarantees) that we
use in preliminary numerical experiments.

###  1.2 Related work

The vast body of work on asymptotic results for bilevel programming, starting
with classical works such as [49, 50, 51, 52, 53, 54], typically fall into two
categories: those based on approximate implicit differentiation: [55, 56, 57,
17, 58, 18, 59, 60] and those via iterative differentiation: [55, 61, 62, 9,
10, 59]. Another recent line of work in this category includes [63, 64, 37,
65], which use various smoothing techniques.

The first non-asymptotic result for bilevel programming was provided by [16],
which was followed by a flurry of work: for example, algorithms that are
single-loop stochastic: [66, 67, 68], algorithms that are projection-free:
[69, 70, 71, 72], those incorporating in them variance-reduction and momentum:
[73, 74, 75, 76], those for single-variable bilevel programs: [77, 78, 79, 70,
80], and for bilevel programs with special constraints: [81, 37, 82, 71].

The most direct predecessors of our work are those of [63, 40, 38, 39, 37]. As
alluded to earlier, [28] recently made a significant contribution by providing
for bilevel programming a first-order algorithm with finite-time stationarity
guarantees. Their work was extended to the stochastic setting by [40] (which
we build upon for our results on the linear equality setting), simplified and
improved by [83], and extended to the constrained setting by [38, 39, 37].

The works of [38, 39] study the more general problem of bilevel programming
with general convex constraints. Their notion of stationarity is that of KKT
stationarity (over both upper and lower-level variables). The guarantee we
provide in Theorem 1.4 holds for the much more restricted setting of linear
inequality constraints; however, our stationarity guarantees are expressed
directly in terms of the hyperobjective of interest, which we believe is more
useful. Extending this notion of stationarity to bilevel programs with
_general_ convex constraints (as studied in [38, 39]) would be a very
interesting direction of study. Moreover, the work of [38] assumes joint
convexity of the lower-level constraints in upper and lower variables to allow
for efficient projections, whereas our guarantee requires convexity only in
the lower-level variable.

The current best result for the linearly constrained setting is that of [37].
However, their algorithm is not first-order method work due to requiring
Hessian computations and also imposes some regularity assumptions on the
hyperobjective $F$, which are, in general, impossible to verify. In contrast,
our Theorem 1.1 imposes assumptions solely on the constituent functions $f$
and $g$ — and none directly on $F$ — thus making substantial progress on these
two fronts.

##  2 Preliminaries

We use $\langle\cdot{},{}\cdot\rangle$ to denote inner products and
$\|{}\cdot{}\|$ for the Euclidean norm. Unless transposed, all vectors are
column vectors. For $f:\mathbb{R}^{d_{2}}\to\mathbb{R}^{d_{1}}$ its Jacobian
with respect to $x\in\mathbb{R}^{d_{2}}$ is $\nabla
f\in\mathbb{R}^{d_{1}\times d_{2}}$. For $f:\mathbb{R}^{d}\to\mathbb{R}$, we
overload $\nabla f$ to refer to its gradient (the transposed Jacobian), a
column vector. We use $\nabla_{x}$ to denote partial derivatives with respect
to $x$.

A function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is $L$-Lipschitz if for any
$x,y$, we have $\|f(x)-f(y)\|\leq L\|x-y\|$. A differentiable function
$f:\mathbb{R}^{n}\to\mathbb{R}$ is convex if for any $x,y\in\mathbb{R}^{n}$ we
have $f(y)\geq f(x)+\nabla f(x)^{\top}(y-x)$; it is $\mu$-strongly convex if
$f-\tfrac{\mu}{2}\|{}\cdot{}\|^{2}$ is convex; it is $\beta$-smooth if $\nabla
f$ is $\beta$-Lipschitz.

######  Definition 2.1.

Consider a locally Lipschitz function $f:\mathbb{R}^{d}\to\mathbb{R}$, a point
$x\in\mathbb{R}^{d}$, and a parameter $\delta>0$. The Goldstein
subdifferential [41] of $f$ at $x$ is the set
$\partial_{\delta}f(x):=\mathrm{conv}(\cup_{y\in\mathbb{B}_{\delta}(x)}\partial
f(y)),$ where $\partial f(x)=\mathrm{conv}\left\\{\lim_{n\to\infty}\nabla
f(x_{n}):x_{n}\to x,~{}x_{n}\in\mathrm{dom}(\nabla f)\right\\}$ is the Clarke
subdifferential [84] of $f$ and $\mathbb{B}_{\delta}(x)$ denotes the Euclidean
ball of radius $\delta$ around $x$. A point $x$ is called
$(\delta,\epsilon)$-stationary if
$\mathrm{dist}(0,\partial_{\delta}f(x))\leq\epsilon$, where
$\mathrm{dist}(x,S):=\inf_{y\in S}\|x-y\|$.

Thus, for a Lipschitz function $f$, a point $x$ is
$(\delta,\epsilon)$-stationary if within a $\delta$-ball around $x$, there
exists a convex combination of subgradients of $f$ with norm at most
$\epsilon$. For a differentiable function $f$, we say that $x$ is
$\epsilon$-stationary if $\|\nabla f(x)\|\leq\epsilon$.

###  2.1 Assumptions

We consider LABEL:{prob:general-constraints} with linear equality constraints
(Section 3) under Assumptions 2.2 and 2.3 and linear inequality constraints
(Sections 5 and 4) under Assumptions 2.2, 2.4 and 2.5. We assume the upper-
level (UL) variable $x\in\mathbb{R}^{d_{x}}$, lower-level (LL) variable
$y\in\mathbb{R}^{d_{y}}$, and $A\in\mathbb{R}^{d_{h}\times d_{x}}$.

######  Assumption 2.2.

For Eq. 1.1, we assume the following for both settings we study:

  1. (i)

Upper-level: The objective $f$ is $C_{f}$-smooth and $L_{f}$-Lipschitz
continuous in $(x,y)$.

  2. (ii)

Lower-level: The objective $g$ is $C_{g}$-smooth. Fixing any
$x\in\mathcal{X}$, $g(x,\cdot)$ is $\mu_{g}$-strongly convex.

  3. (iii)

We assume that the linear independence constraint qualification (LICQ)
condition holds for the LL problem at every $x$ and $y$, i.e., the constraint
$h(x,y)\coloneqq Ax-By-b$ has a full row rank $B$.

######  Assumption 2.3.

For Eq. 3.1 (with linear equality constraints), we additionally assume that
the set $\mathcal{X}$ is convex and compact, and that the objective $g$ is
$S_{g}$-Hessian smooth, that is,
$\left\lVert\nabla^{2}g(x,y)-\nabla^{2}g(\bar{x},\bar{y})\right\rVert\leq
S_{g}\left\lVert(x,y)-(\bar{x},\bar{y})\right\rVert\forall
x,\bar{x}\in\mathcal{X},\text{ and }y,\bar{y}\in\mathbb{R}^{d_{y}}.$

######  Assumption 2.4.

For Eq. 4.1 (with linear inequality constraints), we additionally assume that
$y^{*}$ is $L_{y}$-Lipschitz in $x$, where $y^{*}$ is the LL primal solution
$y^{*}(x),\lambda^{*}(x)=\arg\max\nolimits_{y}\min\nolimits_{\beta\geq
0}g(x,y)+\beta^{\top}h(x,y)$, where $h(x,y)\coloneqq Ax-By-b$.

######  Assumption 2.5.

We provide additional results for Eq. 4.1 under additional stronger
assumptions stated here: Denote the LL primal and dual solution
$y^{*}(x),\lambda^{*}(x)=\arg\max\nolimits_{y}\min\nolimits_{\beta\geq
0}g(x,y)+\beta^{\top}h(x,y)$, where $h(x,y)\coloneqq Ax-By-b$; then, we assume
that we have exact access to $\lambda^{*}$ and that
$\left\lVert\lambda^{*}(x)\right\rVert\leq R$.

2.2(i) and 2.2(ii) are standard in bilevel optimization.
LABEL:{item:assumption_tangen_space} is the same as the complete recourse
assumption in stochastic programming [85], that is, the LL problem is feasible
$y$ for every $x\in\mathbb{R}^{d_{x}}$. Assumption 2.3 is used only in the
equality case and guarantees smoothness of $F$. Assumption 2.4 is used in the
inequality case and implies Lipschitzness of $F$. We need the stronger
assumption detailed in Assumption 2.5 for our dimension-free result for the
linear inequality case.

##  3 Lower-level problem with linear equality constraint

We first obtain improved results for the setting of bilevel programs with
linear equality constraints in the lower-level problem. Our formal problem
statement is:

$\displaystyle\mbox{minimize}_{x\in\mathcal{X}}~{}F(x)\coloneqq
f(x,y^{*}(x))~{}\text{ subject to
}~{}y^{*}(x)\in\arg\min\nolimits_{y:h(x,y)=0}g(x,y),$ (3.1)

where $f$, $g$, $h(x,y)\coloneqq Ax-By-b$, and $\mathcal{X}$ satisfy
Assumptions 2.2 and 2.3. The previous best result on LABEL:{prob:lin-eq}
providing finite-time $\epsilon$-stationarity guarantees, by [37], required
certain regularity assumptions on $F$ as well as Hessian computations. In
contrast, our finite-time guarantees _require assumptions only on $f$ and $g$,
not on $F$_; indeed, in our work, these desirable properties of $F$ are
naturally implied by our analysis. Specifically, our key insight is that the
hypergradient $\nabla
F(x):=\nabla_{x}f(x,y^{*})+\left(\tfrac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*})$
for LABEL:{prob:lin-eq} is Lipschitz-continuous and admits an easily
computable — yet highly accurate — finite-difference approximation. Therefore,
$O(\epsilon^{-2})$ iterations of gradient descent on $F$ with this finite-
difference gradient proxy yield an $\epsilon$-stationary point.

Specifically, for any fixed $x\in\mathcal{X}$, our proposed finite-difference
gradient proxy approximating the non-trivial-to-compute component
$\left(\tfrac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*})$ of the
hypergradient is given by

$v_{x}:={\frac{\nabla_{x}[g(x,y_{\delta}^{*})+\langle\lambda_{\delta}^{*},h(x,y^{*})\rangle]-\nabla_{x}[g(x,y^{*})+\langle\lambda^{*},h(x,y^{*})\rangle]}{\delta},}$

where $(y_{\delta}^{*},\lambda_{\delta}^{*})$ are the primal and dual
solutions to the perturbed lower-level problem:

$\displaystyle
y_{\delta}^{*}\coloneqq\arg\min\nolimits_{y:h(x,y)=0}~{}g(x,y)+\delta f(x,y).$
(3.3)

We show in LABEL:{lem:lineq-finitediff-equals-gradf} that $v$ in Section 3
approximates $\left(\tfrac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*})$
up to an $O(\delta)$-additive error, implying the gradient oracle construction
outlined in the pseudocode presented in Algorithm 1. Our full implementable
algorithm for solving Eq. 3.1 is displayed in LABEL:{alg:LE-full-alg}.

Algorithm 1 Inexact Gradient Oracle for Bilevel Program with Linear Equality
Constraint

1:Input: Current $x$, accuracy $\epsilon$, perturbation $\delta=\epsilon^{2}$.

2:Compute $y^{*}$ (as in Eq. 3.1) and corresponding optimal dual $\lambda^{*}$
(as in Eq. A.1)

3:Compute $y_{\delta}^{*}$ (as in Eq. 3.3) and and corresponding optimal dual
$\lambda_{\delta}^{*}$ (as in Eq. A.6)

4:Compute $v_{x}$ as in Section 3 $\triangleright$ Approximates
$\left({dy^{*}(x)}/{dx}\right)^{\top}\nabla_{y}f(x,y^{*})$

5:Output: $\widetilde{\nabla}F=v_{x}+\nabla_{x}f(x,y^{*})$

Notice that the finite-difference term in LABEL:{eq:part-hypergrad-approx-lin-
eq} avoids differentiating through the implicit function $y^{*}$. Instead, all
we need to evaluate it are the values of
$(y^{*},\lambda^{*},y_{\delta}^{*},\lambda_{\delta}^{*})$ (and gradients of
$g$ and $h$). Since $(y^{*},\lambda^{*})$ are solutions to a smooth strongly
convex linearly constrained problem, they can be approximated at a linear
rate. Similarly, since the approximation error in Section 3 is proportional to
$\delta$ (cf. LABEL:{lem:lineq-finitediff-equals-gradf}), a small enough
$\delta$ in the perturbed objective $g+\delta f$ in Eq. 3.3 ensures that it is
dominated by the strongly convex and smooth $g$, whereby accurate approximates
to $(y_{\delta}^{*},\lambda_{\delta}^{*})$ can also be readily obtained.
Putting it all together, the proposed finite-difference hypergradient proxy in
Section 3 is efficiently computable, yielding the following guarantee.

######  Theorem 3.1.

Consider Eq. 3.1 under Assumption 2.2, and let $\kappa=C_{g}/\mu_{g}$ be the
condition number of $g$. Then Algorithm 6 finds an $\epsilon$-stationary point
(in terms of gradient mapping, see Eq. A.13) after
$T=\widetilde{O}(C_{F}(F(x_{0})-\inf F)\sqrt{\kappa}\epsilon^{-2})$ oracle
calls to $f$ and $g$, where
$C_{F}:=2(L_{f}+C_{f}+C_{g})C_{H}^{3}S_{g}(L_{g}+\|A\|)^{2}$ is the smoothness
constant of the hyperobjective $F$.

We now sketch the proofs of the key components that together imply Theorem
3.1. The complete proofs may be found in LABEL:{sec:appendix_linear_equality}.

###  3.1 Main technical ideas

We briefly outline the two key technical building blocks alluded to above,
that together give us LABEL:{thm:lineq-cost}: the approximation guarantee of
our finite-difference gradient proxy (LABEL:{eq:part-hypergrad-approx-lin-eq})
and the smoothness of hyperobjective $F$ (for LABEL:{prob:lin-eq}). The
starting point for both these results is the following simple observation
obtained by implicitly differentiating, with respect to $x$, the KKT system
associated with $y^{*}=\arg\min_{h(x,y)=0}g(x,y)$ and optimal dual variable
$\lambda^{*}$:

$\begin{bmatrix}\frac{dy^{*}(x)}{dx}\\\
\frac{d\lambda^{*}(x)}{dx}\end{bmatrix}={\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}}\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})\\\
-\nabla_{x}h(x,y^{*})\end{bmatrix}$

The invertibility of the matrix in the preceding equation is proved in
LABEL:{cor:nonsingularH}: essentially, this invertibility is implied by strong
convexity of $g$ and $\nabla_{y}h(x,y^{*})=B$ having full row rank. This in
conjunction with the compactness of $\mathcal{X}$ implies that the inverse of
the matrix is bounded by some constant $C_{H}$ (cf. Corollary A.3 for
details). Our hypergradient approximation guarantee follows:

######  Lemma 3.2.

For Eq. 3.1 under Assumption 2.2, with $v_{x}$ as in Section 3, the following
holds:

$\textstyle{\left\|v_{x}-\left(\frac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*})\right\|\leq
O(C_{F}\delta).}$

###### Proof sketch; see LABEL:{sec:appendix_linear_equality} for the complete
proof.

The main idea for this proof is that the two terms being compared are
essentially the same by an application of the implicit function theorem.
First, we can use the expression for $\tfrac{dy^{*}(x)}{dx}$ from Section 3.1
as follows:

$\displaystyle\left(\frac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*})$
$\displaystyle=\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})\\\
-\nabla_{x}h(x,y^{*})\end{bmatrix}^{\top}\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}\begin{bmatrix}\nabla_{y}f(x,y^{*})\\\
0\end{bmatrix}.$

We now examine our finite-difference gradient proxy. For simplicity of
exposition, we instead consider $\lim_{\delta\rightarrow
0}\frac{\nabla_{x}[g(x,y_{\delta}^{*})+\langle\lambda_{\delta}^{*},h(x,y^{*})\rangle]-\nabla_{x}[g(x,y^{*})+\langle\lambda^{*},h(x,y^{*})\rangle]}{\delta},$
which, by an application of the fundamental theorem of calculus and Assumption
2.3 , equals our finite-difference gradient proxy up to an
$O(\delta)$-additive error. Note that this expression is exactly

$\nabla^{2}_{xy}g(x,y^{*})\frac{dy_{\delta}^{*}}{d\delta}+\nabla_{x}h(x,y^{*})^{\top}\frac{d\lambda_{\delta}^{*}}{d\delta}.$

Since $y^{*}_{\delta}$ is also a minimizer of a strongly convex function over
a linear equality constraint LABEL:{eq:lower_perturb}, the same reasoning that
yields Section 3.1 can be also used to obtain

$\begin{bmatrix}\frac{dy_{\delta}^{*}(x)}{d\delta}\\\
\frac{d\lambda_{\delta}^{*}(x)}{d\delta}\end{bmatrix}\Bigg{|}_{\delta=0}=\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}\begin{bmatrix}-\nabla_{y}f(x,y^{*})\\\
0\end{bmatrix}.$

Combining Section 3.1 and Section 3.1 immediately yields the claim. ∎

Having shown our hypergradient approximation, we now turn to showing
smoothness of the hyperobjective $F$, crucial to getting our claimed rate.

######  Lemma 3.3.

The solution $y^{*}$ (as defined in LABEL:{prob:lin-eq}) is
$O(C_{H}\cdot(C_{g}+\|A\|))$-Lipschitz continuous and $O(C_{H}^{3}\cdot
S_{g}\cdot(C_{g}+\|A\|)^{2})$-smooth as a function of $x$. Thus the hyper-
objective $F$ is gradient-Lipschitz with a smoothness constant of
$C_{F}:=O\\{(L_{f}+C_{f}+C_{g})C_{H}^{3}S_{g}(L_{g}+\|A\|)^{2}\\}.$

###### Proof sketch; see LABEL:{sec:appendix_linear_equality}.

The Lipschitz bound follows by combining LABEL:{eq:def_matrix_H} with the
smoothness of $g$ and bound $C_{H}$ on the matrix inverse therein.
Differentiating LABEL:{eq:def_matrix_H} by $x$ yields a linear system with the
same matrix; we repeat this approach to get the smoothness bound. ∎

##  4 Lower-level problem with linear inequality constraint: nonsmooth
nonconvex optimization with inexact oracles

We now shift gears from the case of linear _equality_ constraints to that of
linear _inequality_ constraints. Specifically, defining $h(x,y)=Ax-By-b$, the
problem we now consider is

$\displaystyle\mathop{\text{minimize}}\nolimits_{x}~{}F(x)\coloneqq
f(x,y^{*}(x))\quad\text{ subject
to}~{}y^{*}(x)\in\operatorname*{arg\,min}\nolimits_{y:h(x,y)\leq 0}g(x,y).$
(4.1)

As noted earlier, for this larger class of problems, the hyperobjective $F$
can be nonsmooth nonconvex, thereby necessitating our measure of convergence
to be the now popular notion of Goldstein stationarity [43].

Our first algorithm for solving Eq. 4.1 is presented in Algorithm 2, with its
convergence guarantee in Theorem 4.1. At a high level, this algorithm first
obtains access to an inexact zeroth-order oracle to $F$ (we shortly explain
how this is done) and uses this oracle to construct a (biased) gradient
estimate of $F$. It then uses this gradient estimate to update the iterates
with an update rule motivated by recent works reducing nonconvex optimization
to online optimization [47]. We explain these components in Section 4.1.

1:Input: Initialization $x_{0}\in\mathbb{R}^{d}$, clipping parameter $D>0$,
step size $\eta>0$, smoothing parameter $\rho>0$, accuracy parameter $\nu>0$,
iteration budget $T\in\mathbb{N}$, inexact zero-order oracle
$\widetilde{F}:\mathbb{R}^{d}\to\mathbb{R}$.

2:Initialize: ${\Delta}_{1}=\mathbf{0}$

3:for $t=1,\dots,T$ do

4: Sample $s_{t}\sim\mathrm{Unif}[0,1]$,
$w_{t}\sim\mathrm{Unif}(\mathbb{S}^{d-1})$

5: $x_{t}=x_{t-1}+{\Delta}_{t}$, $z_{t}=x_{t-1}+s_{t}{\Delta}_{t}$

6: $\widetilde{{g}}_{t}=\tfrac{d}{2\rho}(\widetilde{F}(z_{t}+\rho
w_{t})-\widetilde{F}(z_{t}-\rho w_{t}))w_{t}$

7:
${\Delta}_{t+1}=\mathrm{clip}_{D}\left({\Delta}_{t}-\eta\widetilde{{g}}_{t}\right)$
$\triangleright$ $\mathrm{clip}_{D}(z):=\min\\{1,\tfrac{D}{\left\lVert
z\right\rVert}\\}\cdot z$

8:$M=\lfloor\frac{\nu}{D}\rfloor,~{}K=\lfloor\frac{T}{M}\rfloor$

9:for $k=1,\dots,K$ do

10: $\overline{x}_{k}=\frac{1}{M}\sum_{m=1}^{M}z_{(k-1)M+m}$

11:Sample
$x^{\mathrm{out}}\sim\mathrm{Unif}\\{\overline{x}_{1},\dots,\overline{x}_{K}\\}$

12:Output: $x^{\mathrm{out}}$.

Algorithm 2 Nonsmooth Nonconvex Algorithm with Inexact Zero-Order oracle

######  Theorem 4.1.

Consider LABEL:{prob:ineq} under Assumption 2.2 and Assumption 2.4. Let
$\kappa=C_{g}/\mu_{g}$ be the condition number of $g$. Then combining the
procedure for LABEL:{lem:ZeroOrderApprox} with Algorithm 2 run with parameters
$\rho=\min\left\\{\tfrac{\delta}{2},\tfrac{F(x_{0})-\inf
F}{L_{f}L_{y}}\right\\},\nu=\delta-\rho,~{}D=\Theta\left(\frac{\nu\epsilon^{2}\rho^{2}}{d_{x}\rho^{2}L_{f}^{2}L_{y}^{2}+\alpha^{2}d_{x}^{2}}\right)$,
and
$\eta=\Theta\left(\frac{\nu\epsilon^{3}\rho^{4}}{(d_{x}\rho^{2}L_{f}^{2}L_{y}^{2}+\alpha^{2}d_{x}^{2})^{2}}\right)$
outputs a point $x^{\mathrm{out}}$ such that
$\mathbb{E}[\mathrm{dist}(0,{\partial}_{\delta}F(x^{\mathrm{out}}))]\leq\epsilon+\alpha$
with

$T=O\left(\frac{\sqrt{\kappa}d_{x}(F(x_{0})-\inf
F)}{\delta\epsilon^{3}}\cdot\left(L_{f}^{2}L_{y}^{2}+\alpha^{2}\left(\frac{d_{x}}{\delta^{2}}+\frac{d_{x}L_{f}^{2}L_{y}^{2}}{(F(x_{0})-\inf
F)^{2}}\right)\right)\cdot\log(L_{f}/\alpha)\right)\text{ oracle calls to
}f\text{ and }g.$

###  4.1 Nonsmooth nonconvex optimization with inexact zeroth-order oracle

We can obtain inexact zeroth-order oracle access to $F$ because (as formalized
in Lemma 4.2) despite potential nonsmoothness and nonconvexity of the
hyperobjective $F$ in Eq. 4.1, estimating its _value_ $F(x)$ at any given
point $x$ amounts to solving a single smooth and strongly convex optimization
problem, which can be done can be done using $\widetilde{O}(1)$ oracle calls
to $f$ and $g$ by appealing to a result by [45].

######  Lemma 4.2.

Given any $x$, we can return $\widetilde{F}(x)$ such that
$|F(x)-\widetilde{F}(x)|\leq\alpha$ using
$O(\sqrt{C_{g}/\mu_{g}}\log(L_{f}/\alpha))$ first-order oracle calls to $f$
and $g$.

###### Proof of Lemma 4.2.

We note that it suffices to find $\tilde{y}^{*}$ such that
$\left\lVert\tilde{y}^{*}-y^{*}(x)\right\rVert\leq\alpha/L_{f}$, since setting
$\widetilde{F}(x):=f(x,\tilde{y}^{*})$ will then satisfy
$|\widetilde{F}(x)-F(x)|=|f(x,\tilde{y}^{*})-f(x,y^{*}(x))|\leq
L_{f}\cdot\tfrac{\alpha}{L_{f}}=\alpha$ by Lispchitzness of $f$, as required.
Noting that $y^{*}(x)=\arg\min_{h(x,y)\leq 0}g(x,y)$ is the solution to a
constrained smooth, strongly-convex problem with condition number
$C_{g}/\mu_{g}$, it is possible to approximate it up to $\alpha/L_{f}$ with
$O(\sqrt{C_{g}/\mu_{g}}\log(L_{f}/\alpha))$ first-order oracle calls using the
result of [45]. ∎

Having computed an inexact value of the hyperobjective $F$ (in Theorem 5.3),
we now show how to use it to develop an algorithm for solving Eq. 1.1. To this
end, we first note that $F$, despite being possibly nonsmooth and nonconvex,
is Lipschitz and hence amenable to the use of recent algorithmic developments
in nonsmooth nonconvex optimization pertaining to Goldstein stationarity.

######  Lemma 4.3.

Under Assumption 2.2 and 2.5, $F$ in Eq. 1.1 is $O(L_{f}L_{y})$-Lipschitz in
$x$.

With this guarantee on the Lipschitzness of $F$, we prove
LABEL:{thm:_Lipschitz-min-with-inexact-zero-oracle} for attaining Goldstein
stationarity using the inexact zeroth-order oracle of a Lipschitz function.
Our proof of LABEL:{thm:_Lipschitz-min-with-inexact-zero-oracle} crucially
uses the recent online-to-nonconvex framework of [47]. Combining Lemma 4.2 and
Theorem 4.4 then immediately implies Theorem 4.1.

######  Theorem 4.4.

Suppose $F:\mathbb{R}^{d}\to\mathbb{R}$ is $L$-Lipschitz, and that
$|\widetilde{F}(\cdot)-F(\cdot)|\leq\alpha$. Then running Algorithm 2 with
$\rho=\min\left\\{\tfrac{\delta}{2},\tfrac{F(x_{0})-\inf
F}{L}\right\\},\nu=\delta-\rho,~{}D=\Theta\left(\frac{\nu\epsilon^{2}\rho^{2}}{d\rho^{2}L^{2}+\alpha^{2}d^{2}}\right),\eta=\Theta\left(\frac{\nu\epsilon^{3}\rho^{4}}{(d\rho^{2}L^{2}+\alpha^{2}d^{2})^{2}}\right)$,
outputs a point $x^{\mathrm{out}}$ such that
$\mathbb{E}[\mathrm{dist}(0,{\partial}_{\delta}F(x^{\mathrm{out}}))]\leq\epsilon+\alpha$
with

$T=O\left(\frac{d(F(x_{0})-\inf
F)}{\delta\epsilon^{3}}\cdot\left(L^{2}+\alpha^{2}(\frac{d}{\delta^{2}}+\frac{dL^{2}}{(F(x_{0})-\inf
F)^{2}})\right)\right)\text{ calls to }\widetilde{F}(\cdot).$

###  4.2 Nonsmooth nonconvex optimization with inexact _gradient_ oracle

In the next section (LABEL:{sec:inequality-bilevel}), we provide a way to
generate approximate gradients of $F$. Here, we present an algorithm that
attains Goldstein stationarity of Eq. 4.1 using this inexact gradient oracle.
While there has been a long line of recent work on algorithms for nonsmooth
nonconvex optimization with convergence to Goldstein stationarity [43, 48, 86,
87, 88], these results necessarily require exact gradients. This brittleness
to any error in gradients renders them ineffective in our setting, where our
computed (hyper)gradient necessarily suffers from an additive error. While
inexact oracles are known to be effective for smooth or convex objectives
[89], utilizing inexact gradients in the nonsmooth nonconvex regime presents a
nontrivial challenge. Indeed, without any local bound on gradient variation
due to smoothness, or convexity that ensures that gradients are everywhere
correlated with the direction pointing at the minimum, it is not clear a
priori how to control the accumulating price of inexactness throughout the run
of an algorithm. To derive such results, we again use the recently proposed
connection between online learning and nonsmooth nonconvex optimization by
[47]. In particular, by controlling the accumulated error suffered by online
gradient descent for _linear_ losses (cf. Lemma B.2), we derive guarantees for
our setting of interest, providing Lipschitz optimization algorithms that
converge to Goldstein stationary points even with inexact gradients.

This algorithm matches the best known complexity in first-order nonsmooth
nonconvex optimization [43, 48, 47], merely replacing the convergence to a
$(\delta,\epsilon)$-stationary point by
$(\delta,\epsilon+\alpha)$-stationarity, where $\alpha$ is the inexactness of
the gradient oracle.

1:Input: Initialization $x_{0}\in\mathbb{R}^{d}$, clipping parameter $D>0$,
step size $\eta>0$, accuracy parameter $\delta>0$, iteration budget
$T\in\mathbb{N}$, inexact gradient oracle
$\widetilde{\nabla}F:\mathbb{R}^{d}\to\mathbb{R}^{d}$.

2:Initialize: ${\Delta}_{1}=\mathbf{0}$

3:for $t=1,\dots,T$ do

4: Sample $s_{t}\sim\mathrm{Unif}[0,1]$

5: $x_{t}=x_{t-1}+{\Delta}_{t}$, $z_{t}=x_{t-1}+s_{t}{\Delta}_{t}$

6: $\widetilde{{g}}_{t}=\widetilde{\nabla}F(z_{t})$

7:
${\Delta}_{t+1}=\mathrm{clip}_{D}\left({\Delta}_{t}-\eta\widetilde{{g}}_{t}\right)$
$\triangleright$ $\mathrm{clip}_{D}(z):=\min\\{1,\tfrac{D}{\left\lVert
z\right\rVert}\\}\cdot z$

8:$M=\lfloor\frac{\delta}{D}\rfloor,~{}K=\lfloor\frac{T}{M}\rfloor$

9:for $k=1,\dots,K$ do

10: $\overline{x}_{k}=\frac{1}{M}\sum_{m=1}^{M}z_{(k-1)M+m}$

11:Sample
$x^{\mathrm{out}}\sim\mathrm{Unif}\\{\overline{x}_{1},\dots,\overline{x}_{K}\\}$

12:Output: $x^{\mathrm{out}}$.

Algorithm 3 Nonsmooth Nonconvex Algorithm with Inexact Gradient Oracle

######  Theorem 4.5.

Suppose $F:\mathbb{R}^{d}\to\mathbb{R}$ is $L$-Lipschitz and that
$\|\widetilde{\nabla}F(\cdot)-\nabla F(\cdot)\|\leq\alpha$. Then running
Algorithm 3 with
$D=\Theta(\frac{\delta\epsilon^{2}}{L^{2}}),\eta=\Theta(\frac{\delta\epsilon^{3}}{L^{4}})$,
outputs a point $x^{\mathrm{out}}$ such that
$\mathbb{E}[\mathrm{dist}(0,{\partial}_{\delta}F(x^{\mathrm{out}}))]\leq\epsilon+\alpha$,
with $T=O\left(\frac{(F(x_{0})-\inf F)L^{2}}{\delta\epsilon^{3}}\right)$ calls
to $\widetilde{\nabla}F(\cdot)$.

We defer the proof of Theorem 4.5 to LABEL:{sec:zeroth-order-algs}. Plugging
the complexity of computing inexact gradients, as given by Theorem 5.3, into
the result above, we immediately obtain convergence to a
$(\delta,\epsilon)$-stationary point of Eq. 1.1 with
$\widetilde{O}(\delta^{-1}\epsilon^{-4})$ gradient calls overall.

#### Implementation-friendly algorithm.

While LABEL:{alg:_OIGRM} matches the best-known results in nonsmooth nonconvex
optimization, it could be impractical due to several hyperparameters which
need tuning to guarantee optimal results. Arguably, a more natural application
of the hypergradient estimates would be simply plugging them into gradient
descent, which requires tuning only the stepsize. Since the objective $F$ is
neither smooth nor convex, perturbations are required to guarantee
differentiability along the trajectory. We therefore complement Theorem 4.5 by
analyzing perturbed (inexact) gradient descent in the nonsmooth nonconvex
setting (LABEL:{alg:_PIGD}). Despite its suboptimal worst-case theoretical
guarantees, we find this algorithm easier to implement in practice. We defer
the proof of LABEL:{thm:practical_Lipschitz-min-with-inexact-grad-oracle} to
LABEL:{sec:zeroth-order-algs}.

1:Input: Inexact gradient oracle
$\widetilde{\nabla}F:\mathbb{R}^{d}\to\mathbb{R}^{d}$, initialization
$x_{0}\in\mathbb{R}^{d}$, spatial parameter $\delta>0$, step size $\eta>0$,
iteration budget $T\in\mathbb{N}$.

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

3: Sample $w_{t}\sim\mathrm{Unif}(\mathbb{S}^{d-1})$

4: $\widetilde{{g}}_{t}=\widetilde{\nabla}F(x_{t}+\delta\cdot w_{t})$

5: $x_{t+1}=x_{t}-\eta\widetilde{{g}}_{t}$

6:Output: $x^{\mathrm{out}}\sim\mathrm{Unif}\\{x_{0},\dots,x_{T-1}\\}$.

Algorithm 4 Perturbed Inexact GD

######  Theorem 4.6.

Suppose $F:\mathbb{R}^{d}\to\mathbb{R}$ is $L$-Lipschitz, and that
$\|\widetilde{\nabla}F(\cdot)-\nabla F(\cdot)\|\leq\alpha$. Then running
Algorithm 4 with $\eta=\Theta\left(\frac{{((F(x_{0})-\inf F)+\delta
L)^{1/2}\delta^{1/2}}}{{T^{1/2}L^{1/2}{d}^{1/4}(\alpha+L)}}\right)$ outputs a
point $x^{\mathrm{out}}$ such that
$\mathbb{E}[\mathrm{dist}(0,{\partial}_{\delta}F(x^{\mathrm{out}}))]\leq\epsilon+\sqrt{\alpha
L}$, with

$T=O\left(\frac{(F(x_{0})-\inf F+\delta
L)L^{3}\sqrt{d}}{\delta\epsilon^{4}}\right)\text{ calls to
}\widetilde{\nabla}F(\cdot).$

##  5 Lower-level problem with linear inequality constraint: constructing the
inexact gradient oracle

In this section, we present a method to approximate the hypergradient $\nabla
F(x)$ at any given point $x$, where:

$\displaystyle\nabla
F(x)=\nabla_{x}f(x,y^{*})+\left(\tfrac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*}).$
(5.1)

The key computational challenge here lies in computing $\tfrac{dy^{*}(x)}{dx}$
since it requires differentiating through an argmin operator, which typically
requires second-order derivatives. Instead, here we differentiate (using the
implicit function theorem) through the KKT conditions describing $y^{*}(x)$
and get:

$\displaystyle\begin{bmatrix}\nabla^{2}_{yy}g+(\lambda^{*})^{\top}\nabla_{yy}^{2}h&\nabla_{y}h_{\mathcal{I}}^{\top}\\\
\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{y}h_{\mathcal{I}}&0\end{bmatrix}\begin{bmatrix}\frac{dy^{*}(x)}{dx}\\\
\frac{d\lambda^{*}_{\mathcal{I}}(x)}{dx}\end{bmatrix}=-\begin{bmatrix}\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h\\\
\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{x}h_{\mathcal{I}}\end{bmatrix}$
(5.2)

where given $x$, we assume efficient access to the optimal dual solution
$\lambda^{*}(x)$ of the LL problem in Eq. 4.1. In Eq. 5.2, we use
$\mathcal{I}\coloneqq\\{i\in[d_{h}]:h_{i}(x,y)=0,\lambda^{*}_{i}>0\\}$ to
denote the set of active constraints with non-zero dual solution, with
$h_{\mathcal{I}}\coloneqq[h_{i}]_{i\in\mathcal{I}}$ and
$\lambda^{*}_{\mathcal{I}}\coloneqq[\lambda^{*}_{i}]_{i\in\mathcal{I}}$ being
the constraints and dual variables, respectively, corresponding to
$\mathcal{I}$.

Observe that as currently stated, Eq. 5.2 leads to a second-order computation
of ${dy^{*}(x)}/{dx}$. In the rest of the section, we provide a fully first-
order approximate hypergradient oracle by constructing an equivalent
reformulation of Eq. 4.1 using a penalty function.

###  5.1 Reformulation via the penalty method

We begin by reformulating Eq. 4.1 into a single level constrained optimization
problem:

$\displaystyle\text{minimize}_{x,y}~{}f(x,y)\text{ subject
to}~{}\begin{cases}g(x,y)+(\lambda^{*}(x))^{\top}h(x,y)\leq g^{*}(x)\\\
h(x,y)\leq 0\end{cases},$ (5.3)

where $g^{*}(x)\coloneqq\min_{y:h(x,y)\leq 0}g(x,y)=g(x,y^{*}(x))$ and
$\lambda^{*}(x)$ is the optimal dual solution. The equivalence of this
reformulation to Eq. 4.1 is spelled out in Section C.1. From Eq. 5.3, we
define the following penalty function, crucial to our analysis:

$\displaystyle\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y)=f(x,y)+\alpha_{1}\left(g(x,y)+(\lambda^{*})^{\top}h(x,y)-g^{*}(x)\right)+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y)\right\rVert^{2},$ (5.4)

where $\bm{\alpha}=[\alpha_{1},\alpha_{2}]\geq 0$ are the penalty parameters.
Notably, we can compute its derivative with respect to $x$ of Eq. 5.4 in a
fully first-order fashion by the following expression:

$\displaystyle\nabla_{x}\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y)=\nabla_{x}f(x,y)\\!+\\!\alpha_{1}(\nabla_{x}g(x,y)\\!+\\!\nabla_{x}h(x,y)^{\top}\lambda^{*}\\!-\\!\nabla_{x}g^{*}(x))\\!+\\!\alpha_{2}\nabla_{x}h_{\mathcal{I}}(x,y)^{\top}h_{\mathcal{I}}(x,y).$

To give some intuition for our choice of penalties in LABEL:{eqn:penalty-
lagrangian}, we note that the two constraints in Eq. 5.3 behave quite
differently. The first constraint $g(x,y)+\lambda^{*}(x)^{\top}h(x,y)\leq
g^{*}(x)$ is one-sided, i.e., can only be violated or met, and hence just
needs a penalty parameter $\alpha_{1}$ to weight the “violation”. As to the
second constraint $h(x,y)\leq 0$, it can be arbitrary. To allow for such a
“two-sided” constraint, we penalize only the active constraints $\mathcal{I}$,
i.e., we use $\left\lVert h_{\mathcal{I}}(x,y)\right\rVert^{2}$ to penalize
deviation.

Next, we define the optimal solutions to the penalty function optimization by:

$\displaystyle
y_{\lambda^{*},\bm{\alpha}}^{*}(x):=\arg\min\nolimits_{y}\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y).$
(5.5)

We now show that this minimizer is close to the optimal solution of the LL
problem, while suffering only a small constraint violation.

######  Lemma 5.1.

Given any $x$, the corresponding dual solution $\lambda^{*}(x)$, primal
solution $y^{*}(x)$ of the lower optimization problem in Eq. 4.1, and
$y_{\lambda^{*},\bm{\alpha}}^{*}(x)$ as in LABEL:{eq:def_y_lambda_star},
satisfy:

$\displaystyle\left\lVert
y_{\lambda^{*},\bm{\alpha}}^{*}(x)-y^{*}(x)\right\rVert\leq
O(\alpha_{1}^{-1})\text{~{}~{}and~{}~{}}\left\lVert
h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*}(x))\right\rVert\leq
O(\alpha_{1}^{-1/2}\alpha_{2}^{-1/2}).$ (5.6)

The proof of Lemma 5.1 is based on the Lipschitzness of $f$ and strong
convexity of $g$ for sufficiently large $\alpha_{1}$. The aforementioned
constraint violation bound on $h_{\mathcal{I}}(x,y)$ is later used in Lemma
5.2 to bound the inexactness of our proposed gradient oracle.

###  5.2 Main result: approximating the hypergradient

The main export of this section is the following bound on the approximation of
the hypergradient. This, together with the bounds in Lemma 5.1, validate our
use of the penalty function in Eq. 5.4.

######  Lemma 5.2.

Consider $F$ as in LABEL:{prob:ineq}, $\mathcal{L}$ as in LABEL:{eqn:penalty-
lagrangian}, a fixed $x$, and $y_{\lambda^{*},\bm{\alpha}}^{*}$ as in
LABEL:{eq:def_y_lambda_star}. Then under Assumptions 2.2 and 2.5, we have:

$\displaystyle\left\lVert\nabla
F(x)-\nabla_{x}\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert\leq
O({\alpha_{1}^{-1}})+O({\alpha_{1}^{-1/2}\alpha_{2}^{-1/2}})+O({\alpha_{1}^{1/2}}{\alpha_{2}^{-1/2}})+O({\alpha_{1}^{-3/2}}{\alpha_{2}^{1/2}}).$

The proof can be found in Section C.4. With this hypergradient approximation
guarantee, we design Algorithm 5 to compute an inexact gradient oracle for the
hyperobjective $F$.

Algorithm 5 Inexact Gradient Oracle for General Inequality Constraints

1:Input: Upper level variable $x$, accuracy $\alpha$, penalty parameters
$\alpha_{1}={\alpha^{-2}},\alpha_{2}={\alpha^{-4}}$.

2:Compute $y^{*}$, $\lambda^{*}$, and active constraints $\mathcal{I}$ of the
constrained LL problem $\min_{y:h(x,y)\leq 0}g(x,y)$.

3:Define penalty function $\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y)$ by Eq.
5.4

4:Compute the minimizer
$y^{*}_{\lambda^{*},\bm{\alpha}}=\arg\min\nolimits_{y}\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y)$
(as in Eq. 5.5).

5:Output:
$\widetilde{\nabla}F\coloneqq\nabla_{x}\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y^{*}_{\lambda^{*},\bm{\alpha}})$.

######  Theorem 5.3.

Given any accuracy parameter $\alpha>0$, Algorithm 5 outputs
$\widetilde{\nabla}_{x}F(x)$ such that $\|\widetilde{\nabla}F(x)-\nabla
F(x)\|\leq\alpha$ within $\widetilde{O}({\alpha^{-1}})$ gradient oracle
evaluations.

###### Proof sketch (full proof in Section C.5).

By choosing the penalty parameters $\alpha_{1}=\alpha^{-2}$ and
$\alpha_{2}=\alpha^{-4}$, Lemma 5.2 guarantees the inexactness of the gradient
oracle is upper bounded by $O(\alpha)$.

For the penalty minimization, the objective $\mathcal{L}$ is strongly convex
but with a condition number $\kappa_{\mathcal{L}}=O(\alpha^{-2})$ due to
strong convexity $O(\alpha_{1}\mu_{g})=O({\alpha^{-2}\mu_{g}})$ and smoothness
$O(\alpha_{2}C_{h})=O(\alpha^{-4})$. This leads to
$O(\sqrt{\kappa_{\mathcal{L}}}\log\|\alpha^{-1}\|)=\widetilde{O}({\alpha^{-1}})$
iterations to converge to $\alpha$ accuracy using standard accelerated
methods.

Using these error bounds, we can bound the error propagation of inexact
solution $y^{*}_{\lambda,\alpha}$ and inexact gradient oracle using Lemma 5.2
by total error $O(\alpha)$. Therefore, putting all together, we need only
$O(\alpha^{-1})$ oracle calls and can maintain $O(\alpha)$ inexactness. ∎

Theorem 5.3 shows that we can compute an inexact gradient from Algorithm 5
with $\tilde{O}(\alpha^{-1})$ gradient oracles. This concludes the proof of
finding $(\delta,\epsilon)$-stationary point of Eq. 1.1 with linear inequality
constraints using $\tilde{O}(\delta^{-1}\epsilon^{-4})$ gradient oracle calls.

##  6 Experiments

(a) Convergence comparison and gradient error plot.

(b) Convergence analysis with varying gradient inexactness $\alpha$.

(c) Computation cost per gradient step of varying problem size $d_{y}$.

Figure 1: We run Algorithm 4 using Algorithm 5 on the bilevel optimization in
the toy example in Eq. 6.1 with $d_{x}=100$, $d_{y}=200$,
$n_{\text{const}}=d_{y}/5$, and accuracy $\alpha=0.1$. Figure 1(a), Figure
1(b), Figure 1(c) vary # of iterations, gradient exactness $\alpha$, and
$d_{y}$, respectively, to compare the performance under different settings.

We generate instances of the following constrained bilevel optimization
problem with linear inequality constraints:

$\displaystyle\mbox{minimize}_{x}~{}c^{\top}y^{*}+0.01\left\lVert
x\right\rVert^{2}+0.01\left\lVert y^{*}\right\rVert^{2}~{}\text{ subject to
}~{}y^{*}=\operatorname*{arg\,min}_{y:h(x,y)\leq
0}\frac{1}{2}y^{\top}Qy+x^{\top}Py,$ (6.1)

where $h(x,y)=Ay-b$ is a $d_{h}$-dim linear constraint. The PSD matrix
$Q\in\mathbb{R}^{d_{y}\times d_{y}}$, $c\in\mathbb{R}^{d_{y}}$,
$P\in\mathbb{R}^{d_{x}\times d_{y}}$, and constraints
$A\in\mathbb{R}^{d_{h}\times d_{y}}$, $b\in\mathbb{R}^{d_{h}}$ are randomly
generated from normal distributions (cf. Appendix D). We compare our Algorithm
4 with a non-fully first-order method implemented using cvxpyLayer [18]. Both
algorithms use Adam [90] to control the learning rate in gradient descent. All
the experiments are averaged over ten random seeds.

Figure 1(a) shows that both the algorithms converge to the same optimal
solution at the same rate. Simultaneously, the colorful bars represent the
gradient differences between two methods, showing the inexactness of our
gradients. Figure 1(b) additionally varies this inexactness to demonstrate its
impact on convergence with standard deviation plotted. Figure 1(c) compares
the computation costs for different lower-level problem sizes. Our fully
first-order method significantly outperforms, in computation cost, the non-
fully first-order method implemented using differentiable optimization method.

The implementation can be found in https://github.com/guaguakai/constrained-
bilevel-optimization.

##  7 Limitations and future directions

One limitation to our approach is that the inexact gradient oracle we
constructed in Section 5 requires access to the exact dual multiplier
$\lambda^{*}$. For a first-order algorithm, the closest proxy one could get to
this would be a highly accurate approximation (which could be computed up to
$\epsilon$ error in $O(\log(1/\epsilon))$ gradient oracle evaluations).
Removing this “exact dual access” assumption
(LABEL:{item:assumption_safe_constraints}) would be an important result.

Another important direction for future work would be to extend the
hypergradient stationarity guarantee of Theorem 1.4 to bilevel programs with
_general convex_ constraints. To this end, we conjecture that the use of a
primal-only gradient approximation oracle could be potentially effective.

Finally, our current rate of $\widetilde{O}(\delta^{-1}\epsilon^{-4})$ oracle
calls for reaching $(\delta,\epsilon)$-Goldstein stationarity is not
necessarily inherent to the problem; indeed, it might be the case that an
alternate approach could improve it to the best known rate of
${O}(\delta^{-1}\epsilon^{-3})$ for nonsmooth nonconvex optimization.

## Acknowledgments

GK is supported by an Azrieli Foundation graduate fellowship. KW is supported
by Schmidt Sciences AI2050 Fellowship. SS acknowledges generous support from
the Alexander von Humboldt Foundation. SP is supported by NSF CCF-2112665
(TILOS AI Research Institute) and gratefully acknowledges Ali Jadbabaie for
many useful discussions. ZZ is supported by NSF FODSI Fellowship.

## References

  * [1] Jerome Bracken and James T McGill  “Mathematical programs with optimization problems in the constraints”  In _Operations research_ 21.1 INFORMS, 1973, pp. 37–44
  * [2] Benoît Colson, Patrice Marcotte and Gilles Savard  “An overview of bilevel optimization”  In _Annals of operations research_ 153 Springer, 2007, pp. 235–256
  * [3] Jonathan F Bard  “Practical bilevel optimization: algorithms and applications”  Springer Science & Business Media, 2013
  * [4] Ankur Sinha, Pekka Malo and Kalyanmoy Deb  “A review on bilevel optimization: From classical to evolutionary approaches and applications”  In _IEEE transactions on evolutionary computation_ 22.2 IEEE, 2017, pp. 276–295
  * [5] Jake Snell, Kevin Swersky and Richard Zemel  “Prototypical networks for few-shot learning”  In _Advances in neural information processing systems_ 30, 2017
  * [6] Luca Bertinetto, Joao F Henriques, Philip HS Torr and Andrea Vedaldi  “Meta-learning with differentiable closed-form solvers”  In _arXiv preprint arXiv:1805.08136_ , 2018
  * [7] Aravind Rajeswaran, Chelsea Finn, Sham M Kakade and Sergey Levine  “Meta-learning with implicit gradients”  In _Advances in neural information processing systems_ 32, 2019
  * [8] Kaiyi Ji, Jason D Lee, Yingbin Liang and H Vincent Poor  “Convergence of meta-learning with task-specific adaptation over partial parameters”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 11490–11500
  * [9] Luca Franceschi et al.  “Bilevel programming for hyperparameter optimization and meta-learning”  In _International conference on machine learning_ , 2018, pp. 1568–1577  PMLR
  * [10] Amirreza Shaban, Ching-An Cheng, Nathan Hatch and Byron Boots  “Truncated back-propagation for bilevel optimization”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 1723–1732  PMLR
  * [11] Matthias Feurer and Frank Hutter  “Hyperparameter optimization”  In _Automated machine learning: Methods, systems, challenges_ Springer International Publishing, 2019, pp. 3–33
  * [12] Vijay Konda and John Tsitsiklis  “Actor-critic algorithms”  In _Advances in neural information processing systems_ 12, 1999
  * [13] Richard S Sutton and Andrew G Barto  “Reinforcement learning: An introduction”  MIT press, 2018
  * [14] M Hong, H Wai, Z Wang and Z Yang  “A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic. arXiv e-prints, art”  In _arXiv preprint arXiv:2007.05170_ , 2020
  * [15] Haifeng Zhang et al.  “Bi-level actor-critic for multi-agent coordination”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 34.05, 2020, pp. 7325–7332
  * [16] Saeed Ghadimi and Mengdi Wang  “Approximation methods for bilevel programming”  In _arXiv preprint arXiv:1802.02246_ , 2018
  * [17] Brandon Amos and J Zico Kolter  “Optnet: Differentiable optimization as a layer in neural networks”  In _International Conference on Machine Learning_ , 2017, pp. 136–145  PMLR
  * [18] Akshay Agrawal et al.  “Differentiable convex optimization layers”  In _Advances in neural information processing systems_ 32, 2019
  * [19] Priya Donti, Brandon Amos and J Zico Kolter  “Task-based end-to-end model learning in stochastic optimization”  In _Advances in neural information processing systems_ 30, 2017
  * [20] Bryan Wilder, Bistra Dilkina and Milind Tambe  “Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 33.01, 2019, pp. 1658–1665
  * [21] James Kotary, Ferdinando Fioretto, Pascal Van Hentenryck and Bryan Wilder  “End-to-end constrained optimization learning: A survey”  In _arXiv preprint arXiv:2103.16378_ , 2021
  * [22] Kwonjoon Lee, Subhransu Maji, Avinash Ravichandran and Stefano Soatto  “Meta-learning with differentiable convex optimization”  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2019, pp. 10657–10665
  * [23] Bo Tang and Elias B Khalil  “Pyepo: A pytorch-based end-to-end predict-then-optimize library for linear and integer programming”  In _arXiv preprint arXiv:2206.14234_ , 2022
  * [24] Shaojie Bai, J Zico Kolter and Vladlen Koltun  “Deep equilibrium models”  In _Advances in neural information processing systems_ 32, 2019
  * [25] Akshay Mehra and Jihun Hamm  “Penalty method for inversion-free deep bilevel optimization”  In _Asian conference on machine learning_ , 2021, pp. 347–362  PMLR
  * [26] Kaiyi Ji, Junjie Yang and Yingbin Liang  “Bilevel optimization: Convergence analysis and enhanced design”  In _International conference on machine learning_ , 2021, pp. 4882–4892  PMLR
  * [27] Kai Wang et al.  “Learning mdps from features: Predict-then-optimize for sequential decision making by reinforcement learning”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 8795–8806
  * [28] Bo Liu et al.  “Bome! bilevel optimization made easy: A simple first-order approach”  In _Advances in neural information processing systems_ 35, 2022, pp. 17248–17262
  * [29] Kai Wang et al.  “Coordinating followers to reach better equilibria: End-to-end gradient descent for stackelberg games”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 36.5, 2022, pp. 5219–5227
  * [30] Nicola Basilico, Stefano Coniglio, Nicola Gatti and Alberto Marchesi  “Bilevel programming methods for computing single-leader-multi-follower equilibria in normal-form and polymatrix games”  In _EURO Journal on Computational Optimization_ 8.1 Springer, 2020, pp. 3–31
  * [31] Jiuping Xu, Yan Tu and Ziqiang Zeng  “Bilevel optimization of regional water resources allocation problem under fuzzy random environment”  In _Journal of Water Resources Planning and Management_ 139.3 American Society of Civil Engineers, 2013, pp. 246–264
  * [32] Walter J Gutjahr and Nada Dzubur  “Bi-objective bilevel optimization of distribution center locations considering user equilibria”  In _Transportation Research Part E: Logistics and Transportation Review_ 85 Elsevier, 2016, pp. 1–22
  * [33] Yue Zhang, Oded Berman, Patrice Marcotte and Vedat Verter  “A bilevel model for preventive healthcare facility network design with congestion”  In _IIE Transactions_ 42.12 Taylor & Francis, 2010, pp. 865–880
  * [34] Amir M Fathollahi-Fard, Mostafa Hajiaghaei-Keshteli, Reza Tavakkoli-Moghaddam and Neale R Smith  “Bi-level programming for home health care supply chain considering outsourcing”  In _Journal of Industrial Information Integration_ 25 Elsevier, 2022, pp. 100246
  * [35] Adam N Elmachtoub and Paul Grigas  “Smart “predict, then optimize””  In _Management Science_ 68.1 INFORMS, 2022, pp. 9–26
  * [36] Miguel Angel Muñoz, Salvador Pineda and Juan Miguel Morales  “A bilevel framework for decision-making under uncertainty with contextual information”  In _Omega_ 108 Elsevier, 2022, pp. 102575
  * [37] Prashant Khanduri et al.  “Linearly constrained bilevel optimization: A smoothed implicit gradient approach”  In _International Conference on Machine Learning_ , 2023, pp. 16291–16325  PMLR
  * [38] Wei Yao, Chengming Yu, Shangzhi Zeng and Jin Zhang  “Constrained Bi-Level Optimization: Proximal Lagrangian Value function Approach and Hessian-free Algorithm”, 2024  arXiv:2401.16164 [cs.LG]
  * [39] Zhaosong Lu and Sanyou Mei  “First-order penalty methods for bilevel optimization”, 2024  arXiv:2301.01716 [math.OC]
  * [40] Jeongyeol Kwon, Dohyun Kwon, Stephen Wright and Robert D Nowak  “A fully first-order method for stochastic bilevel optimization”  In _International Conference on Machine Learning_ , 2023, pp. 18083–18113  PMLR
  * [41] AA Goldstein  “Optimization of Lipschitz continuous functions”  In _Mathematical Programming_ 13 Springer, 1977, pp. 14–22
  * [42] Arkadij Semenovič Nemirovskij and David Borisovich Yudin  “Problem complexity and method efficiency in optimization”  Wiley-Interscience, 1983
  * [43] Jingzhao Zhang et al.  “Complexity of finding stationary points of nonconvex nonsmooth functions”  In _International Conference on Machine Learning_ , 2020, pp. 11173–11182  PMLR
  * [44] Guy Kornowski and Ohad Shamir  “Oracle complexity in nonsmooth nonconvex optimization”  In _Journal of Machine Learning Research_ 23.314, 2022, pp. 1–44
  * [45] Zhe Zhang and Guanghui Lan  “Solving Convex Smooth Function Constrained Optimization Is Almost As Easy As Unconstrained Optimization”  In _arXiv preprint arXiv:2210.05807_ , 2022
  * [46] Guy Kornowski and Ohad Shamir  “An algorithm with optimal dimension-dependence for zero-order nonsmooth nonconvex stochastic optimization”  In _Journal of Machine Learning Research_ 25.122, 2024, pp. 1–14
  * [47] Ashok Cutkosky, Harsh Mehta and Francesco Orabona  “Optimal stochastic non-smooth non-convex optimization through online-to-non-convex conversion”  In _International Conference on Machine Learning_ , 2023, pp. 6643–6670  PMLR
  * [48] Damek Davis et al.  “A gradient sampling method with complexity guarantees for Lipschitz functions in high and low dimensions”  In _Advances in neural information processing systems_ 35, 2022, pp. 6692–6703
  * [49] G Anandalingam and DJ White  “A solution method for the linear static Stackelberg problem using penalty functions”  In _IEEE Transactions on automatic control_ 35.10 IEEE, 1990, pp. 1170–1173
  * [50] Yo Ishizuka and Eitaro Aiyoshi  “Double penalty method for bilevel optimization problems”  In _Annals of Operations Research_ 34.1 Springer, 1992, pp. 73–88
  * [51] Douglas J White and G Anandalingam  “A penalty function approach for solving bi-level linear programs”  In _Journal of Global Optimization_ 3 Springer, 1993, pp. 397–419
  * [52] Luis Vicente, Gilles Savard and Joaquim Júdice  “Descent approaches for quadratic bilevel programming”  In _Journal of Optimization theory and applications_ 81.2 Springer, 1994, pp. 379–399
  * [53] DL Zhu  “Optimality conditions for bilevel programming problems”  In _Optimization_ 33.1 Taylor & Francis, 1995, pp. 9–27
  * [54] JJ Ye and DL Zhu  “Exact penalization and necessary optimality conditions for generalized bilevel programming problems”  In _SIAM Journal on optimization_ 7.2 SIAM, 1997, pp. 481–507
  * [55] Justin Domke  “Generic methods for optimization-based modeling”  In _Artificial Intelligence and Statistics_ , 2012, pp. 318–326  PMLR
  * [56] Fabian Pedregosa  “Hyperparameter optimization with approximate gradient”  In _International conference on machine learning_ , 2016, pp. 737–746  PMLR
  * [57] Stephen Gould et al.  “On differentiating parameterized argmin and argmax problems with application to bi-level optimization”  In _arXiv preprint arXiv:1607.05447_ , 2016
  * [58] Renjie Liao et al.  “Reviving and improving recurrent back-propagation”  In _International Conference on Machine Learning_ , 2018, pp. 3082–3091  PMLR
  * [59] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil and Saverio Salzo  “On the iteration complexity of hypergradient computation”  In _International Conference on Machine Learning_ , 2020, pp. 3748–3758  PMLR
  * [60] Jonathan Lorraine, Paul Vicol and David Duvenaud  “Optimizing millions of hyperparameters by implicit differentiation”  In _International conference on artificial intelligence and statistics_ , 2020, pp. 1540–1552  PMLR
  * [61] Dougal Maclaurin, David Duvenaud and Ryan Adams  “Gradient-based hyperparameter optimization through reversible learning”  In _International conference on machine learning_ , 2015, pp. 2113–2122  PMLR
  * [62] Luca Franceschi, Michele Donini, Paolo Frasconi and Massimiliano Pontil  “Forward and reverse gradient-based hyperparameter optimization”  In _International Conference on Machine Learning_ , 2017, pp. 1165–1173  PMLR
  * [63] Risheng Liu et al.  “A value-function-based interior-point method for non-convex bi-level optimization”  In _International Conference on Machine Learning_ , 2021, pp. 6882–6892  PMLR
  * [64] Jane J Ye, Xiaoming Yuan, Shangzhi Zeng and Jin Zhang  “Difference of convex algorithms for bilevel programs with applications in hyperparameter selection”  In _Mathematical Programming_ 198.2 Springer Berlin Heidelberg Berlin/Heidelberg, 2023, pp. 1583–1616
  * [65] Lucy L. Gao et al.  “Moreau Envelope Based Difference-of-weakly-Convex Reformulation and Algorithm for Bilevel Programs”, 2024  arXiv:2306.16761 [math.OC]
  * [66] Tianyi Chen, Yuejiao Sun and Wotao Yin  “Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 25294–25307
  * [67] Tianyi Chen, Yuejiao Sun, Quan Xiao and Wotao Yin  “A single-timescale method for stochastic bilevel optimization”  In _International Conference on Artificial Intelligence and Statistics_ , 2022, pp. 2466–2488  PMLR
  * [68] Mingyi Hong, Hoi-To Wai, Zhaoran Wang and Zhuoran Yang  “A two-timescale stochastic algorithm framework for bilevel optimization: Complexity analysis and application to actor-critic”  In _SIAM Journal on Optimization_ 33.1 SIAM, 2023, pp. 147–180
  * [69] Zeeshan Akhtar, Amrit Singh Bedi, Srujan Teja Thomdapu and Ketan Rajawat  “Projection-free stochastic bi-level optimization”  In _IEEE Transactions on Signal Processing_ 70 IEEE, 2022, pp. 6332–6347
  * [70] Ruichen Jiang, Nazanin Abolfazli, Aryan Mokhtari and Erfan Yazdandoost Hamedani  “A conditional gradient-based method for simple bilevel optimization with convex lower-level problem”  In _International Conference on Artificial Intelligence and Statistics_ , 2023, pp. 10305–10323  PMLR
  * [71] Nazanin Abolfazli, Ruichen Jiang, Aryan Mokhtari and Erfan Yazdandoost Hamedani  “An Inexact Conditional Gradient Method for Constrained Bilevel Optimization”  In _arXiv preprint arXiv:2306.02429_ , 2023
  * [72] Jincheng Cao et al.  “Projection-free methods for stochastic simple bilevel optimization with convex lower-level problem”  In _Advances in Neural Information Processing Systems_ 36, 2024
  * [73] Prashant Khanduri et al.  “A near-optimal algorithm for stochastic bilevel optimization via double-momentum”  In _Advances in neural information processing systems_ 34, 2021, pp. 30271–30283
  * [74] Zhishuai Guo, Quanqi Hu, Lijun Zhang and Tianbao Yang  “Randomized stochastic variance-reduced methods for multi-task stochastic bilevel optimization”  In _arXiv preprint arXiv:2105.02266_ , 2021
  * [75] Junjie Yang, Kaiyi Ji and Yingbin Liang  “Provably faster algorithms for bilevel optimization”  In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 13670–13682
  * [76] Mathieu Dagréou, Pierre Ablin, Samuel Vaiter and Thomas Moreau  “A framework for bilevel optimization that enables stochastic and global variance reduction algorithms”  In _Advances in Neural Information Processing Systems_ 35, 2022, pp. 26698–26710
  * [77] Shoham Sabach and Shimrit Shtern  “A first order method for solving convex bilevel optimization problems”  In _SIAM Journal on Optimization_ 27.2 SIAM, 2017, pp. 640–660
  * [78] Mostafa Amini and Farzad Yousefian  “An iterative regularized incremental projected subgradient method for a class of bilevel optimization problems”  In _2019 American Control Conference (ACC)_ , 2019, pp. 4069–4074  IEEE
  * [79] Mostafa Amini and Farzad Yousefian  “An iterative regularized mirror descent method for ill-posed nondifferentiable stochastic optimization”  In _arXiv preprint arXiv:1901.09506_ , 2019
  * [80] Roey Merchav and Shoham Sabach  “Convex Bi-level Optimization Problems with Nonsmooth Outer Objective Function”  In _SIAM Journal on Optimization_ 33.4 SIAM, 2023, pp. 3114–3142
  * [81] Ioannis Tsaknakis, Prashant Khanduri and Mingyi Hong  “An implicit gradient-type method for linearly constrained bilevel problems”  In _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 5438–5442  IEEE
  * [82] Siyuan Xu and Minghui Zhu  “Efficient gradient approximation method for constrained bilevel optimization”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 37.10, 2023, pp. 12509–12517
  * [83] Lesi Chen, Yaohua Ma and Jingzhao Zhang  “Near-Optimal Fully First-Order Algorithms for Finding Stationary Points in Bilevel Optimization”  In _arXiv preprint arXiv:2306.14853_ , 2023
  * [84] Frank H Clarke  “Generalized gradients of Lipschitz functionals”  In _Advances in Mathematics_ 40.1 Elsevier, 1981, pp. 52–67
  * [85] Alexander Shapiro, Darinka Dentcheva and Andrzej Ruszczynski  “Lectures on stochastic programming: modeling and theory”  SIAM, 2021
  * [86] Michael Jordan et al.  “Deterministic nonsmooth nonconvex optimization”  In _The Thirty Sixth Annual Conference on Learning Theory_ , 2023, pp. 4570–4597  PMLR
  * [87] Siyu Kong and AS Lewis  “The cost of nonconvexity in deterministic nonsmooth optimization”  In _Mathematics of Operations Research_ INFORMS, 2023
  * [88] Benjamin Grimmer and Zhichao Jia  “Goldstein Stationarity in Lipschitz Constrained Optimization”  In _arXiv preprint arXiv:2310.03690_ , 2023
  * [89] Olivier Devolder, François Glineur and Yurii Nesterov  “First-order methods of smooth convex optimization with inexact oracle”  In _Mathematical Programming_ 146 Springer, 2014, pp. 37–75
  * [90] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _arXiv preprint arXiv:1412.6980_ , 2014
  * [91] Guanghui Lan  “First-order and stochastic optimization methods for machine learning”  Springer, 2020
  * [92] Ohad Shamir  “An optimal algorithm for bandit and zero-order convex optimization with two-point feedback”  In _The Journal of Machine Learning Research_ 18.1 JMLR. org, 2017, pp. 1703–1713
  * [93] Simon S Du and Wei Hu  “Linear convergence of the primal-dual gradient method for convex-concave saddle point problems without strong convexity”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 196–205  PMLR
  * [94] Adam Paszke et al.  “Pytorch: An imperative style, high-performance deep learning library”  In _Advances in neural information processing systems_ 32, 2019
  * [95] Steven Diamond and Stephen Boyd  “CVXPY: A Python-embedded modeling language for convex optimization”  In _Journal of Machine Learning Research_ 17.83, 2016, pp. 1–5

## Appendix

##  Appendix A Proofs from Section 3

In this section, we provide the full proofs of claims for bilevel programs
with linear equality constraints, as stated in Section 3. We first state a few
technical results using the implicit function theorem that we repeatedly
invoke in our results for this setting.

######  Lemma A.1.

Fix a point $x$. Given $y^{*}=\arg\min_{y:h(x,y)=0}g(x,y)$ where $g$ is
strongly convex in $y$ and $\lambda^{*}$ is the dual optimal variable for this
problem, define
$\mathcal{L}_{\mathrm{eq}}(x,y,\lambda)=g(x,y)+\langle\lambda,h(x,y)\rangle$.
Then, we have

$\underbrace{\begin{bmatrix}\nabla^{2}_{yy}\mathcal{L}_{\mathrm{eq}}(x,y^{*},\lambda^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}}_{H~{}\text{for linear equality
constraints}}\begin{bmatrix}\frac{dy^{*}}{dx}\\\
\frac{d\lambda^{*}}{dx}\end{bmatrix}=\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})-\nabla^{2}_{yx}\langle\lambda^{*},h(x,y^{*})\rangle\\\
-\nabla_{x}h(x,y^{*})\end{bmatrix}.$

###### Proof.

Since $g$ is strongly convex, by the linear constraint qualification
condition, the KKT condition is both sufficient and necessary for optimality.
Hence, consider the following KKT system obtained via first order optimality
of $y^{*}$, with dual optimal variable $\lambda^{*}$:

$\displaystyle\nabla_{y}g(x,y^{*})+\nabla_{y}\langle\lambda^{*},h(x,y^{*})\rangle=0,\text{
and }h(x,y^{*})=0.$ (A.1)

Differentiating the system of equations in Eq. A.1 with respect to $x$ and
rearranging terms yields:

$\displaystyle\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})+\nabla^{2}_{yy}\langle\lambda^{*},h(x,y^{*})\rangle&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}\begin{bmatrix}\frac{dy^{*}}{dx}\\\
\frac{d\lambda^{*}}{dx}\end{bmatrix}=\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})-\nabla^{2}_{yx}\langle\lambda^{*},h(x,y^{*})\rangle\\\
-\nabla_{x}h(x,y^{*})\end{bmatrix}$ (A.2)

Finally, noting that
$\nabla_{yy}^{2}\mathcal{L}_{\mathrm{eq}}(x,y,\lambda)=\nabla^{2}_{yy}g(x,y)+\nabla^{2}_{yy}\langle\lambda,h(x,y)\rangle$,
we can write Eq. A.2 as claimed. ∎

######  Lemma A.2.

Consider the setup in LABEL:{lem:dystarDxLinEq}. The matrix $H$ defined in
Section 3.1 is invertible if the Hessian
$\nabla_{yy}^{2}\mathcal{L}_{\mathrm{eq}}(x,y^{*},\lambda^{*}):=\nabla^{2}_{yy}g(x,y^{*})+\nabla^{2}_{yy}\langle\lambda^{*},h(x,y^{*})\rangle$
satisfies $\nabla_{yy}^{2}\mathcal{L}_{\mathrm{eq}}(x,y^{*},\lambda^{*})\succ
0$ over the tangent plane $T:=\\{y:\nabla_{y}h(x,y^{*})y=0\\}$ and
$\nabla_{y}h$ has full rank.

###### Proof.

Let $u=[y,\lambda].$ We show that $Hu=0$ implies $u=0$, which in turn implies
invertibility of $H$. If $\nabla_{y}h(x,y^{*})y\neq 0,$ then by construction
of $u$ and $H$, we must also have $Hu\neq 0$. Otherwise if
$\nabla_{y}h(x,y^{*})y=0$ and $y\neq 0$, the quadratic form $u^{\top}Hu$ is
positive, as seen by

$u^{\top}Hu=y^{\top}\nabla^{2}_{yy}\mathcal{L}_{\mathrm{eq}}(x,y^{*},\lambda^{*})y>0,$

where the final step is by the assumption of $\mathcal{L}_{\mathrm{eq}}$ being
positive definite over the defined tangent plane
$T=\\{y:\nabla_{y}h(x,y^{*})y=0\\}$. If $y=0$ while $Hu=0$, then $\nabla_{y}h$
having full rank implies $\lambda=0$. Combined with $y=0$, this means $u=0$,
as required when $Hu=0$. This concludes the proof. ∎

######  Corollary A.3.

For Eq. 3.1 under Assumption 2.2 and Assumption 2.3, the matrix $H$ (as
defined in Section 3.1) is non-singular. Further, there exists a finite
$C_{H}$ such that $\|H^{-1}\|\leq C_{H}$.

###### Proof.

Since we are assuming strong convexity of $g$, LABEL:{lem:non-singular-req}
applies, yielding the claimed invertibility of $H$. Combined with the
boundedness of variables $x$ (per Assumption 2.3) and continuity of the
inverse implies a bound on $\|H^{-1}\|$. ∎

###  A.1 Construction of the inexact gradient oracle

We now show how to construct the inexact gradient oracle for the objective $F$
in Eq. 3.1. As sketched in Section 3, we then use this oracle in a projected
gradient descent algorithm to get the claimed guarantee.

######  Lemma A.4.

Consider Eq. 3.1 under Assumption 2.2 and Assumption 2.3. Let $y_{\delta}^{*}$
be as defined in Eq. 3.3. Then, for any $\delta\in[0,\Delta]$ with
$\Delta\leq\mu_{g}/2C_{f}$, the following relation is valid:

$\|y_{\delta}^{*}-y^{*}\|\leq M(x)\delta,\textrm{ with
}M(x):=\frac{2}{\mu_{g}}\|\nabla_{y}f(x,y^{*})\|\leq\frac{2L_{f}}{\mu_{g}}.$

###### Proof.

The first-order optimality condition applied to $g(x,y)+\delta f(x,y)$ at
$y^{*}$ and $y_{\delta}^{*}$ gives

$\langle\nabla_{y}g(x,y_{\delta}^{*})+\delta\nabla_{y}f(x,y_{\delta}^{*}),y^{*}-y_{\delta}^{*}\rangle\geq
0,$

which upon adding and subtracting $\nabla_{y}f(x,y^{*})$ transforms into

$\langle\nabla_{y}g(x,y_{\delta}^{*})+\delta[\nabla_{y}f(x,y_{\delta}^{*})-\nabla_{y}f(x,y^{*})]+\delta\nabla_{y}f(x,y^{*}),y^{*}-y_{\delta}^{*}\rangle\geq
0.$

Similarly, the first-order optimality condition applied to $g$ at $y^{*}$ and
$y_{\delta}^{*}$ gives

$\langle\nabla_{y}g(x,y^{*}),y_{\delta}^{*}-y^{*}\rangle\geq 0.$

Adding Section A.1 and Section A.1 and rearranging yields

$\displaystyle\langle\nabla_{y}g(x,y_{\delta}^{*})-\nabla_{y}g(x,y^{*})+\delta[\nabla_{y}f(x,y_{\delta}^{*})-\nabla_{y}f(x,y^{*})],y_{\delta}^{*}-y^{*}\rangle$
$\displaystyle\leq\langle\delta\nabla_{y}f(x,y^{*}),y^{*}-y_{\delta}^{*}\rangle.$

Applying to the left side above a lower bound via strong convexity of
$g+\delta f$ and to the right hand side an upper bound via Cauchy-Schwarz
inequality, we have

$s\|y_{\delta}^{*}-y^{*}\|\leq\delta\|\nabla_{y}f(x,y^{*})\|,$

where $s$ is the strong convexity of $g+\delta f$. Since $f$ is
$C_{f}$-smooth, the worst case value of this is $s=\mu_{g}-\delta
C_{f}=\mu_{g}-\frac{\mu_{g}}{2C_{f}}C_{f}=\mu_{g}/2$, which when plugged in
Section A.1 then gives the claimed bound. ∎

######  Lemma A.5.

Consider Eq. 3.1 under Assumption 2.2 and Assumption 2.3. Then the following
relation is valid.

$\lim_{\delta\rightarrow
0}\frac{\nabla_{x}[g(x,y_{\delta}^{*}(x))+\lambda_{\delta}^{*}h(x,y^{*})]-\nabla_{x}[g(x,y^{*}(x))+\lambda^{*}h(x,y^{*})]}{\delta}=\left(\frac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*}(x)).$

###### Proof.

Recall that by definition, $g$ is strongly convex and
$y^{*}=\arg\min_{y:h(x,y)=0}g(x,y)$. Hence, we can apply Lemma A.1. Combining
this with Lemma A.2 and further applying that linearity of $h$ implies
$\nabla^{2}_{yy}h=0$ and $\nabla^{2}_{xy}h=0$, we obtain the following:

$\begin{bmatrix}\frac{dy^{*}}{dx}\\\
\frac{d\lambda^{*}}{dx}\end{bmatrix}=\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})\\\
-\nabla_{x}h(x,y^{*})\end{bmatrix}.$

So we can express the right-hand side of the claimed equation in the lemma
statement by

$\displaystyle\left(\frac{dy^{*}(x)}{dx}\right)^{\top}\nabla_{y}f(x,y^{*}(x))$
$\displaystyle=\begin{bmatrix}\left(\frac{dy^{*}}{dx}\right)^{\top}&\left(\frac{d\lambda^{*}}{dx}\right)^{\top}\end{bmatrix}\begin{bmatrix}\nabla_{y}f(x,y^{*}(x))\\\
0\end{bmatrix},$

which can be further simplified to

$\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})^{\top}&-\nabla_{x}h(x,y^{*})^{\top}\end{bmatrix}\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}\begin{bmatrix}\nabla_{y}f(x,y^{*}(x))\\\
0\end{bmatrix}.$

We now apply Lemma A.1 to the perturbed problem defined in Eq. 3.3. We know
from Lemma A.4 that $\lim_{\delta\rightarrow 0}y_{\delta}^{*}=y^{*}$. The
associated KKT system is given by

$\displaystyle\delta\nabla_{y}f(x,y_{\delta}^{*})+\nabla_{y}g(x,y_{\delta}^{*})+\nabla_{y}\langle\lambda_{\delta}^{*},h(x,y_{\delta}^{*})\rangle=0\text{
and }h(x,y_{\delta}^{*})=0.$ (A.6)

Taking the derivative with respect of Eq. A.6 gives the following implicit
system, where we used the fact that $h$ is linear and hence
$\nabla^{2}_{yy}h=0$:

$\underbrace{\begin{bmatrix}\delta\nabla^{2}_{yy}f(x,y_{\delta}^{*})+\nabla^{2}_{yy}g(x,y_{\delta}^{*})&\nabla_{y}h(x,y_{\delta}^{*})^{\top}\\\
\nabla_{y}h(x,y_{\delta}^{*})&0\end{bmatrix}}_{H_{\delta}}\begin{bmatrix}\frac{dy_{\delta}^{*}}{d\delta}\\\
\frac{d\lambda_{\delta}^{*}}{d\delta}\end{bmatrix}=\begin{bmatrix}-\nabla_{y}f(x,y_{\delta}^{*})^{\top}\\\
0\end{bmatrix}.$ (A.7)

For a sufficiently small $\delta$, we have
$\nabla^{2}_{yy}g(x,y_{\delta}^{*})+\delta\nabla^{2}_{yy}f(x,y_{\delta}^{*})\succeq\tfrac{\mu_{g}}{2}I$,
which implies invertibility of $H_{\delta}$ by an application of
LABEL:{lem:non-singular-req}. Since Lemma A.4 implies $\lim_{\delta\rightarrow
0}y_{\delta}^{*}=y^{*}$, we get

$\begin{bmatrix}\frac{dy_{\delta}^{*}}{d\delta}\\\
\frac{d\lambda_{\delta}^{*}}{d\delta}\end{bmatrix}|_{\delta=0}=\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}\begin{bmatrix}-\nabla_{y}f(x,y^{*})\\\
0\end{bmatrix}.$

So we can express the left-hand side of the expression in the lemma statement
by

$\displaystyle\lim_{\delta\rightarrow
0}\frac{\nabla_{x}[g(x,y_{\delta}^{*}(x))+\langle\lambda_{\delta}^{*},h(x,y^{*})\rangle]-\nabla_{x}[g(x,y^{*}(x))+\langle\lambda^{*},h(x,y^{*})\rangle]}{\delta}$
$\displaystyle=\nabla^{2}_{xy}g(x,y^{*})\frac{dy_{\delta}^{*}}{d\delta}+\nabla_{x}h(x,y^{*})^{\top}\frac{d\lambda_{\delta}^{*}}{d\delta}$
$\displaystyle=\begin{bmatrix}\nabla^{2}_{xy}g(x,y^{*})&\nabla_{x}h(x,y^{*})^{\top}\end{bmatrix}\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&\nabla_{y}h(x,y^{*})^{\top}\\\
\nabla_{y}h(x,y^{*})&0\end{bmatrix}^{-1}\begin{bmatrix}-\nabla_{y}f(x,y^{*})\\\
0\end{bmatrix},$

which matches Section A.1 (since
$(\nabla^{2}_{yx}g)^{\top}=\nabla^{2}_{xy}g$), thus concluding the proof. ∎

See 3.3

###### Proof.

Rearranging LABEL:{{eq:def_matrix_H}} and applying Corollary A.3, we have

$\begin{bmatrix}\frac{dy^{*}}{dx}\\\
\frac{d\lambda^{*}}{dx}\end{bmatrix}=\begin{bmatrix}\nabla^{2}_{yy}g(x,y^{*})&B^{\top}\\\
B&0\end{bmatrix}^{-1}\begin{bmatrix}-\nabla^{2}_{yx}g(x,y^{*})\\\
-\nabla_{x}h(x,y^{*})\end{bmatrix}.$

This implies a Lipschitz bound of $C_{H}\cdot(C_{g}+\|A\|)$. Next, note that
in the case with linear equality constraints, the terms in
LABEL:{{eqs:differentiated-kkt-lin-eq}} involving second-order derivatives of
$h$ are all zero; differentiating LABEL:{{eqs:differentiated-kkt-lin-eq}} with
respect to $x$, we notice that the linear system we get again has the same
matrix $H$ from before. We can therefore again perform the same inversion and
apply the bound on $\|H^{-1}\|$ and on the third-order derivatives of $g$
(Assumption 2.3) to observe that $\|\frac{d^{2}y^{*}}{dx^{2}}\|\leq
O(C_{H}\cdot S_{g}\|\frac{dy^{*}}{dx}\|^{2})=O(C_{H}^{3}\cdot
S_{g}\cdot(C_{g}+\|A\|)^{2})$, where we are hiding numerical constants in the
Big-Oh notation.

As a result, we can calculate the Lipschitz smoothness constant associated
with the hyper-objective $F$ by

$\displaystyle\|\nabla F(x)-\nabla F(\bar{x})\|$
$\displaystyle\leq\|\frac{dy^{*}(x)}{dx}\nabla_{y}f(x,y^{*}(x))-\frac{dy^{*}(\bar{x})}{dx}\nabla_{y}f(\bar{x},y^{*}(\bar{x}))\|+\|\nabla_{x}f(x,y^{*}(x))-\nabla_{x}f(\bar{x},y^{*}(\bar{x}))\|$
$\displaystyle\leq[C_{f}C_{H}(L_{g}+\|A\|)+C_{f}C^{2}_{H}(L_{g}+\|A\|)^{2}+L_{f}C_{H}^{3}S_{g}(L_{g}+\|A\|)^{2}]\|x-\bar{x}\|$
$\displaystyle\ \ +[C_{f}+C_{f}C_{H}(L_{g}+\|A\|)]\|x-\bar{x}\|$
$\displaystyle\leq\underbrace{2(L_{f}+C_{f}+C_{g})C_{H}^{3}S_{g}(L_{g}+\|A\|)^{2}}_{C_{F}}\|x-\bar{x}\|.$

∎

See 3.2

###### Proof.

For simplicity, we adopt the following notation throughout this proof:
$g_{xy}(x,y)=\nabla^{2}_{xy}g,$ and $g_{xyy}$ denotes the tensor such that its
$ijk$ entry is given by $\frac{\partial^{3}g}{\partial x_{i}\partial
y_{j}\partial y_{k}}$. We first consider the terms involving $g$. By the
fundamental theorem of calculus, we have

$\nabla_{x}g(x,y_{\delta}^{*}(x))-\nabla_{x}g(x,y^{*}(x))=\int_{t=0}^{\delta}g_{xy}(x,y_{t}^{*}(x))\frac{dy_{t}^{*}(x)}{dt}dt.$

As a result, we have

$\displaystyle\frac{\nabla_{x}g(x,y_{\delta}^{*}(x))-\nabla_{x}g(x,y^{*}(x))}{\delta}-g_{xy}(x,y^{*}(x))\frac{dy_{t}^{*}(x)}{dt}|_{t=0}$
$\displaystyle=\frac{1}{\delta}\int_{t=0}^{\delta}\left(g_{xy}(x,y_{t}^{*}(x))\frac{dy_{t}^{*}(x)}{dt}-g_{xy}(x,y^{*}(x))\frac{dy_{t}^{*}(x)}{dt}|_{t=0}\right)dt$
$\displaystyle=\frac{1}{\delta}\int_{t=0}^{\delta}\left(g_{xy}(x,y_{t}^{*}(x))\frac{dy_{t}^{*}(x)}{dt}-g_{xy}(x,y^{*}(x))\frac{dy_{t}^{*}(x)}{dt}|_{t=0}\right)dt$
$\displaystyle=\frac{1}{\delta}\int_{t=0}^{\delta}\left(g_{xy}(x,y_{t}^{*}(x))-g_{xy}(x,y^{*}(x))\right)\frac{dy_{t}^{*}(x)}{dt}dt+\frac{1}{\delta}\int_{t=0}^{\delta}g_{xy}(x,y^{*}(x))\cdot\left(\frac{dy_{t}^{*}(x)}{dt}-\frac{dy_{t}^{*}(x)}{dt}|_{t=0}\right)dt.$
(A.8)

We now bound each of the terms on the right-hand side of Eq. A.8. For the
first term, we have

$\displaystyle\|\frac{1}{\delta}\int_{t=0}^{\delta}\left(g_{xy}(x,y_{t}^{*}(x))-g_{xy}(x,y^{*}(x))dt\right)\frac{dy_{t}^{*}(x)}{dt}\|$
$\displaystyle\leq\frac{1}{\delta}\int_{t=0}^{\delta}\|\frac{dy_{t}^{*}(x)}{dt}\|\cdot\int_{s=0}^{t}\|g_{xyy}(x,y_{s}^{*}(x))\|\|\frac{dy_{s}^{*}(x)}{ds}\|ds\cdot
dt$
$\displaystyle\leq\frac{1}{\delta}\int_{t=0}^{\delta}\|\frac{dy_{t}^{*}(x)}{dt}\|\cdot\max_{s\in[0,\delta]}\|g_{xyy}(x,y_{s}^{*}(x))\|\cdot\|\frac{dy_{s}^{*}(x)}{ds}\|tdt$
$\displaystyle\leq\frac{1}{\delta}\cdot\max_{u\in[0,\delta]}\|g_{xyy}(x,y_{u}^{*}(x))\|\cdot\delta^{2}\cdot\max_{t\in[0,\delta]}\|\frac{dy_{t}^{*}(x)}{dt}\|^{2}$
$\displaystyle\leq\delta\cdot\max_{u\in[0,\delta]}\|g_{xyy}(x,y_{u}^{*}(x))\|\cdot\max_{t\in[0,\delta]}\|\frac{dy_{t}^{*}(x)}{dt}\|^{2}$
$\displaystyle=\delta\cdot S_{g}\cdot M_{y}^{2},$ (A.9)

where $M_{y}$ is the Lipschitz bound on $y^{*}$ as shown in Lemma 3.3, and
$S_{g}$ is the smoothness of $g$ from Assumption 2.3. For the second term on
the right-hand side of Eq. A.8, we have

$\displaystyle\|\frac{1}{\delta}\int_{t=0}^{\delta}g_{xy}(x,y^{*}(x))\cdot\left(\frac{dy_{t}^{*}(x)}{dt}-\frac{dy^{*}(x)}{dt}\right)\|$
$\displaystyle\leq\frac{1}{\delta}\cdot\|g_{xy}(x,y^{*}(x))\|\cdot\int_{t=0}^{\delta}\left(\int_{s=0}^{t}\|\frac{d^{2}}{ds^{2}}y_{s}^{*}(x)\|ds\right)dt$
$\displaystyle\leq\frac{1}{\delta}\cdot\|g_{xy}(x,y^{*}(x))\|\cdot\max_{s\in[0,\delta]}\|\frac{d^{2}}{ds^{2}}y_{s}^{*}(x)\|\cdot\delta^{2}$
$\displaystyle\leq\delta\cdot\|g_{xy}(x,y^{*}(x))\|\cdot\max_{s\in[0,\delta]}\|\frac{d^{2}}{ds^{2}}y_{s}^{*}(x)\|$
$\displaystyle=\delta\cdot C_{g}\cdot C_{y},$ (A.10)

where $C_{g}$ is the bound on smoothness of $g$ as in Assumption 2.3, and
$C_{y}$ is the bound on $\|\frac{d^{2}y^{*}}{dx^{2}}\|$ from Lemma 3.3. For
the terms involving the function $h$, we have

$\displaystyle\|\frac{\lambda_{\delta}^{*}-\lambda^{*}}{\delta}-\frac{d\lambda_{\delta}^{*}}{d\delta}|_{\delta=0}\|$
$\displaystyle=\frac{1}{\delta}\int_{t=0}^{\delta}\|\frac{d\lambda_{t}^{*}}{dt}-\frac{d\lambda_{\delta}^{*}}{d\delta}|_{\delta=0}\|dt$
$\displaystyle=\frac{1}{\delta}\int_{t=0}^{\delta}\int_{s=0}^{t}\|\frac{d^{2}}{ds^{2}}\lambda_{s}^{*}\|ds\cdot
dt$
$\displaystyle\leq\frac{1}{\delta}\max_{s\in[0,\delta]}\|\frac{d^{2}}{ds^{2}}\lambda_{s}^{*}\|\cdot\delta^{2}\leq\delta\cdot\max_{s\in[0,\delta]}\|\frac{d^{2}}{ds^{2}}\lambda_{s}^{*}\|$
$\displaystyle=\delta\cdot C_{\ell},$ (A.11)

where $C_{\ell}$ is the bound on $\|\frac{d^{2}\lambda^{*}}{ds^{2}}\|$ from
Lemma 3.3. Combining Eq. A.8, Eq. A.9, Eq. A.10, and Eq. A.11, along with
LABEL:{lem:lineq-in-limit-finitediff-equals-gradf}, Corollary A.3, and
LABEL:{lem:smoothness_of_ydelstar_lamdelstar}, we have that overall bound is

$\displaystyle\delta\cdot(S_{g}M_{y}^{2}+C_{g}C_{y}+C_{\ell})$
$\displaystyle\leq O(\delta\cdot(S_{g}\cdot
C_{H}^{3}\cdot(C_{g}+\|A\|)^{2}\cdot(C_{g}+C_{f}+L_{f}))).$

∎

###  A.2 The cost of linear equality constrained bilevel program

We present here our full (implementable) algorithm for LABEL:{prob:lin-eq} and
restate and prove its convergence guarantee.

Algorithm 6 The Fully First-Order Method for Bilevel Equality Constrained
Problem

1:Input: Current $x_{0}$, accuracy $\epsilon$, perturbation
$\delta=\epsilon^{2}/8C^{2}_{F}R_{\mathcal{X}}$ with
$C_{F}=2(L_{f}+C_{f}+C_{g})C_{H}^{3}S_{g}(L_{g}+\|A\|)^{2}$, accuracy for the
lower level problem $\tilde{\delta}=2(C_{g}+\|A\|)\delta^{2}$.

2:for t=0,1,2,… do

3: Run Algorithm 7 to generate $\tilde{\delta}$-accurate primal and dual
solutions $(\hat{y}^{*},\hat{\lambda}^{*})$ for
$\min_{y:Ax_{t}+By=b}g(x_{t},y)$

4: Run Algorithm 7 to generate $\tilde{\delta}$-accurate primal and dual
solutions $(\hat{y}_{\delta}^{*},\hat{\lambda}_{\delta}^{*})$ for
$\min_{y:Ax_{t}+By=b}g(x_{t},y)+\delta f(x_{t},y)$

5: Compute
$\hat{v}_{t}:=\frac{\nabla_{x}[g(x_{t},\hat{y}_{\delta}^{*})+\hat{\lambda}_{\delta}^{*}h(x,\hat{y}^{*})]-\nabla_{x}[g(x_{t},\hat{y}^{*})+\hat{\lambda}^{*}h(x,\hat{y}^{*})]}{\delta}$,
set  $\widetilde{\nabla}F(x_{t}):=\hat{v}^{t}+\nabla_{x}f(x,\hat{y}^{*}(x)).$

6: Set
$x_{t+1}\leftarrow\arg\min_{z\in\mathcal{X}}\|z-(x_{t}-\frac{1}{C_{F}}\widetilde{\nabla}F(x_{t}))\|^{2}.$

See 3.1

###### Proof.

We first show the inexact gradient $\widetilde{\nabla}F(x_{t})$ generated in
Algorithm 6 is an $\delta$-accurate approximation to the hyper-gradient
$\nabla F(x_{t}).$ Consider the inexact gradient defined in (3)

$\displaystyle\|v_{t}-\hat{v}_{t}\|$
$\displaystyle\leq\frac{1}{\delta}\\{\|[\nabla_{x}g(x_{t},\hat{y}_{\delta}^{*})-\nabla_{x}[g(x_{t},\hat{y}^{*})]-[\nabla_{x}g(x_{t},y_{\delta}^{*})-\nabla_{x}[g(x_{t},y^{*})\|$
$\displaystyle\,\,\,+\|\hat{\lambda}_{\delta}^{*}-\hat{\lambda}^{*}-[\lambda_{\delta}^{*}-\lambda^{*}\|\|A\|\\}$
$\displaystyle\leq\frac{2}{\delta}[C_{g}+\|A\|]\tilde{\delta}.$

Thus we get

$\displaystyle\|\widetilde{\nabla}F(x_{t})-\nabla F(x_{t})\|$
$\displaystyle\leq\|\nabla_{x}f(x_{t},y^{*})-\nabla_{x}f(x_{t},\hat{y}^{*})\|+\left\lVert\hat{v}^{t}-v^{t}\right\rVert+\|v^{t}-\frac{dy^{*}(x^{t})}{dx}\nabla_{y}f(x_{t},y^{*}(x_{t}))\|$
$\displaystyle\leq C_{f}\tilde{\delta}+\frac{2}{\delta}[C_{g}+\left\lVert
A\right\rVert]\tilde{\delta}+C_{F}\delta$
$\displaystyle\leq\frac{2\tilde{\delta}}{\delta}[C_{f}+C_{g}+\|A\|]+C_{F}\delta$
$\displaystyle\leq\frac{\epsilon^{2}}{4C_{F}R_{\mathcal{X}}}.$

Applied to the $C_{F}$-smooth hyper-objective $F$, such an inexact gradient
oracle satisfies the requirement for Proposition A.6. Thus an
$\epsilon$-stationary point with $\|\mathcal{G}_{F}(x^{t})\|\leq\epsilon$ (see
Eq. (A.13) for the definition of gradient mapping) must be found in
$N=O(\frac{C_{F}(F(x^{0})-F^{*})}{\epsilon^{2}})$ iterations. Noting the
evaluation of inexact solutions
$(\hat{y}^{*},\hat{\lambda}^{*},\hat{y}_{\delta}^{*},\hat{\lambda}_{\delta}^{*})$
requires $\tilde{O}(\sqrt{C_{g}/\mu_{g}})$ first order oracle evaluations, we
arrive at the total oracle complexity of
$\tilde{O}(\sqrt{C_{g}/\mu_{g}}\frac{C_{F}(F(x^{0})-F^{*})}{\epsilon^{2}})$
for finding an $\epsilon$-stationary point. ∎

###  A.3 The cost of inexact projected gradient descent method

In this subsection, we state the number of iterations required by projected
gradient descent method to find an $\epsilon$-stationary point using inexact
gradient oracles. Specifically, we consider the following non-convex smooth
problem where the objective $F$ is assumed to be $C_{F}$-Lipschitz smooth:

$\mbox{minimize}_{x\in\mathcal{X}}F(x).$ (A.12)

Since the feasible region $\mathcal{X}$ is compact, we use the norm of the
following gradient mapping $\mathcal{G}_{F}(x)$ as the stationarity criterion

$\mathcal{G}_{F}(x):={C_{F}}(x-x^{+})\text{ where
}x^{+}=\arg\min_{z\in\mathcal{X}}\left\|z-\left(x-\frac{1}{C_{F}}\nabla
F(x)\right)\right\|^{2}.$ (A.13)

Initialized to some $x_{0}$ and the inexact gradient oracle
$\widetilde{\nabla}F$, the updates of the inexact projected gradient descent
method is given by

$\begin{split}\text{{For}}&\text{ t=1,2,..., N {do}:}\\\ &\text{Set
}x_{t}\leftarrow\arg\min_{z\in\mathcal{X}}\left\|z-\left(x_{t-1}-\frac{1}{C_{F}}\widetilde{\nabla}F(x_{t-1})\right)\right\|^{2}.\end{split}$
(A.14)

The next proposition calculates the complexity result.

######  Proposition A.6.

Consider the constrained optimization problem in (A.12) with $F$ being
$C_{F}$-Lipschitz smooth and $\mathcal{X}$ having a radius of $R$. When
supplied with a $\delta=\epsilon^{2}/4C_{F}R$ -inexact gradient oracle
$\widetilde{\nabla}F$, that is, $\|\nabla
F(x)-\widetilde{\nabla}F(x)\|\leq\delta$, the solution generated by the
projected gradient descent method (A.14) satisfies

$\min_{t\in[N]}\|\mathcal{G}_{F}(x_{t})\|^{2}\leq\frac{C_{F}(F(x_{0})-F^{*})}{N}+\delta
C_{F}R,$

that is, it takes at most $O(\frac{C_{F}(F(x^{0})-F^{*})}{\epsilon^{2}})$
iterations to generate some $\bar{x}$ with
$\|\mathcal{G}_{F}(x)\|\leq\epsilon$.

###### Proof.

By $C_{F}$-smoothness of $F$, we have

$\displaystyle
f(x_{t+1})=f(x_{t}-\frac{1}{C_{F}}\widetilde{\mathcal{G}_{F}}(x_{t}))$
$\displaystyle\leq
f(x_{t})-\frac{1}{C_{F}}\widetilde{\mathcal{G}_{F}}(x_{t})^{\top}\nabla
f(x_{t})+\frac{1}{2C_{F}}\|\widetilde{\mathcal{G}_{F}}(x_{t})\|^{2}$
$\displaystyle=f(x_{t})-\frac{1}{2C_{F}}\|\widetilde{\mathcal{G}_{F}}(x_{t})(x_{t})\|^{2}+\frac{1}{C_{F}}\widetilde{\mathcal{G}_{F}}(x_{t})^{\top}(\widetilde{\mathcal{G}_{F}}(x_{t})-\nabla
f(x_{t})).$ (A.15)

We now show that
$\frac{1}{\beta}\widetilde{\mathcal{G}_{F}}(x_{t})^{\top}(\widetilde{\mathcal{G}_{F}}(x_{t})-\nabla
f(x_{t}))\leq 0.$ Let
$\widetilde{y}_{t}=x_{t}-\frac{1}{C_{F}}\widetilde{\nabla}F(x_{t})$, and let
$y_{t}=x_{t}-\frac{1}{C_{F}}{\nabla}f(x_{t})$. Then have that

$\displaystyle\frac{1}{C_{F}}\widetilde{\mathcal{G}_{F}}(x_{t})^{\top}(\frac{1}{C_{F}}\widetilde{\mathcal{G}_{F}}(x_{t})-\nabla
f(x_{t}))$
$\displaystyle=C_{F}(x_{t}-\textrm{proj}_{\mathcal{X}}(\widetilde{y}_{t}))^{\top}(y_{t}-\textrm{proj}_{\mathcal{X}}(\widetilde{y}_{t}))$
$\displaystyle=C_{F}(x_{t}-\textrm{proj}_{\mathcal{X}}(\widetilde{y}_{t}))^{\top}(\widetilde{y}_{t}-\textrm{proj}_{\mathcal{X}}(\widetilde{y}_{t}))$
$\displaystyle\quad+C_{F}(x_{t}-\textrm{proj}_{\mathcal{X}}(\widetilde{y}_{t}))^{\top}(y_{t}-\widetilde{y}_{t})$
$\displaystyle\leq
C_{F}(x_{t}-\textrm{proj}_{\mathcal{X}}(\widetilde{y}_{t}))^{\top}(y_{t}-\widetilde{y}_{t})$
$\displaystyle\leq\delta C_{F}R,$

where the penultimate inequality uses the fact that $\mathcal{X}$ is a convex
set, and $R$ is the diameter of the set $X$. Combining this with Eq. A.15, we
have that the function decrease per iteration is

$F(x_{t+1})\leq
F(x_{t})-\frac{1}{2C_{F}}\|\widetilde{\mathcal{G}_{F}}(x_{t})\|^{2}+\delta
C_{F}R.$

Summing over $N$ iterations telescopes the terms, we get

$\min_{t\in[N]}\|\widetilde{\mathcal{G}_{F}}(x_{t})\|^{2}\leq\frac{1}{N}C_{F}(F(x^{0})-F^{*})+\delta
C_{F}R.$

Substituting in $N=\frac{4}{\epsilon^{2}}C_{F}(F(x^{0})-F^{*})$ and the choice
of $\delta=\epsilon^{2}/4C_{F}R$, we get

$\min_{t\in[N]}\|\widetilde{\mathcal{G}_{F}}(x_{t})\|^{2}\leq\frac{\epsilon^{2}}{2}.$

Taking into account the fact that
$\|\widetilde{\mathcal{G}_{F}}(x_{t})-\mathcal{G}_{F}(x_{t})\|\leq\|\nabla
F(x^{t})-\widetilde{\nabla}F(x^{t})\|\leq\delta$, we obtain the desired
result.

∎

###  A.4 The cost of generating approximate solutions of LABEL:{prob:lin-eq}

In this subsection, we address the issue of generating approximations to the
primal and dual solutions $(y^{*},\lambda^{*})$ associated with the lower-
level problem in LABEL:{{prob:lin-eq}}. These approximations are required for
computing the approximate hypergradient in Algorithm 1. For notational
simplicity, we are going to consider the following constrained strongly convex
problem:

$\begin{array}[]{ll}\mbox{minimize}_{y\in\mathbb{R}^{d}}&g(y)\\\ \mbox{subject
to }&By=b.\end{array}$

We propose the following simple scheme to generate approximate solutions to
LABEL:{eq:simple-linear-cosntrained-problem}.

Compute a feasible $\hat{y}$ such that $\|\hat{y}-y^{*}\|\leq\delta$. Then
solve
$\hat{\lambda}=\arg\min_{\lambda\in\mathbb{R}^{m}}\|\nabla_{y}g(\hat{y})-B^{\top}\lambda\|^{2}.$
(A.17)

The following lemma tells us that $\hat{\lambda}$ is close to $\lambda^{*}$ if
$B$ has full row rank.

######  Lemma A.7.

Suppose $g$ in Section A.4 is a $C_{g}$-Lipschitz smooth, and the matrix $B$
has full row rank such that the following matrix $M_{B}$ is invertible

$M_{B}=\begin{bmatrix}I&B^{\top}\\\ B&0\end{bmatrix}.$

Then the approximate solution $(\hat{\lambda},\hat{y})$ from Eq. A.17
satisfies $\|\hat{\lambda}-\lambda^{*}\|\leq\|M_{B}^{-1}\|(1+C_{g})\delta$.

###### Proof.

Since $(\lambda^{*},y^{*})$ satisfy the KKT conditions, they are the solution
to the following linear system

$\underbrace{\begin{bmatrix}I&B^{\top}\\\
B&0\end{bmatrix}}_{=M_{B}}\begin{bmatrix}y^{*}\\\
\lambda^{*}\end{bmatrix}=\begin{bmatrix}-\nabla_{y}g(y^{*})+Iy^{*}\\\
b\end{bmatrix}.$ (A.18)

That is

$\begin{bmatrix}y^{*}\\\
\lambda^{*}\end{bmatrix}=M_{B}^{-1}\begin{bmatrix}-\nabla_{y}g(y^{*})+Iy^{*}\\\
b\end{bmatrix}.$

On the other hand, the approximate solutions $(\hat{y},\hat{\lambda})$ in Eq.
A.17 satisfies

$\begin{bmatrix}I&B^{\top}\\\ B&0\end{bmatrix}\begin{bmatrix}\hat{y}\\\
\hat{\lambda}\end{bmatrix}=\begin{bmatrix}B^{\top}\hat{\lambda}+I\hat{y}\\\
b\end{bmatrix}.$

We show the right hand side (r.h.s) of the above equation to be close to the
r.h.s of Eq. A.18. Let $S:=\\{B^{\top}\lambda:\lambda\in\mathbb{R}^{m}\\}$
denote the subspace spanned by the rows of $B$. We can rewrite
$B^{\top}\hat{\lambda}$ as the projection of $\nabla g(\hat{y})$ onto $S$,
that is,

$\displaystyle B^{\top}\hat{\lambda}=\arg\min_{s\in
S}\|\nabla_{y}g(\hat{y})-s\|^{2}$ $\displaystyle-\nabla_{y}g(y^{*})=$
$\displaystyle B^{\top}\lambda^{*}=\arg\min_{s\in
S}\|\nabla_{y}g(y^{*})-s\|^{2},$

where the second relation follows from the KKT conditon associated with
$(\lambda^{*},y^{*})$. Since the projection is an non-expansive operation, we
have

$\|B^{\top}\hat{\lambda}-(-\nabla_{y}g(y^{*}))\|=\|B^{\top}\hat{\lambda}-B^{\top}\lambda^{*}\|\leq\|\nabla_{y}g(\hat{y})-\nabla
g(y^{*})\|\leq C_{g}\|\hat{y}-y^{*}\|\leq C_{g}\delta.$

We can rewrite $(\hat{y},\hat{\lambda})$ as solutions to the following linear
system with some $\|\tau\|\leq(1+C_{g})\delta$,

$\begin{bmatrix}\hat{y}\\\
\hat{\lambda}\end{bmatrix}=M_{B}^{-1}\begin{bmatrix}-\nabla_{y}g(y^{*})+Iy^{*}+\tau\\\
b\end{bmatrix}.$

Thus we get

$\|\begin{bmatrix}\hat{y}\\\
\hat{\lambda}\end{bmatrix}-\begin{bmatrix}y^{*}\\\
\lambda^{*}\end{bmatrix}\|=\|M_{B}^{-1}\|\|\begin{bmatrix}\tau\\\
0\end{bmatrix}\leq\|M_{B}^{-1}\|(1+C_{g})\delta.$

∎

Now we can just use the AGD method to generate a close enough approximate
solution $\hat{y}$ and call up the Subroutine in Eq. A.17 to generate the
approximate dual solution $\hat{\lambda}.$

Algorithm 7 The Projected Gradient Method to Generate Primal and Dual
Solutions for a Linearly Constrained Problem

1:Input: accuracy requirement $\epsilon>0$ and linearly constrained problem
$\min_{y:By=b}g(y)$.

2:Starting from $y^{0}=0$ and using $Y:=\\{y\in\mathbb{R}^{d}:By=b\\}$ as the
simple feasible region.

3:Run the Accelerated Gradient Descent (AGD) Method (Section 3.3 in [91]) for
$N=\lceil
4\sqrt{C_{g}/\mu_{g}}\log\frac{\|y^{*}\|\|M^{-1}_{B}\|(C_{g}+1)}{\mu_{g}\epsilon}\rceil$
iterations.

4:Use the $y^{N}$ as the approximate solution $\hat{y}$ to generate
$\hat{\lambda}$ according to Eq. A.17.

5:return $(\hat{y},\hat{\lambda})$

######  Proposition A.8.

Suppose the objective function $g$ is both $L_{g}$-smooth and
$\mu_{g}$-strongly convex, and that the constraint satisfies the assumption in
Lemma A.7. Fix an $\epsilon>0$, the solution $(\hat{y},\hat{\lambda})$
returned by the above procedure satisfies $\|y^{*}-\hat{y}\|\leq\epsilon$ and
$\|\hat{\lambda}-\lambda^{*}\|\leq\epsilon$. In another words, the cost of
generating $\epsilon$-close primal and dual solutions are bounded by
$O(\sqrt{\frac{C_{g}}{\mu_{g}}}\log{\frac{1}{\epsilon}}).$

###### Proof.

With $N:=\lceil
4\sqrt{C_{g}/\mu_{g}}\log\frac{\|y^{*}\|\|M^{-1}_{B}\|(L_{g}+1)}{\mu_{g}\epsilon}\rceil$,
Theorem 3.7 in [91] shows that
$\|y^{N}-\hat{y}\|\leq\epsilon/\|M_{B}^{-1}\|(1+L_{g})$. Then we can apply
Lemma A.7 to obtain the desired bound. ∎

##  Appendix B Proofs from LABEL:{sec:nonsmooth}

In this section we present a variant of Algorithm 3, which as stated in
LABEL:{sec:nonsmooth}, requires evaluating only the hyperobjective $F$ (as
opposed to estimating the hypergradient in Algorithm 3). The motivation for
designing such an algorithm, is that while evaluating $\nabla F$ up to
$\alpha$ accuracy requires $O(\alpha^{-1})$ gradient evaluations, the
hyperobjective value can be estimated at a linear rate, as stated in
LABEL:{lem:ZeroOrderApprox}. Our algorithms are based on the Lipschitzness of
$F$, which we prove below. See 4.3

###### Proof.

By Lemma $2.1$ of [16], the hypergradient of $F$ computed with respect to the
variable $x$ may be expressed as
$\nabla_{x}F(x)=\nabla_{x}f(x,y^{\ast}(x))+\left(\frac{dy^{\ast}(x)}{dx}\right)^{\top}\cdot\nabla_{y}f(x,y^{\ast}(x))$.
Since we impose Lipschitzness on $f$ and $y^{*}$, we can bound each of the
terms of $\nabla_{x}F(x)$ by the claimed bound. ∎

###  B.1 Proof of Theorem 4.4

Denoting the uniform randomized smoothing
$F_{\rho}(x):=\mathbb{E}_{\left\lVert z\right\rVert\leq 1}[F(x+\rho\cdot z)]$
where the expectation, here and in what follows, is taken with respect to the
uniform measure, it is well known [92, Lemma 10] that

$\displaystyle\mathbb{E}_{\left\lVert
w\right\rVert=1}\left[\tfrac{d}{2\rho}(F(x+\rho w)-F(x-\rho w))w\right]$
$\displaystyle=\nabla F_{\rho}(x)~{},$ $\displaystyle\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert\nabla F_{\rho}(x)-\tfrac{d}{2\rho}(F(x+\rho
w)-F(x-\rho w))w\right\rVert^{2}$ $\displaystyle\lesssim dL^{2}~{}.$ (B.1)

We first show that replacing the gradient estimator with the inexact
evaluations $\widetilde{F}(\cdot)$ leads to a biased gradient estimator of
$F$.

######  Lemma B.1.

Suppose $|F(\cdot)-\widetilde{F}(\cdot)|\leq\alpha$. Denoting

$\displaystyle g_{x}$ $\displaystyle=\tfrac{d}{2\rho}(F(x+\rho w)-F(x-\rho
w))w~{},$ $\displaystyle\widetilde{{g}}_{x}$
$\displaystyle=\tfrac{d}{2\rho}(\widetilde{F}(x+\rho w)-\widetilde{F}(x-\rho
w))w~{},$

it holds that

$\mathbb{E}_{\left\lVert w\right\rVert=1}\left\lVert
g_{x}-\widetilde{{g}}_{x}\right\rVert\leq\frac{\alpha
d}{\rho}~{},~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert\widetilde{{g}}_{x}\right\rVert^{2}\lesssim\frac{\alpha^{2}d^{2}}{\rho^{2}}+dL^{2}~{}.$

###### Proof.

For the first bound, we have

$\mathbb{E}_{\left\lVert w\right\rVert=1}\left\lVert
g_{x}-\widetilde{{g}}_{x}\right\rVert\leq\frac{d}{2\rho}(2\alpha)\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert w\right\rVert=\frac{\alpha d}{\rho}~{},$

while for the second bound

$\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert\widetilde{{g}}_{x}\right\rVert^{2}=\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert\widetilde{{g}}_{x}-g_{x}+g_{x}\right\rVert^{2}\leq
2\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert\widetilde{{g}}_{x}-g_{x}\right\rVert^{2}+2\mathbb{E}_{\left\lVert
w\right\rVert=1}\left\lVert
g_{x}\right\rVert^{2}\lesssim\frac{d^{2}}{\rho^{2}}\cdot\alpha^{2}+dL^{2}~{},$

where the last step invoked Eq. B.1. ∎

###### Proof of Theorem 4.4.

We denote $\alpha^{\prime}=\frac{\alpha
d}{\rho},~{}\widetilde{G}=\sqrt{\frac{\alpha^{2}d^{2}}{\rho^{2}}+dL^{2}}$.
Since $x_{t}=x_{t-1}+\Delta_{t}$, we have

$\displaystyle F_{\rho}(x_{t})-F_{\rho}(x_{t-1})$
$\displaystyle=\int_{0}^{1}\left\langle\nabla
F_{\rho}(x_{t-1}+s\Delta_{t}),\Delta_{t}\right\rangle ds$
$\displaystyle=\mathbb{E}_{s_{t}\sim\mathrm{Unif}[0,1]}\left[\nabla
F_{\rho}(x_{t-1}+s_{t}{\Delta}_{t}),\Delta_{t}\right]$
$\displaystyle=\mathbb{E}\left[\left\langle\nabla
F_{\rho}(z_{t}),\Delta_{t}\right\rangle\right]~{}.$

By summing over $t\in[T]=[K\times M]$, we get for any fixed sequence
$u_{1},\dots,u_{K}\in\mathbb{R}^{d}:$

$\displaystyle\inf F_{\rho}\leq F_{\rho}(x_{T})$ $\displaystyle\leq
F_{\rho}(x_{0})+\sum_{t=1}^{T}\mathbb{E}\left[\left\langle\nabla
F_{\rho}(z_{t}),\Delta_{t}\right\rangle\right]$
$\displaystyle=F_{\rho}(x_{0})+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F_{\rho}(z_{(k-1)M+m}),\Delta_{(k-1)M+m}-u_{k}\right\rangle\right]$
$\displaystyle~{}~{}~{}~{}+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F_{\rho}(z_{(k-1)M+m}),u_{k}\right\rangle\right]$ $\displaystyle\leq
F_{\rho}(x_{0})+\sum_{k=1}^{K}\mathrm{Reg}_{M}(u_{k})+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F_{\rho}(z_{(k-1)M+m}),u_{k}\right\rangle\right]$ $\displaystyle\leq
F_{\rho}(x_{0})+KD\widetilde{G}\sqrt{M}+K\alpha^{\prime}DM+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F_{\rho}(z_{(k-1)M+m}),u_{k}\right\rangle\right]$

where the last inequality follows by combining Lemma B.1 and Lemma B.2. By
setting $u_{k}:=-D\frac{\sum_{m=1}^{M}\nabla
F_{\rho}(z_{(k-1)M+m})}{\left\lVert\sum_{m=1}^{M}\nabla
F_{\rho}(z_{(k-1)M+m})\right\rVert}$, rearranging and dividing by $DT=DKM$ we
obtain

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left\lVert\frac{1}{M}\sum_{m=1}^{M}\nabla
F_{\rho}(z_{(k-1)M+m})\right\rVert$
$\displaystyle\leq\frac{F_{\rho}(x_{0})-\inf
F_{\rho}}{DT}+\frac{\widetilde{G}}{\sqrt{M}}+\alpha^{\prime}$
$\displaystyle=\frac{F_{\rho}(x_{0})-\inf
F_{\rho}}{K\nu}+\frac{\sqrt{\frac{\alpha^{2}d^{2}}{\rho^{2}}+L^{2}d}}{\sqrt{M}}+\frac{\alpha
d}{\rho}$ $\displaystyle\leq\frac{F_{\rho}(x_{0})-\inf
F_{\rho}}{K\nu}+\frac{{\frac{\alpha
d}{\rho}}}{\sqrt{M}}+\frac{L\sqrt{d}}{\sqrt{M}}+\frac{\alpha d}{\rho}~{}.$
(B.2)

Finally, note that for all $m\in[M]:\left\lVert
z_{(k-1)M+m}-\overline{x}_{k}\right\rVert\leq MD\leq\nu$, therefore $\nabla
F_{\rho}(z_{(k-1)M+m})\in\partial_{\nu}F_{\rho}(\overline{x}_{k})\subset\partial_{\delta}F(\overline{x}_{k})$,
where the last containment is due to [46, Lemma 4] by using our assignment
$\rho+\nu=\delta$. Invoking the convexity of the Goldstein subdifferential,
this implies that

$\frac{1}{M}\sum_{m=1}^{M}\nabla
F_{\rho}(z_{(k-1)M+m})\in\partial_{\delta}F(\overline{x}_{k})~{},$

thus it suffices to bound the first three summands in (B.2) by $\epsilon$ in
order to finish the proof. This happens as long as $\frac{F_{\rho}(x_{0})-\inf
F_{\rho}}{K\nu}\leq\frac{\epsilon}{3}$, $\frac{{\frac{\alpha
d}{\rho}}}{\sqrt{M}}\leq\frac{\epsilon}{3}$, and
$\frac{L\sqrt{d}}{\sqrt{M}}\leq\frac{\epsilon}{3}$, which imply
$K\gtrsim\frac{F_{\rho}(x_{0})-\inf F_{\rho}}{\nu\epsilon}$,
$M\gtrsim\frac{\alpha^{2}d^{2}}{\rho^{2}\epsilon^{2}}$, and
$M\gtrsim\frac{L^{2}d}{\epsilon^{2}}~{}$. By our assignments of $\rho$ and
$\nu$, these result in

$\displaystyle T=KM$ $\displaystyle=O\left(\frac{F_{\rho}(x_{0})-\inf
F_{\rho}}{\nu\epsilon}\cdot\left(\frac{\alpha^{2}d^{2}}{\rho^{2}\epsilon^{2}}+\frac{L^{2}d}{\epsilon^{2}}\right)\right)$
$\displaystyle=O\left(\frac{(F(x_{0})-\inf
F)d}{\delta\epsilon^{3}}\cdot\left(\frac{\alpha^{2}d}{\rho^{2}}+L^{2}\right)\right)$
$\displaystyle=O\left(\frac{(F(x_{0})-\inf
F)d}{\delta\epsilon^{3}}\cdot\left(\alpha^{2}d\cdot\max\left\\{\frac{1}{\delta^{2}},\frac{L^{2}}{(F(x_{0})-\inf
F)^{2}}\right\\}+L^{2}\right)\right)~{},$

completing the proof of convergence guarantee for Algorithm 2.

∎

###  B.2 Proof of Theorem 4.5

We recall Theorem 4.5 below to keep this section self-contained.

See 4.5

Our analysis is inspired by the reduction from online learning to nonconvex
optimization given by [47]. To that end, we start by proving a seemingly
unrelated result, asserting that online gradient descent minimizes the regret
with respect to inexact evaluations. Recalling standard definitions from
online learning, given a sequence of linear losses
$\ell_{m}(\cdot)=\left\langle g_{m},\cdot\right\rangle$, if an algorithm
chooses ${\Delta}_{1},\dots,{\Delta}_{M}$ we denote the regret with respect to
$u$ as

$\mathrm{Reg}_{M}(u):=\sum_{m=1}^{M}\left\langle
g_{m},{\Delta}_{m}-u\right\rangle.$

Consider an update rule according to online projected _inexact_ gradient
descent:

$\Delta_{m+1}:=\mathrm{clip}_{D}(\Delta_{m}-\eta_{m}\widetilde{{g}}_{m}).$

######  Lemma B.2 (Inexact Online Gradient Descent).

In the setting above, suppose that $(\widetilde{{g}}_{m})_{m=1}^{M}$ are
possibly randomized vectors, such that
$\mathbb{E}\left\lVert\widetilde{{g}}_{m}-g_{m}\right\rVert\leq\alpha$ and
$\mathbb{E}\left\lVert\widetilde{{g}}_{m}\right\rVert^{2}\leq\widetilde{G}^{2}$
for all $m\in[M]$. Then for any $\left\lVert u\right\rVert\leq D$ it holds
that

$\mathbb{E}\left[\mathrm{Reg}_{M}(u)\right]\leq\frac{D^{2}}{\eta_{M}}+\widetilde{G}^{2}\sum_{m=1}^{M}\eta_{m}+\alpha
DM~{},$

where the expectation is with respect to the (possible) randomness of
$(\widetilde{{g}}_{m})_{m=1}^{M}$. In particular, setting
$\eta_{m}\equiv\frac{D}{\widetilde{G}\sqrt{M}}$ yields

$\mathbb{E}\left[\mathrm{Reg}_{M}(u)\right]\lesssim
D\widetilde{G}\sqrt{M}+\alpha DM~{}.$

###### Proof.

For any $m\in[M]:$

$\displaystyle\left\lVert\Delta_{m+1}-u\right\rVert^{2}$
$\displaystyle=\left\lVert\mathrm{clip}_{D}(\Delta_{m}-\eta_{m}\widetilde{{g}}_{m})-u\right\rVert^{2}$
$\displaystyle\leq\left\lVert\Delta_{m}-\eta_{m}\widetilde{{g}}_{m}-u\right\rVert^{2}=\left\lVert\Delta_{m}-u\right\rVert^{2}+\eta_{m}^{2}\left\lVert\widetilde{{g}}_{m}\right\rVert^{2}-2\eta_{m}\left\langle\Delta_{m}-u,\widetilde{{g}}_{m}\right\rangle~{},$

thus

$\left\langle\widetilde{{g}}_{m},\Delta_{m}-u\right\rangle\leq\frac{\left\lVert\Delta_{m}-u\right\rVert^{2}-\left\lVert\Delta_{m+1}-u\right\rVert^{2}}{2\eta_{m}}+\frac{\eta_{m}}{2}\left\lVert\widetilde{{g}}_{m}\right\rVert^{2}~{},$

from which we get that

$\displaystyle\mathbb{E}\left\langle g_{m},\Delta_{m}-u\right\rangle$
$\displaystyle=\mathbb{E}\left\langle\widetilde{{g}}_{m},\Delta_{m}-u\right\rangle+\mathbb{E}\left\langle
g_{m}-\widetilde{{g}}_{m},\Delta_{m}-u\right\rangle$
$\displaystyle\leq\frac{\left\lVert\Delta_{m}-u\right\rVert^{2}-\left\lVert\Delta_{m+1}-u\right\rVert^{2}}{2\eta_{m}}+\frac{\eta_{m}}{2}\mathbb{E}\left\lVert\widetilde{{g}}_{m}\right\rVert^{2}+\mathbb{E}\left\lVert
g_{m}-\widetilde{{g}}_{m}\right\rVert\cdot\left\lVert\Delta_{m}-u\right\rVert$
$\displaystyle\leq\frac{\left\lVert\Delta_{m}-u\right\rVert^{2}-\left\lVert\Delta_{m+1}-u\right\rVert^{2}}{2\eta_{m}}+\frac{\eta_{m}}{2}\widetilde{G}^{2}+\alpha
D~{}.$

Summing over $m\in[M]$, we see that

$\displaystyle\mathbb{E}\left[\mathrm{Reg}_{M}(u)\right]$
$\displaystyle\leq\sum_{m=1}^{M}\left\lVert\Delta_{m}-u\right\rVert^{2}\left(\frac{1}{\eta_{m}}-\frac{1}{\eta_{m-1}}\right)+\frac{\widetilde{G}^{2}}{2}\sum_{m=1}^{M}\eta_{m}+M\alpha
D$
$\displaystyle\leq\frac{D^{2}}{\eta_{M}}+\widetilde{G}^{2}\sum_{m=1}^{M}\eta_{m}+\alpha
DM~{}.$

The simplification for $\eta_{m}\equiv\frac{D}{\widetilde{G}\sqrt{M}}$ readily
follows. ∎

We are now ready to analyze Algorithm 3 in the inexact gradient setting.

###### Proof of Theorem 4.5.

Since Algorithm 3 has $x_{t}=x_{t-1}+\Delta_{t}$, we have

$\displaystyle F(x_{t})-F(x_{t-1})$
$\displaystyle=\int_{0}^{1}\left\langle\nabla
F(x_{t-1}+s\Delta_{t}),\Delta_{t}\right\rangle ds$
$\displaystyle=\mathbb{E}_{s_{t}\sim\mathrm{Unif}[0,1]}\left[\langle\nabla
F(x_{t-1}+s_{t}{\Delta}_{t}),\Delta_{t}\rangle\right]$
$\displaystyle=\mathbb{E}\left[\left\langle\nabla
F(z_{t}),\Delta_{t}\right\rangle\right]~{}.$

By summing over $t\in[T]=[K\times M]$, we get for any fixed sequence
$u_{1},\dots,u_{K}\in\mathbb{R}^{d}:$

$\displaystyle\inf F\leq F(x_{T})$ $\displaystyle\leq
F(x_{0})+\sum_{t=1}^{T}\mathbb{E}\left[\left\langle\nabla
F(z_{t}),\Delta_{t}\right\rangle\right]$
$\displaystyle=F(x_{0})+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F(z_{(k-1)M+m}),\Delta_{(k-1)M+m}-u_{k}\right\rangle\right]$
$\displaystyle~{}~{}~{}~{}+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F(z_{(k-1)M+m}),u_{k}\right\rangle\right]$ $\displaystyle\leq
F(x_{0})+\sum_{k=1}^{K}\mathrm{Reg}_{M}(u_{k})+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F(z_{(k-1)M+m}),u_{k}\right\rangle\right]$ $\displaystyle\leq
F(x_{0})+KD\widetilde{G}\sqrt{M}+K\alpha
DM+\sum_{k=1}^{K}\sum_{m=1}^{M}\mathbb{E}\left[\left\langle\nabla
F(z_{(k-1)M+m}),u_{k}\right\rangle\right]$

where the last inequality follows from Lemma B.2 for
$\widetilde{G}=\sqrt{L^{2}+\alpha^{2}},~{}\eta=\frac{D}{\widetilde{G}\sqrt{M}}$,
since $\left\lVert\widetilde{{g}}_{t}-\nabla F(z_{t})\right\rVert\leq\alpha$
(deterministically) for all $t\in[T]$ by assumption. Letting
$u_{k}:=-D\frac{\sum_{m=1}^{M}\nabla
F(z_{(k-1)M+m})}{\left\lVert\sum_{m=1}^{M}\nabla
F(z_{(k-1)M+m})\right\rVert}$, rearranging and dividing by $DT=DKM$, we obtain

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}\left\lVert\frac{1}{M}\sum_{m=1}^{M}\nabla
F(z_{(k-1)M+m})\right\rVert$ $\displaystyle\leq\frac{F(x_{0})-\inf
F}{DT}+\frac{\widetilde{G}}{\sqrt{M}}+\alpha$
$\displaystyle=\frac{F(x_{0})-\inf
F}{K\delta}+\frac{\widetilde{G}}{\sqrt{M}}+\alpha~{}.$ (B.3)

Finally, note that for all $k\in[K],m\in[M]:\left\lVert
z_{(k-1)M+m}-\overline{x}_{k}\right\rVert\leq MD\leq\delta$, therefore $\nabla
F(z_{(k-1)M+m})\in\partial_{\delta}F(\overline{x}_{k})$. Invoking the
convexity of the Goldstein subdifferential, we see that

$\frac{1}{M}\sum_{m=1}^{M}\nabla
F(z_{(k-1)M+m})\in\partial_{\delta}F(\overline{x}_{k})~{},$

thus it suffices to bound the first two summands on the right-hand side in Eq.
B.3 by $\epsilon$ in order to finish the proof. This happens as long as
$\frac{F(x_{0})-\inf F}{K\delta}\leq\frac{\epsilon}{2}$ and
$\frac{\widetilde{G}}{\sqrt{M}}\leq\frac{\epsilon}{2}$. These are equivalent
to $K\geq\frac{2(F(x_{0})-\inf F)}{\delta\epsilon}$ and
$M\geq\frac{4\widetilde{G}^{2}}{\epsilon^{2}}$, which results in

$T=KM=O\left(\frac{F(x_{0})-\inf
F}{\delta\epsilon}\cdot\frac{L^{2}+\alpha^{2}}{\epsilon^{2}}\right)=O\left(\frac{(F(x_{0})-\inf
F)L^{2}}{\delta\epsilon^{3}}\right),$

completing the proof. ∎

###  B.3 Proof of Theorem 4.6

Throughout the proof we denote $z_{t}=x_{t}+\delta\cdot w_{t}$. Since $F$ is
$L$-Lipschitz,
$F_{\delta}(x):=\mathbb{E}_{w\sim\mathrm{Unif}(\mathbb{S}^{d-1})}[F(x+\delta\cdot
w)]$ is $L$-Lipschitz and $O(L\sqrt{d}/\delta)$-smooth. By smoothness we get

$\displaystyle F_{\delta}(x_{t+1})-F_{\delta}(x_{t})$
$\displaystyle\leq\left\langle\nabla
F_{\delta}(x_{t}),x_{t+1}-x_{t}\right\rangle+O\left(\frac{L\sqrt{d}}{\delta}\right)\cdot\left\lVert
x_{t+1}-x_{t}\right\rVert^{2}$ $\displaystyle=-\eta\left\langle\nabla
F_{\delta}(x_{t}),\widetilde{{g}}_{t}\right\rangle+O\left(\frac{\eta^{2}L\sqrt{d}}{\delta}\right)\cdot\left\lVert\widetilde{{g}}_{t}\right\rVert^{2}$
$\displaystyle=-\eta\left\langle\nabla F_{\delta}(x_{t}),\nabla
F(z_{t})\right\rangle-\eta\left\langle\nabla
F_{\delta}(x_{t}),\widetilde{{g}}_{t}-\nabla
F(z_{t})\right\rangle+O\left(\frac{\eta^{2}L\sqrt{d}}{\delta}\right)\cdot\left\lVert\widetilde{{g}}_{t}\right\rVert^{2}~{}.$

Noting that $\mathbb{E}[\nabla F(z_{t})]=\nabla F_{\delta}(x_{t})$ and that
$\left\lVert\widetilde{{g}}_{t}\right\rVert\leq\left\lVert\widetilde{{g}}_{t}-\nabla
F(z_{t})\right\rVert+\left\lVert\nabla F(z_{t})\right\rVert\leq\alpha+L$, we
see that

$\displaystyle\mathbb{E}[F_{\delta}(x_{t+1})-F_{\delta}(x_{t})]\leq-\eta\mathbb{E}\left\lVert\nabla
F_{\delta}(x_{t})\right\rVert^{2}+\eta
L\alpha+O\left(\frac{\eta^{2}L\sqrt{d}}{\delta}(\alpha+L)^{2}\right)~{},$

which implies

$\mathbb{E}\left\lVert\nabla
F_{\delta}(x_{t})\right\rVert^{2}\leq\frac{\mathbb{E}[F_{\delta}(x_{t})]-\mathbb{E}[F_{\delta}(x_{t+1})]}{\eta}+L\alpha+O\left(\frac{\eta
L\sqrt{d}(\alpha+L)^{2}}{\delta}\right)~{}.$

Averaging over $t=0,\dots,T-1$ and noting that $F_{\delta}(x_{0})-\inf
F_{\delta}\leq(F(x_{0})-\inf F)+\delta L$ results in

$\mathbb{E}\left\lVert\nabla
F_{\delta}(x^{\mathrm{out}})\right\rVert^{2}=\frac{1}{T}\sum_{t=0}^{T-1}E\left\lVert\nabla
F_{\delta}(x_{t})\right\rVert^{2}\leq\frac{(F(x_{0})-\inf F)+\delta L}{\eta
T}+L\alpha+O\left(\frac{\eta L\sqrt{d}(\alpha+L)^{2}}{\delta}\right)~{}.$

By Jensen’s inequality and the sub-additivity of the square root,

$\mathbb{E}\left\lVert\nabla
F_{\delta}(x^{\mathrm{out}})\right\rVert\leq\sqrt{\frac{(F(x_{0})-\inf
F)+\delta L}{\eta T}}+\sqrt{L\alpha}+O\left(\sqrt{\frac{\eta
L\sqrt{d}(\alpha+L)^{2}}{\delta}}\right)~{}.$

Setting $\eta=\frac{\sqrt{((F(x_{0})-\inf F)+\delta
L)\delta}}{\sqrt{TL\sqrt{d}(\alpha+L)^{2}}}$ yields the final bound

$\mathbb{E}\left\lVert\nabla
F_{\delta}(x^{\mathrm{out}})\right\rVert\lesssim\frac{((F(x_{0})-\inf
F)+\delta
L)^{1/4}L^{1/4}d^{1/8}(\alpha+L)^{1/2}}{\delta^{1/4}T^{1/4}}+\sqrt{L\alpha}~{},$

and the first summand is bounded by $\epsilon$ for
$T=O\left(\frac{((F(x_{0})-\inf F)+\delta
L)L\sqrt{d}(L+\alpha)^{2}}{\delta\epsilon^{4}}\right)$.

##  Appendix C Proofs of Section 5

###  C.1 Reformulation equivalence

######  Theorem C.1 (Reformulation equivalence).

When $\lambda^{*}$ matches to an optimal dual solution to the lower level
problem $y^{*}=\arg\min_{y}g(x,y)~{}\text{s.t.}~{}h(x,y)\leq 0$, we show that
for each $x$, the reformulation has the same feasible region of $y$.

###### Proof.

We first show that lower-level feasibility implies feasibility of the
reformulated problem. Let
$y^{*},\lambda^{*}=\min\limits_{y}\max\limits_{\beta\geq
0}g(x,y)+\beta^{\top}h(x,y)$ be the primal and the dual solution to the lower
level problem with parameter $x$. We can verify that $y^{*}$ satisfies all the
constraints in the reformulation problem. The feasibility condition
$h(x,y^{*})$ is automatically satisfied. We just need to check:

$\displaystyle g^{*}(x)$
$\displaystyle\coloneqq\min\limits_{\theta}g(x,\theta)+(\lambda^{*})^{\top}h(x,\theta)$
$\displaystyle=g(x,y^{*})+(\lambda^{*})^{\top}h(x,y^{*}).$ (C.1)

Therefore, $x,y^{*}$ is a feasible point to the reformulation problem.

We now show the other direction, i.e., that feasibility of the reformulaed
problem implies that of the lower-level problem. Given $\lambda^{*}$, let us
assume $y$ satisfies $g(x,y)\leq g^{*}_{\lambda^{*}}(x)$ and $h(x,y)\leq 0$.
On the other hand, assume
$y^{*},\lambda^{*}=\min\limits_{y}\max\limits_{\beta\geq
0}g(x,y)+\beta^{\top}h(x,y)$ be the primal and the dual solution. We can show
that:

$\displaystyle g(x,y)+(\lambda^{*})^{\top}h(x,y)\leq
g^{*}(x)\coloneqq\min_{\theta}g(x,\theta)+(\lambda^{*})^{\top}h(x,\theta).$
(C.2)

By the strong convexity of $g+(\lambda^{*})^{\top}h$, we know that $y$ matches
to the unique minimum $y^{*}$, which implies that $y=y^{*}$ is also a feasible
point to the original bilevel problem. ∎

###  C.2 Active constraints in differentiable optimization

By computing the derivative of the KKT conditions in Section 2.1, we get:

$\displaystyle(\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h)+(\nabla^{2}_{yy}g+(\lambda^{*})^{\top}\nabla_{yy}^{2}h)\frac{dy^{*}}{dx}+(\nabla_{y}h)^{\top}\frac{d\lambda^{*}}{dx}$
$\displaystyle=0$ (C.3)
$\displaystyle\text{diag}(\lambda^{*})\nabla_{x}h+\text{diag}(\lambda^{*})\nabla_{y}h\frac{dy^{*}}{dx}+\text{diag}(h)\frac{d\lambda^{*}}{dx}$
$\displaystyle=0.$ (C.4)

Let $\mathcal{I}=\\{i\in[d_{h}]|h(x,y^{*})_{i}=0,\lambda^{*}_{i}>0\\}$ be the
set of active constraints with positive dual solution, and
$\mathcal{I}_{1}=\\{i|h(x,y^{*})_{i}\neq 0\\}$ be the set of inactive
constraints and $\mathcal{I}_{2}=\\{i|h(x,y^{*})_{i}=0,\lambda^{*}_{i}=0\\}$.
We know that $\bar{\mathcal{I}}=\mathcal{I}_{1}\cup\mathcal{I}_{2}$. For each
$i\in\mathcal{I}_{1}$, due to complementary slackness, we know that
$\lambda^{*}_{i}=0$.

For $i\in\mathcal{I}_{1}$ in Eq. C.3, we have
$\lambda^{*}_{i}\nabla_{x}h(x,y^{*})_{i}+\lambda^{*}_{i}\nabla_{y}h(x,y^{*})_{i}\frac{dy^{*}}{dx}+h(x,y^{*})_{i}\frac{d\lambda^{*}_{i}}{dx}=0$,
which implies $h(x,y^{*})_{i}\frac{d\lambda^{*}_{i}}{dx}=0$ because
$\lambda^{*}_{i}=0$. This in turn implies $\frac{d\lambda^{*}_{i}}{dx}=0$
because $h(x,y^{*})_{i}<0$. That means the dual variable $\lambda^{*}_{i}=0$
and has zero gradient $\frac{d\lambda^{*}_{i}}{dx}=0$ for any index
$i\in\mathcal{I}_{1}$. Therefore, we can remove row $i\in\mathcal{I}_{1}$ in
Eq. C.4 and obtain $\lambda^{*}_{i}=0$ and $\frac{d\lambda^{*}_{i}}{dx}=0$.

For $i\in\mathcal{I}_{2}$, the KKT condition in Eq. C.4 is degenerate.
Therefore, $\frac{d\lambda^{*}_{i}}{dx}$ can be arbitrary, i.e., non-
differentiable. As a subgradient choice, we can set
$\frac{d\lambda^{*}_{i}}{dx}=0$ for such $i$. This choice will also eliminate
its impact on the KKT condition in Eq. C.3 because
$\frac{d\lambda^{*}_{i}}{dx}$ is set to be $0$. By this choice of subgradient,
we can also remove row $i\in\mathcal{I}_{2}$ Eq. C.4.

Thus Eq. C.4 can be written as the following set of equations, for
$h_{\mathcal{I}}=[h_{i}]_{i\in\mathcal{I}}$ and
$\lambda^{*}_{\mathcal{I}}=[\lambda^{*}_{i}]_{i\in\mathcal{I}}$:

$\displaystyle\text{diag}(\lambda^{*})\nabla_{x}h_{\mathcal{I}}+\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{y}h_{\mathcal{I}}\frac{dy^{*}}{dx}+\text{diag}(h_{\mathcal{I}})\frac{d\lambda^{*}_{\mathcal{I}}}{dx}=0$
$\displaystyle\Longrightarrow\quad$
$\displaystyle\text{diag}(\lambda^{*})\nabla_{x}h_{\mathcal{I}}+\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{y}h_{\mathcal{I}}\frac{dy^{*}}{dx}=0\quad\text{(due
to $h_{\mathcal{I}}(x,y^{*})=0$)}.$ (C.5)

In Eq. C.3, due to $\frac{d\lambda^{*}_{i}}{dx}=0$ for all
$i\in\bar{\mathcal{I}}$, we can remove $\frac{d\lambda^{*}_{i}}{dx}~{}\forall
i\in\bar{\mathcal{I}}$ in Eq. C.3 by:

$\displaystyle 0=$
$\displaystyle~{}(\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h)+(\nabla^{2}_{yy}g+(\lambda^{*})^{\top}\nabla_{yy}^{2}h)\frac{dy^{*}}{dx}+(\nabla_{y}h)^{\top}\frac{d\lambda^{*}}{dx}$
$\displaystyle=$
$\displaystyle~{}(\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h)+(\nabla^{2}_{yy}g+(\lambda^{*})^{\top}\nabla_{yy}^{2}h)\frac{dy^{*}}{dx}+(\nabla_{y}h_{\mathcal{I}})^{\top}\frac{d\lambda^{*}_{\mathcal{I}}}{dx}.$
(C.6)

Combining Eq. C.6 and Eq. C.5, we get:

$\displaystyle(\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h)+(\nabla^{2}_{yy}g+(\lambda^{*})^{\top}\nabla_{yy}^{2}h)\frac{dy^{*}}{dx}+(\nabla_{y}h_{\mathcal{I}})^{\top}\frac{d\lambda^{*}_{\mathcal{I}}}{dx}$
$\displaystyle=0$
$\displaystyle\text{diag}(\lambda^{*})\nabla_{x}h_{\mathcal{I}}+\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{y}h_{\mathcal{I}}\frac{dy^{*}}{dx}$
$\displaystyle=0,$

which can be written in its matrix form:

$\displaystyle\begin{bmatrix}\nabla^{2}_{yy}g+(\lambda^{*})^{\top}\nabla_{yy}^{2}h&\nabla_{y}h_{\mathcal{I}}^{\top}\\\
\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{y}h_{\mathcal{I}}&0\end{bmatrix}\begin{bmatrix}\frac{dy^{*}}{dx}\\\
\frac{d\lambda^{*}_{\mathcal{I}}}{dx}\end{bmatrix}=-\begin{bmatrix}\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h\\\
\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{x}h_{\mathcal{I}}\end{bmatrix}$
(C.7)

This concludes the derivation of the derivative of constrained optimization in
Eq. 5.2.

###  C.3 Inequality case: bounds on primal solution error and constraint
violation

See 5.1

###### Proof.

We first provide the claimed bound on
$\|y^{*}_{\alpha_{1},\alpha_{2}}-y^{*}(x)\|$.

Part 1: Bound on the convergence of $y$.

Since $y_{\lambda^{*},\bm{\alpha}}^{*}$ minimizes
$\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y)$, the first-order condition gives
us:

$\displaystyle
0=\nabla_{y}\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y_{\lambda^{*},\bm{\alpha}}^{*}).$

Similarly, we can compute the gradient of
$\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y)$ at $y^{*}$:

$\displaystyle\nabla_{y}\mathcal{L}_{\alpha}(x,y^{*})$
$\displaystyle=\nabla_{y}f(x,y^{*})+\alpha_{1}(\nabla_{y}g(x,y^{*})+(\lambda^{*})^{\top}\nabla_{y}h(x,y^{*}))+\alpha_{2}\nabla_{y}h_{\mathcal{I}}(x,y^{*})^{\top}h_{\mathcal{I}}(x,y^{*})$
$\displaystyle=\nabla_{y}f(x,y^{*}),$

where the second step is due to the property of the primal and dual solution:
$\nabla_{y}g(x,y^{*})+(\lambda^{*})^{\top}\nabla_{y}h(x,y^{*})=0$ by the
stationarity condition in the KKT conditions, and by definition of the active
constraints $h_{\mathcal{I}}$ where the optimal $y^{*}$ must have
$h_{\mathcal{I}}(x,y^{*})=0$.

Since, for a sufficiently large $\alpha_{1}$, the penalty function is
$\alpha_{1}\mu_{g}-L_{f}\geq\frac{\alpha_{1}\mu_{g}}{2}$ strongly convex in
$y$, we have:

$\displaystyle\frac{\alpha_{1}\mu_{g}}{2}\left\lVert
y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert\leq\left\lVert\nabla_{y}\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y^{*})-\nabla_{y}\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert=\left\lVert\nabla_{y}f(x,y^{*})\right\rVert\leq
L_{f}.$

Therefore, upon rearranging the terms, we obtain the claimed bound:

$\displaystyle\left\lVert
y^{*}-y^{*}_{\bm{\alpha},\lambda^{*}}\right\rVert\leq\frac{2L_{f}}{\alpha_{1}\mu_{g}}.$

Part 2: bound on the constraint violation.

When we plug $y^{*}$ into Eq. 5.4, we get:

$\displaystyle\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y^{*})$
$\displaystyle=f(x,y^{*})+\alpha_{1}(g(x,y^{*})+(\lambda^{*})^{\top}h(x,y^{*})-g^{*}_{\lambda^{*}}(x))+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y^{*})\right\rVert^{2}=f(x,y^{*}).$

Plugging in $y^{*}_{\bm{\alpha},\lambda^{*}}$, we may obtain:

$\displaystyle\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})$
$\displaystyle=f(x,y_{\lambda^{*},\bm{\alpha}}^{*})+\alpha_{1}(g(x,y_{\lambda^{*},\bm{\alpha}}^{*})+(\lambda^{*})^{\top}h(x,y_{\lambda^{*},\bm{\alpha}}^{*})-g^{*}(x))+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert^{2}$
$\displaystyle=f(x,y_{\lambda^{*},\bm{\alpha}}^{*})+\alpha_{1}(g(x,y_{\lambda^{*},\bm{\alpha}}^{*})+(\lambda^{*})^{\top}h(x,y_{\lambda^{*},\bm{\alpha}}^{*})-g(x,y^{*})-(\lambda^{*})^{\top}h(x,y^{*}))$
$\displaystyle\qquad+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert^{2}$
$\displaystyle\geq
f(x,y_{\lambda^{*},\bm{\alpha}}^{*})+\alpha_{1}\frac{\mu_{g}}{2}\left\lVert
y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert^{2}+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert^{2},$

where we used the strong convexity (with respect to $y$) of
$g(x,y)+(\lambda^{*})^{\top}h(x,y)$ and the optimality of $y^{*}$ for
$g(x,y)+(\lambda^{*})^{\top}h(x,y)$. By the optimality of
$y_{\lambda^{*},\bm{\alpha}}^{*}$ for $\mathcal{L}_{\bm{\alpha},\lambda^{*}}$,
we know that

$\displaystyle
f(x,y^{*})=\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y^{*})\geq\mathcal{L}_{\bm{\alpha},\lambda^{*}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\geq
f(x,y_{\lambda^{*},\bm{\alpha}}^{*})+\alpha_{1}\frac{\mu_{g}}{2}\left\lVert
y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert^{2}+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert^{2}.$

Therefore, by the Lipschitzness of the function $f$ in terms of $y$, and the
bound
$\|y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\|\leq\frac{2L_{f}}{\alpha_{1}\mu_{g}}$,
we know that:

$\displaystyle\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\right\rVert^{2}$
$\displaystyle\leq
f(x,y^{*})-f(x,y_{\lambda^{*},\bm{\alpha}}^{*})-\alpha_{1}\frac{\mu_{g}}{2}\left\lVert
y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert^{2}$ $\displaystyle\leq
L_{f}\left\lVert
y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert-\alpha_{1}\frac{\mu_{g}}{2}\left\lVert
y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert^{2}$ $\displaystyle\leq
L_{f}\left\lVert y^{*}-y_{\lambda^{*},\bm{\alpha}}^{*}\right\rVert$
$\displaystyle=O({\alpha_{1}^{-1}}).$

Rearranging terms then gives the claimed bound. ∎

The bound on the constraint violation in LABEL:{thm:solution-bound} is an
important step in the following theorem.

###  C.4 Proof of Lemma 5.2: gradient approximation for inequality constraints

See 5.2

###### Proof.

First, we recall Eq. 5.4 here:

$\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y)=f(x,y)+\alpha_{1}\left(g(x,y)+(\lambda^{*})^{\top}h(x,y)-g^{*}(x)\right)+\frac{\alpha_{2}}{2}\left\lVert
h_{\mathcal{I}}(x,y)\right\rVert^{2}.$

Next, recall from Eq. C.1, we can express
$g^{*}(x)=g(x,y^{*})+(\lambda^{*})^{\top}h(x,y^{*})$, which we use in the
first step below:

$\displaystyle\nabla_{x}F(x)-\frac{d}{dx}\mathcal{L}_{\lambda^{*},\bm{\alpha}}(x,y_{\lambda^{*},\alpha}^{*})$
$\displaystyle=$
$\displaystyle\left(\nabla_{x}f(x,y^{*})+\frac{dy^{*}}{dx}^{\top}\nabla_{y}f(x,y^{*})\right)-\Biggl{(}\nabla_{x}f(x,y_{\lambda^{*},\bm{\alpha}}^{*})+\alpha_{1}(\nabla_{x}g(x,y_{\lambda^{*},\bm{\alpha}}^{*})+\nabla_{x}h(x,y_{\lambda^{*},\bm{\alpha}}^{*})^{\top}\lambda^{*}$
$\displaystyle-\alpha_{1}(\nabla_{x}g(x,y^{*})+\nabla_{x}h(x,y^{*})^{\top}\lambda^{*})+\alpha_{2}\nabla_{x}h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})^{\top}h_{\mathcal{I}}(x,y_{\lambda^{*},\bm{\alpha}}^{*})\Biggl{)}$
$\displaystyle=$
$\displaystyle\nabla_{x}f(x,y^{*})-\nabla_{x}f(x,y_{\lambda^{*},\bm{\alpha}}^{*})$
(C.8)
$\displaystyle+\frac{dy^{*}}{dx}^{\top}\nabla_{y}f(x,y^{*})-\frac{dy^{*}}{dx}^{\top}\nabla_{y}f(x,y_{\lambda^{*},\bm{\alpha}}^{*})$
(C.9)
$\displaystyle+\frac{dy^{*}}{dx}^{\top}\nabla_{y}f(x,y_{\lambda^{*},\bm{\alpha}}^{*})-\underbrace{\alpha_{1}\begin{bmatrix}\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h\\\
\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{x}h_{\mathcal{I}}\end{bmatrix}^{\top}\begin{bmatrix}y_{\lambda^{*},\bm{\alpha}}^{*}-y^{*}\\\
0\end{bmatrix}}_{\text{added term 1}}$
$\displaystyle\qquad-\underbrace{\alpha_{2}\begin{bmatrix}\nabla^{2}_{yx}g+(\lambda^{*})^{\top}\nabla_{yx}^{2}h\lambda^{*}\\\
\text{diag}(\lambda^{*}_{\mathcal{I}})\nabla_{x}h_{\mathcal{I}}\end{bmatrix}^{\top}\begin{bmatrix}0\\\