Michael Arbel

Julien Mairal

Non-Convex Bilevel Games with Critical Point Selection Maps

Bilevel optimization problems involve two nested objectives, where an upper-level objective depends on a solution to a lower-level problem. When the latter is nonconvex, multiple critical points may be present, leading to an ambiguous definition of the problem. In this paper, we introduce a key ingredient for resolving this ambiguity through the concept of a selection map which allows one to choose a particular solution to the lower-level problem. Using such maps, we define a class of hierarchical games between two agents that resolve the ambiguity in bilevel problems. This new class of games requires introducing new analytical tools in Morse theory to extend implicit differentiation, a technique used in bilevel optimization resulting from the implicit function theorem. In particular, we establish the validity of such a method even when the latter theorem is inapplicable due to degenerate critical points. Finally, we show that algorithms for solving bilevel problems based on unrolled optimization solve these games up to approximation errors due to finite computational power. A simple correction to these algorithms is then proposed for removing these errors.

Introduction

Bilevel optimization has proven to be a major tool for solving machine learning problems that possess a nested structure such as hyper-parameter optimization [START_REF] Feurer | Hyperparameter optimization[END_REF], meta-learning [START_REF] Bertinetto | Meta-learning with differentiable closed-form solvers[END_REF], reinforcement learning [START_REF] Hong | A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic[END_REF][START_REF] Liu | Investigating bi-level optimization for learning and vision from a unified perspective: A survey and beyond[END_REF], or dictionary learning [START_REF] Mairal | Task-driven dictionary learning[END_REF]. Introduced in the field of economic game theory in [START_REF] Von | MarktformundGleichgewicht[END_REF], a bilevel optimization problem can be understood as a game between a leader and a follower each of which optimizes their own objective function but where the leader can anticipate follower's actions. In the context of machine learning, the leader typically optimizes a hyper-parameter over a validation loss while the follower optimizes the model parameter on a training loss [START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF].

Bilevel optimization introduces many challenges. In particular, when multiple optimal solutions are available to the follower, the leader would need to optimize a different objective depending on the follower's strategy to select an optimal solution. As a result, the bilevel problem becomes ambiguously defined without knowing the follower's strategy [START_REF] Liu | Towards gradient-based bilevel optimization with non-convex followers and beyond[END_REF]. A large body of work on bilevel programs for machine learning gets around these considerations by assuming the follower to have a unique optimal choice, a situation that typically occurs when the follower's objective is strongly convex, leading to efficient and scalable algorithms [START_REF] Ablin | Super-efficiency of automatic differentiation for functions defined as a minimum[END_REF][START_REF] Arbel | Amortized implicit differentiation for stochastic bilevel optimization[END_REF][START_REF] Blondel | Efficient and modular implicit differentiation[END_REF][START_REF] Domke | Generic methods for optimization-based modeling[END_REF][START_REF] Gould | On differentiating parameterized argmin and argmax problems with application to bi-level optimization[END_REF][START_REF] Liao | Reviving and improving recurrent back-propagation[END_REF][START_REF] Liu | Investigating bi-level optimization for learning and vision from a unified perspective: A survey and beyond[END_REF][START_REF] Shaban | Truncated backpropagation for bilevel optimization[END_REF]. However, in many machine learning applications, the strong convexity of the follower's objective is an unrealistic assumption. This is particularly the case in the context of deep learning, where the follower's objective, the training loss, can be highly non-convex in the parameters of the model and can have regions of flat optima due to symmetries and other degeneracies [START_REF] Draxler | Essentially no barriers in neural network energy landscape[END_REF][START_REF] Li | Visualizing the loss landscape of neural nets[END_REF].

In the literature on mathematical optimization, the ambiguity in bilevel problems is often resolved by making an additional assumption on the follower's strategy for choosing their optimal solution. In particular, two problems are often considered: optimistic and pessimistic bilevel programs, see [START_REF] Dempe | New necessary optimality conditions in optimistic bilevel programming[END_REF].

36th Conference on Neural Information Processing Systems (NeurIPS 2022).

Both problems rely on two assumptions: (i) the follower is using a strategy for selecting a solution to their problem that is either improving or degrading the leader's objective and (ii) the leader knows exactly what strategy the follower is using. These assumptions are strong from a game-theoretical perspective and often unrealistic for machine learning problems such as hyper-parameter optimization. Still, optimistic/pessimistic bilevel games are well defined and early works have proposed several algorithms to solve them with strong convergence guarantees [START_REF] Jj Ye | Necessary optimality conditions for optimization problems with variational inequality constraints[END_REF][START_REF] Jj Ye | Optimality conditions for bilevel programming problems[END_REF][START_REF] Jj Ye | Exact penalization and necessary optimality conditions for generalized bilevel programming problems[END_REF]]. Yet, these algorithms are often ill-suited to large-scale and high-dimensional problems arising in machine learning applications as they rely on second-order optimization methods such as Newton's method [START_REF] Guo | Solving mathematical programs with equilibrium constraints[END_REF]. For this reason, scalable first-order algorithms for such games have been proposed recently [START_REF] Liu | Value-function-based sequential minimization for bi-level optimization[END_REF][START_REF] Liu | Towards gradient-based bilevel optimization with non-convex followers and beyond[END_REF].

However, many of the best-performing approaches for hyper-parameter optimization rely neither on an optimistic nor a pessimistic formulation of the bilevel problem [START_REF] Vicol | Implicit regularization in overparameterized bilevel optimization[END_REF]. Instead, they often rely on algorithms initially designed for bilevel problems with strongly convex lower objectives even though the convexity assumption does not hold [START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF]. Consequently, these algorithms are solving a seemingly ill-defined bilevel program due to the ambiguity in the way the follower selects their solution. However, their ability to provide models with good empirical performance raises the question of whether these algorithms are solving another class of well-defined hierarchical problems beyond optimistic and pessimistic bilevel programs that are still relevant for machine learning.

In this work, we answer the above question by introducing Bilevel Games with Selection (BGS), a class of games between two agents: a leader and a follower, where the leader uses a mechanism for anticipating the solution of the follower without knowing the exact follower's strategy. We define such a mechanism using the notion of a selection, which is simply a map for selecting a particular solution to the follower's objective given the current state of the game. In particular, BGS recovers a usual bilevel program when the follower's objective admits a unique solution. By playing a BGS, the agents seek an equilibrium point for which each of their objectives ceases to vary. The equilibria are completely determined by the selection thus resulting in a well-defined problem.

When the selection is differentiable, the equilibrium point can be characterized by a first-order optimality condition which enables gradient-based approximations. More precisely, we show that implicit differentiation [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF], which, a priori, is only valid when the critical points of the follower's objective are non-degenerate, remains applicable for solving BGS even when these critical points are degenerate. To this end, we consider a general construction of the selection as the limit of a gradient flow of the follower's objective and prove the differentiability of such a selection near local minimizers, provided the follower's objective satisfies a generalization of the Morse-Bott property [START_REF] David | Morse-bott theory and equivariant cohomology[END_REF][START_REF] Feehan | On the morse-bott property of analytic functions on banach spaces with łojasiewicz exponent one half[END_REF]. We then characterize the differential of the selection as a solution to a linear system thus extending implicit differentiation to degenerate critical points. Finally, we leverage this characterization to show that popular algorithms based on iterative differentiation (ITD) [START_REF] Gunes Baydin | Automatic differentiation in machine learning: a survey[END_REF] find fixed points approximating the BGS's equilibria up to approximation errors. We then introduce a simple corrective term to these algorithms based on implicit differentiation to remove these errors.

Related Work

Iterative/Unrolled optimization (ITD) is a class of methods approximating the lower-level solution map by a differentiable function obtained through successive gradient updates [START_REF] Gunes Baydin | Automatic differentiation in machine learning: a survey[END_REF]. When the lower-level objective is strongly convex, these algorithms solve a well-defined bilevel problem up to an error that is controlled by increasing the computational budget for the approximate solution [START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF]. Our analysis suggests a simple algorithmic correction to these approaches which can result in solutions to a bilevel game with a constant budget for the approximate solution.

Approximate Implicit Differentiation (AID) is a class of methods approximating the variations of the lower-level solution map using the Implicit Function theorem [START_REF] Franceschi | Bilevel programming for hyperparameter optimization and meta-learning[END_REF][START_REF] Ghadimi | Approximation methods for bilevel programming[END_REF][START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF][START_REF] Rajeswaran | Meta-Learning with Implicit Gradients[END_REF]. The nondegeneracy requirement under which the latter theorem holds restricts the applicability of AID to, essentially, strongly convex lower-level objectives. These algorithms admit fixed points that match the solutions to the bilevel problem [START_REF] Ghadimi | Approximation methods for bilevel programming[END_REF][START_REF] Hong | A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic[END_REF][START_REF] Ji | Lower bounds and accelerated algorithms for bilevel optimization[END_REF][START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF]. As such, they typically require a smaller computational budget than ITD [START_REF] Arbel | Amortized implicit differentiation for stochastic bilevel optimization[END_REF][START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF]. Recently, [START_REF] Bolte | Nonsmooth implicit differentiation for machine-learning and optimization[END_REF][START_REF] Bolte | Automatic differentiation of nonsmooth iterative algorithms[END_REF][START_REF] Bolte | Nonsmooth automatic differentiation: a cheap gradient principle and other complexity results[END_REF] extended AID to non-smooth objectives while still requiring non-degenerate critical points. The present work is complementary to these works as it extends AID to smooth objectives that have possibly degenerate critical points.

Optimistic and pessimistic bilevel optimization. When the lower-level objective is non-convex, the ambiguity of the problem arising from the multiplicity of the lower-level solutions can be re-solved by optimizing the upper-level objective over all such possible solutions [START_REF] Wiesemann | Pessimistic bilevel optimization[END_REF][START_REF] Zemkoho | Solving ill-posed bilevel programs[END_REF]. The optimistic and pessimistic problems arise when either minimizing or maximizing the upper-level over all such lower-level solutions. Early works proposed to solve these problems using exact penalization [START_REF] Jj Ye | Exact penalization and necessary optimality conditions for generalized bilevel programming problems[END_REF], second-order optimization [START_REF] Jj Ye | Necessary optimality conditions for optimization problems with variational inequality constraints[END_REF][START_REF] Jj Ye | Optimality conditions for bilevel programming problems[END_REF] or smoothing method [START_REF] Xu | A smoothing augmented lagrangian method for solving simple bilevel programs[END_REF]. However, these approaches are hard to scale to the high dimensional problems arising in machine learning. More recently, [START_REF] Liu | Towards gradient-based bilevel optimization with non-convex followers and beyond[END_REF][START_REF] Liu | Value-function-based sequential minimization for bi-level optimization[END_REF] considered first-order methods based on unrolled optimization or interior-point methods for solving optimistic bilevel problems and provided approximation guarantees. However, as shown in [START_REF] Vicol | Implicit regularization in overparameterized bilevel optimization[END_REF], most practical applications to bilevel optimization rely on a formulation that goes beyond optimistic or pessimistic formulations. The present work departs from these approaches and instead introduces a bilevel game that is more tractable to solve. We show that popular bilevel algorithms, such as unrolled optimization, yield approximations of these games.

Non-Convex Bilevel Optimization with Selection

Notations. Define X = R p and Y = R d for some positive integers p and d. We consider two real valued functions f and g defined on X ×Y and assume g to be twice-continuously differentiable.

Background on Bilevel Optimization

A bilevel program is an optimization problem where an upper-level objective f defined over a set X × Y of variables (x, y) is optimized in the first variable x under the constraint that the second variable y is optimal for a lower-level objective y → g(x, y) depending on the upper-variable x. When g(x, .) admits a unique minimizer denoted by y (x), which is the case if y → g(x, y) is strongly convex, the bilevel problem is well-defined and can be expressed as:

min x∈X f (x, y (x)), y (x) := arg min y∈Y g(x, y). (BP)
When g is non-convex, the set of minimizers T (x):= arg min y g(x, y) may contain more than one element making (BP) ambiguous. A possible approach for resolving the ambiguity is to adopt a game-theoretical point of view, where a lower-level agent uses a particular strategy for selecting a solution in T (x). For instance, in pessimistic bilevel games, the lower agent chooses a minimizer of g(x, .) that maximizes f (x, .) while the upper agent minimizes the resulting worst-case loss F in x:

(UL): min Similarly, an optimistic bilevel game can be obtained by replacing maximization with minimization so that both agents cooperate. While these approaches are highly relevant from a game-theoretical point of view, many machine learning applications do not rely on a pessimistic/optimistic bilevel formulation. For instance, for hyper-parameter optimization, the lower agent may have access to training data, but it should not have access to the validation data processed (used in f ) by the upper agent. Instead, a popular approach consists of applying algorithms designed for bilevel programs that admit unique solutions for the lower problems, even though this assumption may not hold in practice [START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF]. In the next section, we introduce a class of games that allow characterizing the equilibrium points obtained by these popular algorithms while resolving the ambiguity of non-convex bilevel problems and bypassing the limitations of pessimistic/optimistic bilevel formulations.

Bilevel Games with Selection (BGS)

We introduce a new class of nested games for bilevel optimization with two agents, a leader and a follower. The follower minimizes the lower-level objective g w.r.t. a variable y in Y. Similarly, the leader minimizes the upper-level objective f w.r.t. a variable x ∈ X while anticipating the follower's solution. More precisely, the leader has access to a selection map: φ : X ×Y → Y to choose a unique critical point φ(x, y) of y → g(x, y) given the current state of the game (x, y) ∈ X ×Y thus allowing the leader to anticipate the follower's solution. Typically, the selection φ(x, y) represents the critical point that is selected by an optimization process of g(x, .) starting from an initial condition y (e.g., the limit of a gradient flow for a gradient descent algorithm). The Bilevel Game with Selection (BGS) is therefore defined as the following interdependent optimization problems:

(UL): Given a selection map φ, the game (BGS) is well-defined and does not suffer from the ambiguity problem in (BP). The explicit dependence of φ(x, y) on the initialization y might seem unnecessary at first, as one could simply fix y to some value y 0 and consider only the dependence on the variable x. However, such a dependence on the variable y allows performing warm-start [START_REF] Vicol | Implicit regularization in overparameterized bilevel optimization[END_REF], where the lower-level problem is optimized starting from a previous state of the game, thus resulting in computational savings Figure 1. We provide below a formal definition for the selection map. Definition 1 (Selection map). Given a continuously differentiable function g : X × Y → R, the map φ : X ×Y → Y is a selection if it satisfies the following properties for any pair (x, y) ∈ X ×Y:

min x∈X L φ (x, y) := f (x, φ(x, y)), (
1. Criticality: The element y = φ(x, y) is a critical point of g(x, .), i.e. ∂ y g(x, y ) = 0.

Self-consistency:

If y is a critical point of g(x, .) i.e. ∂ y g(x, y) = 0, then φ(x, y) = y.

Criticality ensures the leader possesses a hierarchical advantage in that they know what are the optimal choices accessible to the follower. Self-consistency implies that the leader makes a guess that is not contradicting the current choice y of the follower. Both properties ensure the leader can rationally anticipate the follower's actions from the current state of the game (x, y). We will see in Section 4, under mild assumptions on g, that it is always possible to define a selection φ as the limit of a continuous-time gradient flow of y → g(x, y) initialized at y. Moreover, as we discuss later in Section 5, the selection does not need to be explicitly constructed for solving (BGS) in practice. It can be simply related to the implicit bias of the algorithm used for solving the follower's problem.

Connection to (BP).

When the lower-level objective y → g(x, y) admits a unique minimizer y (x), it is easy to check that there exists a unique selection map φ satisfies φ(x, y)=y (x). Hence, (BGS) recovers the bilevel problem in (BP) as a particular case.

Connection to (pessimistic-BG) or the optimistic variant. Key differences between (BGS) and pessimistic or optimistic games is that (i) the follower has never access to the upper function f with (BGS), which matches practical hyper-parameter optimization applications where f relies on a validation dataset, whereas g relies on a distinct training set; (ii) the leader in (pessimistic-BG) does not take into account the strategy used by the follower, whereas the leader in (BGS) makes more rational choices by guessing the strategy of the follower through the selection map φ.

First-order equilibrium conditions. The agents can play the game (BGS) by successively taking actions (x k , y k ) to improve their own objectives x → L φ (x, y k-1 ) and y → g(x k , y), by hoping the strategy will reach an equilibrium pair (x * , y * ) Figure 1(Right). In the case where f , g and φ are differentiable at (x * , y * ), the equilibrium pair is characterized by a first-order stationary condition:

∂ x L φ (x , y ) = ∂ x f (x , y ) + ∂ x φ(x , y )∂ y f (x , y ) = 0, ∂ y g(x , y ) = 0. (SC)
When g is smooth and strongly convex in y, the implicit function theorem [28, Theorem 5.9] ensures that φ is differentiable and provides an expression of ∂ x φ(x , y ) as a solution to a linear system which key for implicit differentiation. This allows to devise efficient algorithms using estimates of the gradient ∂ x L φ , see, e.g., [START_REF] Arbel | Amortized implicit differentiation for stochastic bilevel optimization[END_REF]. However, extensions of the implicit function theorem, such as the constant rank theorem [29, Theorem 4.12], for cases where g has possibly degenerate critical points require strong assumptions on g which are unrealistic in machine learning. In the next section, we provide new analytical tools for extending implicit differentiation by studying the differentiability of a family of selection maps corresponding to a large class of functions g. The resulting expression will be key for devising first-order methods to solve (BGS), as discussed in Section 5.

Selection Based on Gradient Flows for Parameteric Morse-Bott Functions

In this section, we extend implicit differentiation to a class of functions with possibly degenerate critical points. To this end, we consider a particular selection φ(x, y) obtained as the limit of a gradient flow (φ t (x, y)) t≥0 of g(x, .) initialized at y. We then study the differentiability w.r.t. x of the selection by analyzing the dynamics of such a gradient flow. For general non-convex functions, the selection might be non-differentiable since a small perturbation to the parameter x can change the geometry of the critical points of g, causing the perturbed flow to move away from the non-perturbed one (see Figure 2). We are therefore interested in functions g preserving the local geometry near critical points as x varies. In Section 4.1, we introduce such a class of functions called parametric Morse-Bott functions, which covers many practical machine learning models. We then show, in Section 4.2, that the selection resulting from such a function is differentiable near local minima.

Parameteric Morse-Bott Functions

We introduce parametric Morse-Bott functions, a class of parametric functions g : X × Y → R with parameter x in X extending the more familiar notion of Morse-Bott functions (Appendix A.1, [START_REF] Feehan | On the morse-bott property of analytic functions on banach spaces with łojasiewicz exponent one half[END_REF]) to account for the effect of the parameter x on the geometry of critical points. Definition 2 (Parametric Morse-Bott function.). Let g : X × Y be a real-valued twice continuousely differentiable function and define the set of augmented critical points M as follows:

M := {(x, y) ∈ X × Y | ∂ y g(x, y) = 0} (1)
Let (x 0 , y 0 ) ∈ M. We say that g is Morse-Bott at y 0 w.r.t. x 0 , if there exists an open neighbordhood V of (x 0 , y 0 ) s.t. the intersection M ∩ V is a C 2 -connected sub-manifold of X ×Y of dimension:

dim(M ∩ V) = dim(X ) + dim Ker(∂ 2 yy g(x 0 , y 0 )) . g is a parametric Morse-Bott function if for any (x 0 , y 0 )∈M, g is Morse-Bott at y 0 w.r.t. x 0 .

The functions in Definition 2 satisfy a condition that is stronger than simply satisfying the Morse-Bott property at any parameter value x (Definition 3 of Appendix A.1). Indeed, we show in Proposition 7 of Appendix A.2 that, for any x 0 ∈ X , the function y → g(x 0 , y) is a Morse-Bott function, meaning that the critical set C(x 0 ) of y → g(x 0 , .) near a critical point y 0 is locally a C 2 connected sub-manifold of Y of dimension equal to the dimension of the null-space of the Hessian ∂ 2 yy g(x 0 , y 0 ). For conciseness, we introduce the following assumption which ensures g satisfies the condition of Definition 2 as well as possesses continuous third-order derivatives. Assumption 1 (Parameteric Morse-Bott property). The function g is at least three-times continuously differentiable and is a parameteric Morse-Bott function as defined in Definition 2.

Examples of parametric Morse-Bott function. A notable class of parametric Morse-Bott functions is the one containing all twice-continuously differentiable functions that are strongly convex or, more generally, possess only non-degenerate critical points in the second variable as shown in Proposition 8 of Appendix A.2. Note that parametric Morse-Bott functions need not be convex and can have multiple (possibly degenerate) local minima, saddle-points, and local maxima.

Another class of functions, this time with possibly degenerate critical points, are those that can be expressed as a composition of some Morse-Bott function h and a family (τ x ) x∈X of diffeomorphisms on Y parameterized by x, i.e. g(x, y)=h(τ x (y)). This particular form is relevant in generative modeling where the diffeomorphisms are defined using normalizing flows of parameter x [START_REF] Rezende | Variational inference with normalizing flows[END_REF].

The condition in Definition 2 ensures that the degree of freedom of the augmented critical set M is exactly determined by the degree of freedom of the parameter x and the degree of degeneracy of the Hessian at a critical point y. This condition is precisely what guarantees the stability of the local shape of critical points when the parameter x varies as we formalize through the next theorem. Theorem 1 (Morse-Bott lemma with parameters). Let g be a function satisfying Assumption 1. Let (x 0 , y 0 ) in M be an augmented critical point of g. Denote by K the null space of the Hessian A 0 :=∂ 2 yy g(x 0 , y 0 ) and by K ⊥ its orthogonal complement in Y. Let J 0 be a diagonal matrix with diagonal element given by the sign of the non-zero eigenvalues of A 0 . Then, there exists open neighborhoods U and V of (x 0 , 0 K , 0 K ⊥ ) and (x 0 , y 0 ) in X ×K×K ⊥ and X ×Y, and a diffeomorphism ψ : U → V preserving the first variable, i.e. ψ(x, r, w)=(x, y) for any (x, r, w) ∈ U, with ψ(x 0 , 0 K , 0 K ⊥ )=(x 0 , y 0 ) such that g admits the representation:

g(ψ(x, r, w)) = g(ψ(x, 0 K , 0 K ⊥ )) + 1 2 w J 0 w, ∀(x, r, w) ∈ U.
Theorem 1, which is proven in Appendix A.3, shows that, near an augmented critical point (x 0 , y 0 ), g looks like a quadratic function up to an additive term that depends only on the parameter x. Moreover, slightly varying the parameter x does not change the quadratic function and thus preserves the local shape near critical points. Theorem 1 is an extension of the Morse-Bott lemma [16, Theorem 2.10] to the case when there is a dependence on a parameter x. It can also be seen as an extension of the Morse lemma with parameters [16, Theorem 4] which allows dependence to a parameter x but requires the critical points to be non-degenerate (invertible matrix A 0 ). To our knowledge, Theorem 1 is the first result in the literature providing a decomposition of parametric functions with degenerate critical points into the sum of a quadratic non-degenerate term and a singular term depending only on the parameter x. We present now a corollary of Theorem 1 which is a strengthened version of the standard Łojasiewicz inequality [START_REF] Lojasiewicz | Sur les trajectoires du gradient d'une fonction analytique[END_REF] that will be essential for our subsequent analysis. Proposition 1 (Locally Uniform Łojasiewicz gradient inequality). Let g be a function satisfying Assumption 1 and let (x 0 , y 0 ) be in M the augmented critical set defined in Definition 2. Then, there exists an open neighborhood U of (x 0 , y 0 ) and a positive number µ > 0 such that y → g(x, y) is constant on the set M ∩ U with some common value G(x) := g(x, y) and the following holds:

µ|g(x, y) -G(x)| ≤ 1 2 ∂ y g(x, y) 2 , ∀(x, y) ∈ U.
Proposition 1, which is proven in Appendix A.3, ensures that the Łojasiewicz gradient inequality holds uniformly on (x, y) near any augmented critical point (x 0 , y 0 ). This result will be essential in Section 4.2 for defining a selection φ obtained as limits of gradient flows and to obtain a locally uniform control of these flows in the parameter x. This in turn will allow us to obtain the differentiability of the selection in the parameter x whenever φ(x, y) is a local minimum.

Smoothness of Selections Based on Gradient Flows of a Parametric Morse-Bott Function

We consider a construction for the selection φ in Definition 1 as a limit of a continuous-time gradient flow of g. More precisely, we define a continuous-time trajectory (φ t (x, y)) t≥0 in Y initialized at φ 0 (x, y) = y and driven by the differential equation: dφ t (x, y) dt = -∂ y g(x, φ t (x, y)).

(GF) Provided φ t (x, y) converges towards some element φ(x, y) as t→+∞, we can expect such a limit to satisfy both conditions of Definition 1, therefore constituting a valid selection. However, for general non-convex functions, φ t (x, y) might not always converge [START_REF] Lojasiewicz | Sur les trajectoires du gradient d'une fonction analytique[END_REF]. To guarantee the existence and convergence of the flow, we make the following assumptions on the function g. Assumption 2 (Smoothness). There exists L>0 such that y →∂ y g(x, y) is L-Lipschitz for any x∈X . Assumption 3 (Coercivity). For any x ∈ X , it holds that g(x, y) → +∞ as y → +∞.

The smoothness assumption in Assumption 2 is standard and guarantees the existence of the flow by the Cauchy-Lipschitz theorem. The coercivity condition in Assumption 3 guarantees that φ t (x, y) cannot escape to infinity. It can be easily enforced by adding a small 2 -penalty to a non-negative loss (such as cross-entropy or mean-squared loss) which is already a common practice in machine learning. These assumptions, along with Assumption 1 ensure that the limit φ(x, y) always exists as we summarize in the following proposition, which is proven in Appendix B. Proposition 2. Under Assumptions 1 to 3, and for any (x, y) ∈ X × Y, the gradient flow (GF) always converges towards a critical point φ(x, y) of y → g(x, y) and the map (x, y) → φ(x, y) is a selection map as defined in Definition 1. We call φ the flow selection relatively to g.

Proposition 2 is a consequence of a general result that holds for functions satisfying a Łojasiewicz gradient inequality [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF][START_REF] Merlet | Convergence to equilibrium for discretizations of gradient-like flows on riemannian manifolds[END_REF] which is the case here by Proposition 1. From now on, we restrict our attention to the selection φ defined in Proposition 2. Even though φ satisfies the implicit equation ∂ y g(x, φ(x, y)) = 0, we cannot rely anymore on the implicit function theorem for studying the differentiability of φ(x, y) in x since g can have degenerate critical points. Instead, we propose to characterize the differentiability of φ by studying the limit of U t (x, y) := ∂ x φ t (x, y) which is formally driven by a linear differential equation of the form:

- dU t (x, y) dt = ∂ 2 xy g(x, φ t (x, y)) + U t (x, y)∂ 2 yy g(x, φ t (x, y)). (2)
Had we known in advance that φ(x, y) is differentiable in x, the limit U ∞ (x, y) of U t (x, y) as t → +∞, whenever defined, would be a promising candidate for the differential of φ(x, y) in x. Such a limit is indeed expected to satisfy the following linear equation:

0 = ∂ 2 xy g(x, φ(x, y)) + U ∞ (x, y)∂ 2 yy g(x, φ(x, y)). (3)
A first challenge is to ensure that U t does not diverge. For critical points φ(x, y) that are not local minima, it is easy to see that the Hessian ∂ 2 yy g(x, φ t (x, y)) must have a negative eigenvalue for t large enough, therefore causing the system (2) to diverge. Intuitively, unless φ(x, y) is a local minimum, there is no reason to expect φ(x, y) to be differentiable or even continuous in x, simply because φ(x, y) would be an unstable fixed-point of the flow φ t (x, y), so that any change in x might cause a large variation in φ(x, y). The possible non-differentiability of φ(x, y) for critical points that are not local minima is not problematic in practice, since for almost all initial conditions y of the flow φ t (x, y), the limit φ(x, y) is guaranteed to be a local minimizer [START_REF] Panageas | Gradient descent only converges to minimizers: Non-isolated critical points and invariant regions[END_REF]. In addition, we show in Proposition 13 of Appendix B.3 that if φ(x 0 , y) is a local minimum, then φ(x, y) must also be a local minimum in a neighborhood of x 0 .

Nevertheless, even for local minima, if the Hessian ∂ 2 yy g(x, φ(x, y)) is non-invertible, (3) might never hold if ∂ 2 xy g(x, φ(x, y)) does not belong to the image of the Hessian. However, we show in Proposition 6 of Appendix A.2 that, for any pair (x, y) of critical points, ∂ 2 xy g(x, y) must always belong to the span of the Hessian ∂ 2 yy g(x, y) as soon as g satisfies Assumption 1, therefore ensuring that (3) admits a solution. The following theorem, which is proven in Appendix C, establishes the differentiability of φ at local minima and shows that ∂ x φ is exactly given by the limit U ∞ . Theorem 2 (Differentiability of the flow selection). Let g be a function satisfying Assumptions 1 to 3 so that the flow selection φ is well-defined. Let (x 0 , y 0 ) be in X ×Y. If φ(x 0 , y 0 ) is a local minimizer of y → g(x 0 , y), then there exists a neighborhood U of x 0 on which x → φ(x, y 0 ) is differentiable with differential ∂ x φ(x, y 0 )=U ∞ (x, y 0 ). Moreover, if y 0 is a local minimizer of y →g(x 0 , y), then, denoting by † the pseudo inverse operator, ∂ x φ(x 0 , y 0 ) is exactly given by:

∂ x φ(x 0 , y 0 )= -∂ xy g(x 0 , y 0 )(∂ yy g(x 0 , y 0 )) † . (4)
The expression in ( 4) is very similar to the one that would arise by application of the implicit function theorem to a strongly convex function g. However, the proof technique does not rely on such a theorem which would not be applicable here. The key technical challenges in proving the above result are: (i) showing that φ(x, y) must be continuous at x 0 and (ii) controlling the error U t (x, y) -U ∞ (x, y) locally uniformly in x. The result follows by the application of classical uniform convergence results [START_REF] Rudin | Principles of mathematical analysis[END_REF]Theorem 7.17]. The continuity of φ is established in Proposition 12 of Appendix B.3 and relies on a stability analysis of the flow φ t performed in Appendix B.2. The uniform convergence of U t towards U ∞ is shown in Proposition 17 of Appendix C and relies on a local uniform convergence of the flow φ t towards φ which is proven in Proposition 14 of Appendix B.4. It is worth noting that, even though we identified ∂ x φ to be U ∞ , the latter is not fully characterized by (3) as it might contain a non-zero component in the null-space of the Hessian. However, when (x 0 , y 0 ) is an augmented critical pair of g, such a component vanishes, and ∂ x φ(x 0 , y 0 ) is exactly determined by the minimal norm solution in [START_REF] David | Morse-bott theory and equivariant cohomology[END_REF]. The latter fact has practical implications when designing algorithms for solving (BGS) as we discuss next.

Algorithms

Unrolled Optimization for BGS

Unrolled optimization constructs a map ϕ T (x, y) approximating a critical point of the function y → g(x, y) for any fixed x by applying a finite number T > 0 of gradient updates starting from some initial condition y. By convention, we set ϕ 0 (x, y)=y. Hence, ϕ T can be understood as an approximation to the selection map defined in Section 4.2. We emphasize that ϕ T is not a selection (Definition 1) since ϕ T (x, y) is not a critical point of g in general. Nevertheless, it provides a tractable approximation to critical points which is key for constructing practical algorithms for bilevel optimization. The gradient of ϕ T (x, y) w.r.t. x is then obtained by differentiating through the optimization steps and used to optimize the approximate upper-level objective:

L T (x, y) := f (x, ϕ T (x, y)).
Given the k-th upper-level iterate x k and an initial condition ỹk for the unrolled optimization, these approaches compute an approximation y k =ϕ T (x k-1 , ỹk ) and find an update direction d k for the upper-level variable x by differentiating L T (x, ỹk ) in x at the current iterate x k-1 . The following iterate x k is obtained by applying an update procedure, such as x k =x k-1 -γd k for positive small enough step-size γ. In Algorithm 1, we present several variants of these schemes, including a simple correction allowing them to solve (BGS) instead of an approximation.

The initial condition ỹk is often computed using a warm-start procedure ỹk =I M (x k-1 , y k-1 ). The simplest procedure is to set ỹk =y k-1 in which case I 0 (x, y)=y. However, it is not uncommon to perform M >0 optimization steps to minimize the objective y →g(x k-1 , y) starting from y k-1 . By doing so, gradient unrolling stops at ỹk and ignores the dependence of ỹk on y k-1 , resulting in Truncated unrolled optimization [START_REF] Shaban | Truncated backpropagation for bilevel optimization[END_REF]. Algorithm 1 summarizes these approaches when the binary variable AddCorrection is set to False. To characterize the limit points of Algorithm 1, we make the following assumptions on I M , ϕ T .

Algorithm 1 BGS-Opt(x 0 , y 0 )

1: Inputs: x 0 , y 0 , 2: Parameters: K, T , M , γ AddCorrection 3: for k ∈ {1, ..., K + 1} do 4: ỹk ← I M (x k-1 , y k-1 ). # Warm-start. 5: y k ←ϕ T (x k-1 , ỹk ) # Unrolled optimization. 6: d k ← ∂ x L T (x k-1 , ỹk ) 7: if AddCorrection= True then 8: v k ← ∂ y L T (x k-1 , ỹk ) 9: ξ k ≈ -(∂ yy g(x k-1 , y k )) † v k # Approx. solver 10: d k ← d k +∂ xy g(x k-1 , y k )ξ k # Grad. correction 11:
end if 12:

x k ← x k-1 -γd k # Updating x 13: end for 14: Return (x K , y K ). Assumption 4. For any non-negative integers M, T ≥ 0, the maps I M and ϕ T are continuous on X × Y and take values in Y, with ϕ T being continuously differentiable. Moreover, for any (x, y) ∈ X × Y s.t. ∂ y g(x, y)=0 and M, T ≥ 0, there exists a matrix D such that:

I M (x, y) = ϕ T (x, y) = y, ∂ x ϕ T (x, y) = ∂ 2 xy g(x, y)D, ∂ y ϕ T (x, y) = I + ∂ 2 yy g(x, y)D.
Finally, for any (x, y)∈X ×Y, and M, T ≥ 0 s.t. T + M > 0, the equality y=ϕ T (x, I M (x, y)) implies that y is a critical point of g, i.e. ∂ y g(x, y) = 0. Assumption 5. ϕ T converges to a selection φ and ∂ x ϕ T converges uniformly near local minima.

Assumption 4 is satisfied by many mappings used in practice such as T -steps of the gradient descent or proximal point algorithms, whenever g is twice-continuousely differentiable and L-smooth as shown in Proposition 19 of Appendix D. Assumption 5 is a discrete-time version of the uniform convergence result in Proposition 17 of Appendix C but that we directly assume here for simplicity. Under these assumptions we show that Algorithm 1 can find equilibria of (BGS) up to an approximation error resulting from the fact that ϕ T is not an exact selection. Proposition 3. Let M, T be non-negative numbers s.t. M + T > 0 and let (x k , y k ) be the iterates of Algorithm 1 using the maps I M and ϕ T and without any correction, i.e. AddCorrection=False.

If (x k , y k ) converges to a limit point (x T , y T ) then, under Assumption 4:

∂ x L T (x T , y T ) = 0, ∂ y g(x T , y T ) = 0.
Let E be the set of limit points (x T , y T ) for T ≥ 0. If E is bounded and y T is a local minimum of g(x T , .) for any T ≥ 0, then, under Assumptions 4 and 5, the elements of E are approximate equilibria for (BGS):

lim sup T ∂ x L φ (x T , y T ) = 0, ∂ y g(x T , y T ) = 0, (∀T > 0).
Proposition 3 shows that unrolled optimization algorithms approximately solve (BGS) in the limit where the number of unrolling steps T of the ϕ T goes to infinity. This result is consistent with the ones obtained in [START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF] for the case where g is strongly convex and illustrates the high computational cost for solving (BGS) without correcting for the bias introduced by unrolling. Next, we show how to get rid of such a bias in light of Theorem 2.

Implicit Gradient Correction

We propose to correct the bias of unrolling by exploiting the expression of the gradient ∂ x φ provided in Theorem 2. The key idea is to obtain an expression for ∂ x L φ (x, y) in terms of L T and the second-order derivatives of g which holds for any local minimizer y of y → g(x, y) as shown by the proposition below. Proposition 4. Let φ be the selection defined in Section 4.2 and (x, y) ∈ X × Y be s.t. y is a local minimum of y → g(x, y). Then, under Assumptions 1 to 4, ∂ x L φ (x, y) is given by the equation:

∂ x L φ (x, y) := ∂ x L T (x, y) -∂ 2 xy g(x, y) ∂ 2 yy g(x, y) † ∂ y L T (x, y).
Proposition 4, which is proven in Appendix D, suggests a simple correction for the gradient estimate d k in Algorithm 1. By doing so, the corrected algorithm would be performing an approximate gradient descent on each of the upper-level and lower-level objectives, suggesting that the algorithm may recover equilibrium points of (BGS) without having to increase the computation budget for the unrolling as we show later in Proposition 5. A simple way to proceed would to compute c k satisfying the approximate equation

c k ≈ -B k (A k ) † v k , where A k =∂ 2 yy g(x k-1 , y k ), B k :=∂ 2 xy g(x k-1 , y k ) and v k =∂ y L T (x k-1 , ỹk ).
More concretely, c k can be computed by setting c k =B k ξ k where ξ k approximates the minimum norm solution to the least squares problem:

ξ k ≈ arg min ξ ξ 2 , s.t. ξ ∈ arg min ξ A k ξ + v k 2 , (5)
Approximate solution to [START_REF] Gunes Baydin | Automatic differentiation in machine learning: a survey[END_REF]. It is possible to solve (5) approximately using an iterative procedure by constructing N iterates ξ t starting from ξ 0 = 0 and performing (conjugate) gradient descent on the quadratic objective. This can be implemented efficiently using only Hessian vector products with the Hessian A k [START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF]. The constrained problem (5) can also be expressed as an unconstrained one by re-parametrizing ξ = A k z:

ξ k ≈ A k z k , s.t. z k ∈ arg min z A 2 k z + v k 2 . (6)
Eq. ( 6) has the advantage that z k solves an unconstrained problem. As such, it is more amenable to applying a warm-start strategy, which can yield efficient approximation z k to z k by exploiting previously computed approximation z k-1 to z k-1 [START_REF] Arbel | Amortized implicit differentiation for stochastic bilevel optimization[END_REF]. This strategy can be achieved using a standard iterative algorithm P for approximately solving the least-squares problems, such as a fixed number of conjugate gradient iterations, that takes as input the matrix A k , vector v k and initialization z k-1 ≈ z k-1 and returns the next iterate z k ≈ z k . More formally we view P as a continuous map of (A, v, z) → P(A, v, z) returning a vector z and such that the only fixed points are exact solutions to the least square problem min z A 2 z + v 2 . We refer to Appendix D.1 for examples of such maps.

We can then define the iterates z k and ξ k as follows:

ξ k = A k z k , z k = P(A k , v k , z k-1 ). (7)
The corrected algorithm is obtained by setting the variable AddCorrection=True in Algorithm 1 and computing the ξ k using any approximate solver including, in particular, the ones based on a warm-start strategy as in [START_REF] Blondel | Efficient and modular implicit differentiation[END_REF]. The following proposition, with proof in Appendix D, shows that the proposed correction indeed yields equilibrium points of (BGS). Proposition 5. Let (x k , y k ) be the iterates obtained using Algorithm 1 with AddCorrection=True and T + M > 0 and assume that ξ k are computed using [START_REF] Blondel | Efficient and modular implicit differentiation[END_REF]. If (x k , y k , z k ) k≥0 converges to a limit point (x , y , z ), then y is a critical point of y → g(x , y) and if, in addition, y is a local minimizer, then (x * , y * ) must be an equilibrium of (BGS) satisfying (SC):

∂ x L φ (x , y ) = 0 and ∂ y g(x , y ) = 0
Proposition 5 shows that the proposed correction allows to recover equilibria of (BGS) without having to increase the number of iterations T of the unrolled algorithm. This is by contrast with Proposition 3 where T must increase to infinity, which would be impractical. We discuss in Appendix D.2 how different choices for the parameters T and M recover known algorithms. In particular, that Algorithm 1 with correction allows interpolating between two families of algorithms: (ITD) and (AID) while still recovering the correct equilibria. Numerical results illustrating the benefits of the correction are presented in Appendix E.

Discussion

We have introduced a bilevel game that resolves the ambiguity in bilevel optimization with nonconvex objectives using the notion of selection maps. We have shown that many algorithms for bilevel optimization approximately solve these games up to a bias due to finite computational power.

Our study of the differentiability properties of the selection maps has resulted in practical procedures for correcting such a bias and required the development of new analytical tools. This study opens the way for several avenues of research to understand the tradeoff between unrolling and implicit gradient correction for designing efficient algorithms. In future work, studying these algorithms in a non-smooth and stochastic setting would also be of great theoretical and practical interest.

Funding This project was supported by ANR 3IA MIAI@Grenoble Alpes (ANR-19-P3IA-0003).

Checklist

The checklist follows the references. Please read the checklist guidelines carefully for information on how to answer these questions. For each question, change the default [TODO] to [Yes] , [No] , or [N/A] . You are strongly encouraged to include a justification to your answer, either by referencing the appropriate section of your paper or providing a brief inline description. Please do not modify the questions and only use the provided macros for your answers. Note that the Checklist section does not count towards the page limit. In your paper, please delete this instructions block and only keep the Checklist section heading above along with the questions/answers below.

1. For all authors... We recall the definition of classical Morse-Bott functions [START_REF] David | Morse-bott theory and equivariant cohomology[END_REF][START_REF] Feehan | On the morse-bott property of analytic functions on banach spaces with łojasiewicz exponent one half[END_REF], which we extend in Section 4.1 to the case where there is a dependence on some additional parameter x in X . Definition 3 (Morse-Bott function). Let h : Y → R be a real-valued twice continuousely differentiable function. Define C h to be the set of critical points of h and consider y 0 ∈ C h . We say that h is Morse-Bott at y 0 , if there exists a open neighbordhood V of y 0 such that C h ∩ V is a connected sub-manifold of Y of dimension dim Ker(∂ 2 yy h(y 0 )) . We say that h is a Morse-Bott function if for any y 0 ∈ C h , h is Morse-Bott at y 0 .

Morse-Bott functions were introduced in the context of differential topology to analyze the geometry of a manifold by studying the properties of differentiable functions defined on that manifold [START_REF] David | Morse-bott theory and equivariant cohomology[END_REF]. Their main property is that all their critical points that are connected have the same type (same number of positive and negative eigenvalues for the Hessian), a fact expressed by the Morse-Bott lemma [16, Theorem 2.10] that we generalize to the parametric setting in Theorem 1. Morse-Bott functions form a generic class of functions [START_REF] Martınez-Alfaro | Topological classification of simple morse bott functions on surfaces[END_REF], meaning that any smooth function can always be slightly perturbed to become a smooth Morse-Bott function. Hence, in principle, requiring that y → g(x, y) is a Morse-Bott function for any parameter x ∈ X is essentially a mild assumption. The Morse-Bott property allows characterizing the geometry of critical points of g(x, .) for any x and ensures that the selection map φ is well-defined [START_REF] Ralph L Cohen | Topics in Morse theory[END_REF]Chapter 15]. However, this condition does not provide any information about how the set of critical points evolves as the parameter x varies, which is crucial for the study of smoothness of the selection φ. This is precisely why we introduced parametric Morse-Bott functions in Section 4.1.

A.2 Properties of Parameteric Morse-Bott Functions.

In this section, we describe some elementary properties of parametric Morse-Bott functions. In particular, Proposition 6 shows that ∂ 2 xy g(x, y) belongs to the range of ∂ 2 yy g(x, y) whenever (x, y) is an augmented critical point of g, i.e. ∂ y g(x, y) = 0. Proposition 7 shows that any parametric Morse-Bott function g satisfies a pointwise Morse-Bott property in the sense of Definition 3. Finally, Proposition 8 and Lemma 1 provide examples of functions that satisfy the parametric Morse-Bott property. Recall M the set of augmented critical points of g:

M = {(x, y) ∈ X × Y|∂ y g(x, y) = 0}. ( 8
)
Proposition 6 (Exact least square solution). Let g be a parametric Morse-Bott function. Let (x 0 , y 0 ) be an element in M defined in [START_REF] Bolte | Nonsmooth implicit differentiation for machine-learning and optimization[END_REF] and define the matrices A := ∂ 2 yy g(x 0 , y 0 ) and B := ∂ 2 xy g(x 0 , y 0 ). Then, B is in the range of A, i.e. there exists a matrix U such that B=U A.

Proof. Recall that M is the set of augmented critical points of g. Since g is a parametric Morse-Bott function, there exists a neighborhood U of (x 0 , y 0 ) such that the augmented critical set

M ∩ U is a C 2 manifold of dimension d M = dim(X ) + dim(Ker(∂ 2 yy g(x 0 , y 0 ))). We know that M ∩ U is characterized locally by the equation ∂ y g(x 0 , y 0 ) = 0, hence the tangent space T M (x0,y0) of M∩U at point (x 0 , y 0 ) consist of the set of directions (u, v) ∈ X × Y for which ∂ y g(x 0 + u, y 0 + v) = O( 2 ). In other words T M (x0,y0) is the set of vectors (u, v) ∈ X × Y of M ∩ U satisfying the equation: u ∂ 2 xy g(x 0 , y 0 ) + v ∂ 2 yy g(x 0 , y 0 )=u B + v A=0. Since M ∩ U is of dimension d M , the tangent space T M (x0,y0) must also have dimension d M .
Therefore, by the rank theorem, it must hold that the matrix D = (B, A) has a rank equal to dim(X ) + dim(Y) -d M = rank(A). On the other hand, we know that 0 B + v A = v A ∈ Range(A) for any v ∈ Y, so that Range(A) ⊂ Range(D). The two subspaces having the same dimension, the inclusion implies equality (Range(A) = Range(D)). Henceforth, there must exist a matrix U such that B can be written as B=U A.

Proposition 7 (Pointwise Morse-Bott property). Let g be a parametric Morse-Bott function. Then for any x ∈ X , the function y → g(x, y) is a Morse-Bott function in the following sense: For any x 0 and any critical point y 0 of g(x 0 , .), there exists an open neighborhood V of y 0 so that C x,y0 := {y ∈ Y|∂ y g(x 0 , y) = 0} ∩ V is a connected sub-manifold of dimension equal to the dimension of the null space of the Hessian ∂ 2 yy g(x, y 0 ).

Proof. Let (x 0 , y 0 ) be in X × Y such that ∂ y g(x 0 , y 0 ) = 0. Then, since g is a parameteric Morse-Bott function, there exists a neighborhood U of (x 0 , y 0 ) such that the augmented critical set M ∩ U is a C 2 manifold of dimension d M = dim(X ) + dim(Ker(∂ 2 yy g(x 0 , y 0 ))). On the other hand, we know that M ∩ U is characterized locally by the equation ∂ y g(x, y) = 0, hence the tangent vectors (u, v) ∈ X × Y of M ∩ U at (x 0 , y 0 ) must satisfy the equation:

u ∂ 2
xy g(x 0 , y 0 ) + v ∂ 2 yy g(x 0 , y 0 )=0. For simplicity, we denote by B = ∂ 2 xy g(x 0 , y 0 ) and A = ∂ 2 yy g(x 0 , y 0 ). By Proposition 6, we know that B can be written in the form B = U A for some matrix. Hence, the tangent space of M at

(x 0 , y 0 ) consists in vectors (u, v) ∈ X × Y satisfying u U + v A = 0.
In particular, for any u ∈ X , we can set v = -u U which ensures that (u, v) is in the tangent space of M at (x 0 , y 0 ). Now consider the sub-manifold {x 0 } × Y, its tangent space at (x 0 , y 0 ) is {0} × Y. For any element (x, y) ∈ X × Y, we have the decomposition (x, y)=(x, -x U ) + (0, y + x U ) where the first tuple belongs to the tangent space of M and the second one belongs to the tangent space of {x 0 } × Y at (x 0 , y 0 ). Hence, the tangent space of X × Y is generated by the both separate tangent spaces which means that both manifolds intersect transversally and that

{x 0 } × C := (M ∩ U) ∩ ({x 0 } × Y) is a sub-manifold of dimension dim Ker ∂ 2
yy g(x 0 , y 0 ) [29, Theorem 6.30]. For a small enough open connected neighborhood V of y 0 , we can ensure that C ∩ V is a connected sub-manifold of Y. This precisely means that y → g(x 0 , y) is Morse-Bott at the point y 0 which concludes the proof.

Proposition 8 (Morse functions with parameters). Let g : X × Y be a three-times continuously differentiable function such that for any (x, y) ∈ X ×Y for which ∂ y g(x, y) = 0, the Hessian matrix ∂ 2 yy g(x, y) is invertible. Then g is a parametric Morse-Bott function.

Proof. Let (x 0 , y 0 ) ∈ X × Y be such that y 0 is a critical point of g(x 0 , .) ( i.e. ∂ y g(x 0 , y 0 )=0). Since, by assumption, the Hessian is invertible, we can apply the implicit function theorem which guarantees the existence of a function x → y(x) defined in a neighborhood U of x 0 and taking values in a neighborhood V of y 0 , such that y(x 0 ) = y 0 and y(x) is the unique critical point of g(x, .) on V, i.e.:

∂ y g(x, y(x)) = 0, ∀x ∈ U.
Moreover, x → y(x) is twice continuously differentiable. This ensures that M the set of augmented critical points of g satisfies:

M ∩ (U × V) = {(x, y(x)) ∈ X × Y|x ∈ U} := S.
We only need to show that M∩(U × V) is a manifold of dimension dim(X ). For this, we will apply the regular level set theorem [START_REF] Lee | Introduction to Smooth Manifolds[END_REF]Corollary 5.14] to the function G : (x, y) → ∂ y g(x, y) defined on U × V. The pre-image of 0 by G is exactly equal to M ∩ (U × V). Moreover, for any (x, y) ∈ M ∩ (U × V), we have that dG(x, y) is of maximal rank since ∂ 2 yy g(x, y) is invertible. Hence, by application of the regular level set theorem theorem to the twice continuously differentiable ( Proof. Consider the function G : (x, y) → ∂ y g(x, y). We have the following equivalence

C 2 ) function G, it follows that M ∩ (U × V) = G -1 ({0}) is a C 2 sub-manifold of X × Y of dimension dim(ker(dG(x, y))) = dim(X ). We have shown that M ∩ (U × V) is sub-manifold of dimension dim(X ),
(x, y) ∈ G -1 ({0}) ⇐⇒ ∂ y h(τ x (h))∂ y τ x (y) = 0 ⇐⇒ τ x (y) ∈ ∂ y h -1 ({0}).
Consider the map T : (x, y) → τ x (y), then we have shown that G -1 ({0}) = T -1 ∂ y h -1 ({0}) . Let (x 0 , y 0 ) ∈ X ×Y be an augmented critical point of g. Set ỹ = T (x 0 , y 0 ) which is a critical point of h. Since h is, by assumption, Morse-Bott at ỹ, then there exists an open neighborhood Ṽ of ỹ such that ∂ y h -1 ({0}) ∩ Ṽ is a sub-manifold of dimension dim(Ker(∂ 2 yy h(ỹ))). By continuity of T , we can always find open connected neighborhoods U and V of (x 0 , y 0 ) so that V :=T (U × V) ⊂ Ṽ. Moreover, since for any x T (x, .) is a diffeomorphism, it must be that V is an open set. Therefore, S := ∂ y h -1 ({0}) ∩ V must be a sub-manifold as of dimension dim(Ker(∂ 2 yy h(ỹ))). It remains to show that T -1 (S) is a sub-manifold. To see this, it suffice to note that the differential of T is surjective which ensures that T is transverse to S and that T -1 (S) is a sub-manifold [START_REF] Lee | Introduction to Smooth Manifolds[END_REF]Theorem 6.30]. Moreover, the dimension of such manifold is equal to dim(X ) + dim(ker(∂ 2 yy h(ỹ))) = dim(X ) + dim(ker(∂ 2 yy g(x 0 , y 0 ))).

A.3 Proof of the Morse-Bott Lemma with Parameters

In this section, we provide a proof of the Morse-Bott lemma with parameters introduced in Theorem 1. We then introduces two results in Proposition 9 and Corollary 1 which are consequences of Theorem 1. Proposition 9 shows that near an augmented critical point (x 0 , y 0 ), the Hessian matrices of nearby augmented critical points are all similar. This result illustrates that the geometry near a critical point is preserved when the parameter x is perturbed. Proposition 9 will be used later in Proposition 16 of Appendix C to show that the pseudo-inverse of the Hessian matrices of critical points near a local minimum are uniformly bounded. Finally, Corollary 1 shows that near any augmented critical point (x 0 , y 0 ) the function y → g(x, y) can be expressed as a slight deformation of y → g(x 0 , y). This result, along with the stability result in Appendix B.2 of the gradient flow to deformations will be key to prove the continuity of the selection map x → φ(x, y) near local minima.

Proof of Theorem 1. Let x 0 ∈ X and y 0 be a critical point of g(x 0 , .). Denote by K the null space of the Hessian A 0 = ∂ 2 yy g(x 0 , y 0 ) and by K ⊥ its orthogonal complement in Y. where J 0 is an invertible diagonal matrix whose diagonal elements are equal to the sign of the non-zero eigenvalues of the Hessian ∂ 2 yy g(x 0 , y 0 ). By convention J 0 =0 in case the Hessian ∂ 2 yy g(x 0 , y 0 ) = 0. Since, the function h(x, r, w) := g(x, s(r, w)) is such that ∂ w h(x 0 , 0, 0) = 0 and the partial Hessian ∂ 2 ww h(x 0 , 0, 0) = J 0 is invertible, we are in position to apply the Morse lemma with parameters [16, Theorem 4]. The lemma ensures that O and O ⊥ can be chosen small enough so that there exits open neighborhoods B and O ⊥ 1 of x 0 ∈ X and 0 ∈ K ⊥ and a diffeomorphism τ from B × O × O ⊥ 1 to B × O × O ⊥ such that τ (x 0 , 0, 0)=(x 0 , 0, 0) and decomposing h locally into a quadratic component and a singular one. More precisely, for any (x, r, w) ∈ B × O × O ⊥ 1 , the map r satisfies τ (x, r, w)=(x, r, w ) for some w ∈ O ⊥ and the following equation holds:

h(τ (x, r, w)) = h(τ (x, r, 0)) + 1 2 w J 0 w. (9)
It remains to show that ξ → h(τ (x, r, 0)) is in fact constant for (x, r) in an open neighborhood of (x 0 , 0) ∈ X × O. To this end, define the sets A, B and C as follows:

A := {(x, y) ∈ B × V | ∂ y g(x, y) = 0}, B := {(x, r, w) | (x, r) ∈ B × O, ∂ r h(τ (x, r, w)) = 0}, C := {(x, r, 0) | (x, r) ∈ B × O, ∂ r h(τ (x, r, 0)) = 0}.
Then by [START_REF] Bolte | Nonsmooth automatic differentiation: a cheap gradient principle and other complexity results[END_REF], it holds that B = C. Moreover, A and B are homeomorphic. Indeed to see this, we introduce the notation s(x, r, w) := (x, s(r, w)) which defines a diffeomorphism from

B × O × O ⊥ to B × V. Hence, g • s • τ = h • τ . This ensures s • τ (B)=A,

which means precisely that

A and B are homeomorphic since s • τ is a homeomorphism. Moreover, by definition of g as a parametric Morse-Bott function, we also know that A is a sub-manifold of X × Y of dimension dim(X )+dim Ker ∂ 2 yy g(x 0 , y 0 ) provided the neighborhoods B and V are small enough. Hence, we can deduce that B and C must also be sub-manifolds of the same dimension. In particular, C is a sub-manifold of B × O × {0} which is of dimension dim(X ) + dim Ker ∂ 2 yy g(x 0 , y 0 ) . Therefore, C is an open sub-manifold of B × O × {0}. Hence, since (x 0 , 0, 0) ∈ C, there must exists an open connected neighborhood B 1 × O 1 × {0} of (x 0 , 0, 0) in B × O × {0} that is contained in C. Hence, we deduce that for any (x, r) ∈ B 1 × O 1 , the function h satisfies ∂ r h(τ (x, r, 0)) = 0 so that h(τ (x, r, 0)) = h(τ (x, 0, 0)) on such neighborhood. Finally, we have shown that there exits

g • s • τ (x, r, w) = g • s • τ (x, 0, 0) + 1 2 w J 0 w.
We conclude the proof by setting ψ(x, r, w) = s • τ (x, r, w) which is the desired diffeomorphism.

Proposition 9. Let g be a real-valued function such that Assumption 1 holds. Consider an augmented critical point (x 0 , y 0 )∈M, with M defined in [START_REF] Bolte | Nonsmooth implicit differentiation for machine-learning and optimization[END_REF]. Then there exists a neighborhood V of (x 0 , y 0 ) and a continuous map (x, y) → P (x, y) defined on V with values in R d×d such that:

• P (x, y) is invertible for any (x, y) ∈ V with singular values contained in an interval [σ min , σ max ] for some positive constants σ min and σ max .

• For any augmented critical point (x, y) ∈ V, the Hessian of g is given by: ∂ 2 yy g(x, y) = P (x, y) ∂ 2 yy g(x 0 , y 0 )P (x, y).

Proof. Denote by K the null space of the Hessian A 0 = ∂ 2 yy g(x 0 , y 0 ) and by K ⊥ its orthogonal complement in Y. Let J 0 be a diagonal matrix with diagonal elements given by the sign of the non-zero eigenvalues of A 0 . Since g satisfies Assumption 1, we apply Theorem 1 which ensures the existence of a diffeomorphism ψ defined on an open neighborhood U of (x 0 , 0, 0) ∈ X × K × K ⊥ with values in an open neighborhood V of (x 0 , y 0 ) in X × Y, s.t. ψ(x 0 , 0, 0) = (x 0 , y 0 ) and for all (x, r, w) ∈ U, ψ satisfies ψ(x, r, w)=(x, y) and g(ψ(x, r, w)) = g(ψ(x, 0, 0)) + 1 2 w J 0 w,

= g(ψ(x, 0, 0)) + 1 2 (r w ) J0 r w ,

where we defined J0 to be the matrix of dimension d × d given by:

J0 = 0 0 0 J 0 .
Since ψ is a diffeomorphism satisfying ψ(x, r, w) = (x, y), we can equivalently write (10) as:

g(x, y) = g(ψ(x, 0, 0)) + 1 2 ψ -1 2,3 (x, y) J0 ψ -1 2,3 (x, y), ∀(x, y) ∈ V, (11)
where ψ -1 2,3 (x, y) are last two components of ψ -1 (x, y) (i.e. ψ -1 (x, y) = (x, ψ 1 2,3 (x, y))). By differentiating (11) w.r.t. y we obtain:

∂ y g(x, y) = ∂ y ψ -1 2,3 (x, y) J0 ψ -1 2,3 (x, y). (12)
∂ y ψ -1 2,3 (x, y) must be invertible since (∂ x ψ -1 , ∂ y ψ -1 ) is invertible and of the form:

∂ x ψ -1 ∂ y ψ -1 = I ∂ x ψ -1 2,3 0 ∂ y ψ -1 2,3
.

Therefore, if y is a critical point of g(x, .), then [START_REF] Daneri | Lecture notes on gradient flows and optimal transport[END_REF] implies that J0 ψ -1 2,3 (x, y) = 0. Let (x, y) be an augmented critical point of g, > 0 and u be vector in Y, then the following holds:

1 ∂ y g(x, y + u) = ∂ y ψ -1 2,3 (x, y + u) J0 1 ψ -1 3 (x, y + u) -ψ -1 2,3 (x, y) .
Hence, by taking the limit when approaches 0, it follows that:

∂ 2 yy g(x, y) = ∂ y ψ -1 2,3 (x, y) J0 ∂ y ψ -1 2,3 (x, y).
Define P 0 := ∂ y ψ -1 2,3 (x 0 , y 0 )∂ y ψ -1 2,3 (x 0 , y 0 ) which is invertible. Then, we can write:

∂ 2 yy g(x, y) = ∂ y ψ -1 2,3 (x, y) P -1 0 P 0 J0 P 0 P -1 0 ∂ y ψ -1 2,3 (x, y) = ∂ y ψ -1 2,3 (x, y) P -1 0 ∂ y ψ -1 2,3 (x 0 , y 0 )A 0 ∂ y ψ -1 2,3 (x 0 , y 0 ) P -1 0 ∂ y ψ -1 2,3 (x, y) = P (x, y) ∂ 2
yy g(x 0 , y 0 )P (x, y),

where we defined P (x, y)=∂ y ψ -1 2,3 (x 0 , y 0 ) P -1 0 ∂ y ψ -1 2,3 (x, y). The matrix P (x, y) is invertible for any (x, y) ∈ V and the map (x, y) → P (x, y) is continuous. Hence, by considering compact neighborhood of (x 0 , y 0 ) contained in V, we can ensure that the singular values of P (x, y) are contained in an interval [σ min , σ max ] where σ min and σ max are positive numbers. Further considering the restriction of such map on an open neighborhood V ⊂ K of (x 0 , y 0 ) yields the desired result.

Corollary 1. Let g be a real-valued function such that Assumption 1 holds. Consider an augmented critical point (x 0 , y 0 )∈M, with M defined in [START_REF] Bolte | Nonsmooth implicit differentiation for machine-learning and optimization[END_REF]. Then, there exists a open neighborhoods B and V of x 0 and y 0 in X and Y and a continuously differentiable map τ from B × V to V such that:

• For any x ∈ B, the map τ x : y → τ (x, y) is a diffeomorphism from V to itself satisfying τ x0 (y) = y for any y ∈ V. Moreover, (x, y) → τ -1 x (y) is continuous. • For any (x, y) ∈ B × V, the function g satisfies g(x, y)=g(x 0 , τ (x, y)) + C(x), where

x → C(x) is a function independent of y.

• There exists positive numbers and L s.t for any (x, y) ∈ B × V:

2 I ≤ ∂ y τ (x, y) ∂ y τ (x, y) ≤ (L ) 2 I. (13)
Proof. We use the notations of Theorem 1 where K is the null subspace of the Hessian ∂2 yy g(x 0 , y 0 ) and K ⊥ its orthogonal complement in Y. By Theorem 1 g satisfies: g(ψ(x, r, w)) = g(ψ(x, 0, 0)) + 1 2 ȳ J 0 ȳ, with ψ and J 0 being the diffeomorphism and matrix defined in Theorem 1. Recall that ψ is defined on an open neighborhood B × O × O ⊥ of (x 0 , 0, 0) ∈ X × K × K ⊥ and whose image by ψ is an open neighborhood B × V of (x 0 , y 0 ). Hence, we can write:

g(ψ(x, r, w)) = C(x) + g(ψ(x 0 , r, w)),
with C(x) := g(ψ(x, 0, 0)) -g(ψ(x 0 , 0, 0)). We also know that ψ preserves x, meaning that ψ(x, r, w)=(x, y). Hence, we can define (x, r, w) → τx (r, w) ∈ Y, s.t. ψ(x, r, w)=(x, τx (r, w)).

For any x ∈ B, (r, w) → τx (r, w) defines a diffeomorphism from O×O ⊥ onto its image. Moreover, its image must be equal to V. Indeed, since ψ(B ×O ×O ⊥ ) = B ×V, it follows that for any (x, y) ∈ B × V, there exists (r, w) ∈ O × O ⊥ such that ψ(x, r, w) = (x, τx (r, w)) = (x, y). In particular, if (x, y) ∈ B × V and (r, w) = τ -1 x (y), we can write ψ(x 0 , r, w)=(x 0 , τx0 (r, w))=(x 0 , τx0 τ -1 x (y)). Therefore, the following expression holds for any (x, y) ∈ B × V:

g(x, y) = C(x) + g(x 0 , τ (x, y)),
where we defined τ (x, y)=τ x0 • τ -1

x (y). For any x ∈ B, the map τ x : y → τ (x, y) is a diffeomorphism satisfying τ (x 0 , y)=y. Moreover, (x, r, w) → τx (r, w) and (x, y) → τ -1

x (y) are continuously differentiable since ψ is a diffeomorphism. As a result, τ is continuously differentiable as well and (x, y) → τ -1

x (y) is continuously differentiable. Finally, since ∂ y τ (x, y) is jointly continuous in x and y and ∂ y ψ x (y) is invertible, then, provided that B and V are small enough, there must exist two positive numbers and L such that for any (x, y) ∈ B × V:

Proof of Proposition 1 . Recall M = {(x, y) ∈ X × Y|∂ y g(x, y) = 0} the set of augmented critical points of g and let (x 0 , y 0 ) be in M. First, since Assumption 1 holds, we know by Proposition 7 that g(x 0 , .) is a Morse-Bott function. Hence, by [START_REF] Feehan | On the morse-bott property of analytic functions on banach spaces with łojasiewicz exponent one half[END_REF]Theorem 1], it follows that g(x 0 , .) satisfies a Łojasiewicz inequality near y 0 . In other words, there exists a neighborhood V of y 0 and a positive constant µ > 0 such that:

µ |g(x 0 , y) -g(x 0 , y 0 )| ≤ 1 2 ∂ y g(x 0 , y) 2 , ∀y ∈ V.
By Corollary 1, there exists a continuous function τ defined on an open neighborhood B × V of (x 0 , y 0 ) whose image is V and for which g(x, y)=g(x 0 , τ (x, y))+C(x) for any (x, y) ∈ B × V, where C(x) is a function of x independent of y. Moreover, for any x ∈ B, y → τ (x, y) is a diffeomorphism from V to itself whose inverse is written as τ -1 (x, y) by an abuse of notion. In particular, for y=τ -1 (x, y 0 ) we set G(x):=g(x, τ -1 (x, y 0 )) = g(x 0 , y 0 ) + C(x). Note that τ -1 (x, y 0 ) is critical point of g(x, .) since ∂ y g(x, τ -1 (x, y 0 ))∂ y τ -1 (x, y 0 )=∂ y g(x 0 , y 0 ) = 0 and ∂ y τ -1 (x, y 0 ) is invertible. Hence, the following holds for any (x, y) ∈ B × V.

µ |g(x, y) -G(x)| = µ |g(x 0 , τ (x, y)) -g(x 0 , y 0 )| ≤ 1 2 ∂ y g(x 0 , τ (x, y) 2 . ( 14
)
Moreover, by construction of τ , we know that ∂ y τ (x, y) satisfies ( 13) for any (x, y) ∈ B × V. Therefore, we deduce that:

∂ y g(x, y) 2 = ∂ y g(x 0 , τ (x, y))∂ y τ (x, y) 2 ≥ 2 ∂ y g(x 0 , τ (x, y)) 2 ,
Finally, combining the above inequality with ( 14), we get that, for any (x, y) ∈ B × V:

2 µ |g(x, y) -G(x)| ≤ 1 2 ∂ y g(x, y) 2 , ∀(x, y) ∈ U.
The result follows by setting µ = 2 µ > 0 and U = B × V.

B Asymptotic Properties of Gradient Flows

B.1 Convergence of the gradient flow.

Recall that the gradient flow φ t (x, y) satisfies the differential equation

dφ t (x, y) dt = -∂ y g(x, φ t (x, y)), φ 0 (x, y) = y.
The next proposition shows that the gradient flow φ t (x, y) converges towards a well-defined selection map φ(x, y). Proposition 10 (Convergence of φ t .). Let x, y be in X × Y. Under Assumptions 1 to 3, (t, x, y) → φ t (x, y) is continuous and for any (x, y) ∈ X × Y, φ t (x, y) converges towards a unique critical point φ(x, y) of y → g(x, y) as t goes to +∞.

Proof. First, Assumption 2 ensures that the gradient flow φ t (x, y) is uniquely defined at all times t [START_REF] Daneri | Lecture notes on gradient flows and optimal transport[END_REF]. φ t (x, y) is jointly continuous in (t, x, y) by Cauchy-Lipschitz theorem. Moreover, t → φ t (x, y) remains bounded thanks to Assumption 3. Otherwise, there exists a subsequence φ tn (x, y) such that g(x, φ tn (x, y)) diverges to +∞. This contradicts the fact that g(x, φ tn (x, y)) is decreasing since φ t (x, y) is a gradient flow of g. Hence, we deduce that φ t (x, y) must have at least one accumulation point y . Moreover, y must be a critical point of g(x, .). To see this, note that g(x, φ t (x, y)) is a decreasing function in time and is lower-bounded. Hence, it admits a finite limit l. Moreover, by differentiating g(x, φ t (x, y)) is time, it follows that:

d dt g(x, φ t (x, y)) = -∂ y g(x, φ t (x, y)) 2
This implies that +∞ 0 ∂ y g(x, φ s (x, y)) 2 ds = g(x, y) -l is finite. Since, g is L-smooth by Assumption 2, this is only possible if ∂ y g(x, φ s (x, y)) converges to 0. In particular, by continuity of ∂ y g(x, y), it follows that ∂ y g(x, y ) = 0. We only need to show that y is the unique accumulation point of φ t (x, y). To show this, we apply Proposition 1, which implies, in particular, that g satisfies a Łojasiewicz inequality in a neighborhood V of y :

µ|g(x, y) -G(x)| ≤ ∂ y g(x, y) 2 , ∀y ∈ V.
We can therefore apply [START_REF] Merlet | Convergence to equilibrium for discretizations of gradient-like flows on riemannian manifolds[END_REF]Theorem 2.7] which ensure that y is the unique accumulation point of φ t (x, y) and that φ t (x, y) converges towards y . We can therefore defined the map φ(x, y) = lim t→∞ φ t (x, y) which constitues a selection.

B.2 Stability of the Gradient Flow Near Local Minima

In this section, we provide a general result establishing the stability of gradient flows to perturbations. This result shows that deforming a gradient flow by a family of diffeomorphisms yields trajectories that are not too far from the unperturbed flow. We will use this result later in Appendix B.3 in conjunction with the formulation of y → g(x, y) as a perturbation of y → g(x 0 , y) provided in Corollary 1 to prove that the gradient flow φ t (x, y) remain stable as the parameter x varies. Proposition 11 (Stability near local minima). Let h be a real valued differentiable function defined on Y and y 0 be a local minimizer of h. We assume that h satisfies the Łojasiewicz inequality near y 0 , meaning that there exists µ > 0 and R > 0 s.t.:

µ(h(y) -h(y 0 )) ≤ 1 2 ∂ y h(y) 2 , ∀y ∈ B(y 0 , R). ( 15
)
Let V be an open neighborhood of y 0 , R > 0 such that B(y 0 , 2R ) ⊂ V and P a family of diffeomorphisms defined from V to itself and satisfying:

1. For any τ ∈P, the pre-image y ψ :=ψ -1 (y 0 ) of y 0 by ψ belongs to B(y 0 , R ).

2. There exists positive numbers and L s.t. for any τ ∈ P and any y ∈ V:

2 I ≤ ∂ y τ (y) ∂ y τ (y) ≤ (L ) 2 I. ( 16
)
For some τ ∈ P, consider a maximal solution (z t ) of the following ODE:

żt = -∂ y h(τ (z t ))∂ z τ (z t ), z 0 ∈ B(y τ , R ). ( 17
)
Then, there exists 0 < C ≤ R , such that for any 0 < ≤ C, there exists 0 < η ≤ 2 with the following property:

For any τ ∈ P and any z 0 s.t. z 0 -y τ ≤ η:

1. The solution z t to ( 17) is well-defined at all times t ≥ 0.

2. For all t ≥ 0, it holds that z t -y τ ≤ .

Proof. The proof is inspired from the the abstract stability result in [START_REF] Li | Convergence and stability of generalized gradient systems by łojasiewicz inequality with application in continuum kuramoto model[END_REF]. We know that y 0 is a local minimizer of h, therefore there exists R" > 0 such that for any y satisfying y -y 0 ≤ R", it holds that L(y) := h(y) -h(y 0 ) ≥ 0. Moreover, by [START_REF] Draxler | Essentially no barriers in neural network energy landscape[END_REF], we also have that:

2µL(y) ≤ ∂ y L(y) 2 , ∀y ∈ B(y 0 , R). ( 18
)
Take < 1 L min(R, R", L R ) := C. To simplify subsequent calculations, we will choose y close enough to y 0 so that 2 -1 2 µ L(y)

1 2
≤ , where µ is the positive constant appearing in [START_REF] Gould | On differentiating parameterized argmin and argmax problems with application to bi-level optimization[END_REF] and is the positive constant in [START_REF] Feehan | On the morse-bott property of analytic functions on banach spaces with łojasiewicz exponent one half[END_REF]. This is possible by continuity of L which that there exists 0 < η ≤ 2 for which any y ∈ B(y 0 , L η) satisfies:

L(y) 1 2 ≤ 1 2 µ 2 . ( 19
)

B.3 Continuity of the Flow Selection

Proposition 12 shows that x → φ(x, y) is continuous at x 0 whenever φ(x 0 , y) is a local minimum of g(x 0 , .). Proposition 13 shows that, near x 0 , φ(x, y) are local minima as well provided φ(x 0 , y) is a local minimum of g(x 0 , .). Proposition 12 (Continuity near local minima). Let x 0 ∈ X and y ∈ Y. Let g be such that Assumptions 1 to 3 hold. Assume that y 0 = φ(x 0 , y) is a local minimizer of y → g(x 0 , y). Then, for any > 0 small enough, there exists T > 0 and η > 0, s.t.:

φ t (x, y) -φ(x 0 , y) ≤ , ∀t ≥ T, ∀x ∈ B(x 0 , η).
In particular, x → φ(x, y) is continuous at φ(x 0 , y).

Proof. We will apply Proposition 11 to the function h(y) = g(x 0 , y) and the well-chosen family P of local diffeomorphisms on Y. By application of Corollary 1, there exists a open neighborhoods B and V of x 0 and y 0 in X and Y and a continuously differentiable map τ from B × V to V such that y → τ (x, y) is a diffeomorphism from V onto itself and for which g satisfies for any (x, y) ∈ B × V:

g(x, y) = g(x 0 , τ (x, y)) + C(x).
For simplicity, we write τ x : y → τ (x, y) by an abuse of notations. We know, by Corollary 1, that x → τ -1 x (y 0 ) is continuous and converges to τ -1 x0 (y 0 ) = y 0 . Hence, by restricting x to a smaller neighborhood B ⊂ B, we can ensure that τ -1

x (y 0 ) belongs to B(y 0 , R ) with R small enough so that B(y 0 , 2R ) ⊂ V. Consider now the family of diffeomorphisms

P P = {V y → τ (x, y) ∈ V|x ∈ B }.
We have constructed P satisfying the conditions of Proposition 11. Moreover, by Proposition 1, the function h(y):=g(x 0 , y) satisfies a Łojasiewicz inequality in an open neighborhood V of y 0 :

µ|h(y) -G(x 0 )| ≤ ∂ y h(y) 2 , ∀y ∈ V .
We can always choose the neighborhood V to be an open ball B(y 0 , R) of radius R > 0 centered in y 0 . Therefore, we have shown so far that h and P satisfy the conditions of Proposition 11.

For any τ ∈ P, consider the ODE:

z t = -∂ y g(x, τ (z t )), z 0 ∈ B(y 0 , R ).
Following the notation in Proposition 11, we define y τ := τ -1 (y 0 ) for any τ ∈ P. We apply Proposition 11 which ensures stability of z t . More precisely, there exists a positive constant C smaller than R so that for any 0 < < C, the solution z t is well-defined at all times and satisfies z t -y τ ≤ for any t ≥ 0, provided that the initial condition z 0 satisfies z 0 -y τ ≤ η for some positive η < 2 that is independent of the choice of τ is P:

∀τ ∈ P : z 0 -y τ ≤ η =⇒ z t -y τ ≤ . (22)
We will apply this result to a particular choice for z 0 . From now on, we fix 0 < < C and let 0 < η ≤ 2 be as in Proposition 11. Using Proposition 10, we know that φ t (x 0 , y) converges to y 0 = φ(x 0 , y), hence there exits T > 0 s.t. φ T (x 0 , y) -y 0 ≤ η 3 . Moreover, since the maps x → φ T (x, y) and x → y τx are continuous at x 0 with y τx 0 =y 0 , there exits η satisfying 0 < η such that B(x 0 , η ) ⊂ B and φ T (x, y) -φ T (x 0 , y) ≤ η 3 and y τx -y 0 ≤ η 3 for any x ∈ B(x 0 , η ). Therefore:

φ T (x, y) -y τx ≤ φ T (x, y) -φ T (x 0 , y) + φ T (x 0 , y) -y 0 + y 0 -y τx ≤ η.
For any x ∈ B(x 0 , η ), by choosing z 0 = φ T (x, y), we have that z 0 -y τx ≤ η. Therefore, we deduce by ( 22) that z t -y τx ≤ and subsequently that

z t -y 0 ≤ z t -y τx + y τx -y 0 ≤ + η 3 ≤ 7 6 ,
since we imposed that η < 2 . Recall now that z t satisfies the ODE:

żt = -∂ y g(x 0 , τ x (z t ))∂ z τ x (z t ).
By definition of τ x , we have g(x, y) = g(x 0 , τ x (y)) for any y ∈ B(y 0 , 2R ) ⊂ V. In particular, as we have shown that z t -y 0 ≤ 7 6 < 2R , it follows that z t satisfies the ODE:

żt = -∂ y g(x, z t ) = -∂ y g(x 0 , τ x (z t ))∂ z τ x (z t ).
By Cauchy-Lipschtz theorem, the solution of the above ODE is unique. Moreover, since we know that φ T +t (x, y) is a solution to the above ODE, then we deduce that z t = φ T +t (x, y). We have shown that for any < C, there exists T > 0 and η such that:

φ t (x, y) -y 0 ≤ 7 6 , ∀t ≥ T, ∀x ∈ B(x 0 , η ). (23)
Since φ t (x, y) converges towards φ(x, y) by Proposition 10, taking the limit t→∞ in [START_REF] Hong | A two-timescale framework for bilevel optimization: Complexity analysis and application to actor-critic[END_REF], we obtain:

φ(x, y) -y 0 ≤ 7 6 , ∀x ∈ B(x 0 , η ).
The above inequality imply in particular that x → φ(x, y) is continuous at x 0 .

Proposition 13 (Stability of local minimizers). Let x 0 ∈ X and y ∈ Y and g be such that Assumptions 1 to 3 hold. Assume that y 0 = φ(x 0 , y) is a local minimizer of y → g(x 0 , y). Then, for any x 1 in a neighborhood of x 0 , φ(x 1 , y) is a local minimizer of g(x 1 , .).

Proof. By assumption, y 0 := φ(x 0 , y) is a local minimizer of g(x 0 , .) ensuring that ∂ 2 yy g(x 0 , y 0 ) is positive semi-definite. Moreover, by Corollary 1, there exists a neighborhood B(x 0 , η)×B(y 0 , 2R ) of (x 0 , y 0 ) such that for any augmented critical point (x 1 , y 1 ) ∈ M ∩ B(x 0 , η) × B(y 0 , 2R ), the Hessian ∂ 2 yy g(x 1 , y 1 ) is similar to ∂ 2 yy g(x 0 , y 0 ). Hence, for any (x 1 , y 1 ) ∈ M ∩ B(x 0 , η) × B(y 0 , 2R ), ∂ 2 yy g(x 1 , y 1 ) must be positive semi-definite so that y 1 is a local minimizer of g(x 1 , .). We can then apply Proposition 12 which ensures that x → φ(x, y) is continuous at x 0 . Therefore, there exists η < η so that, for any x 1 ∈ B(x 0 , η ), y 1 := φ(x 1 , y) belongs to B(y 0 , 2R ). As a result, y 1 must be a local minimizer of g(x 1 , .) since the augmented critical point (x 1 , y 1 ) belongs to M ∩ B(x 0 , η) × B(y 0 , 2R ).

B.4 Uniform Convergence of the Gradient Flow

The result bellow shows that the gradient flow φ t (x, y) converges locally uniformly in x near x 0 at an exponential rate, whenever φ(x 0 , y) is a local minimum. It relies on the locally uniform convergence result in Proposition 12 and the locally uniform Łojasiewicz inequality in Proposition 1. Proposition 14. Let x 0 ∈ X and y ∈ Y and g be such that Assumptions 1 to 3 hold. Assume that y 0 := φ(x 0 , y) is a local minimum. Then there exists positive constants η, T , µ and C such that:

φ t (x, y) -φ(x, y) ≤ Ce -tµ , ∀t ≥ T, x ∈ B(x 0 , η).
A fortiori, φ(x, y) is continuous on B(x 0 , η).

Proof. Proposition 1 ensures the existence of > 0 and η > 0 be such that the following inequality holds:

µ|g(x, y ) -G(x)| ≤ 1 2 ∂ y g(x, y ) 2 , ∀x, y ∈ B(x 0 , η ) × B(y 0 , ). (24)
By Proposition 12 and for >0 small enough, there exists T >0 and η > η>0 for which:

φ t (x, y) -y 0 ≤ , ∀t ≥ T, ∀x ∈ B(x 0 , η ).
Therefore, choosing y =φ t (x, y) in (24) implies:

µ|g(x, φ t (x, y)) -G(x)| ≤ 1 2 ∂ y g(x, φ t (x, y)) 2 , ∀x ∈ B(x 0 , η). (25)
Note that G(x) is the common value of g(x, y) when y is a critical point of g(x, .) in B(y 0 , ). In particular, since φ(x, y) ∈ B(y 0 , ), it holds that G(x) = g(x, φ(x, y)). Moreover, by Proposition 13, φ(x, y) is a local minimum for y → g(x, y). Hence, we must have g(x, φ t (x, y)) -G(x) ≥ 0. We may assume that the inequality is strict otherwise the φ t (x, y) would be a fixed point and we would have φ t (x, y) = φ(x, y). The following inequality holds for any t ≥ T :

φ t (x, y) -φ(x, y) ≤ +∞ t ∂ y g(x, φ s (x, y)) ds ≤ - 2 µ +∞ t Ḣ(s) ds = 2 µ H(t).
where we introduced H(t) = (g(x, φ t (x, y)) -G(x))

1 2 . Thus, we only need to study the evolution of H(t) in time. Computing the derivatives of H(t) and using the inequality in [START_REF] Ji | Bilevel optimization: Convergence analysis and enhanced design[END_REF] yields

Ḣ(t) = - 1 2 H(t) -1 ∂ y g(x, φ t (x, y)) 2 ≤ -µH(t).
By integrating the above inequality, it follows that H(t) ≤ H(T )e -µ(t-T ) . Moreover, using the smoothness of y → g(x, y), we know that

H(T ) ≤ L 2 φ T (x, y) -φ(x, y) ≤ L 2 , ∀x ∈ B(x 0 , η).
Finally, we have shown that φ t (x, y)-φ(x, y) ≤ L µ e -(t-T )µ for any x ∈ B(x 0 , η) and t≥T . Since φ t are continuous in x and converge uniformly in x on B(x 0 , η), then their limit must be continuous on B(x 0 , η).

C Differentiability of the Flow Selection

In this section, we study the differentiability of x →φ(x, y) through the evolution of ∂ x φ t (x, y). The following result establishes that ∂ x φ t (x, y) is well-defined and satisfies a linear differential equation. Proposition 15. Assume g is twice continuously differentiable and satisfies Assumption 2. Then, (x, t) → φ t (x, y) is continuously differentiable with ∂ x φ t (x, y) := U t (x, y) satisfying the differential equation:

Ut (x, y) = -B t (x, y) -A t (x, y)U t (x, y), (26)
where B t and A t are given by:

B t (x, y) = ∂ 2 xy g(x, φ t (x, y)), A t (x, y) = ∂ 2 yy g(x, φ t (x, y)).
Proof. The differentiability of the flow φ t (x, y) in x follows by the application of Cauchy-Lipschitz theorem. It suffices to differentiate the equation defining the flow w.r.t. to obtain [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF].

Note that, by Proposition 10 and continuity of ∂ 2 xy g(x, y) and ∂ 2 yy g(x, y), the matrices A t (x, y) and B t (x, y) must converge to the following matrices A ∞ and B ∞ for any (x, y) ∈ X × Y:

A ∞ (x, y) := ∂ 2 yy g(x, φ(x, y)), B ∞ (x, y) := ∂ 2 xy g(x, φ(x, y))
The following proposition shows that the pseudo-inverse of A ∞ (x, y) remain bounded near x 0 provided that φ(x 0 , y) is a local minimum. Proposition 16. Let (x 0 , y) be in X × Y and set y 0 and g be such that Assumptions 1 to 3 hold. Assume that y 0 :=φ(x 0 , y) is a local minimum of g(x 0 , .). Then there exists an open neighborhood U of x 0 and a positive constant λ > 0 such that:

λ A ∞ (x, y) op ≤ 1, ∀x ∈ U,
where A ∞ (x, y) = ∂ 2 yy g(x, φ(x, y)).

Proof. We apply Proposition 9 which ensures the existence of an open neighborhood V of (x 0 , y 0 ) for which:

∂ 2 yy g(x, y) = P (x, y) ∂ 2 yy g(x 0 , y 0 )P (x, y).

where (x, y) → P (x, y) is continuous map with values in R d×d , and P (x, y) is invertible for any (x, y) ∈ V with singular values in [σ min , σ max ] for σ min > 0 and σ max < +∞. Moreover, since y 0 := φ(x 0 , y) is a local minimum of g(x 0 , .), we know, by Proposition 12, that x → φ(x, y) is continuous at x 0 . Hence, there exists a neighborhood U of x 0 for which (x, φ(x, y)) ∈ V for any x ∈ U. Therefore, it follows that:

A ∞ (x, y) = P (x, φ(x, y)) ∂ 2 yy g(x 0 , y 0 )P (x, φ(x, y)), ∀x ∈ U.

In particular, it follows that:

A ∞ (x, y) † = P (x, φ(x, y)) -1 ∂ 2 yy g(x 0 , y 0 ) † P (x, φ(x, y)) -.
Hence, we easily deduce that the operator norm of A ∞ (x, y) † satisfies:

A ∞ (x, y) † op ≤ σ -2 min ∂ 2 yy g(x 0 , y 0 ) † op .
The result follows by setting λ = σ 2 min ∂ 2 yy g(x 0 , y 0 ) † -1 op .

We will need to introduce the following matrix U (x, y) defined as:

U (x, y) := -(A ∞ (x, y)) † B ∞ (x, y).
The following proposition shows, under mild conditions, that U t (x, y) converges towards a limiting element U ∞ (x, y) satisfying the equation:

A ∞ U ∞ = A ∞ U .
Proposition 17. Let (x 0 , y) be in X × Y and set y 0 and g be such that Assumptions 1 to 3 hold. Assume that y 0 :=φ(x 0 , y) is a local minimum of g(x 0 , .). Then there exists η > 0 such that, for any x ∈ B(x 0 , η), U t (x, y) converges towards an element U ∞ (x, y) satisfying

A ∞ (x, y)U ∞ (x, y)=A ∞ (x, y)U (x, y).
In paticular, if y is a critical point of y → g(x, .) then U ∞ (x, y) := U (x, y). Moreover, there exists a time T > 0 and constants C > 0, µ such that for any x ∈ B(x 0 , η) and t ≥ T :

U t (x, y) -U ∞ (x, y) ≤Ce -µt ,
Proof. For simplicity, we omit the dependence on (x, y) as they remain fixed. Let P be a projection matrix that commutes with A ∞ , i.e. : P A ∞ = A ∞ P . We will choose P to be either P ∞ = A ∞ A † ∞ or P = I -P ∞ . Define V t = P (U t -U ). By differentiating in time, it is easy to see that V t satisfies:

Vt = Bt -P A t V t . (27)
where Bt := P (B ∞ -B t + (A ∞ -A t )U ). Denote by (s, t) → R t s the resolvant of the linear system [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF], i.e. the squared matrix satisfying dR t s dt = -P A t R t s for t ≥ s and R s s = I . Standard results for linear differential equations [START_REF] Clark | An introduction to dynamical systems: continuous and discrete[END_REF]Chapter 2] ensure that R t s is always invertible at any time and that V t can be expressed in terms of R t s as follows:

V t = -R t 0 P U + t 0 R t s Bs ds.
Controlling R t s op : We will show the following inequality:

log R t s op ≤ t s (-λ P + A ∞ -A u ) du, (28)
where . op refers to the operator norm and λ P is the smallest eigenvalue of P A ∞ P . To achieve this, we define L t = 1 2 R t s u 2 for t ≥ s and u a vector in Y. We then differentiate L t in time to get:

Lt = -R t s u, P A t R t s u = -R t s u, P A ∞ R t s u + R t s u, P (A ∞ -A t )R t s u , = -R t s u, P A ∞ P R t s u + R t s u, P (A ∞ -A t )R t s u , ≤ 2 -λ P + A ∞ -A t op L t ,
where we used that P A ∞ P =P 2 A ∞ =P A ∞ since P and A ∞ commute. We also used elementary properties of the trace of product of matrices to get the last inequality. By integrating the above inequality, we obtain:

1 2 R t s u 2 L t ≤ 1 2 R s s u 2 e 2 t s (-λP + A∞-Au op ) du , ≤ 1 2 u 2 e 2 t s (-λP + A∞-Au op ) du ,
where we used that R s s = I. The desired bound on R t s op follows by taking the supremum over u in the unit ball.

Controlling B ∞ -B t op and A ∞ -A t op : By Proposition 14, there exists η > 0 and T > 0 such that:

φ t (x, y) -φ(x, y) ≤ Ce -tµ , ∀t ≥ T, x ∈ B(x 0 , η).
Moreover, since φ(x, y) is continuous at x 0 by Proposition 12, we can always choose η small enough so that φ(x, y) remains bounded. Hence, there exists a compact set K containing φ t (x, y) for any t ≥ T and x ∈ B(x 0 , η). Denote by |K| its diameter. By continuity of φ t (x, y), we can also take K large enough so that φ t (x, y) ∈ K for any 0 ≤ t ≤ T and x ∈ B(x 0 , η). Since g is three-times continuity differentiable by Assumption 1, there exists a positive constant L s.t. for all x ∈ B(x 0 , η) and y, y ∈ K:

∂ 2 xy g(x, y) -∂ 2 xy g(x, y ) op , ≤ L y -y , ∂ 2 yy g(x, y) -∂ 2 yy g(x, y ) op ≤ L y -y . As a result, we can write max B ∞ -B t op , A ∞ -A t op ≤ L φ(x, y) -φ t (x, y) ≤ c t . (29)
where, we defined c t to be:

c t = LCe -tµ , t ≥ T, 2L|K|, t < T.
Controlling V t : For simplicity define C t = t 0 c u du ≤ C ∞ := 2L|K|T + LCe -T µ /µ. We will first control the error term t 0 R t s Bs ds. For t > T , the following holds: where we used elementary linear algebra inequalities for the first line and ( 28) and ( 29) for the second line. We need to control U (x, y) op = A ∞ (x, y) † B ∞ (x, y) op . To achieve this, we use Proposition 16 which ensures that A ∞ (x, y) † op ≤ λ -1 for some positive λ provided x is close enough to x 0 . Thus, we can choose η small enough so that A ∞ (x, y) † op ≤ λ -1 for any x ∈ B(x 0 , η). Moreover, by Proposition 14, we know that x → φ(x, y) is continuous on B(x 0 , η) provided η is small enough. Therefore, we can ensure that B ∞ (x, y) is bounded by some value B max on B(x 0 , η). Hence, we deduce that U (x, y) op ≤ M = λ -1 B max for any x ∈ B(x 0 , η). We can finally write the upper-bound bellow:

t 0 R t s Bs op ds ≤ t 0 R t s op ( B ∞ -B t + A ∞ -A t U ∞ ) ds, (30)
t 0 R t s Bs op ds ≤e C∞ (1 + M ) t 0 c s e -(t-s)λ P ds Et . (31)
We distinguish two cases depending on the choice of P :

• Case P = A ∞ A † ∞ .
In the case where A ∞ (x 0 , y)=0, then by Proposition 9 and for η > 0 small enough, it holds that A ∞ (x, y) = 0 for any x ∈ B(x 0 , η). In this case, the dynamics is trivial. Instead, if A ∞ (x 0 , y) =0, then by Proposition 9 and for η > 0 small enough, A ∞ (x, y) = 0 for any x ∈ B(x 0 , η). In this case, we know that A ∞ (x, y) op is the inverse of the smallest positive eigenvalue of A ∞ (x, y) which is also equal to λ P by definition. Moreover, by Proposition 16, there exists η > 0 small enough and λ > 0 such that λ A ∞ (x, y) op ≤ 1 for any x ∈ B(x 0 , η). We then deduce that λ < λ P . Hence, for t ≥ T , we have:

E t =c T T 0 e -λ(t-s) + LC t T e -λ(t-s)-(s-T )µ , = c T λ e -λ(t-T ) + LC λ -µ e -µ(t-T ) -e -λ(t-T ) .
By abuse of notation, we still write 1 λ-µ e -µ(t-T ) -e -λ(t-T ) even when when λ = µ, to refer to the limit (t -T )e -λ(t-T ) when µ approaches λ. By introducing μ = 1 2 min(λ, µ), we get the simpler bound:

E t ≤ (c T + LC) μ e -μ(t-T ) .
On the other hand, recalling the upper-bound on R t s op we deduce that R t 0 op ≤ e C∞-λ P t . Hence, we can write for any t ≥ T :

V t ≤e C∞ (1 + M ) e -λ P t + c T + LC μ e -μ(t-T ) , ≤e C∞ (1 + M ) 1 + c T + LC μ e μT e -μt .
Hence, V t converges towards 0 at an exponential rate.

•

Case P = I -A ∞ A † ∞ .
In this case, λ P = 0 and P U = -P A † ∞ B ∞ = 0. Therefore, V t simplifies to V t = t 0 R t s Bs ds. We will simply show that such integral is absolutely convergent. To achieve this, we consider t ≥ T and compute E t :

E t = t 0 c s ds = T 0 c s ds + t T c s ds, =T c T + LC t T e -µ(s-T ) ds, =T c T + LC µ 1 -e -µ(t-T ) ≤ T c T + LC µ := E ∞ .
Hence, E t converges to a finite quantity E ∞ . Using [START_REF] Li | Convergence and stability of generalized gradient systems by łojasiewicz inequality with application in continuum kuramoto model[END_REF], we deduce that t 0 R t s Bs ds is absolutely convergent so that V t converges to an element V ∞ . Moreover, we have:

V t -V ∞ ≤ ∞ t R t s Bs ds ≤e C∞ (1 + M )(E ∞ -E t ), ≤Ce C∞ L µ (1 + M )e -µ(t-T ) .
Hence, we have shown that there exists η > 0 small enough such that for any x ∈ B(x 0 , η), U t (x, y) converges to an element U ∞ (x, y) satisfying A ∞ (x, y)U ∞ (x, y) = A ∞ (x, y)U (x, y). Moreover, the there exists a time T and positive constants C and µ such that:

U t (x, y) -U ∞ (x, y) ≤ C e -µ t , ∀t ≥ T, ∀x ∈ B(x 0 , η).
Proof of Theorem 2. By Proposition 10, we have that φ t (x, y) converges to φ(x, y). Moreover, since φ(x 0 , y) is a local minimizer of g(x 0 , .), Proposition 12 ensures that φ(x, y) is continuous at x 0 . Finally, we know by Proposition 15 that φ t (x, y) is differentiable in x and by Proposition 17 that ∂ x φ t (x, y) := U t (x, y) converges uniformly towards U ∞ (x, y). Therefore, by [START_REF] Rudin | Principles of mathematical analysis[END_REF]Theorem 7.17], we conclude that φ(x, y) is differentiable in a neighborhood of x 0 with differential given by ∂ x φ(x, y) = U ∞ (x, y). If in addition, y is a local minimizer, then, by Proposition 17, ∂ x φ(x, y)= -∂ xy g(x, y)(∂ yy g(x, y)) † .

D Limits Points of Bilevel Optimization Algorithms

Proposition 18. Let g be a real-valued function on X × Y such that Assumption 2 holds. Consider the maps ϕ T and I M defined in [START_REF] Liu | Investigating bi-level optimization for learning and vision from a unified perspective: A survey and beyond[END_REF] and let T and M be non-negative integers, such that T +M > 0.

Let (x, y) ∈ X × Y such that ϕ T (x, I M (x, y)) = y. Then, ∂ y g(x, y) = 0.

Proof. Let us fix (x, y) ∈ X × Y and consider the iterates y T = ϕ T (x, y). We will show that y T satisfy a sufficient decrease condition for some positive constant a:

g(x, y T +1 ) + L 2 y T -y T +1 2 ≤ g(x, y T ). (32)
To see this, we can use the smoothness of g to write:

g(x, y T +1 ) -g(x, y T ) ≤ -d H T d + L 2 H T d 2 ,
where d = ∂ y g(x, y T ) and we write H T =H T (x, y) by abuse of notation. Hence, it follows that:

g(x, y T +1 ) -g(x, y T ) + L 2 y T +1 -y T 2 ≤ -d H T -LH 2 T d ≤ 0,
where we used that H T ≤ 1 L I. Similarly, we obtain a sufficient decrease condition for the iterates defined by I M . Consider now T and M , such that T + M > 0, and let (x, y) be such that ϕ T (x, I M (x, y))=y. Consider the iterates y k =I k (x, y) for m ≤ M , and y k =ϕ t (x, y M ) for t ≤ T . Then the iterates y k define a non-increasing sequence g(x, y k ). Moreover, since y T +M = y 0 = y, it must be that g(x, y k ) = g(x, y). The sufficient decrease condition in [START_REF] Liao | Reviving and improving recurrent back-propagation[END_REF] implies that the iterates are all constant y k = y 0 . In particular, if M > 0, this implies that H M (x, y)∂ y g(x, y) = 0 so that ∂ y g(x, y) = 0 since H M (x, y) is invertible. On the other hand, if M = 0, then the condition T +M > 0 implies that T > 0, so that y = y 1 = y -H T (x, y)∂ y g(x, y). Similarly, since H T (x, y), we deduce that ∂ y g(x, y) = 0.

Proposition 19 (Properties of the maps ϕ T and I M ). Let g be a function satisfying Assumption 2 with a smoothness constant L. Consider ϕ T (x, y) and I M (x, y) defined by the following recursion which holds for any x, y ∈ X × Y:

ϕ T +1 (x, y) = ϕ T (x, y) -H T (x, y)(∂ y g(x, ϕ T (x, y))), ϕ 0 (x, y) = y (33) I M +1 (x, y) = I M (x, y) -H M (x, y)(∂ y g(x, I M (x, y))), I 0 (x, y) = y,
where H T (x, y) and H M (x, y) are positive symmetric matrices satisfying H M (x, y) ≤ 1 L I and H T (x, y) ≤ 1 L I for any (x, y) ∈ X × Y and non-negative integers T, M . Moreover, assume that H T (x, y) is continuously differentiable. Then ϕ T and I M satisfy Assumption 4.

Proof. It is clear that for any (x, y) ∈ X × Y, s.t. y is a critical point of g(x, .), we have that I M (x, y) = ϕ T (x, y) = y. Moreover, if T, M are such that T + M > 0 and (x, y) ∈ X × Y satisfy ϕ T (x, I M (x, y)) = y, then Proposition 18 ensures that ∂ y g(x, y)=0. It remains to obtain an expression for ∂ x ϕ T and ∂ y ϕ T in terms of second-order derivatives of g. We proceed by recursion. For T = 0, by setting D = 0, we have that:

∂ x ϕ 0 (x, y) = 0 = ∂ 2 xy g(x, y)D, ∂ y ϕ 0 (x, y) = I = I + ∂ 2
xy g(x, y)D. Let (x, y) be an augmented critical point of g. Assume now that for some T ≥ 0, there exists a matrix D T , such that:

∂ x ϕ T (x, y) = ∂ 2 xy g(x, y)D T , ∂ y ϕ T (x, y) = I + ∂ 2 yy g(x, y)D T . (34)
Differentiating the expression of ϕ T +1 (x, y) w.r.t. x and y yields:

∂ x ϕ T +1 (x, y) =∂ x ϕ T (x, y) -∂ xy g(x, ϕ T (x, y)) + ∂ x ϕ T (x, y)∂ 2 yy g(x, ϕ T (x, y)) H T (x, y) -∂ x H T (x, y)∂ y g(x, ϕ T (x, y)). =∂ xy g(x, y) D T -I + D T ∂ 2 yy g(x, y) H T (x, y) = ∂ xy g(x, y)D T +1 ,
Where we defined D T +1 (x, y)=D T -I + D T ∂ 2 yy g(x, y) H T (x, y). In the above expression, the last line follows by recalling that ϕ T (x, y)=0 and ∂ x H T (x, y)∂ y g(x, ϕ T (x, y)) = 0 since (x, y) is an augmented critical point of g and by using the recursion assumption on ∂ x ϕ T (x, y).

Similarly, for ∂ y ϕ T +1 (x, y), the following holds:

∂ y ϕ T +1 (x, y) =∂ y ϕ T (x, y) -∂ y ϕ T (x, y)∂ 2 yy g(x, ϕ T (x, y))H T (x, y) -∂ y H T (x, y)∂ y g(x, ϕ T (x, y)), =I + ∂ 2 yy g(x, y) D T -I + D T ∂ 2 yy g(x, y) H T (x, y) = I + ∂ 2 yy g(x, y)D T +1 .
Hence, by recursion, ϕ T (x, y) satisfies the equation ( 34) for any T ≥ 0. We have shown that ϕ T and I M satisfy Assumption 4.

Proof of Proposition 3. Fix T ≥ and consider the iterates (x k , y k ) of Algorithm 1 using ϕ T . By assumption (x k , y k ) k≥0 converges to an element (x T , y T ) in X × Y. By continuity of the maps ϕ T , I M and ∂ x L T , we have that:

y T = lim k y k = lim k ϕ T (x k-1 , I M (x k-1 , y k-1 )) = ϕ T (x T , I M (x T , y T )), lim k ∂ x L T (x k-1 , I M (x k-1 , y k-1 )) = ∂ x L T (x T , I M (x T , y T )) := d .
By Assumption 4, the first equation implies that y T is a critical point of g(x T , .) (i.e. ∂ y g(x T , y T ) = 0). Moreover, taking the limit in the update equation x k = x k-1 -γd k yields d = 0. Hence, we also have that ∂ x L T (x T , I M (x T , y T )) = 0. Finally, recall that I M (x T , y T )=y T by Assumption 4 since (x T , y T ) is an augmented critical point of g. Thus we have shown that:

∂ x L T (x T , y T ) = 0, ∂ y g(x T , y T ).
Assume now that y T is a local minimum of g(x T , .) and that (x T , y T ) T ≥0 is bounded. Hence, there exists a subsequence of (x T , y T ) T ≥0 converging towards an accumulation point (x , y ). By abuse of notation, we denote (x T , y T ) T ≥0 such subsequence. By continuity of the Hessian of g, it follows that y must also be a local minimum of g(x , .). We can now use Assumption 5 which ensures that ϕ T converges to a selection φ. Moreover, since ∂ x ϕ T converges uniformly near local minima, it follows by [START_REF] Rudin | Principles of mathematical analysis[END_REF]Theorem 7.17] that φ(x, y) is differentiable w.r.t. x near (x , y ) and that ∂ x ϕ T (x, y) converges uniformly near (x , y ) towards ∂ x φ(x, y). Hence, we can write for T large enough:

∂ x L φ (x T , y T ) =∂ x L T (x T , y T ) + (∂ x φ(x T , y T ) -∂ x ϕ T (x T , y T ))∂ y f (x T , y T ), =(∂ x φ(x T , y T ) -∂ x ϕ T (x T , y T ))∂ y f (x T , y T ).
By uniform convergence of ∂ x ϕ T (x, y) to ∂ x φ(x, y) and recalling that (x T , y T ) is bounded, we deduce that ∂ x L φ (x T , y T ) converges to 0. In particular, this holds true for a subsequence satisfying lim sup T ∂ x L φ (x T , y T ) = lim T ∂ x L φ (x T , y T ) , which proves the desired result.

Proof of Proposition 4. Let (x, y) ∈ X × Y be such that y is a local minimum of g(x, .). Define d to be:

d = ∂ x L T (x, y) -∂ 2 xy g(x, y) ∂ 2 yy g(x, y) † ∂ y L T (x, y).
By Theorem 2, x → φ(x, y) is differentiable at x and since y is a critical point of g(x, .), the differential of φ(x, y) is given by ∂ x φ(x, y) = -∂ 2 xy g(x, y) ∂ 2 yy g(x, y) † . Hence, d is equal to:

d = ∂ x L T (x, y) + ∂ x φ(x, y)∂ y L T (x, y).
Using the definition of L T and recalling that ϕ T satisfies Assumption 4, the following holds: Since f and ϕ T are continuously differentiable we get that u

d = ∂ x f (x, ϕ T (x, y)) + ∂ x ϕ T (x, y)∂ y f (x, ϕ T (x, y)) + ∂ x φ(x, y)∂ y ϕ T (x, y)∂ y f (x, ϕ T (x, y)), = ∂ x f (x, y) + (∂ x ϕ T (x, y) + ∂ x φ(x, y)∂ y ϕ T (x, y))∂ y f (x, y), = ∂ x f (x, y) + ∂ 2 xy g(x, y)D + ∂ x φ(x, y) I + ∂ 2 yy g(x, y)D ∂ y f (x, y), = ∂ x f (x, y) + ∂ x φ(x, y)∂ y f (x, y) + ∂ 2 xy g(x, y) + ∂ x φ(x, y)∂ 2 yy g(x, y) D∂ y f (x, y), = ∂ x L φ (x, y) + ∂ 2 xy g(x, y) + ∂ x φ(x,
k , v k -→ k ∂ x L T (x , y ), ∂ y L T (x , y ).
Moreover, recalling that L T (x, y)=f (x, ϕ T (x, y)), by application of the chain rule and using that ϕ T (x , y )=y it follows that:

u k -→ k u := ∂ x f (x , y ) + ∂ x ϕ T (x , y )∂ y f (x , y ), v k -→ k v := ∂ y ϕ T (x , y )∂ y f (x , y ).
By continuity of the higher-order derivatives of g, it holds that:

∂ 2 yy g(x k , y k+1 ) -→ k A * := ∂ 2 yy g(x , y ), ∂ 2 xy g(x k , y k+1 ) -→ k B * := ∂ 2 xy g(x , y ).
Recall that z k is given by the update equation z k = P(A k , v k , z k-1 ), where P is a continuous map for which z = P(A, v, z) if and only if z ∈ arg min z A 2 z + v 2 . By continuity of P, it follows that z satisfies:

z = P(A , v , z ).
Therefore, z minimizes z → (A ) 2 z + v 2 and satisfies the fixed point equation

(A ) 3 z + A v =0 so that A z = -(A ) † v . Moreover, recall that ξ k =A k z k , hence ξ k converges towards ξ := A z . Therefore, ξ = -(A ) † v .
Taking the limit as k goes to +∞, we get that d k defined in Algorithm 1 converges towards d defined by:

d : = u + B ξ , = u -B (A ) † v .
By Proposition 4, it is easy to see that d = ∂ x L φ (x , y ). Finally, recalling the update equation

x k+1 =x k -γd k and that x k -→ k
x , we directly deduce that d k -→ k 0, so that d = 0. This shows that (x , y ) is an equilibrium point of (BGS) and satisfies (SC).

D.1 Warm-start Strategy

In this section, we provide simple examples for the map P(A, v, z) to find approximate solutions minimizing Q(z) := 1 2 A 2 z + v 2 , where A is a symmetric matrix in R d×d satisfying A ≤ LI, with L being the smoothness constant of g in Assumption 2. The algorithm P can be as simple as N -step of conjugate gradient descent on Q with a step-size α ≤ 1 L 4 where L is the smoothness constant of g in Assumption 2. More formally, P(A, v, z)=z N where z N is the N iterate of the following recursion:

z n+1 = z n -α∂ z Q(z n ), z 0 = z. (35)
Algorithm T M Correction ITD [START_REF] Gunes Baydin | Automatic differentiation in machine learning: a survey[END_REF] T > 0 M = 0 False Corrected ITD T > 0 M = 0 True Truncated ITD [START_REF] Shaban | Truncated backpropagation for bilevel optimization[END_REF] T > 0 M > 0 False Corrected Truncated ITD T > 0 M > 0 True AID [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF] T = 0 T > 0 True Table 1: Recovering bilevel optimization algorithms from Algorithm 1.

It is clear that P(A, v, z) is continuous in its arguments. Moreover, using a similar argument as in Proposition 18, one can prove that whenever z is a fixed point of P(A, v, z), then z must be a critical point of Q and therefore satisfies the equation A 3 z + Av = 0. The update equation in [START_REF] Liu | Towards gradient-based bilevel optimization with non-convex followers and beyond[END_REF] depends however on the step-size α which needs to be smaller than 1 L 4 . to avoid the dependence on such step-size, A more efficient choice for the map P which does not require using a step-size, is to perform N conjugate gradient iterations on Q starting from an initial condition z.

D.2 Recovering Existing Algorithms

Table 1 below summarizes how to recover well-known gradient-based algorithms for bilevel optimization from Algorithm 1. Hence, Algorithm 1 recovers the most popular bilevel optimization algorithms but also introduces a corrected version to them to ensure that they recover the equilibria of (BGS).

E Experiments

To illustrate the effect of the corrective term introduced in Section 5 , we consider two sets of experiments: a synthetic problem for which the optimal solutions can be computed in closed form and a dataset distillation task on Cifar10 [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF] using a ResNet18 architecture [START_REF] He | Deep Residual Learning for Image Recognition[END_REF].

E.1 Synthetic Problem

Motivated by the instrumental variable regression problem [START_REF] Singh | Kernel Instrumental Variable Regression[END_REF] which solves a bilevel problem with quadratic objectives for both levels, we consider lower and upper-level objectives of the form: f (x, y) := 1 2

x A f x + C f y g(x, y) := 1 2 y A g y + y B g x where A f and A g are symmetric positive matrices of size d x ×d x and d y ×d y , B g is a d y ×d x matrix and C f is a d y vector with d x =2000 and d y =1000. To allow for multiple solutions to the LL objective, we choose A g to be non-invertible with a null-space of dimension 100 while we choose A f to be invertible for simplicity. Furthermore, to ensure that f admits a finite minimum value we choose B g to be of the form A g U for some randomly sampled matrix U . We construct the matrices A f and A g so that the highest eigenvalues of A f and A g are smaller than 1 and their conditioning is equal to 10. Here, we define the conditioning of a matrix to be the ratio between the highest and smallest non-zero eigenvalues. For a given x, the minimizers of g are of the form:

y = -A † g B g x + (I -A g A † g )y 0 ,
where y 0 is any vector in R dy . Replacing the optimal y in the UL objective results in the expression which holds for any y 0 ∈ R dy .

1 2 x A f x -C f A † g B g x + C f (I -A g A † g )y 0 .
At this point, it is easy to check that either maximizing or minimizing the above objective over y 0 results in an infinite value of the objective whenever C f (I -A g A † g ) is non-zero. This implies that the optimistic and pessimistic formulations of the bilevel problem result in an infinite optimal loss. However, the (BGS) has a well-defined solution. To see this, it is possible to define a selection of the form φ(x, y)= -A † g B g x + (I -A g A † g )y which corresponds to the limit of a gradient flow of g initialized at y. The upper objective of (BGS) is therefore given by:

L φ (x, y) = 1 2 x A f x -C f A † g B g x + C f (I -A g A † g )y.
Instead of optimizing L φ (x, y) over x and y which would result in an infinite loss, (BGS) optimizes L φ (x, y) over x only, while y is optimized for f (x, y), thus seeking an equilibrium (x , y ) satisfying (SC) which can be expressed in closed form as

x := A -1 f B g A † g C f , y := -A † g B g x + (I -A g A † g )y 0
where y 0 is any vector in R dy . Hence, while there exist multiple equilibria, they all have the same value for x and yield a finite objective.

We solve the above problem using Algorithm 1 either using the correction or not. When using the correction, we compute the approximate solution ξ k to the linear system (5) using the following update rule:

ξ k = ξ k-1 -β(∂ yy g(x k-1 , y k )ξ k-1 + v k ) (36)
where β = 0.9 is a positive step-size. For the lower-level problem, we use T steps of gradient descent with a step-size α = 0.9 while we set the upper-level step-size to γ = 1.. We then set the warm-start parameter value M to 0 and vary T .

Results. We consider the distance of the iterate x k to the optimal equilibrium x as measured by the metric induced by A f : A f as a function of time (in seconds) for different algorithmic choices, while Figure 3(right) shows the evolution of the approximate upper-level gradient d k used in Algorithm 1. We first observe that, without correction, and when using a small number of unrolled iterations (T ≤ 10), the algorithm does not converge towards x , (the distance to the iterate is larger than 10 3 ). Instead, the algorithm reaches a different equilibrium as suggested by the evolution of the gradient approximation d k towards 0 (Figure 3-(right)). As the number of unrolling steps T increases, the algorithm takes more time to converge as suggested by Figure 3-(right) (green trace T = 1000). However, the limit gets closer to the equilibrium x (Figure 3-(right), green trace). This confirms our first convergence result in Proposition 3 stating that unrolled optimization finds an approximate solution to (BGS).

x k -x 2 A f := 1 2 (x k -x ) A f (x k -x )
When using the correction, Algorithm 1 is able to recover the equilibrium x while still using a small number of unrolling steps T ≤ 10 and requiring less time to converge. This observation supports the result in Proposition 5.

E.2 Dataset Distillation on Cifar10

We consider the task of learning a small synthetic dataset so that a classifier trained on such a dataset achieves a small error on a training set. More formally, we consider a classification problem with C classes using a model with parameters y and a training dataset D tr = {(ξ i , c i )} consisting of N i.i.d. samples ξ i and corresponding labels c i . The goal is to learn a synthetic dataset of F C points, where F is a positive integer, such that each class c contains F representative samples. We can collect the synthetic points into a vector x to be learned and denote by D x the synthetic dataset. For a given dataset D, denote by L D (y) the cross-entropy loss of a model with parameters y evaluated on D. The bi-level formulation of the distillation task consists in optimizing a lower-level objective g(x, y) = L Dx (y) to learn the model parameters y that best predicts the classes of the synthetic dataset. The upper-level objective g(x, y) = L Dtr (y) evaluates the optimal model on the training set and optimizes the synthetic samples.

Setup . We consider a setup similar to [START_REF] Wang | Improving mmd-gan training with repulsive loss function[END_REF] for distilling Cifar10 [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF] on 100 synthetic points. We set F =10, thus requiring 10 synthetic points for each of the C=10 classes of Cifar10. We then use ResNet18 [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] as a classifier and apply Algorithm 1 to learn the optimal synthetic points. For the lower level, we use gradient descent with 1 unrolled iteration (i.e. T = 1, M = 0) and a step-size of α=0.001. For the upper level, we use Adam optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF], with the default parameters, a step-size of γ = 0.01 and a batch-size of 1024. When using the corrective term, we use the update equation ( 36) with a step-size β = 0.0001.

Results. Figure 4 shows the evolution of the training and test accuracy of the model as a function of time in two settings, either with or without correction. While the training accuracy for both versions of the algorithm is similar, the corrective term yields an improved final test accuracy (54.19% vs 48.6%). Note that these accuracies are of the same order as those obtained in [START_REF] Wang | Dataset distillation[END_REF] suggesting that distilling Cifar10 in only 100 samples is not sufficient to capture all variability in the dataset. While the additional correction increases the computational cost per iteration, it provides a better gradient estimate which results in a faster/better performance overall.

x∈XF

  (x), and (LL): F (x) := max y∈Y f (x, y) s.t. y ∈ T (x). (pessimistic-BG)

Figure 1 :

 1 Figure 1: Left: Heatmap of the lower-level objective g(x, y). The local minimizers of y → g(x, y) are represented by the 'critical lines' in blue. The selection map φ(x, y) is defined by following the vector field ∂yg(x, y), in black. Right: Iterates (x k , y k ) (in red) obtained by playing a BGS. The follower finds the next update y k by optimizing y → g(x k , y) starting from previous iterate y k-1 . The leader finds the next update x k by optimizing the upper-level objective f along the 'critical lines' (iterates in green).

Figure 2 :

 2 Figure 2: Two examples of functions g with different behaviors of the gradient flow under perturbations of x. In both figures, the green surface represents a function y → g(x0, y) with y ∈ R 2 resembling a Mexican hat which has a manifold of (degenerate) local minimizers (in dark green). The blue surfaces represent deformed versions of the Mexican hat function when the parameter x is slightly perturbed x1≈x0. Depending on the deformation, the resulting function y →g(x1, y) can either preserve the same type of critical points as the unperturbed function, i.e. local minimizers remain local minimizers (Left), or change their type, i.e.: local minimizers can become saddle-points (Right). Left: the selection behaves smoothly as a function of the deformation. Right: the selection is discontinuous since the gradient flow is pushed away from φ(x0, y0) which is deformed into a saddle point.

  (a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] (b) Did you describe the limitations of your work? [Yes] (c) Did you discuss any potential negative societal impacts of your work? [N/A] (d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] 2. If you are including theoretical results... (a) Did you state the full set of assumptions of all theoretical results? [Yes] (b) Did you include complete proofs of all theoretical results? [Yes] 3. If you ran experiments... (a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [N/A] (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [N/A] (c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [N/A] (d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [N/A] 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets... (a) If your work uses existing assets, did you cite the creators? [N/A] (b) Did you mention the license of the assets? [N/A] (c) Did you include any new assets either in the supplemental material or as a URL? [N/A] (d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] (e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] 5. If you used crowdsourcing or conducted research with human subjects... (a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] (b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A] (c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A] A Morse-Bott Lemma with Parameters A.1 Background on Morse-Bott Functions

Lemma 1 .

 1 which proves the result. Let h be a smooth Morse-Bott function defined on Y. Let T : X × Y → Y be a smooth function, such that y → T (x, y) = τ x (y) is a diffeomorphism on Y for any x ∈ X . Then the function g(x, y) = h(τ x (y)) is a parametric Morse-Bott function.

  The function g(x 0 , .) is a Morse-Bott function by Proposition 7, therefore by the Morse-Bott lemma [16, Theorem 2.10], there exists three open neighborhoods O, O ⊥ and V of 0 ∈ K, 0 ∈ K ⊥ and y 0 ∈ Y and a diffeomorphism s : O × O ⊥ → V s.t. s(0, 0) = y 0 and for any r, w ∈ O × O ⊥ it holds that: g(x 0 , s(r, w)) = g(x 0 , y 0 ) + 1 2 w J 0 w, ∀r, w ∈ O × O ⊥ .

  P +cu du 1 + U op c s ds, ≤e C∞ 1 + U op t 0 c s e -(t-s)λ P ds,

  y)∂ 2 yy g(x, y) D∂ y f (x, y). The last term of the above equation vanishes, since by definition of ∂ x φ(x, y), it holds that ∂ 2 xy g(x, y) + ∂ x φ(x, y)∂ 2 yy g(x, y) = 0. Therefore, we have shown that d = ∂ x L φ (x, y), which concludes the proof.Proof of Proposition 5. By continuity of the maps ϕ T and I M and since (x k , y k , z k ) → (x , y , z ), it holds that y = ϕ T (x , I M (x , y )). Hence, by Assumption 4, it follows that y must be a critical point of y → g(x , y ), i.e. ∂ y g(x , y )=0. Moreover, we have that ỹk = I M (x k , y k ) -→ k I M (x , y ) = y by continuity of I M and the condition in Assumption 4.

  vs time With Corr, T = 1 With Corr, T = 10 Without Corr, T = 1 Without Corr, T = 10 Without Corr, T = 100 Without Corr, T = 1000

Figure 3 :

 3 Figure 3: (left) Evolution of the distance of the UL iterate x k to the equilibrium x vs time (in seconds). (right) evolution of the norm of approximate gradient d k vs time in seconds. In all cases, algorithms are run until convergence, i.e. d k converges to 0.

Figure 3 (

 3 Figure 3 (left) shows the evolution of x k -x 2

Figure 4 :

 4 Figure 4: Evolution in time of the training and test accuracy of a ResNet18 model on Cifar10 dataset. Each iteration corresponds to the accuracy of the model with parameter y k trained on a synthetic dataset of 100 points x k to minimize the LL objective. The synthetic points x k are learned by minimizing the training error when using the running model y k .

I ≤ ∂ y ψ(y) ∂ y ψ(y) ≤ (L ) 2 I.

Consider now τ ∈ P. Equation [START_REF] Feehan | On the morse-bott property of analytic functions on banach spaces with łojasiewicz exponent one half[END_REF] implies that τ is L -Lipschitz on B(y 0 , 2R ). Moreover, for any z in B(y τ , η), it holds that z ∈ B(y 0 , 2R ) since η ≤ R and y τ ∈ B(y 0 , R ) by definition of y τ . Therefore, we can write the following inequality:

We have shown that τ (z) ∈ B(y 0 , L η) for any z ∈ B(y τ , η), so that [START_REF] Ghadimi | Approximation methods for bilevel programming[END_REF] holds for τ (z):

Additionally, by [START_REF] Gould | On differentiating parameterized argmin and argmax problems with application to bi-level optimization[END_REF] and using that 2 ∂ y L(τ (z))

2 by ( 16), it holds for any < C that:

From now on, we fix τ , and consider z t to the ODE [START_REF] Feurer | Hyperparameter optimization[END_REF] with initial condition z 0 ∈ B(y τ , η). Define T = {t ∈ R + |. z s -y τ < C ∀s ∈ [0, t)} which is not empty by construction since z 0 -y τ < C s → z s is continuous. Hence, t 1 := sup T is positive. We will show that t 1 = +∞. We will also consider the time until which L(τ (z t )) remains positive:

We may assume that L(τ (z 0 )) > 0 so that t + > 0 by continuity of the solution z t . The case where L(τ (z 0 )) = 0 will be treated separately. Denote by t + 1 := min(t 1 , t + ) so that, for any t ∈ [0, t + 1 ) the following holds:

-dL(τ (z t ))

where the first equality follows by differentiating z t in time and using the ODE equation ( 17), while the last inequality uses the inequality (20) which holds since z t -y τ < C. Integrating between 0 and t ∈ [0, t + 1 ), we get:

Since z 0 -y τ ≤ η and using [START_REF] Ghadimi | Approximation methods for bilevel programming[END_REF], it holds that L(z 0 )

We can therefore deduce that t 0 ∇L • τ (z s ) ds ≤ 2 . This allows to write for all t ∈ [0, t + 1 )

We distinguish two cases depending on whether t + < t 1 or t 1 ≤ t + .

Case 1: t + < t 1 or. In this case we have t + 1 = t + < +∞. This case also accounts for when L • τ (z 0 )=0 which implies that t + = 0 < t 1 . If t + =0, then z t + -y τ ≤ by construction. Otherwise, we still have that z t + -y τ ≤ by [START_REF] Guo | Solving mathematical programs with equilibrium constraints[END_REF] and the continuity of z t at t + . Moreover, by definition of t + , it must also hold that L • τ (z t + ) = 0. We only need to show that ∇L • τ (z t + ) = 0. By contradiction, if ∇L•τ (z t + ) = 0, then we would have L•τ (z t + +s ) < 0 for s > 0 small enough. However, since t + < t 1 , then t + + s < t 1 for s small enough, so that z t + +s -y τ < C. The latter means that L • τ (z t + +s ) ≥ 0 since y 0 is a local minimizer of L. This contradicts L • (z t + +s ) < 0. Therefore ∇L • τ (z t + )=0 which implies that τ (z t + ) is a critical point of y → L(y) so that z t = z t + for any t ≥ t + . This directly means that z t -y τ ≤ for any t ≥ 0, hence t 1 = +∞.

Case 2: t + ≥ t 1 . In this case, t + 1 = t 1 . If by contradiction we had t 1 < +∞, then we would directly get z t1 -y τ ≤ by continuity of t at t 1 and maximality of the solution z t . However, by defintion of t 1 , we also have z t1 -y τ = C. This contradicts the condition < C and therefore means that t 1 = +∞. Hence, it holds that z t -y τ ≤ for any t ≥ 0 and that the solution z t is well-defined at all times.