<figure><img src="image_1.png" /></figure>

HAL
open science

# Bilevel gradient methods and the Morse parametric qualification condition

Jérôme Bolte, Quoc-Tung Le, Edouard Pauwels, Samuel Vaiter

► To cite this version:

Jérôme Bolte, Quoc-Tung Le, Edouard Pauwels, Samuel Vaiter. Bilevel gradient methods and the Morse parametric qualification condition. 2025. hal-04942322v2

HAL Id: hal-04942322

https://inria.hal.science/hal-04942322v2

Preprint submitted on 17 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bilevel gradient methods and the Morse parametric qualification condition

Jérôme Bolte* Tùng Lê * Edouard Pauwels* Samuel Vaiter †

## Abstract

We introduce the *Morse parametric qualification condition* for bilevel programming. Generic semi-algebraic functions are Morse parametric in a piecewise sense. Thus, bilevel programs with a Morse parametric lower level constitute a relevant intermediate class between strongly convex and fully generic lower levels. In this framework, we study bilevel gradient algorithms with two strategies: the *single-step multi-step strategy*, which involves a sequence of steps on the lower-level problems followed by one step on the upper-level problem, and a *differentiable programming strategy* that optimizes a smooth approximation of the bilevel problem. While the first is shown to be a biased gradient method on the problem with rich properties, the second, inspired by meta-learning applications, is less stable but offers simplicity and ease of implementation.

## 1 Introduction

Bilevel optimization provides a versatile formalism which encompasses a diversity of application settings [19, 21]. It has recently attracted significant attention in machine learning as it allows to formalize a broad spectrum of problem situations, including hyper-parameter tuning [52], meta-learning [22, 23], data augmentation [55], deep equilibrium networks [7] or neural architecture search [41]. The present work investigates solution methods for bilevel optimization. We focus on bilevel gradient methods: first-order algorithms combining gradient steps for both levels of the problem.

We consider bilevel optimization problems of the following form

$$ \min_{x \in \mathbb{R}^n, y \in \mathbb{R}^m} f(x, y) \quad \text{s.t.} \quad y \in \operatorname*{argmin}_{x \in \mathbb{R}^n} g(x, \cdot). \qquad (\text{BL}) $$

where $f, g: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ are respectively the *upper level objective* and the *lower level objective*. Throughout this paper, we assume that the lower level minimum is achieved.

### 1.1 The Morse qualification condition and two algorithms

**Parametric Morse qualification condition:** A delicate aspect of bilevel programming is the gap between qualification regimes: it jumps from very strong conditions, when the lower level is strongly convex, to a great generality which must be handled by rather complex KKT machineries, see e.g., [1, 19, 21]. In the first case, applications are restrictive, in particular for machine learning, whereas in the general case one needs to resort to complex qualification conditions to rule out pathologies even in the polynomial case [10]<fnref n="1" />. To partly remedy this difficulty, we introduce and study parametric Morse qualification conditions (see Assumption Morse QC)

<footnote>*Toulouse School of Economics, University of Toulouse Capitole, Toulouse, France.</footnote>

<footnote>†CNRS & Université Côte d'Azur, Laboratoire J. A. Dieudonné. Nice, France.</footnote>

<footnote n="1">¹The value function is typically discontinuous; any closed set can be expressed with smooth bilevel constraints; and polynomial bilevel optimization is $\Sigma_2^p$-hard, above NP in the polynomial hierarchy. For example, there is no numerically efficient certificate of stationarity or approximate stationarity for general bilevel problems.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Structures of crit₁g (see the definition in §2.1) under Assumption Morse QC; see also Proposition 3.6. The decomposition of the sets of critical points and local minima of g as a finite union of manifolds is represented by several branches y⁽¹⁾(x), y⁽²⁾(x), y⁽³⁾(x). Saddle-point manifolds (blue) and local-minima manifolds (red) illustrate the stratification of the critical set.</figcaption></figure>

for the lower-level function $g$. An intuitive picture is that as the parameters $x$ move, the lower-level landscape stays essentially the same, with an invariant “Morse profile”: the number and type of critical points remain constant, and each one traces a smooth branch as the parameter varies.

A first obvious interest is that this condition provides a meaningful intermediate between strongly convex and general nonconvex lower-level problems. There are moreover two other features that makes the setting both appealing and workable: (i) generically, a semi-algebraic lower-level objective $g$ satisfies a piecewise parametric Morse property (Proposition 3.4); (ii) under this qualification condition, the lower-level critical and local-minima sets split into a finite union of $C^2$ manifolds (Proposition 3.6), as depicted in Figure 1.

A consequence of the latter is that we obtain the following mixed-integer nonlinear programming relaxation for the bilevel problem (see also Table 1):

$$ \begin{array}{ccc} \min_{y \in \argmin} f(x, y) & \xrightarrow{\text{relaxation}} & \min_{y \in \argmin\text{-loc}} f(x, y) \\ y \in \argmin\text{-loc} g(x, \cdot) & & i \in \{1, \dots, N\} \end{array} \xrightarrow{\text{equivalent}} \min_{y \in \argmin\text{-loc}} f(x, y^{(i)}(x)) $$

where $y^{(i)}$ are $C^2$ functions, $i = 1, \dots, N$, highlighting the disconnected nature of the constraints in bilevel programming. Let us draw the attention of the reader on the fact that the relation $y = y^{(i)}(x)$ is implicitly defined through a nonconvex inner minimization, which must be approximated algorithmically, and the integer $N$ is typically unknown. Throughout this paper, we approximate lower-level solutions using a simple gradient-descent scheme (Algorithm LLD). Our main results (cf. Theorems 4.2, 5.3 and 5.6) relate to this local minima relaxation (see Table 1).

**Convergence of a single-step multi-step bilevel gradient method:** Our first algorithmic result analyzes a simple alternating scheme: one outer step followed by multiple inner steps that approximately minimize the lower-level problem (see Algorithm SMBG). The scheme can be recast as an *inexact* method on the upper-level value function. Such strategies are common in bilevel gradient methods, with notable extensions [4, 24, 13, 17, 18, 33, 48].

Our analysis departs from previous work in two ways: (i) we allow general nonconvex lower-level problems with a possibly discontinuous argmin set-valued map; (ii) our assumptions are explicit and verifiable a priori. In Theorem 4.2, under the parametric Morse qualification condition and natural regularity and stepsize conditions, we prove that the method provides approximate solutions of the bilevel problem. This strengthens [4] in several directions and leverages the theory of *inexact gradient methods* for semi-algebraic optimization [11].
**Pseudo stability of a differentiable programming bilevel gradient method:** Motivated by the analysis of Model-Agnostic Meta-Learning (MAML) models [22, 54], we consider a fully differentiable programming approximation of the lower level constraints in Algorithm DPBG, where the initialization of the lower level problem is treated as an upper level parameter. Such approaches were promoted in several works [23, 35, 32, 34, 43]. Although conceptually simple, the resulting minimization problems exhibit a fuzzy connection with the original bilevel problem. First, the differentiable programming approximation essentially erases the bilevel constraint as its stationary points are the same as the unconstrained upper level, see Proposition 5.2. This negative result is balanced by a stability analysis that may explain the method's partial empirical success. Indeed, well-qualified bilevel solutions are pseudo-stable in the sense that, if the algorithm is to leave them, it generally needs many iterations to exit their neighborhood (Theorem 5.3). Moreover, unconstrained local minimizers, which are irrelevant to the bilevel program, are, in practice, hard to reach (Theorem 5.6). Taken together, these two facts are rather positive and point in the same direction. Yet the picture is imperfect and instabilities can still arise, consistent with meta-learning evidence [3]. This suggests either adjusting differentiable-programming pipelines with appropriate safeguards or modifications or turning to classical gradient methods (as in the previous approach).

## 1.2 Related work

Gradient methods and steepest descent methods are natural candidate solutions for bilevel optimization. Early works on the topic include differentiation of second-order local minima [37], bilevel optimality conditions, and the design of steepest descent methods [56]. These studies focus on unique, well-posed lower level solutions, a simplifying hypothesis which allows to reduce the bilevel problem to that of the minimization of a locally Lipschitz objective for which general purpose methods could be considered [51]. Overviews of bilevel programming contributions, prior to the latter machine learning developments, can be found in [61, 14].

The motivations for machine learning applications [52, 23, 45] have triggered the investigation of dedicated solution methods, the vast majority of them being gradient techniques [4, 24, 13, 17, 18, 33], most of them with a stochastic flavour. These more recent works have specificities compared to the broader bilevel optimization literature. First, the proposed theoretical analyses take explicitly into account the fact that the lower level is solved by an algorithm, hence producing only an approximate solution to the lower level problem. Second, the fact that the lower level solution algorithm is explicitly considered allows us to leverage all variations of algorithmic differentiation [28] to obtain not only an approximate solution to the lower level, but also its derivative with respect to the upper level parameters, with limited overheads. This concept is known as the cheap gradient principle and the overall algorithmic differentiation framework sometimes called *differentiable programming*. This represents a conceptual advantage and has a lot of practical benefits. The bilevel gradient algorithms which we study follow these principles: the lower level solution algorithm is explicit and we consider that it is possible to propagate derivatives through the algorithm. Another specificity of this literature is the possibility of considering the initialization of a lower level gradient scheme as an upper level parameter to be optimized. This is common in variants of Model-Agnostic Meta-Learning (MAML) models [22, 54] and was explicitly considered by several authors [23, 35, 32, 34, 43]. We will see in the study of our differentiable programming bilevel gradient method that this has benefits, but also drawbacks.

In terms of theoretical guarantees, a majority of existing studies make simplifying assumptions ensuring the uniqueness and smoothness of the lower level solution, such as strong convexity. The bilevel problem then reduces to a smooth or at least Lipschitz composite optimization problem [51, 62, 49]. Beyond strong convexity, there is no complete general treatment of non-convex lower level problems, with possibly multiple lower level solutions. A specificity of our analysis is the generality of the lower-level problem, as we do not assume neither convex-
ity nor the existence of a unique minimizer. This lack of structure significantly makes algorithm design and theoretical guarantees harder. Some approaches introduce auxiliary initializations to guide optimization [43], while others rely on interior-point methods to mitigate computational costs [42]. Certain works shift their focus on finding approximated solutions and combine the single-level penalized functions uniform growth conditions on the lower level problem to establish convergence properties [57]. The ambiguity arising from multiple lower-level solutions can be addressed using selection maps (to allow to prove an implicit function theorem even in degenerate cases) [4]. Reformulating bilevel problems as single-level ones introduces nonsmooth constraints requiring specialized algorithms [40]. From an analytical perspective, variational analysis and generalized differentiation provide optimality conditions when classical constraint qualifications fail [49], often requiring calmness conditions [62] and sensitivity analysis to properly characterize solutions [20]. The assumption of a unique lower-level solution is typically unrealistic and some work [44] propose to reformulate bilevel optimization as an alternative kind of optimization. Finally, bilevel optimization is provably inherently hard, with bilevel smooth problems being as difficult as general lower semicontinuous minimization and bilevel polynomial problems reaching complexity beyond NP-hard problems [10].

# 2 Presentation of the problem and main algorithms

## 2.1 The bilevel problem and the lower-level gradient routine

Given $x \in \mathbb{R}^n$, the global minimizer of $g_x := g(x, \cdot)$ does not have a closed-form expression and one often resorts to algorithms to solve approximately the lower level problem. In this work, we assume to have access to an iterative algorithm. Given $x \in \mathbb{R}^n$, one initializes $y_0 = z$ and run an iterative method $y_{i+1} = A(x, y_i)$ for a fixed number of iterations $i = 1, \dots, k$. Here, our algorithm $A$ is chosen to be the classical Cauchy's gradient descent (GD) algorithm with a fixed stepsize. Denoting by $A^k$, $k$ iterations of the algorithm, we consider $y = A^k(x, z)$ where $z$ is the initialization of the algorithm, see Algorithm LLD.

**Algorithm LLD** Lower level descent $A^k$ ($k$-step fixed-stepsize gradient descent)

**Require:** $x \in \mathbb{R}^n, z \in \mathbb{R}^m, \alpha_g > 0, k \in \mathbb{N}$.

1: $y_0 \leftarrow z$.

2: **for** $j = 1, \dots, k$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp;$y_j \leftarrow y_{j-1} - \alpha_g \nabla_y g(x, y_{j-1})$

4: **return** $y_k$

▷{gradient descent update}

We will work under the following assumptions which ensure that for any $k \in \mathbb{N}$, $A^k$ is a continuously differentiable mapping.

**Assumption 2.1** (Regularity of the upper and lower functions). Let $m, n > 0$.

(i) The functions $f, g: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ are $C^2$ and $C^3$ respectively, and semi-algebraic.

(ii) For all $x \in \mathbb{R}^n$, $g_x$ has a $L_g$-Lipschitz continuous gradient and is coercive.

(iii) The lower level step size parameter $\alpha_g$ satisfies $0 < \alpha_g L_g < 1$.

Since $g_x(\cdot) = g(x, \cdot)$ is not assumed to be convex, Algorithm LLD merely approaches critical points of $g_x$. On the other hand, it is known that, generically the algorithm finds local minima [53, 26].

This leads us to consider two fundamental relaxations of our original bilevel program: the *critical set relaxation* and the *local-min relaxation*. The latter is already considered and studied with the so-called Karush-Kuhn-Tucker (KKT) approach [1].

Denote by $\text{crit}_1 g$ the (set-valued) mapping $x \mapsto \text{crit}_1 g_x$, and $\text{argmin-loc}_1 g$ the (set-valued) mapping $x \mapsto \text{argmin-loc}_1 g_x$, our problems can actually be rewritten:
<table><thead><tr><th>Critical points relaxation</th><th>Local minima relaxation</th></tr></thead><tbody><tr><td>min <i>f</i>(<i>x</i>, <i>y</i>)<br>(<i>x</i>, <i>y</i>) ∈ graph <i>crit</i><sub>1</sub><i>g</i></td><td>min <i>f</i>(<i>x</i>, <i>y</i>)<br>(<i>x</i>, <i>y</i>) ∈ graph <i>argmin-loc</i><sub>1</sub><i>g</i></td></tr></tbody></table>

## 2.2 Two bilevel gradient algorithms

We consider two bilevel gradient methods with different strategies. For a given $k \in \mathbb{N}$, we consider the approximate value function objective,

$$ \varphi^k(x, z) := f(x, \mathcal{A}^k(x, z)) \qquad (\text{DP-BL}) $$

which may be seen as an unconstrained smooth optimization problem approximating (BL). Indeed, due to the convergence properties of the gradient method for semi-algebraic functions [6], under Assumption 2.1, the corresponding function has the property that the limit below exists and defines a function $\varphi^\infty$ as follows:

$$ \lim_{k \to \infty} \varphi^k(x, z) := \varphi^\infty(x, z) \\ \text{graph } \varphi^\infty = \left\{ ((x, z), f(x, y(x, z))) : y(x, z) = \lim_{k \to \infty} \mathcal{A}^k(x, z) \in \text{crit} g_x \right\}. $$

**Single-step Multi-step strategy.** We first consider a two-stage method: one approximates the lower-level critical points (or local minima) for a fixed $x$, in our case, by Algorithm LLD with $\mathcal{A}^k$, then performs a gradient step on the upper level. This is what we call, in this work, the *Single-step Multi-step* strategy.

<caption>**Algorithm SMBG** Single-step Multi-step bilevel Gradient Algorithm</caption>

**Require:** $\alpha_f > 0, \alpha_g > 0, k \in \mathbb{N}, L \in \mathbb{N}$.
1. Initialize $(x_0, y_0)$
2. **for** $\ell = 1, \dots, L$ **do**
3.   $y_\ell \leftarrow \mathcal{A}^k(x_{\ell-1}, y_{\ell-1})$ &rarr; $\{k \text{ GD steps for lower-level } g(x_{\ell-1}, \cdot)\}$
4.   $x_\ell \leftarrow x_{\ell-1} - \alpha_f \nabla_x \varphi^k(x_{\ell-1}, y_{\ell-1})$ &rarr; $\{\text{one GD step for upper-level}\}$
5.   **return** $(x_L, y_L)$

The alternation of update steps between each level is classical in bilevel programming. Variations close to the proposed description were explored recently motivated by machine learning applications [4, 24, 13, 17, 18, 33]. These works, however, harness the implicit function theorem and approximate implicit differentiation (AID), we consider the most basic version of this idea.

The discrepancy between the type of update and number of iterations dedicated to the two levels is natural: in the first place the lower level constraint must be (approximately) enforced, then and only then, progress can be made on the upper level.

**Differentiable programming strategy** Since $\varphi^k$ approximates $\varphi^\infty$ which corresponds to the critical point relaxation (BL-crit) of problem (BL), it is tempting to minimize directly the unconstrained, smooth, $\varphi^k$ for a given fixed $k$. This is the main idea of the *differentiable programming* strategy.
**Require:** $\alpha_f > 0, \alpha_g > 0, k \in \mathbb{N}, L \in \mathbb{N}$.

1: Initialize $(x_0, z_0)$.

2: for $\ell = 1, \dots, L$ do

3: $\begin{pmatrix} z_\ell \\ x_\ell \end{pmatrix} = \begin{pmatrix} z_{\ell-1} - \alpha_f \nabla_z \varphi^k(x_{\ell-1}, z_{\ell-1}) \\ x_{\ell-1} - \alpha_f \nabla_x \varphi^k(x_{\ell-1}, z_{\ell-1}) \end{pmatrix}$

4: return $(x_L, z_L)$

Such strategies are widespread in MAML [22] which formulate bilevel problems to find a common algorithmic initialization, good for several learning tasks. They combine the minimization of the smooth approximation $\varphi^k$ [23] with the idea that derivatives of $\varphi^k$ are accessible at a moderate cost using algorithmic differentiation [46, 47, 27]. In a broader context, the idea of minimizing the surrogate $\varphi^k$ jointly in $x$ and $z$ appears in [43, Algorithm 1]. In particular, they argue that when the lower-level problem $g_x$ admits multiple optimal solutions (or critical points), it is necessary to also update $z$ w.r.t. $\varphi^k$ to avoid initializing in a “bad” region, resulting in a sub-optimal $y^*(x)$ (remind that we want to minimize $f$ w.r.t. both $x$ and $y$ in (BL)).

**Remark 2.1** (Machine learning perspective). For both algorithms, the gradient estimation for $\mathcal{A}^k$ (cf. Algorithm LLD) and its compositions is called iterative differentiation (ITD) in the context in machine learning. This relies on an intensive usage of algorithmic differentiation routines [25]. In addition, Algorithm SMBG and Algorithm DPBG combine $k$ gradient steps on the inner problem with an algorithmic step on the outer problem, a feature which is described as a *double loop* strategy in the machine learning literature. This literature focuses on complexity estimates and stochastic approximation aspects. In contrast, the main specificity of our work is to consider a simpler algorithmic setting, but a broader class of problems with a focus on the difficulties caused by the absence of lower level convexity.

# 3 Morse-parametric constraint qualifications and the lower level gradient method

## 3.1 Morse-parametric constraint qualifications

Let us now consider the fundamental subject of constraint qualification. To establish necessary conditions for local minimizers, works on this question [62, 20, 49, 31] requires various assumptions such as the so-called *partial calmness* [62] and the inner-semicontinuity of the set-valued mappings related to the mapping $x \Rightarrow \operatorname{argmin} g_x$ [19, 49, 31]. Despite the mathematical insights of previous works, finding a class of optimization problems that satisfy these requirements is non-trivial.

In this section, we introduce a family of nonconvex bilevel problems whose constraints qualification is much simpler. In particular, we require the lower-level function $g$ to satisfy the so-called *parametric Morse* property, whose definition is given below.

**Definition 3.1** (Morse and parametric Morse functions). Let $\phi: \mathbb{R}^m \to \mathbb{R}$ be a $C^2$ function. We say that $\phi$ is a *Morse function* if its Hessian matrix $\nabla^2\phi(y)$ is invertible for any critical point $y$ of $\phi$, i.e., $\nabla\phi(y) = 0$. Let $g: U \times \mathbb{R}^m \to \mathbb{R}$ be a $C^2$ function where $U \subset \mathbb{R}^n$ is open and connected. The function $g$ is called *parametric Morse* on $U$ if, for all $x \in U$, $g_x := g(x, \cdot)$ is a Morse function. If $U = \mathbb{R}^n$, the function is simply called parametric Morse.

This is a particular case of the Morse-Bott parametric property, proposed in [4]. Further variants in bilevel optimization and machine learning were proposed in [12, 38]. Here are two examples of parametric Morse functions:
**Example 3.2.** (a) *Strongly convex functions:* If $g_x$ is a $C^2$ and $\mu$-strongly convex for all $x \in \mathbb{R}^n$, then $g$ is parametric Morse since $\partial_{yy}f(x,y) \supseteq \mu \mathbf{I}$. This class is usually assumed in many works on lower-level convex bilevel optimization such as [17, 18].

(b) *Composition of Morse functions and diffeomorphisms:* If $g = F \circ h$ where $F : \mathbb{R}^m \to \mathbb{R}$ is a $C^2$, Morse function and $h : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^m$ is $C^2$, each mapping $h(x, \cdot)$ is a diffeomorphism from $\mathbb{R}^m$ to itself, then $g$ is parametric Morse.

We introduce the constraint qualification assumptions for the bilevel optimization in Morse QC.

**Assumption Morse QC (Morse-parametric qualification conditions).** The function $g$ is parametric Morse and the set-valued mapping $\text{crit}_{1}g : x \mapsto \text{crit}_{1}g_x$ is non-empty and locally bounded.

**3.2 Relevance of the parametric Morse property**

In mathematical programming, constraint qualifications have two distinctive features. First, they allow to obtain an algebraic description of problem solutions, and second, they hold true for the typical, or generic, problems of the class. For example, Slater's interiority condition for convex programs can always be obtained by a slight perturbation of inequality constraints and is sufficient to ensure that all solutions are described using KKT conditions. This has numerous extensions in nonlinear programming, see for example [8] and references therein. In the bilevel programming context, qualification conditions satisfying these two qualitative features are missing.

We will see that, while Assumption Morse QC allows for a useful description of problem solutions, parametric Morse functions are not dense in the class of $C^2$ functions, and hence cannot be considered as typical in the $C^2$ topology (see Example 3.5). However, we argue that the parametric Morse property possibly has a central role in bilevel optimization. Indeed, a generic semi-algebraic bilevel problem satisfies the following *piecewise parametric Morse property*.

**Definition 3.3 (Piecewise parametric Morse functions).** A $C^2$ function $g : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ is piecewise parametric Morse if $g_x := g(x, \cdot)$ is Morse, for every $x$ in an open dense set.

The main differences with parametric Morse property are 1/ that the property is required to hold on a dense and open subset, and 2/ this set is not required to be connected. The piecewise parametric Morse property turns out to be generic for semi-algebraic functions.

**Proposition 3.4 (Genericity of piecewise parametric Morse functions).** Given a $C^2$, semi-algebraic function $g : \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$, the set of vectors $a \in \mathbb{R}^m$ such that $g(x,y) - \langle a, y \rangle$ is piecewise parametric Morse is semi-algebraic and dense in $\mathbb{R}^m$ (hence residual and of full measure).

Before proving the result, we recall that a semi-algebraic set is the finite union of solution sets of systems of polynomial inequalities. A semi-algebraic piecewise Morse parametric function is, therefore, Morse parametric in the sense of Definition 3.1 on each of the finitely many connected components of semi-algebraic dense and open set required by Definition 3.3. This result suggests that a general description of the solution set of a generic bilevel program requires a deep understanding of the simpler Morse parametric assumption. Therefore, this property appears as a relevant intermediate between widely studied strongly convex lower levels, to the most general nonconvex lower levels.

*Proof of Proposition 3.4.* The fact that the resulting set is semi-algebraic is a direct consequence of the Tarski-Seidenberg principle, see for example [15].

Define $A := \{(x, z) | z = \nabla_y g(x, y) \text{ where } \det(\nabla_{yy}^2 g(x, y)) = 0\}$. Since $g$ is semi-algebraic, so are $A$ and its fibers i.e., $A(x) = \{z | (x, z) \in A\}$. Moreover, for any $x$, $A(x)$ is also the set
of critical values of the mapping $\nabla_y g_x = \nabla_y g(x, \cdot)$, hence of zero Lebesgue measure (thanks to Sard theorem). By Fubini-Tonelli theorem, the Lebesgue measure of A is zero as well. Combine with the fact that A is semi-algebraic, A is of dimension at most $(m+n-1)$.

Consider the projection mapping (onto the last $m$ coordinate) $\pi: A \to \mathbb{R}^m$. As a consequence of Hardt's triviality theorem [16, Corollary 4.2], the set

$$P := \{z \in \mathbb{R}^m \mid \dim \pi^{-1}(z) = n\}$$

is semi-algebraic of dimension at most $\dim A - n \le m - 1$. This concludes the proof since for a generic $z \in \mathbb{R}^m$, $z \notin P$, therefore $\dim \pi^{-1}(z) < n$ and $g_x$ is Morse for any $x \notin \pi^{-1}(z)$. $\square$

Finally, we also illustrate why parametric Morse functions property are not dense among $C^2$ functions in the following example.

**Example 3.5.** Consider the function $g: \mathbb{R}^2 \to \mathbb{R}: (x, y) \mapsto (x - y^2)^2$. This function is not Morse parametric since its critical sets is given by:

$$\text{crit } g_x := \begin{cases} 0 & x \le 0 \\ \{0, \pm\sqrt{x}\} & x > 0 \end{cases}.$$

Hence, $g_x$ is not Morse at $x = 0$. Moreover, an arbitrarily $C^2$ small perturbation $e$ cannot make the function $g + e$ parametric Morse. Here, the size of a $C^2$ perturbation $e$, can be considered as $\|e\|_{C^2} = \sup_{x,y} |e(x, y)| + \|\nabla e(x, y)\| + \|\nabla^2 e(x, y)\|$. Note that

$$\sup_{x,y} |e_x(y)| + |e'_x(y)| + |e''_x(y)| \le \sup_{x,y} |e(x, y)| + \|\nabla e(x, y)\| + \|\nabla^2 e(x, y)\|, \quad (1)$$

since $\|\nabla e(x, y)\| \ge |e'_x(y)|$ and $\|\nabla^2 e(x, y)\| \ge |e''_x(y)|$ for all $x, y \in \mathbb{R}$. In other words $\|e\|_{C^2} \ge \sup_x \|e_x\|_{C^2}$. Assume that $\|e\|_{C^2} < 1$.

Returning to $g$, we have $g'_x(y) = 4y(y^2 - x)$ and $g''_x(y) = 12y^2 - 4x$. It is clear that the critical points of $g_x + e_x$ are locally bounded as $\lim_{|y| \to +\infty} \inf_{x \in K} |g'_x(y)| = +\infty$ for any compact interval $K$. Let us show that the number of critical points of $g_x + e_x$ cannot be constant in $x$.

For $x = -1$, it is clear that $g_{-1} + e_{-1}$ is strongly convex as $g''_{-1}(y) \ge 4$ and $|e''_{-1}(y)| \le 1$ for all $y$. So it has a single critical point which corresponds to its global minimum. Let us consider $x = 1$, we have the following variations for each interval $I$ as below, chosen to contain none or a single critical point of $g_1$

<table><thead><tr><th>I</th><th>(-∞, -2]</th><th>[-2, -1/√2]</th><th>[1/√2, 1/2]</th><th>[1/2, 1/√2]</th><th>[1/2, 2]</th><th>[1/√2, +∞)</th></tr></thead><tbody><tr><th>g'<sub>1</sub>(I)</th><td>≤ -24</td><td>⊂ [-1, 1]</td><td>≥ 1</td><td>⊂ [-1, 1]</td><td>≤ -1</td><td>⊂ [-1, 1]</td><td>≥ 24</td></tr></tbody></table>

From this table, we see that since $|e'_1(y)| < 1$ for all $y$, for each $I$ as above such that $g_1$ has no critical point, $g_1 + e_1$ also has no critical point, and for each $I$ such that $g_1$ has one critical point, $g_1 + e_1$ has at least one critical point. We conclude that $g_1 + e_1$ has at least three critical points.

Finally, the contradiction can be deduced from Lemma 3.7 (see Section 3.3 for more detail): supposing that $g$ is parametric Morse, since for each critical point of $g_1$, one can find a smooth function defined globally on $\mathbb{R}^2$ that gives the solution to $\nabla_y g(x, y) = 0$, $g_{-1}$ must have at least three distinct critical points (otherwise, these functions are identical due to the uniqueness stated in Lemma 3.7). This is, however, not possible due to our argument about $g_{-1}$.

All following results in Sections 3.3, 3.4 and 4 only hold under the assumption of parametric Morse property (cf. Definition 3.1), and not its piecewise version (cf. Definition 3.3). Nevertheless, as we argued earlier, due to the genericity of the piecewise parametric Morse property,
the parametric Morse property becomes an important intermediate case that is worth studying. Moreover, given a piecewise parametric Morse function $g$, if one supposes that the iterates always stay in the same connected component of the set $x$ where $g_x$ is Morse, then all of our results transfer naturally. An important future extension is to understand better the case where the iterates can “jump” between the connected components and identify the conditions in which interesting theoretical results (such as convergence) can be established.

### 3.3 Structure of the crit and argmin-loc mappings under Assumption Morse QC

The interest of Assumption Morse QC (and also Definition 3.1) lies in the fact that the constraint $y \in \text{crit } g(x, \cdot)$ can be written equivalently as: $y \in \{y^{(i)}(x) | i = 1, \dots, N\}$ where $N$ is a natural number and $y^{(i)}, i = 1, \dots, N$ are $C^1$ functions. This representation can actually be made of constant index, we only distinguish local minima from other critical points.

**Proposition 3.6** (Critical points and local minima are finite union of manifolds). Under Assumption Morse QC and $g$ being semi-algebraic $C^k$ ($k > 1$), there exists integers $M \ge N \ge 0$, and $M$ semi-algebraic $C^{k-1}$ functions $y^{(1)}, \dots, y^{(N)}, \dots, y^{(M)}: \mathbb{R}^n \to \mathbb{R}^m$ with distinct values at each points such that<fnref n="2" />:

$$ \text{graph } \text{crit}_1 g = \bigcup_{i=1}^{M} \text{graph } y^{(i)}. \qquad (2) $$

$$ \text{graph } \text{argmin-loc}_1 g = \bigcup_{i=1}^{N} \text{graph } y^{(i)}. \qquad (3) $$

*Proof. Proof of (2).* The graph of $\text{crit}_1 g$ is a semi-algebraic set (since $g$ is semi-algebraic) such that for each $x \in \mathbb{R}^n$, the fibers $\text{crit } g_x$, if non-empty, contain only isolated points thanks to the implicit function theorem. Hence, the connected components of $\text{crit } g_x$ are singletons. By [60, Properties 4.4], the number of connected components of $\text{crit } g_x$ is uniformly bounded above by an integer. Moreover, consider the following claim:

**Lemma 3.7.** Suppose that Assumption Morse QC holds such that $g$ is $C^k$ ($k > 1$), and let $(\bar{x}, \bar{y})$ be a solution to the equation $\nabla_y g(\bar{x}, \bar{y}) = 0$. There is a unique globally defined $C^{k-1}$ function $y: \mathbb{R}^n \to \mathbb{R}^m$ such that for any $x \in \mathbb{R}^n$, $y(x) \in \mathbb{R}^m$ is the unique solution to the equation $\nabla_y g(x, y) = 0$ locally, and $y(\bar{x}) = \bar{y}$.

Lemma 3.7 is proved at the end of this proof. Denote by $M(x)$ the size of $\text{crit } g_x$, Lemma 3.7 allows to conclude that this is a (finite) constant, i.e., $M(x) = M, \forall x \in \mathbb{R}^n$. Indeed, each solution at a given $\bar{x}$ might be extended to globally defined solution, and if any two of such functions are equal at a given $x$, they are equal on the whole space.

To finish the proof of (2), it remains to show that each function $y^{(i)}, i = 1, \dots, M$ is semi-algebraic. This is because for each $i$, graph $y^{(i)}$ is a connected component of the semi-algebraic set graph $\text{crit}_1 g$ (see [60, Property 4.3]).

**Proof of (3).** It is sufficient to prove that for each $i \in \{1, \dots, M\}$, either $y^{(i)}(x)$ is a local minimum of $g_x$ for all $x \in \mathbb{R}^n$ (positive definite Hessian), or $y^{(i)}(x)$ is a strict saddle point of $g_x$ for all $x$ (indefinite Hessian). Note that these are only two possibilities due to the parametric Morse assumption. By the continuity of the eigenvalues and Morse QC, we can conclude that $\nabla^2_{yy} g(x, y^{(i)}(x))$ is either positive definite for all $x$ or indefinite for all $x$. In fact, the index, i.e., the number of positive eigenvalues, of $\nabla^2_{yy} g(x, y^{(i)}(x))$ is constant for all $x$, as a change would result in a zero eigenvalue which contradicts the parametric Morse assumption. The result follows after permuting $\{1, \dots, M\}$ so that local minima functions take the first $N$ positions. □

<footnote n="2">²We use the convention that for $N = 0$ or $M = 0$ the union is empty</footnote>
We consider an intermediate claim before proving Lemma 3.7.

**Lemma 3.8.** Suppose that Assumption Morse QC holds. Let $B_1, B_2 \subset \mathbb{R}^n$ be two open balls and $y_1: B_1 \to \mathbb{R}^m$ and $y_2: B_2 \to \mathbb{R}^m$ be continuous solutions to $\nabla_y g(x, y) = 0$. If $y_1(x) = y_2(x)$ for some $x \in B_1 \cap B_2$, then $y_1(x) = y_2(x)$ for all $x \in B_1 \cap B_2$.

*Proof of Lemma 3.8.* Set $\mathcal{E} := \{x \in B_1 \cap B_2 | y_1(x) = y_2(x)\}$. This set is closed (in the induced topology on $B_1 \cap B_2$) because it is the pre-image of 0 by the continuous function $y_1 - y_2$. It is also open, thanks to the implicit function theorem. Since $B_1 \cap B_2$ is convex, hence path connected, we have either $\mathcal{E} = \emptyset$ or $\mathcal{E} = B_1 \cap B_2$. The conclusion follows because $\mathcal{E} \neq \emptyset$ by assumption. $\square$

*Proof of Lemma 3.7.* Define $R$ as:

$$R = \sup_{r>0} \left\{ r \mid \exists y_r: B(\bar{x}, r) \to \mathbb{R}^m, C^{k-1} \text{ s.t. } \forall x \in B(\bar{x}, r), \nabla_y g(x, y_r(x)) = 0, y_r(\bar{x}) = \bar{y} \right\}. \quad (4)$$

We first argue that there exists a unique $C^{k-1}$ function $y: B(\bar{x}, R) \to \mathbb{R}^m$ such that $\nabla_y g(x, y(x)) = 0$ and $y(\bar{x}) = \bar{y}$. Indeed, by definition of $R$, for any $0 < a \le b < R$, there exists two functions $y_a: B(\bar{x}, a) \to \mathbb{R}^m$ and $y_b: B(\bar{x}, b) \to \mathbb{R}^m$ satisfying the conditions of Equation (4). Lemma 3.8 ensures that $y_a(x) = y_b(x), \forall x \in B(\bar{x}, a)$. For any $x \in B(\bar{x}, R)$, we set $y(x) = y_r(x)$ for any $r > 0$ such that $\text{dist}(\bar{x}, x) < r < R$. This is well defined since, by the argument above the chosen value does not depend on $r$. In addition, $y$ is trivially $C^{k-1}$ (since $y_r$ is $C^{k-1}$ by the implicit function theorem), and unique by Lemma 3.8.

For the sake of contradiction, suppose that $R$ is finite. Using the boundedness assumption, for any $\tilde{x} \in \partial B(\bar{x}, R)$ – the boundary of $B(\bar{x}, R)$, there is a solution $\tilde{y}$ which is an accumulation point of $y(x)$ as $x \to \tilde{x}$. By the implicit function theorem, there is an open ball neighborhood of $\tilde{x}$, $V_{\tilde{x}}$ and a $C^{k-1}$ function $f: V_{\tilde{x}} \to \mathbb{R}^m$ such that $\nabla_y g(x, f(x)) = 0, \forall x \in V_{\tilde{x}}$, and $f(\tilde{x}) = \tilde{y}$. Since the boundary is compact, it can be covered by finitely many such open balls, says $B_i$, $i = 1, \dots, K$ with centers on the sphere. In each open ball $B_i$, let $x_i$ be its center and $f_i: B_i \to \mathbb{R}^m$ be the $C^{k-1}$ function satisfying $f_i(x_i) = y_i$. Denote $B(\bar{x}, R)$ by $B_0$ and $f_0 := y$, we extend the function $y$ to the function $\tilde{y}$ of a larger domain $\mathcal{B} := \bigcup_{i=0}^K B_i$ as follows:

$$\tilde{y}(x) = f_i(x) \text{ if } x \in B_i. \quad (5)$$

A contradiction results from the following two claims:

1. $\tilde{y}(x)$ in (5) is well defined for all $x \in \mathcal{B}$, and the resulting function is $C^{k-1}$.

2. There exists $\delta > 0$ such that $B(\bar{x}, R + \delta) \subseteq \mathcal{B}$.

In particular, the second claim contradicts the fact that $R < +\infty$ is the maximum, concluding the proof.

1. The well-definedness of $\tilde{y}$: we need to show that the value in (5) does not depend on $i$. It is sufficient to prove that for any $0 \le i < j \le k$, $f_i(x) = f_j(x)$, for all $x \in B_i \cap B_j$. Fix any $i < j$ such that $B_i \cap B_j \neq \emptyset$, by Lemma 3.8, it suffices to show that the set $\{x \in B_i \cap B_j | f_i(x) = f_j(x)\}$ is non-empty. Indeed if $i = 0$, the implicit function theorem ensures local uniqueness of solutions to $\nabla_y g(x, y) = 0$ around $(x_j, y_j)$, and $y_j$ is a point of accumulation of $f_0$ on $B_0$, so that we have $f_i(x) = f_j(x)$ for some $x \in B_0 \cap B_j$ close to $x_j$. If both $i, j > 0$, since $x_i, x_j$ belong to $\partial B_0$, there exists $x \in B_i \cap B_j \cap B_0$ and, by the argument for $i = 0$, we have $f_i(x) = f_0(x) = f_j(x)$.

2. The existence of $\delta$: we claim that $\delta := \text{dist}(\overline{B(\bar{x}, R)}, \mathcal{B}^C) > 0$. Indeed, this is the distance between a compact ball i.e., $\overline{B(\bar{x}, R)}$ and a closed set (the complement of the open set $\mathcal{B}$), so the infimum is attained. If $\delta = 0$, there is a point on the sphere not belonging to any $B_i$, which contradicts the fact that we have a covering. Hence $B(\bar{x}, R + \delta) \subset \mathcal{B}$.
□

**Remark 3.9.** If in Proposition 3.6, one replaces the condition *g* is semi-algebraic by the uniform boundedness of the cardinality of crit *g*<sub>*x*</sub>, the conclusions of Proposition 3.6 still hold, but *y*<sup>(i)</sup> is not necessarily semi-algebraic.

Figure 1 illustrates the claims in Proposition 3.6. We remark that the same claim for global minimizers of $g_x$ is not true in general. Indeed, under Assumption Morse QC, $\argmin_{x=1,...,N} g(x, y^{(i)}(x))$, which is generally only piecewise smooth and can be set-valued for arguments $x$ with multiple solutions.

**Remark 3.10** (Necessity of local-boundedness assumption). Consider the function $g = \tau \circ h$ where $\tau$ is the real-valued function $\tau : x \mapsto \frac{x^2}{2}$ and $h : \mathbb{R} \times \mathbb{R} \to \mathbb{R}$ defined by $h : (x, y) \mapsto x + y/\sqrt{y^2 + 1}$. Note that $h(x, \cdot)$ is a diffeomorphism from $\mathbb{R}$ to the open interval $(x - 1, x + 1)$, for all $x \in \mathbb{R}$. Hence, $g$ is a parametric Morse function (cf. Example 3.2b). Its sets of critical points (and also local/global minima) of a fixed $x$ (given below) is not the graph of a $C^1$ function due to escape at infinity.

$$ \text{crit } g_x = \text{argmin-loc } g_x = \text{argmin } g_x = \begin{cases} \left\{ -\frac{x}{\sqrt{1-x^2}} \right\} & \text{if } x \in (-1, 1), \\ \emptyset & \text{otherwise} \end{cases}. $$

A consequence of the above is that, under Assumption Morse QC, the relaxed bilevel optimization problems (BL-crit) and (BL-loc) can be rewritten as mixed continuous-discrete variable programs in the following ways (Table 1).

<table><caption>**Table 1:** Equivalent reformulation of (BL-crit) and (BL-loc) under Assumption Morse QC. The constants $M, N$ are defined as in Proposition 3.6.</caption><thead><tr><th></th><th>Critical points relaxation</th><th>Local minima relaxation</th></tr></thead><tbody><tr><td>Constrained optimization</td><td>$$ \begin{array}{l} \min_{x \in \mathbb{R}^n, y \in \mathbb{R}^m} f(x, y) \\ \text{s.t. } y \in \{y^{(i)}(x), 0 \le i \le M\} \end{array} $$</td><td>$$ \begin{array}{l} \min_{x \in \mathbb{R}^n, y \in \mathbb{R}^m} f(x, y) \\ \text{s.t. } y \in \{y^{(i)}(x), 0 \le i \le N\} \end{array} $$</td></tr><tr><td>Mixed-integer optimization</td><td>$$ \min_{x \in \mathbb{R}^n, 1 \le i \le M} f(x, y^{(i)}(x)) $$</td><td>$$ \min_{x \in \mathbb{R}^n, 1 \le i \le N} f(x, y^{(i)}(x)) $$</td></tr></tbody></table>

## 3.4 Lower level gradient descent under parametric Morse QC

**Uniform linear convergence** The following lemma states that if $z^* \in \text{argmin-loc } g_{x^*}$ (cf. eq. (3)), then the algorithm $\mathcal{A}$ behaves uniformly in terms of convergence in a neighborhood of $(x^*, z^*)$. These are classical results which we recall for completeness.

**Lemma 3.11** (Uniform linear convergence near a local minimum). Under Assumption 2.1, if $g$ is parametric Morse, for any $z^* \in \text{argmin-loc } g_{x^*}$, denote by $y^*$ the implicit function associated to the equation $\nabla_y g(x, y) = 0$, defined in a neighborhood of $x^*$ and set $\varphi: x \mapsto f(x, y^*(x))$. Then, there exists a neighborhood $U$ of $(x^*, z^*)$ and constants $0 \le \rho < 1$, $C > 0$, such that for all $(x, z) \in U$, we have for all $k \in \mathbb{N}$ and all $(x, z) \in U$:

$$ \begin{align*} \| \mathcal{A}^k(x, z) - y^*(x) \| &\le \rho^k \| z - y^*(x) \|, \\ \| \nabla_z \varphi^k(x, z) \| &\le C \rho^k, \\ \| \nabla_x \varphi^k(x, z) - \nabla_y \varphi(x) \| &\le C \rho^k, \end{align*} $$

where $\mathcal{A}$ is defined in Algorithm LLD and $\varphi^k$ is the function of (DP-BL).
Proof. Since we have a Morse local minimizer, $\partial_{yy}^2 g(x^*, z^*)$ is positive definite. Let $U = B_1 \times B_2$ be the product of two Euclidean balls centered at $x^*$ and $z^*$ such that:

* $\partial_{yy}^2 g$ remains positive definite on $U$.

* For any $x \in B_1$, we have $y^*(x) \in B_2$. In other words, for all $(x, y) \in U$, $y^*(x) \in B_2$.

The first item can be satisfied using the continuity of $\partial_{yy}^2 g$ and the second one can be satisfied by reducing $B_1$ if necessary, using the local continuity of $y^*$. We have the following properties, for all $x \in B_1$:

* The function $g_x$ is strongly convex over $B_2$.

* $\operatorname{argmin}_{B_2} g_x$ is achieved.

It is well known that the gradient recursion is a strict contraction with a factor $\rho < 1$ for strongly convex objectives and that its unique fixed point is the global minimizer. Without loss of generality, we may assume that the same $\rho < 1$ is a uniform bound on the quantity $\|\mathbf{I} - \alpha_g \partial_{yy}^2 g(x, y)\|_{\text{op}}$ on $U$. All algorithmic iterates remain in the chosen neighborhood because of the algorithmic contraction property. This concludes the proof of the first point.

Regarding derivatives, for a given $(x, z) \in U$, we introduce the following notations for all $k \in \mathbb{N}$, and $(x, z) \in U$,

$$D_k(x, z) = \frac{\partial}{\partial x} A^k(x, z)$$

$$A_k(x, z) = \mathbf{I} - \alpha_g \partial_{yy}^2 g(x, A^k(x, z))$$

$$A(x) = \mathbf{I} - \alpha_g \partial_{yy}^2 g(x, y^*(x))$$

$$C_k(x, z) = \frac{\partial}{\partial z} A^k(x, z)$$

$$B_k(x, z) = -\alpha_g \partial_{yx}^2 g(x, A^k(x, z))$$

$$B(x) = -\alpha_g \partial_{yx}^2 g(x, y^*(x))$$

We have the following identities from the chain rule of differential calculus and implicit differentiation:

$$C_{k+1}(x, z) = A_k(x, z)C_k(x, z),$$

$$D_{k+1}(x, z) = A_k(x, z)D_k(x, z) + B_k(x, z),$$

$$\text{Jac } y^*(x) = A(x)\text{Jac } y^*(x) + B(x).$$

We have for all $k$, uniformly in $x, z$, $\|A_k(x, z)\|_{\text{op}} \le \rho$. We conclude that $\|C_k(x, z)\| \le \|C_0(x, z)\| \rho^k = m\rho^k$ which proves the second point. We also have the following:

$$\begin{aligned} D_{k+1}(x, z) - \text{Jac } y^*(x) &= A_k(x, z)D_k(x, z) + B_k(x, z) - A(x)\text{Jac } y^*(x) - B(x) \\ &= A_k(x, z)(D_k(x, z) - \text{Jac } y^*(x)) + (A_k(x, z) - A(x))\text{Jac } y^*(x) + B_k(x, z) - B(x) \end{aligned} \quad (6)$$

Since $g$ is $C^3$, we have that $\partial_{yy}^2 g$ and $\partial_{yx}^2 g$ are both Lipschitz on $U$, say with constant $M/\alpha_g$. By continuity, we may also consider that $H$ a uniform bound on $\|\text{Jac } y^*(x)\|$ and set $K = M(H+1)$. Combining with the first convergence estimate, we have for all $k$:

$$\|A_k(x, z) - A(x)\| \le M\rho^k \|z - y^*(x)\|,$$

$$\|B_k(x, z) - B(x)\| \le M\rho^k \|z - y^*(x)\|.$$

Getting back to (6), setting $\delta_k = \|D_k(x, z) - \text{Jac } y^*(x)\|$ for all $k$, we have the recursion

$$\delta_{k+1} \le \rho\delta_k + M\rho^k \|z - y^*(x)\| \|\text{Jac } y^*(x)\| + M\rho^k \|z - y^*(x)\| \le \rho\delta_k + K\rho^k \|z - y^*(x)\| \quad (7)$$
with $\delta_0 \le H \le K$. We verify using induction that for all $k \in \mathbb{N}$, $\delta_k \le K\rho^k + kK\|z - y^*(x)\|\rho^{k-1}$. This is true for $k = 0$, and by induction

$$ 
\begin{aligned}
\delta_{k+1} &\le \rho\delta_k + K\rho^k\|z - y^*(x)\| \\
&\le K\rho^{k+1} + kK\|z - y^*(x)\|\rho^k + K\rho^k\|z - y^*(x)\| \\
&= K\rho^{k+1} + (k+1)K\|z - y^*(x)\|\rho^k,
\end{aligned}
$$

the induction holds. Increasing $\rho$ slightly and choosing an appropriate constant $C > 0$, these estimates may be upper bounded by $C\rho^k$.

We have, for all $k \in \mathbb{N}$,

$$ 
\begin{aligned}
\nabla_z\varphi^k(x, z) &= C_k(x, z)^\top \nabla_y f(x, \mathcal{A}^k(x, z)) \\
\nabla_x\varphi^k(x, z) &= \nabla_x f(x, \mathcal{A}^k(x, z)) + D_k(x, z)^\top \nabla_y f(x, \mathcal{A}^k(x, z)) \\
\nabla\varphi(x) &= \nabla_x f(x, y^*(x)) + \text{Jac } y^*(x)^\top \nabla_y f(x, y^*(x))
\end{aligned}
\quad (8) $$

and the third estimate follows, eventually increasing $C$, using the first estimate, since the gradient of $f$ is Lipschitz continuous and bounded on $U$ since $f$ is $C^2$. $\square$

**Convergence to local minima with high probability over initialization** Another interesting property of parametric Morse functions is that with high probability over initialization $(x, z)$, Algorithm LLD eventually finds a local minimum of $g_x$ for large $k$. This is of course in the line of classical results for Morse functions [53, 26] and based on the stable-unstable manifold theorems.

**Proposition 3.12** (Uniform convergence to local minima over compact initialization). Let $g$ satisfy all the assumptions in Assumption 2.1 – except semi-algebraicity – and assume in addition that $g$ is parametric Morse 3 (cf. Definition 3.1). For any compact set $C \subseteq \mathbb{R}^n \times \mathbb{R}^m$, any $\delta, \Delta > 0$, there exists a set $A \subseteq \mathbb{R}^n \times \mathbb{R}^m$ (only dependant on $\Delta$) and a natural number $K \in \mathbb{N}$ (dependant on $C, \delta$ and $\Delta$) such that:

i) The Lebesgue measure of $A$ is smaller than $\Delta$, i.e., $\mu(A) \le \Delta$.

ii) We have:

$$ \forall (x, z) \in C \setminus A, \forall k \ge K \quad \text{dist}(\mathcal{A}^k(x, z), \text{argmin-loc}_{g_x}) \le \delta. $$

$$ \text{where } \text{dist}(x, A) := \inf_{y \in A} \|x - y\|_2. $$

*Proof.* We use two following facts:

1. Since $g_x$ is Morse ($g$ is parametric Morse), $g_x$ satisfies the KL inequality [5, Section 4]. By assumptions, $g$ is $L_g$-smooth, coercive, and the learning rate $0 < \alpha_g < 1/L_g$. Using [6, Theorem 3.2], we can prove that for any $(x, z) \in \mathbb{R}^m \times \mathbb{R}^n$, the sequence $\mathcal{A}^k(x, z)$ will converge to a critical point of $g_x$.

2. Due to the result of [53, 26, 39, 50], for a fixed $x$, the set of $z$ such that the sequence $\mathcal{A}^k(x, z)$ converges to a strict saddle point (i.e. a critical point $y$ of $g_x$ such that $\nabla_{yy}^2 g(x, y)$ has at least one negative eigenvalue) is of Lebesgue measure zero. Using the Fubini-Tonelli theorem, we have the set

$$ \mathcal{B} := \{(x, z) \mid \{\mathcal{A}^k(x, z)\}_{k \in \mathbb{N}} \text{ converges to a saddle point of } g_x\}, $$

is of zero Lebesgue measure.

<footnote n="3">³Note that we do not assume local boundedness of $\text{crit}_{g_x}$ as in Assumption Morse QC.</footnote>
Consider $\tilde{B}$ the complement of the set $B$. By combining the first two remarks, we can conclude that $\tilde{B}$ contains $(x, z)$ whose $\{A^k(x, z)\}_{k \in \mathbb{N}}$ converges to a local minimum $z^*$ of $g_x$.

To finish the proof, for any $\Delta > 0$, we take an open set $A$ such that $\mu(A) \le \Delta$ such that $B \subseteq A$. Note that this is possible because $\mu(B) = 0$. Since $A$ is open, $C \setminus A$ is compact.

Consider a point $(x, z) \in C \setminus A$. Let $z^*$ be the local minimizer of $g_x$ that $A^k(x, z)$, $k \in \mathbb{N}$ converges to. By Lemma 3.11, there exists an open neighborhood $U_1(x, z^*)$ of $(x, z^*)$ such that $\forall (x', z') \in U_1(x, z^*)$, $\forall k \in \mathbb{N}$, $\|A^k(x', z') - \bar{y}(x')\| \le \|z' - \bar{y}(x')\| \rho^k$ for some $0 < \rho < 1$ and the function $\bar{y}$ is defined as in Lemma 3.11. Note that $\bar{y}(x')$ is also the local minimizer of $g(x', \cdot)$ (using the same argument as in the proof of Equation (3)).

Let $k(x, z)$ be the smallest number such that: $(x, A^{k(x,z)}(x, z)) \in U_1(x, z^*) \cap B((x, z^*), \delta)$ (which exists due to the convergence assumption). Since $A^{k(x,z)}$ is a continuous function in both variables $x$ and $z$, there exists an open neighborhood $U_2(x, z)$ such that its image under $A^{k(x,z)}$ is included in $U_1(x, z^*)$. By construction, we have in addition:

$$ \|\bar{y}(x') - A^k(x', z')\| \le \|\bar{y}(x') - A^{k(x,z)}(x', z')\| \le \delta, \forall k \ge k(x, z), \forall (x', z') \in U_2(x, z). $$

Note that $\mathcal{U} := \{U_2(x, z) \mid (x, z) \in C \setminus A\}$ is an open coverage of the compact set $C \setminus A$. Therefore, there exists a finite subset $D \subseteq C \setminus A$ of such that:

$$ C \setminus A = \bigcup_{(x,z) \in D} U(x, z). $$

This implies the result immediately because one can take $K = \max_{(x,z) \in D} k(x, z)$. $\square$

**Remark 3.13** (On the necessity of compactness). Let us observe that the compactness of $C$ is necessary: a trivial counterexample is $g(x, y) = (x - y)^2$ and $C = \mathbb{R}^{n+m}$. In this setting, there does not exist any $A, K$ as in Proposition 3.12 (and we leave the verification of this claim for the readers)

**Attractivity of the connected components of graph (argmin-loc$_1$g)** We shall also need a modified version of the attraction properties of graph (argmin-loc$_1$g) on a compact set.

**Lemma 3.14** (Attraction properties of the lower level constraints). Consider a compact set $S \subseteq \mathbb{R}^n$. Under Assumption Morse QC and Assumption 2.1, there exist constants $r := r(S) > 0$, $\rho := \rho(S) \in (0, 1)$, $c := c(S) > 0$ such that for all $(x, z) \in S \times \mathbb{R}^m$:

1. If $\text{dist}((x, z), \text{graph argmin-loc}_1 g) \le r$, there is a unique $C^2$ function $y$ with graph $y \subseteq \text{argmin-loc}_1 g$ such that:

$$ \|A^k(x, z) - y(x)\| \le \rho^k \|z - y(x)\|. $$

2. If $\text{dist}((x, A^l(x, z)), \text{graph argmin-loc}_1 g) \le r$, there is a unique $C^2$ function $y$ with graph $y \subseteq \text{argmin-loc}_1 g$ such that:

$$ \|\nabla_x \varphi^{k+l}(x, z) - \nabla_x \varphi^\infty(x, z)\|_2 \le (\|\text{Jac}_x A^l(x, z) - \text{Jac}_x y(x)\| + k)c\rho^k $$

$$ \text{where } \varphi^k(x, z) := f(x, A^k(x, z)) \text{ and } \varphi^\infty(x, z) := f(x, y(x)). $$

*Proof.* Using Assumption Morse QC and Assumption 2.1), thanks to Proposition 3.6, there exists $N$ functions $y^{(1)}, \dots, y^{(N)}$ which are $C^2$ and have distinct values such that $\text{argmin-loc}_1 g$ is partitioned into their graphs, i.e., eq. (3). Therefore, the condition:

$$ \text{dist}((x, z), \text{graph argmin-loc}_1 g) \le r $$

is equivalent to:

$$ \min_{i=1,\dots,N} \text{dist}((x, z), \text{graph } y^{(i)}) \le r. $$
Using this observation, we construct $(r, \rho, c)$ as follows: For any $(x, z) \in \text{argmin-loc}_1 g \cap \mathcal{S} \times \mathbb{R}^m$, let $y = y^{(i)}$ (for some $i \in \{1, \dots, N\}$) be the implicit function given by $\nabla_y g(x, y) = 0$ in a neighbor of $(x, z)$. By Lemma 3.11, there exists an open ball centered at $(x, z)$ with radius $r(x, z)$ such that for any point $(\bar{x}, \bar{z})$ in this open ball, there are two constants $K(x, z) > 0, \rho(x, z) \in (0, 1)$ depending on $(x, z)$ satisfying:

1. Using the bound in Lemma 3.11, we immediately get:

$$ \|\mathcal{A}^k(x, z) - y(x)\| \le \rho(x, z)^k \|z - y(x)\|. \quad (9) $$

2. By applying (7) $k$ times, we obtain:

$$ \| \text{Jac}_x \mathcal{A}^{k+l}(x, z) - \text{Jac}_x y(x) \| \le (\| \text{Jac}_x \mathcal{A}^l(x, z) - \text{Jac}_x y(x) \| + K(x, z)k) \rho(x, z)^k \quad (10) $$

Combining (8), (9), (10) and the Lipschitz continuity of $\nabla_x f, \nabla_y f$ in the ball centered at $(x, z)$ with radius $r(x, z)$, we conclude that there is a constant $c(x, z)$ such that:

$$ \| \nabla_x \varphi^{k+l}(x, z) - \nabla_x \varphi(x) \|_2 \le (\| \text{Jac}_x \mathcal{A}^l(x, z) - \text{Jac}_x y(x) \| + K(x, z)k) c(x, z) \rho(x, z)^k. $$

By Equation (3), $[\text{graph argmin-loc}_1 g] \cap (\mathcal{S} \times \mathbb{R}^m)$ is the union of $N$ graphs of $C^2$ functions $y^{(1)}, \dots, y^{(N)}$, restricted to $\mathcal{S}$. By compactness, there exists a finite set: $\mathcal{F} \subseteq [\text{graph argmin-loc}_1 g] \cap (\mathcal{S} \times \mathbb{R}^m)$ and a finite covering:

$$ [\text{graph argmin-loc}_1 g] \cap (\mathcal{S} \times \mathbb{R}^m) \subseteq \bigcup_{(x,z) \in \mathcal{F}} B((x,z), r(x,z)). $$

The lemma can be proved by choosing

$$ r := \min_{(x,z) \in \mathcal{F}} \frac{r(x,z)}{2} > 0, \quad c := \max_{(x,z) \in \mathcal{F}} c(x,z) > 0, \quad \rho := \max_{(x,z) \in \mathcal{F}} \rho(x,z) \in (0,1). $$

It remains to show the uniqueness of the function $y$. Note that with this choice of $r$, from each $(x, z)$ whose distance from $\text{argmin-loc}_1 g$ is at most $r(\mathcal{S})$, $y$ is the only $C^2$ function whose $y(x) = \lim_k \mathcal{A}^k(x, z)$. $\square$

# 4 Convergence of the single-step multi-step gradient method

Under Assumption 2.1, Assumption Morse QC), and natural assumptions on $f$, Algorithm SMBG is actually an *inexact gradient descent method* [59, 11, 4] with non-vanishing steps on a function $f(x, y(x))$ where $y(x)$ is a $C^2$ function giving a local minimizers of $g_x$.

Formally, this can be stated as follows. Let $\varphi: x \mapsto f(x, y(x))$ be the value function. For all $\epsilon > 0$, there is a lower bound on the number of lower level gradient steps $k$ in Algorithm LLD, such that Algorithm SMBG can be rewritten as follows

$$ x_{\ell+1} = x_{\ell} - \alpha_f (\nabla_x \varphi(x) + \xi_{\ell}) \quad \text{where} \quad \|\xi_{\ell}\| \le \epsilon \quad \text{for all } \ell \in \mathbb{N}. $$

## 4.1 Main convergence result

Before stating our main result let us recall a definition.

**Definition 4.1** ($\epsilon$-critical points and critical values). Given a differentiable function $h: \mathbb{R}^n \to \mathbb{R}$, the sets of $\epsilon$-critical points $\text{crit}_\epsilon h$ and $\epsilon$-critical values $\text{vcrit}_\epsilon h$ of $h$ are respectively defined defined through:

$$ \text{crit}_\epsilon h := \{x \mid \|\nabla h(x)\| \le \epsilon\}, $$

$$ \text{vcrit}_\epsilon h := h(\text{crit}_\epsilon h). $$
In the following, we consider a random choice of initialization whose distribution is absolutely continuous w.r.t. Lebesgue measure. It covers all practically used random initializations such as those based on Gaussian or uniform distributions, see for example [30].

**Theorem 4.2** (Convergence property of Algorithm SMBG). Suppose that Assumption 2.1 and Assumption Morse QC hold and assume in addition that:

* there exists some $\lambda > 0$ such that $\inf_{x,y} f(x, y) - \lambda \|x\|^2 > -\infty$,

* the initialization $x_0$ is chosen randomly with a distribution absolutely continuous with respect to the Lebesgue measure.

Then, for any $\epsilon, \Delta > 0$, there is a sufficiently large inner iterations $k$ (cf. Algorithm LLD) and constants $C > 0$ and $\rho \in (0, 1)$, such that with probability at least $1 - \Delta$, there exists $y = y^{(i)}$, $i \in \{1, \dots, N\}$ (cf. Equation (3)) satisfying:

$$ \lim_{\alpha_f \to 0^+} \limsup_{\ell \to \infty} \operatorname{dist}(\varphi(x_\ell(\alpha_f)), \operatorname{vcrit}_\epsilon \varphi) = 0 \\ \limsup_{\alpha_f \to 0^+} \limsup_{\ell \to \infty} \operatorname{dist}(x_\ell(\alpha_f), \operatorname{crit} \varphi) \le C \epsilon^\rho. $$

where $\varphi(x) = f(x, y(x))$ is the value function, and $x_\ell(\alpha_f)$ is the sequence generated by Algorithm SMBG with step size $\alpha_f$.

**Remark 4.3.** In Theorem 4.2, the assumption that $f$ is semi-algebraic (cf. Assumption 2.1) and $\inf_{x,y} f(x, y) - \lambda \|x\|^2 > -\infty$ is equivalent to $f(x, y) = h(x, y) + \lambda \|x\|^2$ with $h(x, y)$ is semi-algebraic and bounded from below. This is a coercivity assumption. This condition is satisfied in many practical settings (e.g., $f$ is a regularized version of a semi-algebraic non-negative data fitting term or a distance function $h(x, y)$). In addition, given a semi-algebraic function $f$ where $\inf f > -\infty$, by adding a $\| \cdot \|^2$ perturbation, the assumption in Theorem 4.2 on $f$ will be satisfied.

Except Proposition 3.12 and Lemma 3.11, semi-algebraicity is necessary in the announcement of Theorem 4.2 and intermediate results used in its proof. They ensure that the minimum values of certain finite quantities are strictly bigger than zero and that the inexact gradient method approximates the set of critical points [11]. Relaxing this condition is an interesting direction for future works, given our remark Remark 3.9 on the uniform boundedness of the cardinality of $\operatorname{crit} g_x$ and [11, Remark 3.2].

Before proving Theorem 4.2, we comment on related results in the literature, notably [4, Proposition 3] which studies the same algorithm:

1. Our proof of Theorem 4.2 is based on Lemma 4.4 and Lemma 4.5, which imply that $y_\ell$ stays close to a local minimum of $g_{x_\ell}$, for all $\ell \ge 1$, with high probability over initialization. It allows us to relax the assumption used in [4, Proposition 1]: $(x_\ell, y_\ell)$ converges to $(x^*, y^*)$ where $y^* \in \operatorname{argmin-loc} g(x^*, \cdot)$.

2. The fact that Algorithm SMBG is equivalent to a biased gradient descent on a certain solution function was previously observed, see [36] (strongly convex $g_x$) and [4] (local analysis for non-convex $g_x$). Our improvement on this aspect follows from the coercivity assumption on $f$ (cf. Remark 4.3). This is sufficient for the sequence to be bounded and allows for a global analysis of Algorithm SMBG for non-convex $g_x$ without a priori assumption on the sequence. In addition, by harnessing the semi-algebraicity, we obtain estimates for the distance between the sequence $\{x_\ell\}_{\ell \in \mathbb{N}}$ and the set of critical points of the bilevel problem, even for finite (and sufficiently large) $k$, whereas existing results [4, 43] have to settle with the limit $k \to \infty$ (remind that $k$ is the number of inner iterations, and not to be confused with the number of outer iterations $\ell$).
In conclusion, we think it is possible to relax the parametric Morse assumption to the parametric Morse-Bott one in Theorem 4.2 by combining analysis from this article and [4]. While this relaxation might be non-trivial and mathematically interesting, its significance is likely to be minor since the parametric Morse-Bott property is not generic, as that is the case for the parametric Morse property (see Example 3.5).

## 4.2 Gradient descent on the lower level

We first prove that Algorithm LLD reaches a neighborhood of local minimizer of $g_{x_0}$ for a typical initialization. The higher the number of iterations $k$, the smaller the neighborhood and the error term $||\xi_0||$.

**Lemma 4.4** (First iterations approximate local minima). Suppose that Assumption 2.1 and Assumption Morse QC holds.

For any compact $C \subset \mathbb{R}^n \times \mathbb{R}^m$, any $\Delta, r, \epsilon > 0$, there exists $A \subset \mathbb{R}^n \times \mathbb{R}^m$ with Lebesgue measure at most $\Delta$, an integer $K \in \mathbb{N}$ and $\alpha > 0$ such that if one runs Algorithm SMBG with $k \ge K$, $0 < \alpha_f \le \alpha$ and $(x_0, y_0) \in C \setminus A$, there exists a $C^2$ function $y$ such that graph $y \subseteq$ graph $\text{argmin-loc}_1 g$ (see Proposition 3.6) and that:

$$ \begin{aligned} & ||\nabla_x \varphi^k(x_0, y_0) - \nabla_x \varphi(x_0)|| \le \epsilon, \\ & \text{dist}((x_1, y_1), \text{graph } y) \le r, \end{aligned} $$

where $\varphi(x) = f(x, y(x))$, $\varphi^k(x, z) = f(x, A^k(x, z))$, $(x_0, y_0), (x_1, y_1)$ are the zeroth (initialization) and the first iterates of Algorithm SMBG, respectively.

*Proof.* We use several results and the notation of Lemma 3.14. By Proposition 3.12, we choose $\bar{K} := K(S, r'/2, \Delta)$ where:

1. $S = \Pi(C)$ where $\Pi(\cdot)$ is the projection onto the first $n$ coordinates.

2. $r' = \min(r(S), r)$, $r(\cdot)$ is defined as in Lemma 3.14).

3. $A \subseteq \mathbb{R}^n \times \mathbb{R}^m$, has Lebesgue measure at most $\mu(A) \le \Delta$ as given by Proposition 3.12.

Then for any $(x_0, y_0) \in C \setminus A$ we have:

$$ \text{dist}(A^{\bar{K}}(x_0, y_0), \text{argmin-loc}_{g_{x_0}}) \le r'/2. $$

Therefore, there exists a $C^2$ function $y = y^{(i)}$, $i \in \{1, \dots, N\}$ (where $N$ and $y^{(i)}$ are defined as in Equation (3)) such that:

$$ ||A^{\bar{K}}(x_0, y_0) - y(x_0)|| \le r'/2 \le r/2. $$

Let $\varphi(x) = f(x, y(x))$. Our next step is to bound $||\nabla_x \varphi^k(x_0, y_0) - \nabla_x \varphi(x_0)||_2$ (note that we cannot use Lemma 3.14 yet because we need $||y_0 - y(x_0)|| \le r(\Pi(S))$ to do so). To this end, we only consider iterations where $A^k(x_0, y_0)$ already entered the band of radius $r(S)$ around the graph of $y$. Let $\rho = \rho(S) \in (0, 1)$, $c = c(S) > 0$ as in Lemma 3.14. For $k \ge \bar{K}$, we have:

$$ ||A^k(x_0, y_0) - y(x_0)|| \le r/2\rho^{k-\bar{K}} \le r/2, $$

by a simple induction. Moreover, also by Lemma 3.14,

$$ ||\nabla_x \varphi^k(x_0, y_0) - \nabla_x \varphi(x_0)|| \le (||\text{Jac}_x A^{\bar{K}}(x_0, y_0) - \text{Jac}_x y(x_0)|| + k - \bar{K})c\rho^{k-\bar{K}}, $$

by our assumption on $k$. Increasing $k$ if necessary, the right hand side can be made smaller than $\epsilon$, because $\text{Jac}_x A^{\bar{K}}(x, z)$ and $\text{Jac}_x y(x)$ are continuous, and thus, bounded on the bounded set $S$.
Finally, define:

$$ \alpha = \min_{i=1,...,N, x \in \Pi(S)} \frac{r}{2(\|\nabla_x f \circ (\mathbf{id}, y^{(i)})(x)\| + \epsilon)}, \quad (11) $$

we have:

$$ \begin{aligned} \mathrm{dist}((x_1, y_1), \mathrm{graph}\ y) &\le \|(x_1, y_1) - (x_0, y(x_0))\| \\ &\le \|x_1 - x_0\| + \|y_1 - y(x_0)\| \\ &\le \alpha \|\nabla_x \varphi^k(x_0, y_0)\| + \|\mathcal{A}^k(x_0, y_0) - y(x_0)\| \\ &\le \frac{r}{2(\|\nabla_x (f \circ (\mathbf{id}, y))(x)\| + \epsilon)} \|\nabla_x \varphi^k(x_0, y_0)\| + \|\mathcal{A}^k(x_0, y_0) - y(x_0)\| \\ &\le r/2 + r/2 = r, \end{aligned} \quad (12) $$

by the choice of $\alpha$ in (11).

<figure><img src="image_4.png" /></figure>

## 4.3 Convergence proof of Algorithm SMBG

We establish, in the spirit of “path following methods”, see e.g. [2] and references therein, that the dynamics evolves near the manifold of local minimizers corresponding to the local minimum selected in Lemma 4.4.

**Lemma 4.5 (Dynamics of remaining iterations).** Assume that Assumption 2.1 and Assumption Morse QC hold, then for all $M > 0, \epsilon > 0$, there exists $r > 0, K \in \mathbb{N}, \alpha > 0$ such that: for any $y^{(*)} = y^{(i)}, i \in \{1, ..., N\}$ ($y^{(i)}$ and $N$ are defined as in Equation (3)), for all $\|x\| \le M, k \ge K, \alpha_f \le \alpha$, if $\mathrm{dist}((x, z), \mathrm{graph}\ y^{(*)}) \le r$, then:

$$ \begin{aligned} & \|\nabla_x \varphi(x) - \nabla_x \varphi^k(x, z)\| \le \epsilon, \\ & \mathrm{dist}((x', z'), \mathrm{graph}\ y^{(*)}) \le r \end{aligned} \quad (13) $$

where $\varphi(x) = f(x, y^*(x))$, and $(x', y')$ is the iterate update of $(x, z)$ performed by Algorithm SMBG, i.e., if $(x, z) = (x_k, y_k)$, then $(x', z') = (x_{k+1}, y_{k+1})$.

*Proof.* Let $S = B(0, M) \subseteq \mathbb{R}^n$ and consider $r = r(S), \rho = \rho(S)$ and $c = c(S)$ as in Lemma 3.14. Using Lemma 3.14, we get:

$$ \begin{aligned} & \|\mathcal{A}^k(x, z) - y^{(*)}(x)\| \le \rho^k \|z - y^{(*)}(x)\|, \\ & \|\nabla_x \varphi^k(x, z) - \nabla_x \varphi(x)\|_2 \le (\|\mathrm{Jac}_x y^{(*)}(x)\| + k) c \rho^k. \end{aligned} $$

where $\varphi(x) = f(x, y^{(*)}(x))$. Thus, for $k$ sufficiently large, we get:

$$ \begin{aligned} & \|\mathcal{A}^k(x, z) - y^{(*)}(x)\| \le r/2, \\ & \|\nabla_x \varphi^k(x, z) - \nabla_x \varphi(x)\|_2 \le \epsilon. \end{aligned} $$

Finally, to show that $\mathrm{dist}(x', y') \le r$, we choose $\alpha$ as in (11) and repeat the same argument as in the proof of Lemma 4.4. □

<figure><img src="image_5.png" /></figure>

In the light of Lemma 4.4, Lemma 4.5, the convergence of Algorithm SMBG boils down to the study of an *inexact gradient descent* methods provided that the iterates remain *bounded* (otherwise, Lemma 4.5 would fail). To establish the convergence results for Algorithm SMBG, we first introduce a condition on $f, g$ ensuring that iterates are uniformly bounded. This condition is also used in [11] to study inexact gradient descent methods. We start with a descent mechanism, the arguments are classical.
**Lemma 4.6** (Quasi-descent lemma). Consider a $C^2$, semi-algebraic, and lower-bounded function $f$ whose crit$_{\epsilon}$f is non-empty and bounded for some $\epsilon > 0$. For any $l > \max_{z \in \text{crit}_\epsilon f} f(z)$, there exist $\alpha > 0$ and $M > 0$, such that for any $\alpha_f$, $0 < \alpha_f \le \alpha$, and any $z \in \mathbb{R}^n$ such that $f(z) \le l$, the sequence defined by:

$$x_0 = z, \quad x_{l+1} = x_l - \alpha_f (\nabla f(x_l) + \xi_l), \|\xi_l\| \le \epsilon, \forall l \in \mathbb{N}, \quad (14)$$

satisfies $f(x_l) \le l$ and $\|x_l\| \le M$, for all $l \in \mathbb{N}$.

*Proof.* [11, Lemma 3.4] ensures that $f$ is coercive. Set $S = \{x, f(x) \le l\}$ which is compact and set $B_0$ a ball of radius $M_0$ containing $S$. Denote by $L_0$ a Lipschitz constants of $\nabla f$ on $B_0$ ($L_0$ exists because $f$ is $C^2$). Set $M = M_0 + \frac{1}{L_0}(\max_{\|z\| \le M_0} \|\nabla f(z)\| + \epsilon)$, denote by $B$ the corresponding ball (with the same center as $B_0$), set $L_f$ and $L_{\nabla f} \ge L_0$ the Lipschitz constants of $f$ and $\nabla f$ on $B$ respectively. Set

$$\alpha = \min \left\{ \frac{1}{L_{\nabla f}}, \frac{l - \max_{z \in \text{crit}_\epsilon f} f(z)}{L_f(L_f + \epsilon)} \right\}$$

Fix arbitrary $x \in S$, $\alpha_f \in [0, \alpha]$ and $\|\xi\| \le \epsilon$, let us show that $y = x - \alpha_f(\nabla f(x) + \xi) \in S$. We have $x \in B_0 \subset B$ and $\alpha_f \|\nabla f(x) + \xi\| \le \frac{1}{L_0}(\max_{\|z\| \le M_0} \|\nabla f(z)\| + \epsilon)$ so that $y = x - \alpha_f(\nabla f(x) + \xi) \in B$. First if $x \notin \text{crit}_\epsilon f$, we can use the descent lemma, since $x, y \in B$,

$$\begin{align*}
f(y) - f(x) &\le \langle \nabla f(x), y - x \rangle + \frac{L_{\nabla f}}{2} \|y - x\|^2 \\
&\le \langle \nabla f(x), y - x \rangle + \frac{1}{2\alpha_f} \|y - x\|^2 \\
&= \frac{1}{2\alpha_f} \left( 2 \langle \alpha_f \nabla f(x), y - x \rangle + \|y - x\|^2 \right) \\
&= \frac{1}{2\alpha_f} \left( \|y - x + \alpha_f \nabla f(x)\|^2 - \|\alpha_f \nabla f(x)\|^2 \right) \\
&= \frac{\alpha_f}{2} \left( \|\xi\|^2 - \|\nabla f(x)\|^2 \right) < 0.
\end{align*}$$

We deduce that $y \in S$. Second, if $x \in \text{crit}_\epsilon$, we have $\|y - x\| = \alpha_f \|\nabla f(x) + \xi\| \le \alpha_f(L_f + \epsilon)$, and

$$f(y) \le f(x) + L_f \|y - x\| \le f(x) + \alpha_f L_f (L_f + \epsilon) \le f(x) + l - \max_{z \in \text{crit}_\epsilon f} f(z) \le l$$

and $y \in S$. This proves the desired result. $\square$

Combining Lemmas 4.4 to 4.6, we get the following result:

**Theorem 4.7** (Inexact gradient descent equivalence of Algorithm SMBG). Under Assumption 2.1 and Assumption Morse QC, assume that $\inf_{x,y} f(x,y) - \lambda \|x\|^2 > -\infty$ for some $\lambda > 0$. For any compact set $C \subset \mathbb{R}^n \times \mathbb{R}^m$, any $\epsilon, \Delta > 0$, there exists a set $A \subset \mathbb{R}^n \times \mathbb{R}^m$ with Lebesgue measure at most $\Delta$, $K \in \mathbb{N}$, $\alpha > 0$, $M > 0$ such that if $x_0 \in C \setminus A$, $k \ge K$, $\alpha_f \le \alpha$:

1. The sequence $\{(x_l, y_l) | l \in \mathbb{N}\}$ produced by Algorithm SMBG is bounded in $B(0, M) \subset \mathbb{R}^n \times \mathbb{R}^m$.

2. There is $y = y^{(i)}$, $i \in \{1, \dots, N\}$ (c.f. Proposition 3.6) such that, $\|\nabla_x \varphi(x_l) - \nabla_x \varphi^k(x_l, y_l)\| \le \epsilon, \forall l \in \mathbb{N}$, where $\varphi(x) = f(x, y(x))$.

*Proof.* Since $f$ and $y^{(i)}$, $i = 1, \dots, N$ are semi-algebraic (cf. Proposition 3.6), $f_i: x \mapsto f(x, y^{(i)}(x))$ is semi-algebraic. Moreover, each function $f_i - \lambda \| \cdot \|^2$, $i = 1, \dots, N$ is trivially lower-bounded and each set crit$_{\epsilon}$f$_i$, $i = 1, \dots, N$ is bounded as proved in [11, Proposition 3.1]. The critical values vcrit$_{\epsilon}$f$_i$ are also bounded for $i = 1, \dots, N$ by continuity.
Choose $l > 0$ such that

$$l > \max\{\text{vcrit}_\epsilon f \circ (\text{id}, y^{(i)})\} \qquad l \ge \max_{x,y \in C} \{f(x, y^{(i)}(x))\}, \forall i = 1, \dots, N.$$

We choose $M > 0, \alpha_1 > 0$ as the maximal values given by Lemma 4.6 applied to each $f_i$, $i = 1, \dots, N$. Note that it follows from the statement of Lemma 4.6 that the ball of radius $M$ should contain the sublevel sets of value $l$ for each $f_i$, $i = 1, \dots, N$.

With $\epsilon$ and $M$ thus given, we choose $(r, K_1, \alpha_2)$ equal to the quantities $(r, K, \alpha)$ of Lemma 4.5. In particular, if $\|x_\ell\| \le M$, $k \ge K_1$, $\alpha_f \le \alpha_2$ and $\text{dist}((x_\ell, y_\ell), \text{graph } y^{(i)}) \le r$, we have:

1. $x_{\ell+1}$ is updated by (14) applied to $f_i: x \mapsto f(x, y^{(i)}(x))$.

2. $\text{dist}((x_{\ell+1}, y_{\ell+1}), \text{graph } y^{(i)}) \le r$.

for any index $i \in \{1, \dots, N\}$.

Next, given $C, \Delta > 0, r > 0, \epsilon > 0$, we obtain a subset $A$ of Lebesgue measure at most $\Delta$, $K_2 \in \mathbb{N}, \alpha_3 > 0$ equal to $K$ and $\alpha$ defined as in Lemma 4.4, respectively and a $C^2$ solution $y$ such that for any $(x_0, y_0) \in C \setminus A$, we have for any $k \ge K_2$:

$$\|\nabla_x \varphi^k(x_0, y_0) - \nabla_x \varphi\| \le \epsilon, \qquad (15) \\ \text{dist}((x_1, y_1), \text{graph } y) \le r.$$

where $\varphi(x) := f(x, y(x))$ and $y = y^{(i)}$ for some $i \in \{1, \dots, N\}$.

Finally, we construct $\alpha = \min(\alpha_1, \alpha_2, \alpha_3)$, $K = \max(K_1, K_2)$. The main claim will be proved by induction: for all initializations $(x_0, y_0) \in C \setminus A$, for all iterations $\ell \in \mathbb{N}$, we have:

1. $\|x_\ell\| \le M$.

2. $x_\ell$ is updated to $x_{\ell+1}$ by (14) applied to $\varphi$ given in (15) with the function $y$.

3. $\text{dist}((x_{\ell+1}, y_{\ell+1}), \text{graph } y) \le r$.

* **Base case:** by the construction of $K \ge K_2$, we have (15) satisfied and $f(x_0, y(x_0)) \le l$ (hence, $\|x_0\| \le M$).

* **Induction:** Since $x_{\ell-1}, \ell \ge 1$ is updated by (14) applied to $\varphi$, by the construction of $\alpha \le \alpha_1$, we have $\|x_\ell\| \le M$. Combining with the induction hypothesis $\text{dist}((x_\ell, y_\ell), \text{graph } y) \le r$, using Lemma 4.5 and the construction of $\alpha_f \le \alpha_2$, $K \ge K_1$, $x_{\ell+1}$ is updated by (14) applied to $\varphi$ and $\text{dist}((x_{\ell+1}, y_{\ell+1}), \text{graph } y) \le r$. $\square$

We now turn to the convergence of Algorithm SMBG. We recall the main convergence result for the inexact gradient descent given in [11]:

**Theorem 4.8 (Convergence for inexact gradient method with constant step size [11, Theorem 2]).** Consider a function $f$ L-Lipschitz, lower-bounded, semi-algebraic with $\text{crit}_\epsilon$ bounded for some $\epsilon > 0$, there is $C > 0, \rho > 0$ such that for any $x_0 \in \mathbb{R}^p$ and $x_\ell(\alpha)$ generated by inexact gradient descent method with error bounded by $\epsilon$ and fixed step size $\alpha$, we have:

$$\lim_{\alpha \to 0^+} \limsup_{\ell \to \infty} \text{dist}(f(x_\ell(\alpha)), \text{vcrit}_\epsilon f) = 0 \\ \limsup_{\alpha \to 0^+} \limsup_{\ell \to \infty} \text{dist}(x_\ell(\alpha), \text{crit} f) \le C\epsilon^\rho.$$

**Remark 4.9.** One needs to argue that in [11, Lemma 3.2], the quantity $\bar{\epsilon}$ can be chosen to be the same as in [11, Assumption 1]. This is indeed the case, by coercivity, all arguments take place in a compact set, KL inequality is trivial for non-critical point, and the definable metric subregularity result used is given for an arbitrary compact.
Combining Theorem 4.7 and Theorem 4.8, we obtain the proof for Theorem 4.2.

*Proof of Theorem 4.2.* Call $P$ the initial distribution over $(x_0, y_0)$ which is assumed to be absolutely continuous. Let $C \subset \mathbb{R}^n \times \mathbb{R}^m$ be compact such that $P((x, y) \in C) \ge 1 - \frac{\Delta}{2}$. By absolute continuity of the initialization, there is $\Delta' > 0$ such that $P(A) \le \frac{\Delta}{2}$ for any $A$ of Lebesgue measure at most $\Delta'$. We may now consider $A, K, \alpha, M$ as given by Theorem 4.7 (with measure threshold $\Delta'$), so that with probability $1 - \Delta$ over the random choice of $(x_0, y_0)$, the result of Theorem 4.8 applies.

In the light of Theorem 4.7, we only need to verify the assumptions of Theorem 4.8. Due to the nature of the main claim, we can always assume that $\alpha_f$ is sufficiently small. For each $i = 1, \dots, N$, set $f_i: x \mapsto f(x, y^{(i)}(x))$ as in the proof of Theorem 4.8

1. $f_i, i = 1, \dots, N$ are semi-algebraic: see proof of Theorem 4.7.

2. $f_i, i = 1, \dots, N$ are lower-bounded, with bounded $\text{crit}_e$: see proof of Theorem 4.7.

3. By Theorem 4.7, and with sufficiently small $\alpha_f$, with high probability, the sequence $x_\ell(\alpha_f)$ are uniformly bounded in $B(0, M)$. Since $f_i$ is continuously differentiable, it is Lipchitz when restricted to a compact domain. Thus, with high probability, we can assume that $f_i$ is L-Lipschitz, for all $i \in \{1, \dots, N\}$.

By applying Theorem 4.8, we conclude the proof. $\square$

# 5 Pros and cons of the differentiable programming strategy

In this section, we study differentiable programming strategy for bilevel gradient method (i.e., Algorithm DPBG), where one minimizes $f(x, A^k(x, z))$ as a surrogate for (BL). An investigation of the behavior of this algorithm is of practical interest: for example, standard gradient algorithms to train MAML [22] for meta learning can be interpreted as Algorithm DPBG applied to the MAML bilevel formulation [54]. In addition, replacing the constraint in (BL) by a fixed number of GD steps provides a scalable method to deal with (BL), and was employed and studied in the bilevel literature (e.g., [43, Algorithm 1]).

Certain results in this section are announced with the “genericity” notion, similar to Proposition 3.4. In this section, genericity has the following interpretation. Given $\mathcal{F}$, a class of functions from $\mathbb{R}^n$ to $\mathbb{R}$, if a property $P$ is generic in $\mathcal{F}$, then for any $f$, the set $\{a \in \mathbb{R}, g(\cdot) := f(\cdot) + \langle a, \cdot \rangle \text{ satisfies } P\}$ is open, dense and of full Lebesgue measure. The term “for generic $f$” can therefore be understood as “for any $f$ up to a typical linear perturbation”, see for instance Lemma 5.4.

The purpose of this section is to compare the optimization landscape of bilevel optimization problem (BL) and the minimization problem

$$ \min \varphi^k(x, z) := f(x, A^k(x, z)) \text{ for } k \in \mathbb{N}. $$

As we cannot formally take $k = \infty$, this approach is connected to the unconstrained version of problem (BL) which we refer to as the single-level problem

$$ \min_{x \in \mathbb{R}^n, y \in \mathbb{R}^m} f(x, y). \qquad (\text{SL}) $$

## 5.1 Equivalence with the unconstrained single level problem

A fundamental property of the gradient descent map is that it defines a *global diffeomorphism* of the space. The following result is well known in optimization, we provide a proof for completeness.
**Lemma 5.1** (Gradient descent defines a global diffeomorphism). Under Assumption 2.1 (except semi-algebraicity of f and g), the mapping $(x, z) \mapsto (x, A^k(x, z))$ is a global diffeomorphism and the inverse is given by k applications of

$$A^{-1}(x, \cdot) = \text{prox}_{-\alpha g_x}$$

which is well-defined and single valued.

*Proof.* For $k = 1$, and for any $x \in \mathbb{R}^n$, the Jacobian of $A(x, \cdot)$ at $z \in \mathbb{R}^m$ is given by

$$Jac_z A(x, z) = \mathbf{I} - \alpha_g \nabla_{yy}^2 g(x, z)$$

and by assumption $\alpha_g \|\nabla_{yy}^2 g(x, z)\|_{op} < 1$, so that this matrix is invertible and $A(x, \cdot)$ is a local diffeomorphism. It is actually a global diffeomorphism since

$$\begin{align*}
A(x, y) = A(x, z) &\implies \|y - z\| = \|\alpha_g(\nabla_y g(x, z) - \nabla_y g(x, y))\| \le \alpha_g L_g \|y - z\| \\
&\implies \|y - z\|(1 - \alpha_g L_g) \le 0 \\
&\implies \|y - z\| = 0,
\end{align*}$$

so that it is injective. Regarding the prox characterization, we have

$$0 = z - u - \alpha_g \nabla_y g(x, z) \iff z = \operatorname{argmin}_v \frac{1}{2} \|v - u\|^2 - \alpha_g g(x, v).$$

The result for general k follows by composition. □

**Proposition 5.2** (Equivalence with the single level problem). Under Assumption 2.1 (except semi-algebraicity of f and g), for any $k \in \mathbb{N}$, up to a diffeomorphism, the critical points of f are the same as those of $\varphi^k$, with the same objective value and the same index (local minimum/maximum, second order structure).

*Proof.* Since $\varphi^k(x, z) = f(x, A^k(x, z))$ the equivalence follows from Lemma 5.1. We have

$$\nabla \varphi_k(x, z) = V_k^\top \nabla f(x, A^k(x, z))$$

where

$$V_k = \begin{pmatrix} I & 0 \\ Jac_x A^k(x, z) & Jac_z A^k(x, z) \end{pmatrix},$$

and the left matrix is invertible so that $\nabla \varphi_k(x, z) = 0$ if, and only if, $\nabla f(x, A^k(x, z)) = 0$. The local minimum / maximum structure is preserved by diffeomorphisms. As for the second order structure, for any critical point of $\varphi^k$, we have

$$\nabla^2 \varphi^k(x, z) = V_k^\top \nabla^2 f(x, A^k(x, z)) V_k,$$

because $\nabla f(x, A^k(x, z)) = 0$. Therefore, the Hessian of $\varphi^k$ is congruent to that of f and they have therefore the same number of positive, negative, and null eigenvalues. □

## 5.2 Pseudo-stability of the differentiable programming strategy

Proposition 5.2 is a rather negative result as it entails that there is no qualitative difference between $\varphi^k$ and f in terms of their critical points. In other words, the approximate value function $\varphi^k$ ignores the bilevel constraint. On the other hand, the goal of this section is to establish that near a local minimizer the algorithm has a pseudo-stability property: if the iterates meet a solution neighborhood, they tend to remain close for a long time. Although
we describe an exponentially long stability duration, it is different from genuine stability as we empirically illustrate in Section 5.4, hence the term “pseudo-stability”.

Let us recall that, given a $C^2$ function $h: \mathbb{R}^n \to \mathbb{R}$, the point $x^* \in \mathbb{R}^m$ is a strong local minimizer of $h$, if

$$ \nabla h(x^*) = 0, \text{ and } \nabla^2 h(x^*) > 0. $$

Assume $y^*$ is a strong local minimizer for $g_{x^*}$ for some couple $(x^*, y^*)$. Observe that, by the implicit function theorem, there is a local map $\bar{y}$ such that $\bar{y}(x)$ is a strict local minimum of $g_x$ in a neighborhood of $x^*$.

When $y$ is a $C^2$ curve which satisfies $y(x) \in \text{crit}g_x$, and $x^*$ is strong local minimizer of $\varphi: x \mapsto f(x, y(x))$, the pair $(x^*, y^*)$ is called a strong local minimizer of problem (BL). In particular, if one assumes Morse QC, this describes a strong local minimizer of $f(x, y^{(i)}(x))$ with $i \in \{1, \dots, N\}$ (the quantities $N$, $y^{(i)}$ being defined as in Proposition 3.6).

**Theorem 5.3** (Pseudo-stability around strong local minimizers). Let Assumption 2.1 (except semi-algebraicity of $f$ and $g$) hold and let $(x^*, z^*)$ be a strong local minimum of (BL).

Then, there exists $\delta > 0$ and $\gamma > 0$, as well as constants $C > 0$ and $0 < \rho < 1$, such that for all $0 < \alpha_f \le \gamma$, for all $k \in \mathbb{N}$, Algorithm DPBG, initialized at $(x, z) \in \mathbb{R}^n \times \mathbb{R}^m$ such that $\|x - x^*\| + \|z - z^*\| < \delta$ satisfies $\|x_\ell - x^*\| + \|z_\ell - z^*\| \le \delta$, for all

$$ \ell \le \frac{\delta - \|x - x^*\| - \|z - z^*\|}{\alpha_f C \rho^k}. $$

*Proof of Theorem 5.3.* Let $U$, $\rho$ and $C$ be given by Lemma 3.11. Since $g$ is $C^3$, the implicit function $y$ of the function $\nabla_y f(x, y) = 0$ is $C^2$. As a consequence, $\varphi(x) := f(x, y(x))$ is $C^2$. We may pick $\gamma > 0$ such that $\gamma \|\nabla^2\varphi(x)\| < 2$ uniformly on $U$. Furthermore, since $x^*$ is a strict second order minimizer of $\varphi$ we may pick $\delta > 0$ such that $\nabla^2\varphi(x) \ge 0$ for all $x$ such that $\|x - x^*\| \le \delta$, and such that $B(x^*, \delta) \times B(z^*, \delta) \subset U$.

The function $\varphi$ is strongly convex on $B$ with a global minimizer at $x^*$. This implies that the gradient mapping $x \mapsto x - \alpha \nabla\varphi(x)$ is a contraction on $B$ with global fixed point at $x^*$. From the update of Algorithm DPBG, we obtain for all $\ell \ge 1$,

$$ \|z_{\ell} - z^{*}\| = \|z_{\ell-1} - z^{*} - \alpha_{f} \nabla_{z} \varphi^{k}(x_{\ell-1}, z_{\ell-1})\| \le \|z_{\ell-1} - z^{*}\| + \alpha_{f} C \rho^{k} \le \|z - z^{*}\| + \ell \alpha_{f} C \rho^{k}, $$

since $\|\nabla_z \varphi^k(x_{\ell-1}, z_{\ell-1})\| \le C\rho^k$ for some constants $C, \rho > 0$ in the neighborhood $U$ (cf. Lemma 3.11). Moreover,

$$ \begin{aligned} \|x_{\ell} - x^{*}\| &= \|x_{\ell-1} - x^{*} - \alpha_{f} \nabla_{z} \varphi^{k}(x_{\ell-1}, z_{\ell-1})\| \\ &\le \|x_{\ell-1} - x^{*} - \alpha_{f} \nabla_{z} \varphi(x_{\ell-1})\| + \alpha_{f} \|\nabla_{x} \varphi(x_{\ell-1}) - \nabla_{x} \varphi^{k}(x_{\ell-1}, z_{\ell-1})\| \\ &\le \|x_{\ell-1} - x^{*}\| + \alpha_{f} C \rho^{k} \\ &\le \|x - x^{*}\| + \ell \alpha_{f} C \rho^{k}, \end{aligned} $$

again, by using the estimation $\|\nabla\varphi(x_{\ell-1}) - \nabla_x\varphi^k(x_{\ell-1}, z_{\ell-1})\| \le C\rho^k$ (cf. Lemma 3.11). We conclude that if $\|z_\ell - z^*\| + \|x_\ell - x^*\| > \delta$, one has

$$ \|x - x^*\| + \|z - z^*\| + 2\ell\alpha_f C \rho^k \ge \delta $$

and the result follows (ignoring the constant 2).

<figure><img src="image_6.png" /></figure>

The following is a consequence of the Morse-Sard theorem (see proof in [29, Chapter 1.7]).
**Lemma 5.4** (Linear perturbation of a semialgebraic function). Let $f, g$ be as in Assumption 2.1 (except semi-algebraicity of $f$ and $g$) and let $\bar{y}: \mathbb{R}^n \to \mathbb{R}^m$ be $C^2$. Generically, for $a \in \mathbb{R}^n$, the function $x \mapsto f(x, \bar{y}(x)) + \langle x, a \rangle$ is Morse.

Combining Lemma 5.4 with Theorem 5.3, we obtain the following corollary.

**Corollary 5.5.** *Under Assumption 2.1, for a generic $f$ and a parametric Morse $g$, the result of Theorem 5.3 holds for all local minimizers of the bilevel problem (BL).*

The consequence of corollary 5.5 is that there are neighborhoods of local minimizers of the original bilevel program (BL) which are pseudo-stable for the differentiable programming approximation, minimized by Algorithm DPBG, for large $k$. This means that for very large $k$ if the recursion visits this neighborhood, it will not necessarily remain in it for an infinitely long time, but it will stay in the neighborhood for a time which is exponential in $k$. In particular, the iterates will eventually leave this neighborhood unless this region contains critical points of $\varphi^k$ (indeed, since $\alpha \le \gamma < 2/L$ where $L$ is the Lipchitz constant of gradient in the neighborhood), we are in the small step-size regime where the objective function decreases monotonically and any accumulation point of the iterates is a critical point). To appreciate this result, it is important to remark that in Theorem 5.3, the quantities $\gamma$ which bounds step size $\alpha$ which defines the outer recursion, the radius $\delta$, the constants $C$ and $\rho$ do not depend on $k$.

## 5.3 Repulsivity of unconstrained critical points

The main result of this section states that under favorable assumptions the critical points of $\varphi^k$ can actually be reached under specific and unlikely circumstances.

**Theorem 5.6** (Escape at infinity and sharpness). Under Assumption 2.1, consider $(x^*, y^*) \in \mathbb{R}^n \times \mathbb{R}^m$ satisfying $\nabla f(x^*, y^*) = 0$ and such that $\nabla^2 f(x^*, y^*)$ invertible. Assume $y^*$ is not a local minimum of $g_{x^*}$ and that $g_{x^*}$ is a Morse function. Consider for each $k \in \mathbb{N}$ the unique $y^{-k} \in \mathbb{R}^m$ such that $A^k(x^*, y^{-k}) = y^*$ as given in Lemma 5.1. We have the following alternatives:

1. $\lim_{k \to \infty} \|y^{-k}\| = +\infty$.

2. There exists $C > 0$, $\rho > 1$ and $k_0 \in \mathbb{N}$, such that: $\|\nabla^2 \varphi^k(x^*, y^{-k})\|_{\text{op}} \ge C\rho^{2k}$, for all $k \ge k_0$.

Before proving Theorem 5.6, we explain its implication: it is difficult to converge to the critical points of $\varphi^k$ that are not related to the original bilevel problem. Indeed, we remind the readers that $(x^*, y^*)$ is a critical point of $f$ if and only if $(x^*, y^{-k})$ is a critical point of $\varphi^k$. If $y^*$ is not a local minimum of $g_{x^*}$ (i.e., $(x^*, y^{-k})$ is a “fake” critical point of $\varphi^k$), then Algorithm DPBG (which is equivalent to gradient descent on $\varphi^k$ with step-size $\alpha_f$) either:

1. Takes a very long time to reach these “fake” critical points since $y^{-k}$ diverges to infinity when $k \to \infty$ (the first possibility in Theorem 5.6).

2. Converges to a critical point with exponentially large curvature (w.r.t. $k$) (the second possibility in Theorem 5.6). This requires either an exponentially small step size (i.e., $\alpha_f = o(1/\rho^{2k})$ for some $\rho > 1$), or constitutes a very unlikely event since in this case the Jacobian is not a contraction in this case (see e.g. [58] or more precisely [9, Theorem 2.1 or Proposition 3.1]). Given the usual range of learning rates $\alpha_f$ in practical settings such as deep learning, Algorithm DPBG is unlikely to converge to these points if $k$ is sufficiently large.

We will need the following lemma on operator norms.

**Lemma 5.7.** Given two matrices $A \in \mathbb{R}^{n \times m}$, $B \in \mathbb{R}^{n \times n}$, with $B$ symmetric. We have

$$ \|A^T B A\|_{\text{op}} \ge \|A\|_{\text{op}}^2 \lambda_{\min}(B). $$

where $\lambda_{\min}(B) > 0$ is the smallest absolute value of eigenvalues of $B$.
Proof. We have for any $x$, $\|x^\top A^\top B A x\| \ge \lambda_{\min}(B) \|Ax\|_2^2$, so that

$$ \|A^\top BA\|_{\text{op}} = \max_{\|x\|=1} \|x^\top A^\top B A x\| \ge \max_{\|x\|=1} \lambda_{\min}(B) \|Ax\|_2^2 = \|A\|_{\text{op}}^2 \lambda_{\min}(B). $$

□

Proof of Theorem 5.6. Thanks to Lemma 5.1, one can interpret the sequence ${y^{-k}}_{k \in \mathbb{N}}$ as a proximal point sequence to minimize $-g(x^*, \cdot)$, that is, to maximize $g(x^*, \cdot)$, initialized at $y^*$. By [6], there are exactly two alternatives:

1. $\lim_{k \to \infty} \|y^{-k}\| = +\infty$.

2. $\lim_{k \to \infty} y^{-k} = \bar{y}$ where $\bar{y}$ is a critical point of $g(x^*, \cdot)$.

Let us consider the second case. The point $\bar{y}$ is not a local minimizer of $g(x^*, \cdot)$. Indeed, if $\bar{y}$ is a local minimizer of $g(x^*, \cdot)$, since $-\alpha_g g(x^*, \cdot) + \frac{1}{2}\|y^* - \cdot\|^2$ attains its global minimum at $y^*$ (because $\alpha_g L_g < 1$), the proximal iterates are therefore stationary so that $y^* = \bar{y}$ is also a local minimizer of $g(x^*, \cdot)$. This last case is ruled out by our hypotheses.

We recall that from the proof of Proposition 5.2, since $\nabla f(x^*, y^*) = 0$, we have

$$ \nabla^2 \varphi^k(x^*, y^{-k}) = V_{-k}^\top \nabla^2 f(x^*, y^*) V_{-k} $$

where

$$ V_{-k} = \begin{pmatrix} I & 0 \\ \text{Jac}_x A^k(x^*, y^{-k}) & \text{Jac}_z A^k(x^*, y^{-k}) \end{pmatrix}. $$

We are going to show that, for large $k$, $\text{Jac}_z A^k(x^*, y^{-k})$ has an exponentially large operator norm and the result will follow from Lemma 5.7. We recall from the proof of Lemma 3.11 the following

$$ \partial_z A^k(x^*, y^{-k}) = \prod_{i=1}^{k} \underbrace{(\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, y^{-i}))}_{\mathbf{H}_i}, $$

Let $\mathbf{H}_{\alpha,b} = \prod_{i=a}^b \mathbf{H}_i$, it is sufficient to prove that $\|\mathbf{H}_{k_0,k}\| \ge K \rho^{k-k_0}$ for some $\rho > 1, K > 0$ and $k_0 \in \mathbb{N}$. It is true because $\mathbf{H}_i$ are symmetric, positive definite with $\lambda_{\min}(\mathbf{H}_i) \ge 1 - \alpha_g L_g > 0$ (cf. Assumption 2.1) for all $i \in \mathbb{N}$.

Since $\partial_{yy}^2 g(x^*, \bar{y})$ is symmetric, there exists an orthogonal matrix $\mathbf{P}$ such that $\mathbf{P}^\top \partial_{yy}^2 g(x^*, \bar{y}) \mathbf{P}$ is diagonal. We conclude that the same matrix $\mathbf{P}$ also makes:

$$ \mathbf{D} = \mathbf{P}^\top (\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, \bar{y})) \mathbf{P} $$

diagonal. Moreover, due to the Morse assumption of $g(x^*, \cdot)$ and its Lipchitz gradient (cf. Assumption 2.1), the matrix $\mathbf{D}$ contains two sets of eigenvalues: $\{\lambda | \lambda > 1\}$ and $\{\lambda | 0 < \lambda < 1\}$. Without loss of generality, we can suppose that:

$$ \mathbf{D} = \begin{pmatrix} \mathbf{D}_1 & 0 \\ 0 & \mathbf{D}_2 \end{pmatrix}, $$

where $\mathbf{D}_1$ (resp. $\mathbf{D}_2$) contains only positive elements that are bigger (resp. smaller) than 1. Note that the size of $\mathbf{D}_1$ is at least one because $\partial_{yy}^2 g(x^*, \bar{y})$ has at least one negative eigenvalue.

Since $\lim_{k \to \infty} y^{-k} = \bar{y}$ and $g$ is $C^2$ (cf. Assumption 2.1), for any $\epsilon > 0$, one can choose $k_0 \in \mathbb{N}$ such that for $\ell \ge k_0$, one can write: $\partial_{yy}^2 g(x^*, y^{-\ell}) = \partial_{yy}^2 g(x^*, \bar{y}) + \mathbf{E}_\ell$ where $\|\mathbf{E}_\ell\| < \epsilon/\alpha_g$. Thus,

$$ \mathbf{P}^\top (\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, y^{-\ell})) \mathbf{P} = \mathbf{P}^\top (\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, \bar{y}) - \alpha_g \mathbf{E}_\ell) \mathbf{P} = \mathbf{D} - \underbrace{\alpha_g \mathbf{P}^\top \mathbf{E}_\ell \mathbf{P}}_{\mathbf{E}'_\ell}. $$
Since the matrix $\mathbf{P}$ is orthogonal, $\|\mathbf{E}'_\ell\| = \alpha_g \|\mathbf{E}_\ell\| \le \epsilon$. Thus,

$$
\begin{aligned}
\|\mathbf{H}_{k_0,k}\| &= \left\| \prod_{\ell=k_0}^{k} (\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, y^{-\ell})) \right\| \\
&= \left\| \mathbf{P}^{\top} \left[ \prod_{\ell=k_0}^{k} (\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, y^{-\ell})) \right] \mathbf{P} \right\| \\
&= \left\| \prod_{\ell=k_0}^{k} \mathbf{P}^{\top} (\mathbf{I} - \alpha_g \partial_{yy}^2 g(x^*, y^{-\ell})) \mathbf{P} \right\| \\
&= \left\| \prod_{\ell=k_0}^{k} (\mathbf{D} + \mathbf{E}'_\ell) \right\|
\end{aligned}
$$

In the following, WLOG, we will assume $k_0 = 1$ and prove $\|\mathbf{H}_{k_0,k}\| = \|\mathbf{H}\| \ge K\rho^k$.

Consider the product between $\mathbf{H}$ and $e_1$ - the first canonical basis in $\mathbb{R}^n$. Define the sequence of pairs of vectors $\{(u_n, v_n)\}$ whose sizes are equal to the number of elements bigger than and smaller than 1 of $\mathbf{D}$ respectively as follows:

$$ \begin{pmatrix} u_0 \\ v_0 \end{pmatrix} = e_1, \quad \begin{pmatrix} u_{\ell+1} \\ v_{\ell+1} \end{pmatrix} = (\mathbf{D} + \mathbf{E}'_{\ell}) \begin{pmatrix} u_{\ell} \\ v_{\ell} \end{pmatrix}, \forall \ell \ge 0. $$

In particular, we have:

$$ \mathbf{H}e_1 = \begin{pmatrix} u_k \\ v_k \end{pmatrix} $$

Moreover, if one write $\mathbf{E}'_\ell$ as a block matrix:

$$ \mathbf{E}'_{\ell} = \begin{pmatrix} \mathbf{E}_{\ell}^{1,1} & \mathbf{E}_{\ell}^{1,2} \\ \mathbf{E}_{\ell}^{2,1} & \mathbf{E}_{\ell}^{2,2} \end{pmatrix} $$

we get the following relation:

$$ u_{\ell+1} = (\mathbf{D}_1 + \mathbf{E}_{\ell}^{1,1})u_{\ell} + \mathbf{E}_{\ell}^{1,2}v_{\ell} $$

$$ v_{\ell+1} = \mathbf{E}_{\ell}^{2,1}u_{\ell} + (\mathbf{D}_2 + \mathbf{E}_{\ell}^{2,2})v_{\ell} $$

Denote $\tau_1 > 1 > \tau_2$ the smallest (resp. the largest) non-zero elements of $\mathbf{D}_1$ and $\mathbf{D}_2$, we obtain an estimation on the norm of $u_\ell, v_\ell$, using the fact that $\|\mathbf{E}_\ell^{i,j}\| \le \|\mathbf{E}'_\ell\| \le \epsilon, \forall i, j \in \{1, 2\}$:

$$ \|u_{\ell+1}\| \ge (\tau_1 - \epsilon)\|u_\ell\| - \epsilon\|v_\ell\| $$

$$ \|v_{\ell+1}\| \le \epsilon\|u_\ell\| + (\tau_2 + \epsilon)\|v_\ell\| $$

Therefore, if $\epsilon$ is small enough so that $\tau_1 - \epsilon > 1 > \tau_2 + \epsilon$, we get:

$$
\begin{aligned}
\|u_{\ell+1}\| &\ge (\tau_1 - \epsilon)\|u_\ell\| - \epsilon\|v_\ell\| \\
&\ge (\tau_1 - \epsilon)\|u_\ell\| - \epsilon^2\|u_{\ell-1}\| - \epsilon(\epsilon + \tau_2)\|v_{\ell-1}\| \\
&\ge (\tau_1 - \epsilon)\|u_\ell\| - \epsilon^2\|u_{\ell-1}\| - \epsilon^2(\epsilon + \tau_2)\|u_{\ell-2}\| - \epsilon(\epsilon + \tau_2)^2\|v_{\ell-2}\| \\
&\dots
\end{aligned}
$$

$$
\begin{aligned}
&\ge (\tau_1 - \epsilon)\|u_\ell\| - \epsilon^2 \left( \sum_{j=0}^{\ell-1} (\tau_2 + \epsilon)^{\ell-1-j} \|u_j\| \right) \\
&\ge (\tau_1 - \epsilon)\|u_\ell\| - \frac{\epsilon^2}{1 - \tau_2 - \epsilon} \max_{j=0, \dots, \ell-1} \|u_j\|
\end{aligned}
$$
For $\epsilon$ sufficient small so that $\rho := \tau_1 - \epsilon - \epsilon^2/(1 - \tau_1 - \epsilon) > 1$, we get $\|u_\ell\| \ge \rho^\ell, \forall \ell \le 0$. Hence, $\|\mathbf{H}_{k_0,k}\| = \|\mathbf{H}\| \ge \|u_k\| \ge \rho^k$, and the result follows. $\square$

**Corollary 5.8.** Under Assumption 2.1, for a generic $f$ and $g$, the result of Theorem 5.6 holds for all critical points of the single level problem (SL).

We finish this section with another corollary, which is an application of Theorem 5.6 when we have in addition the parametric Morse property (cf. Definition 3.1) for $g$.

**Corollary 5.9.** Consider bilevel optimization where Assumption 2.1 holds, $g$ is parametric Morse and $f$ is a generic function. We have: for any bounded set $K \subseteq \mathbb{R}^n \times \mathbb{R}^m$, there exists $k_0 \in \mathbb{N}$ such that for any $k \ge k_0$, there are three possibilities regards a critical point $(x, y)$ of $\varphi^k$:

1. $A^k(x, y)$ is a local minimum of $g(x, \cdot)$.

2. $(x, y) \notin K$.

3. $\|\nabla^2\varphi^k(x, y)\| \ge C\rho^{2k}$ for some constants $C > 0, \rho > 1$. In addition, there exists a neighborhood $U$ of $(x, y)$ where $\|\nabla^2\varphi^k(\tilde{x}, \tilde{y})\| \ge C\rho^{2k}/2$ for all $(\tilde{x}, \tilde{y}) \in U$.

*Proof.* This is just an application of Theorem 5.6. Indeed, the fact that $g$ is parametric Morse ensures the condition $g_x$ being Morse as in Theorem 5.6. Three possibilities then follow correspondingly. Lastly, the neighborhood $U$ in the last claim exists because $\varphi^k$ is a $C^2$ function (remember that $g$ and $f$ are $C^2$ and $C^3$ respectively as assumed in Assumption 2.1). $\square$

## 5.4 Illustrations of theoretical findings

We illustrate the findings of this section in Figures 2 and 3 with a simple example. Let $h: \mathbb{R} \to \mathbb{R}$ denotes the Huber loss, $h(t) = t^2/2$ for $t \in [-1, 1]$ and $h(t) = |t| - 1/2$ for $|t| > 1$. This function is $C^1$ and has vanishing curvature outside $[-1, 1]$ which allows us to compose it with a square while preserving the global Lipschitz-continuity of the gradient. The following examples do not contradict [23, Theorem 3.2] which requires both uniform convergence of the lower level recursion to a unique argmin, and a bounded set of initialization.

**Illustrating instabilities:** We consider the bilevel problem as in (BL) with

$$ \min_{x \in \mathbb{R}^n, y \in \mathbb{R}} f(x, y) := (y - 0.5)^2 \quad \text{s.t.} \quad y \in \operatorname*{argmin}_z g(x, z) := h(z^2 - 1) \quad (16) $$

where the two functions $f(x, y) = (y - 0.5)^2$ and $g(x, y) = h(y^2 - 1)$, are displayed in Figure 2. This is a degenerate problem as there is no dependency in $x$, it corresponds to the following trivial problem, yet, it illustrates our results very well.

$$ \min_{y \in \mathbb{R}} (y - 0.5)^2 \quad \text{s.t.} \quad y = \pm 1. \quad (17) $$

The inner gradient descent recursion in (DP-BL) is attracted by $\pm 1$, with a smooth dependency on the initialization around $\pm 1$, and a sharp transition around 0, the local maximizer of $g$. The higher the value of $k$, the steeper this transition. As a consequence, when composing $f$ and $A^k$, the landscape of $f$ is completely modified. It has two very flat regions where the derivative almost vanishes in a large neighborhood around $\pm 1$, and a sharp global minimum, corresponding to the global minimum of $f$. These are illustrated in Figure 2 where we represent $A^k$ and $f \circ A^k$ for $k = 9$. The larger $k$, the flatter the flat areas, and the sharper the global minimizer.
<figure><img src="image_7.png" /><figcaption>**Figure 2:** Illustration of the instability phenomenon of the diagonal method. Left: inner objective *g* and outer objective *f*, with the outcome of *k* = 9 gradient steps on *g*, depending on the initialization *A*^k(*z*) and the corresponding approximation *φ*^k(*z*) = *f*(*A*^k(*z*)). Middle: value of *y*<sub>*l*</sub> *φ*^k(*y*<sub>*l*</sub>) along outer iterations. Many iterations have objective values corresponding to solutions of the bilevel problem (BL) but tends to be attracted by the sharp global minimizer of *f*, which is repulsive for Algorithm DPBG. Right: profile of *φ*^k = *f* ◦ *A*^k for different values of *k*.</figcaption></figure>

Now running the gradient recursion on $f \circ A^k$ as in Algorithm DPBG, we obtain the following
behavior, illustrated in Figure 2, middle. The recursion spends most of its iterations around
the large flat areas. It does not stay there however since they do not contain stationary points.
it is attracted by the sharp minimizer which is a repulsive point of the recursion because of
the high curvature. If one increases the value of $k$, the behavior remains qualitatively similar,
with flatter and more curved areas as illustrated in Figure 2 right. Actually, in this case, $\varphi^k$
converges pointwise to $f(-1) = 2.25$ for $y < 0$, $f(1) = 0.25$ for $y > 0$ and $f(0) = 0.25$ for $y = 0$.
However, $\min \varphi^k = 0$ for all $k$ and the argmin converges to 0, the local maximum of $g$. This is
an example of pointwise, but non uniform convergence, and the limit is indeed discontinuous.

<figure><img src="image_8.png" /><figcaption>**Figure 3:** Same as Figure 2 with $f(x, y) = (y - 2)^2(\text{left})$. We see (middle) that the recursion spends time with value $\varphi^k(y_l) \sim f(-1)$ and then converges to a point corresponding the global minimum of $f$. The specificity of this setting is that when the inner iteration counter $k$ is increasing, the corresponding argminimizer of $f \circ A^k$ is pushed to infinity (right).</figcaption></figure>

**Illustrating escape at infinity:** We consider the same experiment with $f(x, y) = (y - 2)^2$
instead. This time the global minimizer of $f$ is 2 which is on an infinite branch of $g$. The
inner gradient descent steps are still attracted by $\pm 1$. But for a given $k$, it is possible to find
$A^k(y) \approx 2$ by choosing $y$ large enough. For large $k$, the corresponding initialization $y$ diverges.
Furthermore, this corresponds to a region of positive curvature for $g$ and it does not generate
the sharp minimizer behavior for $f \circ A^k$, which we had previously. It modifies the landscape
of $f$, pushing smoothly the global argminimizer to infinity. In this case, the gradient descent
algorithm applied to $f \circ A^k$ finds the global minimum of $f$, ignoring the bilevel constraints as
illustrated in Figure 3. In this case $\varphi^k$ also converges pointwise to a discontinuous function.
## References

[1] Gemayqzel Bouza Allende and Georg Still. Solving bilevel programs with the kkt-approach. *Mathematical Programming*, 138:309–332, 2012.

[2] Eugene L Allgower and Kurt Georg. Numerical path following. *Handbook of numerical analysis*, 5:3–207, 1997.

[3] Andreas Antoniou, Harrison Edwards, and Amos Storkey. How to train your MAML. In *International Conference on Learning Representations*, 2019.

[4] Michael Arbel and Julien Mairal. Non-Convex Bilevel Games with Critical Point Selection Maps. In *NeurIPS 2022 - 36th Conference on Neural Information Processing Systems*, Advances in Neural Information Processing Systems (NeurIPS) 2022, pages 1–34, New Orleans, United States, November 2022.

[5] Hédy Attouch, Jérôme Bolte, Patrick Redont, and Antoine Soubeyran. Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the kurdyka-lojasiewicz inequality. *Mathematics of operations research*, 35(2):438–457, 2010.

[6] Hédy Attouch, Jérôme Bolte, and Benar Fux Svaiter. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods. *Mathematical Programming*, 137:91–129, 2013.

[7] Shaojie Bai, J. Zico Kolter, and Vladlen Koltun. Deep equilibrium models. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc., 2019.

[8] Jérôme Bolte, Antoine Hochart, and Edouard Pauwels. Qualification conditions in semialgebraic programming. *SIAM journal on Optimization*, 28(2):1867–1891, 2018.

[9] Jerome Bolte, Quoc-Tung Le, and Edouard Pauwels. Convergence of optimizers implies eigenvalues filtering at equilibrium. working paper or preprint, September 2025.

[10] Jérôme Bolte, Quoc-Tung Le, Edouard Pauwels, and Samuel Vaiter. Geometric and computational hardness of bilevel programming. *arXiv preprint arXiv:2407.12372*, 2024.

[11] Jérôme Bolte, Tam Le, Éric Moulines, and Edouard Pauwels. Inexact subgradient methods for semialgebraic functions, April 2024. working paper or preprint.

[12] Lesi Chen, Jing Xu, and Jingzhao Zhang. On finding small hyper-gradients in bilevel optimization: Hardness results and improved analysis. In *The Thirty Seventh Annual Conference on Learning Theory*, pages 947–980. PMLR, 2024.

[13] Tianyi Chen, Yuejiao Sun, and Wotao Yin. Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, *Advances in Neural Information Processing Systems*, 2021.

[14] Benoît Colson, Patrice Marcotte, and Gilles Savard. An overview of bilevel optimization. *Annals of operations research*, 153:235–256, 2007.

[15] Michel Coste. *An introduction to o-minimal geometry*. Istituti editoriali e poligrafici internazionali Pisa, 2000.

[16] Michel Coste. An introduction to semialgebraic geometry, 2000.
[17] 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 (NeurIPS)*, 2022.

[18] Mathieu Dagréou, Thomas Moreau, Samuel Vaiter, and Pierre Ablin. A Lower Bound and a Near-Optimal Algorithm for Bilevel Empirical Risk Minimization. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, 2024.

[19] S. Dempe. *Foundations of Bilevel Programming*. Nonconvex Optimization and Its Applications. Springer US, US, 2002.

[20] S. Dempe, B. S. Mordukhovich, and A. B. Zemkoho. Sensitivity analysis for two-level value functions with applications to bilevel programming. *SIAM Journal on Optimization*, 22(4):1309–1343, 2012.

[21] Stephan Dempe and Alain Zemkoho. Bilevel optimization. In *Springer optimization and its applications*, volume 161. Springer, 2020.

[22] Chelsea Finn, Pieter Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In *Proceedings of the 34th International Conference on Machine Learning - Volume 70*, ICML’17, page 1126–1135. JMLR.org, 2017.

[23] Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In *International Conference on Machine Learning*, 2018.

[24] Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. *arXiv: Optimization and Control*, 2018.

[25] Jean Charles Gilbert. Automatic differentiation and iterative processes. *Optimization Methods and Software*, 1(1):13–21, 1992.

[26] Xavier Goudou and Julien Munier. The gradient and heavy ball with friction dynamical systems: the quasiconvex case. math. program. ser. b 116(1-2), 173-191. *Math. Program.*, 116:173-191, 01 2009.

[27] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration complexity of hypergradient computation. In *International Conference on Machine Learning*, pages 3748–3758. PMLR, 2020.

[28] Andreas Griewank and Andrea Walther. *Evaluating derivatives: principles and techniques of algorithmic differentiation*. SIAM, 2008.

[29] V. Guillemin and A. Pollack. *Differential Topology*. AMS Chelsea Publishing. AMS Chelsea Pub., 2010.

[30] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In *2015 IEEE International Conference on Computer Vision (ICCV)*, pages 1026–1034, 2015.

[31] René Henrion and Thomas Surowiec. On calmness conditions in convex bilevel programming. *Applicable Analysis*, 90:951–970, 06 2011.

[32] Kaiyi Ji, Jason D Lee, Yingbin Liang, and H Vincent Poor. Convergence of meta-learning with task-specific adaptation over partial parameters. *Advances in Neural Information Processing Systems*, 33:11490–11500, 2020.
[33] Kaiyi Ji and Yingbin Liang. Lower bounds and accelerated algorithms for bilevel optimization. *J. Mach. Learn. Res.*, 24(1), March 2024.

[34] Kaiyi Ji, Mingrui Liu, Yingbin Liang, and Lei Ying. Will bilevel optimizers benefit from loops. *Advances in Neural Information Processing Systems*, 35:3011–3023, 2022.

[35] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Theoretical convergence of multi-step model-agnostic meta-learning. *J. Mach. Learn. Res.*, 23:29:1–29:41, 2020.

[36] Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In *Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning*, volume 139 of *Proceedings of Machine Learning Research*, pages 4882–4892. PMLR, 18–24 Jul 2021.

[37] Charles D Kolstad and Leon S Lasdon. Derivative evaluation and computational experience with large bilevel mathematical programs. *Journal of optimization theory and applications*, 65:485–499, 1990.

[38] Jeongyeol Kwon, Dohyun Kwon, Steve Wright, and Robert Nowak. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. *arXiv preprint arXiv:2309.01753*, 2023.

[39] Jason D. Lee, Max Simchowitz, Michael I. Jordan, and Benjamin Recht. Gradient descent only converges to minimizers. In *Vitaly Feldman, Alexander Rakhlin, and Ohad Shamir, editors, 29th Annual Conference on Learning Theory*, volume 49 of *Proceedings of Machine Learning Research*, pages 1246–1257, Columbia University, New York, New York, USA, 23–26 Jun 2016. PMLR.

[40] Gui-Hua Lin, Mengwei Xu, and Jane Ye. On solving simple bilevel programs with a nonconvex lower level program. *Mathematical Programming*, 144, 04 2014.

[41] Hanxiao Liu, Karen Simonyan, and Yiming Yang. DARTS: Differentiable Architecture Search. In *International Conference on Learning Representations (ICLR)*, 2019.

[42] Risheng Liu, Xuan Liu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. A value-function-based interior-point method for non-convex bi-level optimization. In *Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event*, volume 139 of *Proceedings of Machine Learning Research*, pages 6882–6892, virtual, 2021. PMLR.

[43] Risheng Liu, Yaohua Liu, Shangzhi Zeng, and Jin Zhang. Towards gradient-based bilevel optimization with non-convex followers and beyond. In *A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems*, 2021.

[44] Risheng Liu, Pan Mu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang. A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton. In *International Conference on Machine Learning*, Vienna, Austria, 2020.

[45] Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing Millions of Hyperparameters by Implicit Differentiation. In *International Conference on Artificial Intelligence and Statistics (AISTAT)*, 2020.

[46] Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient-based hyperparameter optimization through reversible learning. In *International conference on machine learning*, pages 2113–2122. PMLR, 2015.
[47] Sheheryar Mehmood and Peter Ochs. Automatic differentiation of some first-order methods in parametric optimization. In *International Conference on Artificial Intelligence and Statistics*, pages 1584–1594. PMLR, 2020.

[48] Roey Merchav, Shoham Sabach, and Marc Teboulle. A fast algorithm for convex composite bi-level optimization, 2024.

[49] Boris S. Mordukhovich. Bilevel optimization and variational analysis. In Stephan Dempe and Alain Zemkoho, editors, *Bilevel Optimization: Advances and Next Challenges*, pages 197–226. Springer International Publishing, Cham, 2020.

[50] Ioannis Panageas and Georgios Piliouras. Gradient descent only converges to minimizers: Non-isolated critical points and invariant regions. In *8th Innovations in Theoretical Computer Science Conference (ITCS 2017)*. Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2017.

[51] Jong-Shi Pang, Shih-Ping Han, and Narayan Rangaraj. Minimization of locally lipschitzian functions. *SIAM Journal on Optimization*, 1(1):57–82, 1991.

[52] Fabian Pedregosa. Hyperparameter optimization with approximate gradient. In Maria Florina Balcan and Kilian Q. Weinberger, editors, *Proceedings of The 33rd International Conference on Machine Learning*, volume 48 of *Proceedings of Machine Learning Research*, pages 737–746, New York, New York, USA, 20–22 Jun 2016. PMLR.

[53] Robin Pemantle. Nonconvergence to Unstable Points in Urn Models and Stochastic Approximations. *The Annals of Probability*, 18(2):698–712, 1990.

[54] Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine. Meta-learning with implicit gradients. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc., 2019.

[55] Cédric Rommel, Thomas Moreau, Joseph Paillard, and Alexandre Gramfort. CADDA: Class-wise automatic differentiable data augmentation for EEG signals. In *International Conference on Learning Representations*, 2022.

[56] Gilles Savard and Jacques Gauvin. The steepest descent direction for the nonlinear bilevel programming problem. *Operations Research Letters*, 15(5):265–272, 1994.

[57] Han Shen and Tianyi Chen. On penalty-based bilevel gradient descent method. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, *Proceedings of the 40th International Conference on Machine Learning*, volume 202 of *Proceedings of Machine Learning Research*, pages 30992–31015. PMLR, 23–29 Jul 2023.

[58] M. Shub, A. Fathi, and R. Langevin. *Global Stability of Dynamical Systems*. Springer, 1987.

[59] Vladislav B. Tadić and A. Doucet. Asymptotic bias of stochastic gradient search. In *2011 50th IEEE Conference on Decision and Control and European Control Conference*, pages 722–727, 2011.

[60] Lou van den Dries and Chris Miller. Geometric categories and o-minimal structures. *Duke Mathematical Journal*, 84(2):497 – 540, 1996.

[61] Luis N Vicente and Paul H Calamai. Bilevel and multilevel programming: A bibliography review. *Journal of Global optimization*, 5(3):291–306, 1994.
[62] J. J. Ye and D. L. Zhu. Optimality conditions for bilevel programming problems. *Optimization*, 33(1):9-27, 1995.