Functional Bilevel Optimization for Machine Learning

Authors: ieva  petrulionyte
julien  mairal
michael  arbel
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 15:37:11

Abstract: In this paper, we introduce a new functional point of view on bilevel optimization problems for machine learning, where the inner objective is minimized over a function space. These types of problems are most often solved by using methods developed in the parametric setting, where the inner objective is strongly convex with respect to the parameters of the prediction function. The functional point of view does not rely on this assumption and notably allows using over-parameterized neural networks as the inner prediction function. We propose scalable and efficient algorithms for the functional bilevel optimization problem and illustrate the benefits of our approach on instrumental regression and reinforcement learning tasks.

Main paper: 
 introduction: 
 Bilevel optimization methods solve problems with hierarchical structures, optimizing two interdependent objectives: an inner-level objective and an outer-level one. Initially used in machine learning for model selection [Bennett et al., 2006] and sparse feature learning [Mairal et al., 2012] , these methods gained popularity as efficient alternatives to grid search for hyper-parameter tuning [Feurer and Hutter, 2019 , Lorraine et al., 2019 , Franceschi et al., 2017] . Applications of bilevel optimization include meta-learning [Bertinetto et al., 2019] , auxiliary task learning [Navon et al., 2021] , reinforcement learning [Hong et al., 2023 , Liu et al., 2021a , Nikishin et al., 2022] , inverse problems [Holler et al., 2018] and invariant risk minimization [Arjovsky et al., 2019 , Ahuja et al., 2020] .Bilevel problems are challenging to solve, even in the well-defined bilevel setting with a unique inner-level solution. This difficulty stems from approximating both the inner-level solution and its sensitivity to the outer-level variable during gradient-based optimization. Methods like Iterative Differentiation (ITD, Baydin et al., 2017) and Approximate Implicit Differentiation (AID, Ghadimi and Wang, 2018) were designed to address these challenges in the well-defined setting, resulting in scalable algorithms with strong convergence guarantees [Domke, 2012 , Gould et al., 2016 , Ablin et al., 2020 , Arbel and Mairal, 2022a , Blondel et al., 2022 , Liao et al., 2018 , Liu et al., 2022b , Shaban et al., 2019] . These guarantees usually require the inner-level objective to be strongly convex. However, when the inner-level variables are neural network parameters, the lower-level problem becomes non-convex and may have multiple solutions due to over-parameterization. While nonconvexity is considered "benign" in this setting [Allen-Zhu et al., 2019 , Liu et al., 2022a] , multiplicity of inner-level solutions makes their dependence on the outer-level variable ambiguous [Liu et al., 2021b] , posing a major challenge in bilevel optimization for modern machine learning applications.We identify a common functional structure in bilevel machine learning problems to address the ambiguity challenge that arises with flexible models like neural networks. Specifically, we consider a prediction function h optimized by the inner-level problem over a Hilbert space H. This space consists of functions defined over an input space X and taking values in a finite dimensional vector space V. The optimal prediction function is then evaluated in the outer-level to optimize an outer-level parameter ω in a finite dimensional space Ω = R d , resulting in a functional bilevel problem:min ω∈Ω F(ω) : = L out (ω, h ⋆ ω ) s.t. h ⋆ ω = arg min h∈H L in (ω, h) . (FBO)Figure 1 : Parametric vs functional approaches for solving FBO by implicit differentiation.In contrast to classical bilevel formulations involving neural networks, where the inner objective is non-convex with respect to the network parameters, the inner objective in (FBO) defines an optimization problem over a prediction function h in a functional vector space H.A crucial consequence of adopting this new viewpoint is that it renders the strong convexity of the inner objective with respect to h a mild assumption, which ensures the uniqueness of the solution h ⋆ ω for any outer parameter value ω. Strong convexity with respect to the prediction function is indeed much weaker than strong convexity with respect to model parameters and often holds in practice. For instance, a supervised prediction task with pairs of features/labels (x, y) drawn from some training data distribution is formulated as a regularized empirical minimization problem:EQUATIONwhere H is the L 2 space of square integrable functions w.r.t. the distribution of x, and ω > 0 controls the amount of regularization. The inner objective is strongly convex in h, even though the optimal prediction function h ⋆ ω can be highly nonlinear in x. The function h ⋆ ω may then be approximated, e.g., by an overparameterized deep neural network, used here as a function approximation tool.Although appealing, the (FBO) formulation necessitates the development of corresponding theory and algorithms, which is the aim of our paper. To the best of our knowledge, this is the first work to propose a functional bilevel point of view that can leverage deep networks for function approximation. The closest works are either restricted to kernel methods [Rosset, 2008 , Kunapuli et al., 2008] and thus cannot be used for deep learning models, or propose abstract algorithms that can only be implemented for finite Hilbert spaces [Suonperä and Valkonen, 2024] .We introduce in Section 2 a theoretical framework for functional implicit differentiation in an abstract Hilbert space H that allows computing the total gradient ∇F(ω) using a functional version of the implicit function theorem [Ioffe and Tihomirov, 1979] and the adjoint sensitivity method [Pontryagin, 2018] . This involves solving a well-conditioned functional linear system, equivalently formulated as a regression problem in H, to find an adjoint function a ⋆ ω used for computing the total gradient ∇F(ω). We then specialize this framework to the common scenario where H is an L 2 space and objectives are expectations of point-wise losses. This setting covers many machine learning problems (see Sections 4.1, 4.2, and Appendix A). In Section 3, we propose an efficient algorithm where the prediction and adjoint functions can be approximated using parametric models, like neural networks, learned with standard stochastic optimization tools. We further study its convergence using analysis for biased stochastic gradient descent [Demidovich et al., 2024] .The proposed method, called functional implicit differentiation (FuncID), adopts a novel "differentiate implicitly, then parameterize" approach (left Figure 1 ): functional strong convexity is first exploited to derive an unambiguous implicit gradient in function space using a well-defined adjoint function. Then, both the lower-level solution and adjoint function are approximated using neural networks. This contrasts with traditional AID/ITD approaches (right Figure 1 ), which parameterize the inner-level solution as a neural network, leading to a non-convex 'parametric' bilevel problem in the network's parameters. An ambiguous implicit gradient is then computed by approximately solving an unstable or ill-posed linear system [Arbel and Mairal, 2022b] . Consequently, FuncID addresses the ambiguity challenge by exploiting the functional perspective and results in a stable algorithm with reduced time and memory costs. In Section 4, we demonstrate the benefits of our approach in instrumental regression and reinforcement learning tasks, which admit a natural functional bilevel structure.Related work on bilevel optimization with non-convex inner objectives. In principle, considering amended versions of the bilevel problem can resolve the ambiguity arising from non-convex inner objectives. This is the case of optimistic/pessimistic versions of the problem, often considered in the literature on mathematical optimization, where the outer-level objective is optimized over both outer and inner variables, under the optimality constraint of the inner-level variable [Dempe et al., 2007 , Ye and Ye, 1997 , Ye and Zhu, 1995 , Ye et al., 1997] . While tractable methods were recently proposed to solve them [Liu et al., 2021a ,b, 2023 , Kwon et al., 2024 , Shen and Chen, 2023] , it is unclear how well the resulting solutions behave on unseen data in the context of machine learning. For instance, when using over-parameterized models for the inner-level problem, their parameters must be further optimized for the outer-level objective, possibly resulting in over-fitting [Vicol et al., 2022] . More recently, Arbel and Mairal [2022b] proposed a game formulation involving a selection map to deal with multiple inner-level solutions. Such a formulation justifies the use of ITD/AID outside the well-defined bilevel setting, by viewing those methods as approximating the Jacobian of the selection map. However, the justification only hold under rather strong geometric assumptions. Additional related work is discussed in Appendix B on bilevel optimization with strongly-convex inner objectives, the adjoint sensitivity method that is often used in the context of ordinary differential equations, and amortization techniques [Amos et al., 2023] that have been also exploited for approximately solving bilevel optimization problems [MacKay et al., 2019, Bae and Grosse, 2020] . 
a theoretical framework for functional bilevel optimization: 
 The functional bilevel problem (FBO) involves an optimal prediction function h ⋆ ω for each value of the outer-level parameter ω. Solving (FBO) by using a first-order method then requires characterizing the implicit dependence of h ⋆ ω on the outer-level parameter ω to evaluate the total gradient ∇F(ω) in R d . Indeed, assuming that h ⋆ ω and L out are Fréchet differentiable (this assumption will be discussed later), the gradient ∇F(ω) may be obtained by an application of the chain rule:EQUATIONThe Fréchet derivative ∂ ω h ⋆ ω : H → R d is a linear operator acting on functions in H and measures the sensitivity of the optimal solution on the outer variable. We will refer to this quantity as the "Jacobian" in the rest of the paper. While the expression of the gradient in Equation ( 2) might seem intractable in general, we will see in Section 3 a class of practical algorithms to estimate it.Our starting point is to characterize the dependence of h ⋆ ω on the outer variable. To this end, we rely on the following implicit differentiation theorem (proven in Appendix C) which can be seen as a functional version of the one used in AID [Domke, 2012 , Pedregosa, 2016] , albeit, under a much weaker strong convexity assumption that holds in most practical cases of interest. Theorem 2.1 (Functional implicit differentiation). Consider problem (FBO) and assume that:• For any ω ∈ Ω, there exists µ > 0 such that h → L in (ω ′ , h) is µ-strongly convex for any ω ′ near ω.• h → L in (ω, h) has finite values and is Fréchet differentiable on H for all ω ∈ Ω.• ∂ h L in is Hadamard differentiable on Ω × H (in the sense of Definition C.1 in Appendix C).Then, ω → h ⋆ ω is uniquely defined and is Fréchet differentiable with a Jacobian ∂ ω h ⋆ ω given by:EQUATIONC ω := ∂ 2 h L in (ω, h ⋆ ω ). (3)The strong convexity assumption on the inner-level objective ensures the existence and uniqueness of the solution h ⋆ ω , while differentiability assumptions on L in and ∂ h L in ensure Fréchet differentiability of the map ω → h ⋆ ω . Though the implicit function theorem for Banach spaces [Ioffe and Tihomirov, 1979] could yield similar conclusions, it demands the stronger assumption that ∂ h L in is continuously Fréchet differentiable, which is quite restrictive in our setting of interest: for instance, when H is an L 2 -space and L in is an integral functional of the form L in (ω, h) = ℓ in (w, h(x)) dx, with ℓ in defined on Ω × V and satisfying mild smoothness and growth assumptions on ℓ in , then h → ∂ h L in (ω, h) cannot be Fréchet differentiable with uniformly continuous differential on bounded sets, unless v → ℓ in (ω, v) is a polynomial of degree at most 2 (see [Nemirovski and Semenov, 1973, Corollary 2, p 276] and discussions in [Noll, 1993 , Goodman, 1971] ). Instead, Theorem 2.1 employs the weaker notion of Hadamard differentiability for ∂ h L in , widely used in statistics, particularly for deriving the delta-method [van der Vaart and Wellner, 1996, Chapter 3.9 ]. Consequently, Theorem 2.1 allows us to cover a broader class of functional bilevel problems, as we see in Section 2.2.Similarly to AID, only a Jacobian-vector product is needed when computing the total gradient ∇F(ω). The result in Proposition 2.2 below, relies on the adjoint sensitivity method [Pontryagin, 2018] to provide a more convenient expression for ∇F(ω) and is proven in Appendix C.2. Proposition 2.2 (Functional adjoint sensitivity). Under the same assumption on L in as in Theorem 2.1 and further assuming that L out is jointly differentiable in ω and h, the total objective F is differentiable with ∇F(ω) given by:EQUATIONwhere the adjoint function a ⋆ ω := -C -1 ω d ω is an element of H that minimizes the quadratic objective:EQUATIONEquation ( 4) indicates that computing the total gradient requires optimizing the quadratic objective (5) to find the adjoint function a ⋆ ω . The strong convexity of the adjoint objective L adj ensures the existence of a unique minimizer, and stems from the positive definiteness of its Hessian operator C ω due to the inner-objective's strong convexity. Both adjoint and inner-level problems occur in the same function space H and are equivalent in terms of conditioning, as the adjoint Hessian operator equals the inner-level Hessian at optimum. The strong convexity in H of the adjoint objective guarantees well-defined solutions and holds in many practical cases, as opposed to classical parametric bilevel formulations which require the more restrictive strong convexity condition in the model's parameters, and without which instabilities may arise due to ill-conditioned linear systems (see Appendix F.1).Specializing the abstract results from Section 2.1 to a common scenario in machine learning, we consider both inner and outer level objectives of (FBO) as expectations of point-wise functions over observed data. Specifically, we have two data distributions P and Q defined over a product space X × Y ⊂ R dx × R dy , and denote by H the Hilbert space of functions h : X → V, where V = R dv . Given an outer parameter space Ω, we address the following functional bilevel problem:EQUATIONwhere ℓ out , ℓ in are point-wise loss functions defined on Ω × V × X × Y. This setting encompasses various deep learning problems discussed in Sections 4.1 and 4.2, and in Appendix A, representing a specific instance of FBO. The Hilbert space H of square-integrable functions not only models a broad range of prediction functions but also facilitates obtaining concrete expressions for the total gradient ∇F(ω), enabling the derivation of practical algorithms in Section 3.The following proposition, proved in Appendix D, makes mild technical assumptions on probability distributions P, Q and the point-wise losses ℓ in , ℓ out . It gives an expression for the total gradient in the form of expectations under P and Q. Proposition 2.3 (Functional Adjoint sensitivity in L 2 spaces.). Under the technical Assumptions (A) to (L) stated in Appendix D.1, the conditions on L in and L out in Proposition 2.2 hold and the total gradient ∇F(ω) of F is expressed as ∇F(ω) = g ω + B ω a ⋆ ω , with a ⋆ ω ∈ H being the minimizer of the objective L adj in Equation (5). Moreover, L adj , g ω and B ω a ⋆ ω are given by:EQUATIONEQUATIONwhere ∂ ω ℓ out , ∂ v ℓ out , ∂ ω,v ℓ in , and ∂ 2 v ℓ in are partial first and second order derivatives of ℓ out and ℓ in with respect to their first and second arguments ω and v.Assumptions (A) and (B) on P and Q ensure finite second moments and bounded Radon-Nikodym derivatives, maintaining square integrability under both distributions in Equation (6). Assumptions (C) to (L) on ℓ in and ℓ out primarily involve integrability, differentiability, Lipschitz continuity, and strong convexity of ℓ in in its second argument, typically satisfied by objectives like mean squared error or cross entropy (see Proposition D.1 in Appendix D.1). Next, by leveraging Proposition 2.3, we derive practical algorithms for solving FBO using function approximation tools like neural networks.3 Methods for Functional Bilevel Optimization in L 2 SpacesWe propose Functional Implicit Differentiation (FuncID), a flexible class of algorithms for solving the functional bilevel problem in L 2 spaces described in Section 2.2 when samples from distributions P and Q are available. FuncID relies on three main components detailed in the next subsections:1. Empirical objectives. These approximate the objectives L out , L in and L adj as empirical expectations over samples from inner and outer datasets D in and D out . 2. Function approximation. The search space for both the prediction and adjoint functions is restricted to parametric spaces with finite-dimensional parameters θ and ξ. Approximate solutions ĥω and âω to the optimal functions h ⋆ ω and a ⋆ ω are obtained by minimizing the empirical objectives. 3. Total gradient approximation. FuncID estimates the total gradient ∇F(ω) using the empirical objectives, and the approximations ĥω and âω of the prediction and adjoint functions.Algorithm 1 provides an outlines of FuncID which has a nested structure similar to AID: (1) innerlevel optimizations (InnerOpt and AdjointOpt) to update the prediction and adjoint models using scalable algorithms such as stochastic gradient descent [Robbins and Monro, 1951] , and (2) an outer-level optimization to update the parameter ω using a total gradient approximation TotalGrad.An optional warm-start allows initializing the parameters of both the prediction and adjoint models for the current outer-level iteration with those obtained from the previous one.Algorithm 1 FuncID Input: initial outer, inner, and adjoint parameter ω 0 , θ 0 , ξ 0 ; warm-start option WS.for n = 0, . . . , N -1 do # Optional warm-start if WS=True then (θ 0 , ξ 0 ) ← (θ n , ξ n ) end if # Inner-level optimization ĥωn , θ n+1 ← InnerOpt(ω n , θ 0 , D in ) # Adjoint optimization âωn , ξ n+1 ← AdjointOpt(ω n , ξ 0 , ĥωn , D) # Outer gradient estimation Sample a mini-batch B = (B out , B in ) from D = (D out , D in ) g out ← TotalGrad(ω n , ĥωn , âωn , B) ω n+1 ← update ω n using g out ; end for 
from population losses to empirical objectives: 
 We assume access to two datasets D in and D out , comprising i.i.d. samples from P and Q, respectively. This assumption can be relaxed, such as when using samples from a Markov chain or a Markov Decision Process to approximate population objectives. For scalability, we operate in a minibatch setting, where batches B = (B out , B in ) are sub-sampled from datasets D := (D out , D in ). Approximating both inner and outer level objectives in (6) can be achieved using empirical versions:Lout (ω, h, B out ) := 1 |Bout| (x,ỹ)∈Bout ℓ out (ω, h(x), x, ỹ) , Lin (ω, h, B in ) := 1 |Bin| (x,y)∈Bin ℓ in (ω, h(x), x, y) .Adjoint objective. Using the expression of L adj from Proposition 2.3, we derive a finite-sample approximation of the adjoint loss by replacing the population expectations by their empirical counterparts. More precisely, assuming we have access to an approximation ĥω to the inner-level prediction function, we consider the following empirical version of the adjoint objective:Ladj ω, a, ĥω , B := 1 2 1 |Bin| (x,y)∈Bin a(x) ⊤ ∂ 2 v ℓ in (ω, ĥω (x), x, y) a(x) + 1 |Bout| (x,y)∈Bout a(x) ⊤ ∂ v ℓ out ω, ĥω (x), x, y . (9)The adjoint objective in Equation ( 9) requires computing a Hessian-vector product with respect to the output v in R dv of the prediction function ĥω , which is typically of reasonably small dimension, unlike traditional AID methods that necessitate a Hessian-vector product with respect to some model parameters. Importantly, compared to AID, FuncID does not requires differentiating twice w.r.t the model's parameters τ (θ) which results in memory and time savings as discussed in Appendix F.2.To find approximate solutions to the prediction and adjoint functions we rely on three steps: 1) specifying parametric search spaces for both functions, 2) introducing optional regularization to prevent overfitting and, 3) defining a gradient-based optimization procedure on the empirical objectives.Parametric search spaces. We approximate both prediction and adjoint functions using parametric search spaces. A parametric family of functions defined by a map τ : Θ → H over parameters Θ ⊆ R pin constrains the prediction function h as h(x) = τ (θ)(x). We only require τ to be continuous and differentiable almost everywhere such that back-propagation can be applied [Bolte et al., 2021] . Notably, unlike AID, we do not need τ to be twice differentiable, as functional implicit differentiation computes the Hessian w.r.t. the output of τ , not w.r.t. its parameters θ. For flexibility, we can consider a different parameterized model ν : Ξ → H for approximating the adjoint function, defined over parameters Ξ ⊆ R p adj , constraining the adjoint similarly to τ . In practice, we often use the same parameterization, typically a neural network, for both the inner-level and the adjoint models.Regularization. With empirical objectives and parametric search spaces, we can directly optimize parameters of both the inner-level model τ and the adjoint model ν. However, to address finite sample issues, regularization may be introduced to these empirical objectives for better generalization.The method allows flexibility in regularization choice, accommodating functions θ → R in (θ) and ξ → R adj (ξ), such as ridge penalty or other commonly used regularization techniquesOptimization. The function InnerOpt (defined in Algorithm 2) optimizes inner model parameters for a given ω, initialization θ 0 , and data D in , using M gradient updates. It returns optimized parameters θ M and the corresponding inner model ĥω = τ (θ M ), approximating the inner-level solution. Similarly, AdjointOpt (defined in Algorithm 3) optimizes adjoint model parameters with K gradient updates, producing the approximate adjoint function âω = ν(ξ K ). Other optimization procedures may also be used, especially when closed-form solutions are available, as exploited in some experiments in Section 4. Operations requiring differentiation can be implemented using standard optimization procedures with automatic differentiation packages like PyTorch [Paszke et al., 2019] or Jax [Bradbury et al., 2018] .Algorithm 2 InnerOpt(ω, θ 0 , D in ) for m = 0, . . . , M -1 do Sample batch B in from D in g in ← ∇ θ [ Lin (ω, τ (θ m ), B in )+R in (θ m )] θ m+1 ← Update θ m using g in end for Return τ (θ M ), θ M Algorithm 3 AdjointOpt(ω, ξ 0 , ĥω , D) for k = 0, . . . , K -1 do Sample batch B from D g adj ← ∇ ξ [ Ladj (ω, ν(ξ t ), ĥω , B)+R adj (ξ k )] ξ k+1 ← Update ξ k using g adj end for Return ν(ξ K ), ξ KWe exploit Proposition 2.3 to derive Algorithm 4, which allows us to approximate the total gradient ∇F(ω) after computing the approximate solutions ĥω and âω . There, we decompose the gradient into two terms: g Exp , an empirical approximation of g ω in Equation ( 8) representing the explicit dependence of L out on the outer variable ω, and g Imp , an approximation to the implicit gradient term B ω a ⋆ ω in Equation ( 8). Both terms are obtained by replacing the expectations by empirical averages batches B in and B out , and using the approximations ĥω and âω instead of the exact solutions.Algorithm 4 TotalGrad(ω, ĥω , âω , B)g Exp ← ∂ ω Lout (ω, ĥω , B out ) g Imp ← 1 |Bin| (x,y)∈Bin ∂ ω,v ℓ in (ω, ĥω (x), x, y) âω (x) Return g Exp + g ImpConvergence of Algorithm 1 to stationary points of F depends on approximation errors, which result from sub-optimal inner and adjoint solutions, as shown by the convergence result below. Theorem 3.1. Assume that F is L-smooth and admits a finite lower bound F ⋆ . Use an update rule ω n+1 = ω n -ηg out with step size 0 < η ≤ 1 4Lω in Algorithm 1. Under Assumption (a) on sub-optimality of ĥω and âω , Assumptions (b) to (h) on the smoothness of ℓ in and ℓ out , all stated in Appendix E.1, and the assumptions in Proposition 2.3, the iterates {ω n } n≥0 of Algorithm 1 satisfy:min 0≤n≤N -1 E ∥∇F(ω n )∥ 2 ≤ 4 (F(ω 0 ) -F ⋆ ) ηN + 2ηLσ 2 ef f + (c 1 ϵ in + c 2 ϵ adj ),where c 1 , c 2 , σ 2 ef f are positive constants, and ϵ in , ϵ adj are sub-optimality errors that result from the inner and adjoint optimization procedures.Theorem 3.1 is proven in Appendix E and relies on the general convergence result in Demidovich et al. [2024, Theorem 3] for stochastic biased gradient methods. The key idea is to control both bias and variance of the gradient estimator in terms of generalization errors ϵ in and ϵ adj when approximating the inner and adjoint solutions. These generalization errors can, in turn, be made smaller in the case of over-parameterized networks, by increasing network capacity, number of steps and sample size [Allen-Zhu et al., 2019 , Du et al., 2019 , Zou et al., 2020] . 
applications: 
 We consider two applications of the functional bilevel optimization problem: Two-stage least squares regression (2SLS) and Model-based reinforcement learning. To illustrate its effectiveness we compare it with other bilevel optimization approaches like AID or ITD, as well as state-of-the-art methods for each application. We provide a versatile implementation of FuncID (https://github.com/inria-thoth/funcBO) in PyTorch [Paszke et al., 2019] , compatible with standard optimizers (e.g., Adam [Kingma and Ba, 2015] ), and supports common regularization techniques. For the reinforcement learning application, we extend an existing JAX [Bradbury et al., 2018] implementation of model-based RL from Nikishin et al. [2022] to apply FuncID. To ensure fairness, experiments are conducted with comparable computational budgets for hyperparameter tuning using the MLXP experiment management tool [Arbel and Zouaoui, 2024] . Additionally, we maintain consistency by employing identical neural network architectures across all methods and repeating experiments multiple times with different random seeds.Two-stage least squares regression (2SLS) is commonly used in causal representation learning, including instrumental regression or proxy causal learning [Stock and Trebbi, 2003] . Recent studies have applied bilevel optimization approaches to address 2SLS, yielding promising results [Xu et al. , 2021b,a, Hu et al., 2023] . We particularly focus on 2SLS for Instrumental Variable (IV) regression, a widely-used statistical framework for mitigating endogeneity in econometrics [Blundell et al., 2007 [Blundell et al., , 2012]] , medical economics [Cawley and Meyerhoefer, 2012] , sociology [Bollen, 2012] , and more recently, for handling confounders in off-line reinforcement learning [Fu et al., 2022] .Problem formulation. In an IV problem, the objective is to model f ω : t → o that approximates the structural function f struct using independent samples (o, t, x) from a data distribution P, where x is an instrumental variable. The structural function f struct delineates the true effect of a treatment t on an outcome o. A significant challenge in IV is the presence of an unobserved confounder ϵ, which influences both t and o additively, rendering standard regression ineffective for recovering f ω . However, if the instrumental variable x solely impacts the outcome o through the treatment t and is independent from the confounder ϵ, it can be employed to elucidate the direct relationship between the treatment t and the outcome o using the 2SLS framework, under mild assumptions on the confounder [Singh et al., 2019] . This adaptation replaces the regression problem with a variant that averages the effect of the treatment t conditionally on xEQUATION)Directly estimating the conditional expectation E P [f ω (t)|x] is hard in general. Instead, it is easier to express it, equivalently, as the solution of another regression problem predicting f ω (t) from x:EQUATION)Both equations result in the bilevel formulation in Equation ( 6) with y = (t, o), Q = P and the point-wise losses ℓ in and ℓ out given by ℓin (ω, v, x, y) = ℓ in (ω, v, x, (t, o)) = ∥f ω (t) -v∥ 2 and ℓ out (ω, v, x, y) = ℓ out (ω, v, x, (t, o)) = ∥o -v∥ 2 .It is, therefore, possible to directly apply Algorithm 1 to learn f ω as we illustrate below.Experimental setup. We study the IV problem using the dsprites dataset [Matthey et al., 2017] , comprising synthetic images representing single objects generated from five latent parameters: shape, scale, rotation, and posX, posY positions on image coordinates. Here, the treatment variable t is the images, the hidden confounder ϵ is the posY coordinate, and the other four latent variables form the instrumental variable x. The outcome o is an unknown structural function f struct of t, contaminated by confounder ϵ as detailed in Appendix G.1. We follow the setup of the Deep Feature Instrumental Variable Regression (DFIV) dsprites experiment by Xu et al. [2021a, Section 4.2] , which achieves state-of-the-art performance. In this setup, neural networks serve as the prediction function and structural model, optimized to solve the bilevel problem in Equations ( 10) and (11). We explore two versions of our method: FuncID, which optimizes all adjoint network parameters, and FuncID linear, which learns only the last layer in closed-form while inheriting hidden layer parameters from the inner prediction function. We compare our method with DFIV, AID, ITD, and Penalty-based methods: gradient penalty (GD penal.) and value function penalty (Val penal.) [Shen and Chen, 2023] , using identical network architectures and computational budgets for hyperparameter selection. Full details on network architectures, hyperparameters, and training settings are provided in Appendix G.2.Results. Figure 2 compares structural models learned by different methods using 5K training samples (refer to Figure 6 in Appendix G.3 for 10K sample results). The left subplot illustrates out-of-sample mean squared error of learned structural models compared to ground truth outcomes (uncontaminated by noise ϵ), while the middle and right subplots show the evolution of outer and inner objectives over iterations. For the 5K dataset, FuncID outperforms DFIV (p-value=0.003, one-sided paired t-test), while showing comparable performance on the 10K dataset (Figure 6 ). AID and ITD perform notably worse, indicating their parametric approach fails to fully leverage the functional structure. FuncID outperforms the gradient penalty method and performs either better or comparably to the value function penalty method, though the latter shows higher variance with some particularly bad outliers. While all methods achieve similar outer losses, this criterion alone is only reliable as an indicator of convergence when evaluated near the 'exact' inner-level solution corresponding to the lowest inner-loss values. Interestingly, FuncID obtains the lowest inner-loss values, suggesting its outer-loss is a more reliable indicator of convergence.Model-based reinforcement learning (RL) naturally yields bilevel optimization formulations, since several components of an RL agent need to be learned using different objectives. Recently, Nikishin et al. [2022] showed that casting model-based RL as a bilevel problem can result in better tolerance to model-misspecification. Our experiments show that the functional bilevel framework yields improved results even when the model is well-specified, suggesting a broader use of the bilevel formulation.Problem formulation. In model-based RL, the Markov Decision Process (MDP) is approximated by a probabilistic model q ω with parameters ω that can predict the next state s ω (x) and reward r ω (x), given a pair x := (s, a) where s is the current environment state and a is the action of an agent. A second model h can be used to approximate the action-value function h(x) that computes the expected cumulative reward given the current state-action pair. Traditionally, the action-value function is learned using the current MDP model, while the latter is learned independently from the action-value function using Maximum Likelihood Estimation (MLE) [Sutton, 1991] .In the bilevel formulation of model-based RL by Nikishin et al. [2022] , the inner-level problem involves learning the optimal action-value function h ⋆ ω with the current MDP model q ω and minimizing the Bellman error. The inner-level objective can be expressed as an expectation of a point-wise loss f with samples (x, r ′ , s ′ ) ∼ P, derived from the agent-environment interaction:EQUATIONHere, the future state and reward (r ′ , s ′ ) are replaced by the MDP model predictions r ω (x) and s ω (x).In practice, samples from P are obtained using a replay buffer. The buffer accumulates data over several episodes of interactions with the environment, and can therefore be considered independent of the agent's policy. The point-wise loss function f represents the error between the action-value function prediction and the expected cumulative reward given the current state-action pair:f (v, r ′ , s ′ ) := 1 2 v -r ′ -γ log a ′ e h(s ′ ,a ′ ) 2 ,with h a lagged version of h (exponentially averaged network) and γ a discount factor. The MDP model is learned implicitly using the optimal function h ⋆ ω , by minimizing the Bellman error w.r.t. ω:EQUATIONEquations ( 12) and ( 13) define a bilevel problem as in Equation ( 6), where Q = P, y = (r ′ , s ′ ), and the point-wise losses ℓ in and ℓ out are given by:ℓ in (ω, v, x, y) = f (v, r ω (x), s ω (x)) and ℓ out (ω, v, x, y) = f (v, r ′ , s ′ ).Therefore, we can directly apply Algorithm 1 to learn both the MDP model q ω and the optimal action-value function h ⋆ ω .Experimental setup. We apply FuncID to the CartPole control problem, a classic benchmark in reinforcement learning [Brockman et al., 2016 , Nagendra et al., 2017] . The goal is to balance a pole attached to a cart by moving the cart horizontally. Following Nikishin et al. [2022] , we use a modelbased approach and consider two scenarios: one with a well-specified network accurately representing the MDP, and another with a misspecified model having fewer hidden layer units, limiting its capacity. Both plots show mean reward over 10 runs where the shaded region is the 95% confidence interval.Using the bilevel formulation in Equations ( 12) and ( 13), we compare FuncID with the Optimal Model Design (OMD) algorithm [Nikishin et al., 2022] , a variant of AID. Additionally, we compare against a standard single-level model-based RL formulation using Maximum Likelihood Estimation (MLE) [Sutton, 1991] . For the adjoint function in FuncID, we derive a simple closed-form expression based on the structure of the adjoint objective (see Appendix H.1). We follow the experimental setup of Nikishin et al. [2022] , providing full details and hyperparameters in Appendix H.2.Results. Figure 3 illustrates the training reward evolution for FuncID, OMD, and MLE in both well-specified and misspecified scenarios. FuncID consistently performs well across settings. In the well-specified case, where OMD achieves a reward of 4, FuncID reaches the maximum reward of 5, matching MLE (left Figure 3 ). In the misspecified scenario, FuncID performs comparably to OMD and significantly outperforms MLE (right Figure 3 ). Moreover, FuncID tends to converge faster than MLE (see Figure 7 in Appendix H.3) and yields consistently better prediction error than OMD (see Figure 8 in Appendix H.3). These findings align with Nikishin et al. [2022] , suggesting that MLE may prioritize minimizing prediction errors, potentially leading to overfitting irrelevant features.In contrast, OMD and FuncID focus on maximizing expected returns, especially in the presence of model misspecification. Our results highlight the effectiveness of (FBO) even in well-specified settings, suggesting, for future work, further investigations for more general RL tasks. 
discussion and concluding remarks: 
 This paper introduced a functional paradigm for bilevel optimization in machine learning, shifting focus from parameter space to function space. The proposed approach specifically addresses the ambiguity challenge arising from using deep networks in bilevel optimization. The paper establishes the validity of the functional framework by developing a theory of functional implicit differentiation, proving convergence for the proposed FuncID method, and numerically comparing it with other bilevel optimization methods.The theoretical foundations of our work rely on several key assumptions worth examining. While our convergence guarantees assume both inner and adjoint optimization problems are solved to some optimality, this assumption is supported by recent results on global convergence in over-parameterized networks [Allen-Zhu et al., 2019 , Liu et al., 2022a] . However, quantifying these optimality errors more precisely and understanding their relationship to optimization procedures remains an open challenge. Additionally, like other bilevel methods, our approach requires careful hyperparameter selection, which can impact practical implementation.Several promising directions emerge for future research. While we focus on L 2 spaces, exploring alternative function spaces (such as Reproducing Kernel Hilbert Spaces or Sobolev spaces) could reveal additional advantages for specific applications. Furthermore, extending our framework to non-smooth objectives or constrained optimization problems, potentially building on existing work in non-smooth implicit differentiation [Bolte et al., 2022] , would broaden its applicability.The functional bilevel setting applies to various practical bilevel problems where objectives depend solely on model predictions, not their parameterization. Below, we discuss a few examples.Auxiliary task learning. As in Equation ( 1), consider a main prediction task with features x and labels y, equipped with a loss function f (y, h(x)). The goal of auxiliary task learning is to learn how a set of auxiliary tasks represented by a vector f aux (y, h(x)) could help solve the main task. This problem is formulated by Navon et al. [2021] as a bilevel problem, which can be written as (FBO) withL out (ω, h) = E (y,x)∼D val [f (y, h(x))] ,where the loss is evaluated over a validation dataset D val , andL in (ω, h) = E (y,x)∼Dtrain [f (y, h(x)) + g ω (f aux (y, h(x)))] ,where an independent training dataset D train is used, and g ω is a function that combines the auxiliary losses into a scalar value.Task-driven metric learning. Considering now a regression problem with features x and labels y, the goal of task-driven metric learning formulated by Bansal et al. [2023] is to learn a metric parameterized by ω for the regression task such that the corresponding predictor h ⋆ ω performs well on a downstream task L task . This can be formulated as (FBO) with L out (ω, h) = L task (h) andL in (ω, h) = E (y,x) ∥y -h(x)∥ 2 Aω(x) ,where ∥ • ∥ 2 ω is the squared Mahalanobis norm with parameters ω and A ω (x) is a data-dependent metric that allows emphasizing features that are more important for the downstream task.Bilevel optimization in machine learning. Two families of bilevel methods are prevalent in machine learning due to their scalability: iterative (or 'unrolled') differentiation (ITD, Baydin et al., 2017) and Approximate Implicit Differentiation (AID, Ghadimi and Wang, 2018) . ITD approximates the optimal inner-level solution using an 'unrolled' function from a sequence of differentiable optimization steps, optimizing the outer variable via back-propagation [Shaban et al., 2019 , Bolte et al., 2024] . The gradient approximation error decreases linearly with the number of steps when the inner-level is strongly convex, though at increased computational and memory costs [Grazzi et al., 2020, Theorem 2.1] . ITD is popular for its simplicity and availability in major deep learning libraries [Bradbury et al., 2018] , but can be unstable, especially with non-convex inner objectives [Pascanu et al., 2013 , Bengio et al., 1994 , Arbel and Mairal, 2022b] . AID uses the Implicit Function Theorem (IFT) to derive the Jacobian of the inner-level solution with respect to the outer variable [Lorraine et al., 2019 , Pedregosa, 2016] , solving a finite-dimensional linear system for an adjoint vector representing optimality constraints. AID offers strong convergence guarantees for smooth, strongly convex inner objectives [Ji et al., 2021, Arbel and Mairal, 2022a] . However, without strong convexity, the linear system can become ill-posed due to a degenerate Hessian, leading to instabilities, especially with overparameterized deep neural networks. Our proposed approach avoids this issue, even when using deep networks for function approximation.Adjoint sensitivity method. The adjoint sensitivity method [Pontryagin, 2018] efficiently differentiates a controlled variable with respect to a control parameter. In bilevel optimization, AID applies a finite-dimensional version of this method [Margossian and Betancourt, 2021, Section 2] . Infinite-dimensional versions have differentiated solutions of ordinary differential equations (ODEs) with respect to defining parameters [Margossian and Betancourt, 2021, Section 3] , and have been used in machine learning to optimize parameters of a vector field describing an ODE [Chen et al., 2018] . Here, the ODE's vector field, parameterized by a neural network, is optimized to match observations. The adjoint sensitivity method offers an efficient alternative to the unstable process of back-propagation through ODE solvers, requiring only the solution of an adjoint ODE to compute gradient updates, improving performance [Jia and Benson, 2019 , Zhong et al., 2019 , Li et al., 2020] . This method has been adapted to meta-learning [Li et al., 2023] , viewing the inner optimization as an ODE evolution with gradients obtained via the adjoint ODE. Recently, Marion et al. [2024] employ the adjoint sensitivity method for optimizing diffusion models, where an adjoint SDE is solved to compute the total gradient. Unlike these works, which use the adjoint method for ODE/SDE solutions as functions of time, our work applies an infinite-dimensional version of the adjoint sensitivity method to general learning problems, where solutions are functions of input data.Amortization. Recently, several methods have used amortization to approximately solve bilevel problems [MacKay et al., 2019, Bae and Grosse, 2020] . These methods employ a parametric model called a hypernetwork [Ha et al., 2017 , Brock et al., 2018 , Zhang et al., 2019] , optimized to directly predict the inner-level solution given the outer-level parameter ω. Amortized methods treat the two levels as independent optimization problems: (1) learning the hypernetwork for a range of ω, and (2) performing first-order descent on ω using the hypernetwork as a proxy for the inner-level solution. Unlike ITD, AID, or our functional implicit differentiation method, amortized methods do not fully exploit the implicit dependence between the two levels. They are similar to amortized variational inference [Kingma and Welling, 2014, Rezende et al., 2014] , where a parametric model produces approximate samples from a posterior distribution. Amortization methods perform well when the inner solution's dependence on ω is simple but may fail otherwise [Amos et al., 2023, pages 71-72] . In contrast, our functional implicit differentiation framework adapts to complex implicit dependencies between the inner solution and ω.We recall the definition of Hadamard differentiability and provide in Proposition C.2 a general property for Hadamard differentiable maps that we will exploit to prove Theorem 2.1 in Appendix C.2. Definition C.1. Hadamard differentiability. Let A and B be two separable Banach spaces. A function L : A → B is said to be Hadamard differentiable [van der Vaart, 1991, Fang and Santos, 2018] if for any a ∈ A, there exist a continuous linear map d a L(a) : A → B so that for any sequence (u n ) n≥1 in A converging to an element u ∈ A, and any real valued and non-vanishing sequence (t n ) n≥1 converging to 0, it holds that:EQUATION)Proposition C.2. Let A and B be two separable Banach spaces. Let L : A → B be a Hadamard differentiable map with differential d a L at point a. Consider a bounded linear map defined over a euclidean space R n of finite dimension n and taking values in A, i.e J : R n → A. Then, the following holds:L(a + Ju) = L(a) + d a L(a)Ju + o(∥u∥).Proof. Consider a sequence (u k ) k≥1 in R n so that u k converges to 0 with ∥u k ∥ > 0 for all k ≥ 1 and define the first order error E k as follows:E k = 1 ∥u k ∥ ∥L(a + Ju k ) -L(a) -d a L(a)u k ∥ .The goal is to show that E k converges to 0. We can write u k as u k = t k ũk with t k = ∥u k ∥ and ∥ũ k ∥ = 1, so that:E k = 1 ∥t k ∥ (L(a + t k J ũk ) -L(a)) -d a L(a)ũ k .If E k were unbounded, then, by contradiction, there must exist a subsequence (E ϕ(k) ) k≥1 converging to +∞, with ϕ(k) increasing and ϕ(k) → +∞. Moreover, since ũk is bounded, one can further choose the subsequence E ϕ(k) so that ũϕ(k) converges to some element ũ. We can use the following upper-bound:EQUATIONwhere we used that d a L(a) is bounded. Since L is Hadamard differentiable, Ẽϕ(k) converges to 0. Moreover, ũϕ(k) -ũ also converges to 0. Hence, E ϕ(k) converges to 0 which contradicts E ϕ(k) →+∞. Therefore, E k is bounded.Consider now any convergent subsequence of (E k ) k≥1 . Then, it can be written as (E ϕ(k) ) k≥1 with ϕ(k) increasing and ϕ(k) → +∞. We then have E ϕ(k) → e < +∞ by construction. Since ũk is bounded, one can further choose the subsequence E ϕ(k) so that ũϕ(k) converges to some element ũ. Using again Equation ( 15) and the fact that L is Hadamard differentiable, we deduce that Ẽϕ(k) must converge to 0, and by definition of ũϕ(k) , that ũϕ(k) -ũ converges to 0. Therefore, it follows that E ϕ(k) → 0, so that e = 0. We then have shown that (E k ) k≥1 is a bounded sequence and every subsequence of it converges to 0. Therefore, E k must converge to 0, which concludes the proof.Proof of Theorem 2.1. The proof strategy consists in establishing the existence and uniqueness of the solution map ω → h ⋆ ω , deriving a candidate Jacobian for it, then proving that ω → h ⋆ ω is differentiable.Existence and uniqueness of a solution map ω → h ⋆ ω . Let ω in Ω be fixed. The map h → L in (ω, h) is lower semi-continuous since it is Fréchet differentiable by assumption. It is also strongly convex. Therefore, it admits a unique minimier h ⋆ ω [Bauschke and Combettes, 2011, Corollary 11.17 ]. We then conclude that the map ω → h ⋆ ω is well-defined on Ω. Strong convexity inequalities. We provide two inequalities that will be used for proving differentiability of the map ω → h ⋆ ω . The map h → L in (ω, h) is Fréchet differentiable on H and µ-strongly convex (with µ positive by assumption). Hence, for all h 1 , h 2 in H the following quadratic lower-bound holds:EQUATIONFrom the inequality above, we can also deduce that h → ∂ h L in (ω, h) is a µ-strongly monotone operator:⟨∂ h L in (ω, h 1 ) -∂ h L in (ω, h 2 ), h 1 -h 2 ⟩ H ≥ µ ∥h 1 -h 2 ∥ 2 H .(17) Finally, note that, since h → L in (ω, h) is Fréchet differentiable, its gradient must vanish at the optimum h ⋆ ω , i.e :EQUATION)Candidate Jacobian for ω → h ⋆ ω . Let ω be in Ω. Using Equation ( 17) with h 1 = h + tv and h 2 = h for some h, v ∈ H, and a non-zeros real number t we get:EQUATION)By assumption, h → ∂ h L in (ω, h) is Hadamard differentiable and, a fortiori, directionally differentiable. Thus, by taking the limit when t → 0, it follows that:EQUATION)Hence, ∂ 2 h L in (ω, h) : H → H defines a coercive quadratic form. By definition of Hadamard differentiability, it is also bounded. Therefore, it follows from Lax-Milgram's theorem [Debnath and Mikusinski, 2005, Theorem 4 .3.16], that ∂ 2 h L in (ω, h) is invertible with a bounded inverse. Moreover, recalling that B ω = ∂ ω,h L in (ω, h ⋆ ω) is a bounded operator, its adjoint (B ω ) ⋆ is also a bounded operator from Ω to H. Therefore, we can define J = -C -1 ω (B ω ) ⋆ which is a bounded linear map from Ω to H and will be our candidate Jacobian.→ h ⋆ ω .By the strong convexity assumption (locally in ω), there exists an open ball B centered at the origin 0 that is small enough so that we can ensure the existence of µ > 0 for which h → L in (ω + ϵ, h) is µ-strongly convex for all ϵ ∈ B. For a given ϵ ∈ B, we use the µ-strongmonotonicity of h → ∂ h L in (ω + ϵ, h) (17) at points h ⋆ ω + Jϵ and h ⋆ ω+ϵ to get: µ h ⋆ ω+ϵ -h ⋆ ω -Jϵ 2 ≤ ⟨ ∂ h L in ω + ϵ, h ⋆ ω+ϵ -∂ h L in (ω + ϵ, h ⋆ ω + Jϵ) , h ⋆ ω+ϵ -h ⋆ ω -Jϵ ⟩ H = ⟨-∂ h L in (ω + ϵ, h ⋆ ω + Jϵ) , h ⋆ ω+ϵ -h ⋆ ω -Jϵ ⟩ H ≤ ∥∂ h L in (ω + ϵ, h ⋆ ω + Jϵ)∥ H h ⋆ ω+ϵ -h ⋆ ω -Jϵ H ,where the second line follows from optimality of h ⋆ ω+ϵ (Equation ( 18)), and the last line uses Cauchy-Schwarz's inequality. The above inequality allows us to deduce that:EQUATIONMoreover, since ∂ h L in is Hadamard differentiable, by Proposition C.2 it follows that:EQUATIONwhere the first term vanishes as a consequence of Equation ( 18), sinceEQUATIONwhere∂ 2 h L in (ω, h) : H → H and ∂ ω,h L in (ω, h) : H → Ω are bounded operators and (∂ ω,h L in (ω, h))⋆ denotes the adjoint of ∂ ω,h L in (ω, h). By definition of J, and using Equation ( 23), it follows that:d (ω,h) ∂ h L in (ω, h)(ϵ, Jϵ) = C ω Jϵ + B ω ϵ = 0.Therefore, combining Equation ( 22) with the above equality yields:EQUATIONFinally, combining Equation ( 21) with the above equality directly shows that h⋆ ω+ϵ -h ⋆ ω -Jϵ ≤ 1 µ o(∥ϵ∥). We have shown that ω → h ⋆ ω is differentiable with a Jacobian map ∂ ω h ⋆ ω given by J ⋆ = -B ω C -1 ω .C.3 Proof of the functional adjoint sensitivity in Proposition 2.2Proof of Proposition 2.2. We use the assumptions and definitions from Proposition 2.2 and express the gradient ∇F(ω) using the chain rule:∇F(ω) = ∂ ω L out (ω, h ⋆ ω ) + [∂ ω h ⋆ ω ] ∂ h L out (ω, h ⋆ ω ).The Jacobian ∂ ω h ⋆ ω is the solution of a linear system obtained by applying Theorem 2.1 :∂ ω h ⋆ ω = -B ω C -1 ω .We noteg ω = ∂ ω L out (ω, h ⋆ ω ) and d ω = ∂ h L out (ω, h ⋆ ω ).It follows that the gradient ∇F(ω) can be expressed as:∇F(ω) = g ω + [∂ ω h ⋆ ω ] d ω = g ω + B ω a ⋆ ω a ⋆ ω := -C -1 ω d ω .In other words, the implicit gradient ∇F(ω) can be expressed using the adjoint function a ⋆ ω , which is an element of H and can be defined as the solution of the following functional regression problem:a ⋆ ω = arg min a∈H L adj (ω, a) := 1 2 ⟨a, C ω a⟩ H + ⟨a, d ω ⟩ H .In this section we provides full proofs of Proposition 2.3. We start by stating the assumptions needed on the data distributions and the point-wise losses in Appendix D.1, then provide some differentiation results in Appendix D.2 and conclude with the main proofs in Appendix D.3.Assumption on P and Q. Assumptions on ℓ in .(C) For any ω ∈ Ω, there exists a positive constant µ and a neighborhood B of ω for which ℓ in is µ-strongly convex in its second argument for all (ω ′ , x, y) ∈ B × X × Y.(D) For any ω ∈ Ω, E P |ℓ in (ω, 0, x, y)| + ∥∂ v ℓ in (ω, 0, x, y)∥ 2 < +∞. (E) v → ℓ in (ω, v, x, y) is continuously differentiable for all (ω, x, y) ∈ Ω × X × Y.(F) For any fixed ω ∈ Ω, there exists a constant L and a neighborhood B of ω EQUATIONEQUATION(I) For any ω ∈ Ω, there exists a positive constant C and a neighborhood B of ω s.t. for allEQUATIONAssumptions on ℓ out .(J) For any ω ∈ Ω, E Q [|ℓ out (ω, 0, x, y)|] < +∞. (K) (ω, v) → ℓ out (ω, v, x, y) is jointly continuously differentiable on Ω × V for all (x, y) ∈ X × Y.(L) For any ω ∈ Ω, there exits a neighborhood B of ω and a positive constant C s.t. for all(ω ′ , v, v ′ , x, y) ∈ B × V × V × X × Y we have: ∥∂ ω ℓ out (ω ′ , v, x, y) -∂ ω ℓ out (ω ′ , v ′ , x, y)∥ ≤ C (1 + ∥v∥ + ∥v ′ ∥ + ∥x∥ + ∥y∥) ∥v -v ′ ∥ , ∥∂ v ℓ out (ω ′ , v, x, y) -∂ v ℓ out (ω ′ , v ′ , x, y)∥ ≤ C ∥v -v ′ ∥ , ∥∂ v ℓ out (ω ′ , v, x, y)∥ ≤ C (1 + ∥v∥ + ∥x∥ + ∥y∥) , ∥∂ ω ℓ out (ω ′ , v, x, y)∥ ≤ C 1 + ∥v∥ 2 + ∥x∥ 2 + ∥y∥ 2 .Example. Here we consider the squared error between two vectors v, z in V. Given a map (ω, x, y) → f ω (x, y) defined over Ω × X × Y and taking values in V, we define the following point-wise objective:EQUATION)We assume that for any ω ∈ Ω, there exists a constant C > 0 such that for all ω ′ in a neighborhood of ω and all x, y ∈ X × Y, the following growth assumption holds:EQUATION)This growth assumption is weak in the context of neural networks with smooth activations as discussed by Bińkowski et al. [2018, Appendix C.4 ]. Proposition D.1. Assume that the map ω → f ω (x, y) is continuously differentiable for any x, y ∈ X × Y, and that Equation (28) holds. Additionally, assume that P and Q admit finite second order moments. Then the point-wise objective ℓ in Equation (27) satisfies Assumptions (C) to (L).Proof. We show that each of the assumptions are satisfied by the classical squared error objective.• Assumption (C): the squared error is 1-strongly convex in v, since ∂ 2 v ℓ ⪰ I. Hence, the strong convexity assumption holds with µ = 1.• Assumption (D): For any ω ∈ Ω, we haveE P |ℓ(ω, 0, x, y)| + ∥∂ v ℓ (ω, 0, x, y)∥ 2 = E P 1 2 ∥f ω (x, y)∥ 2 + ∥f ω (x, y)∥ 2 < +∞,which holds by the growth assumption on f ω (x, y), and P having finite second moments.• Assumption (E): With a perturbation u ∈ V we have:ℓ(ω, v + u, x, y) = 1 2 ∥v -z∥ 2 + ⟨v -z, u⟩ + o(∥u∥ 2 ), z = f ω (x, y),with o(∥u∥ 2 ) = 1 2 ∥u∥ 2 . The mapping v → v -z is continuous, thus the assumption holds. • Assumption (F): For any two points v 1 , v 2 ∈ V using the expression of ∂ v ℓ(ω, v, x, y) = v -z with z = f ω (x, y) we have:∥∂ v ℓ(ω, v 1 , x, y) -∂ v ℓ(ω, v 2 , x, y)∥ = ∥(v 1 -z) -(v 2 -z)∥ = ∥v 1 -v 2 ∥,We see that ℓ is L-smooth with L = 1 and the assumption holds.• Assumption (K): By the differentiation assumption on f ω (x, y), with a perturbation ϵ ∈ Ω we can write:f ω+ϵ (x, y) = f ω (x, y) + ∂ ω f ω (x, y)ϵ + o(ϵ).With a perturbation ϵ × u ∈ Ω × V and substituting f ω+ϵ (x, y) with the expression above we have:ℓ(ω + ϵ, v + u, x, y) = 1 2 ∥(v + u) -(f ω (x, y) + ∂ ω f ω (x, y)ϵ + o(ϵ))∥ 2 = 1 2 ∥v -f ω (x, y)∥ 2 + ⟨ϵ, ∂ ω f ω (x, y) ⊤ (f ω (x, y) -v)⟩ + ⟨u, v -f ω (x, y)⟩ + o(∥ϵ∥ + ∥u∥),which allows us to conclude that (ω, v) → ℓ(ω, v, x, y) is continuously differentiable on Ω × V for all x, y ∈ X × Y and the assumption holds.• Assumption (G): With a perturbation ϵ × u ∈ Ω × V using the expression of ∂ v ℓ(ω, v, x, y) we can write:∂ v ℓ(ω + ϵ, v + u, x, y) = (v + u) -f ω+ϵ (x, y) = (v + u) -(f ω (x, y) + ∂ ω f ω (x, y)ϵ + o(ϵ)) = (v -f ω (x, y)) + u -∂ ω f ω (x, y)ϵ + o(ϵ),by continuously differentiable f ω (x, y), we have that the assumption holds.• Assumptions (H) and (I): From the expression of ∂ v ℓ(ω, v, x, y):∂ ω,v ℓ (ω, v, x, y) = ∂ ω (v -f ω (x, y)) = ∂ ω f ω (x, y),then using the expression above and the growth assumption on f ω ′ (x, y) we have that the two assumptions hold.• Assumption (J): For any ω ∈ Ω we have:E Q [|ℓ(ω, 0, x, y)|] = E Q 1 2 ∥f ω (x, y)∥ 2 < +∞,by the growth assumption on f ω (x, y), and P having finite second moments, thus the assumption is verified.• Assumption (L): Using the growth assumption on f ω ′ (x, y), we have the following inequalities:∥∂ ω ℓ (ω ′ , v, x, y) -∂ ω ℓ (ω ′ , v ′ , x, y)∥ = f ω ′ (x, y) ⊤ (v ′ -v) ≤ f ω ′ (x, y) ⊤ + ∥v ′ -v∥ ≤ C (1 + ∥v∥ + ∥v ′ ∥ + ∥x∥ + ∥y∥) ∥v -v ′ ∥ , ∥∂ v ℓ (ω ′ , v, x, y)∥ = ∥v -f ω ′ (x, y)∥ ≤ ∥v∥ + ∥f ω ′ (x, y)∥ ≤ ∥v∥ + C (1 + ∥x∥ + ∥y∥) ≤ C (1 + ∥v∥ + ∥x∥ + ∥y∥) ∥∂ ω ℓ (ω ′ , v, x, y)∥ = ∂ ω f ω ′ (x, y) ⊤ (f ω ′ (x, y) -v) ≤ ∥∂ ω f ω ′ (x, y)∥ ∥(f ω ′ (x, y) -v)∥ ≤ ∥∂ ω f ω ′ (x, y)∥ (∥f ω ′ (x, y)∥ + ∥v∥) ≤ C 1 + ∥v∥ 2 + ∥x∥ 2 + ∥y∥ 2 ,combining the above with L-smoothness of ℓ we can conclude that the assumption holds.The next lemmas show differentiability of L out , L in and ∂ h L in and will be used to prove Proposition 2.3.Lemma D.2 (Differentiability of L in in its second argument). Under Assumptions (D) to (F), the function h → L in (ω, h) is differentiable in H with partial derivative vector ∂ h L in (ω, h) ∈ H given by: ∂ h L in (ω, h) : X → V x → E P [∂ v ℓ in (ω, h(x), x, y) |x] .Proof. We decompose the proof into three parts: verifying that L in is well-defined, identifying a bounded map as candidate for the differential and showing that it is the Fréchet differential of L in .Well-defined objective. Consider (ω, h) in Ω × H. To show that L in (ω, h) is well-defined, we need to prove that ℓ in (ω, h(x), x, y) is integrable under P. We use the following inequalities to control ℓ in (ω, h(x), x, y):|ℓ in (ω, h(x), x, y)| ≤ |ℓ in (ω, h(x), x, y) -ℓ in (ω, 0, x, y)| + |ℓ in (ω, 0, x, y)| = 1 0 dt h(x) ⊤ ∂ v ℓ in (ω, th(x), x, y) + |ℓ in (ω, 0, x, y)| ≤ ∥h(x)∥ 1 0 dt ∥∂ v ℓ in (ω, th(x), x, y) -∂ v ℓ in (ω, 0, x, y)∥ + ∥h(x)∥ ∥∂ v ℓ in (ω, 0, x, y)∥ + |ℓ in (ω, 0, x, y)| ≤ L 2 ∥h(x)∥ 2 + 1 2 ∥h(x)∥ 2 + ∥∂ v ℓ in (ω, 0, x, y)∥ 2 + |ℓ in (ω, 0, x, y)| ,where the first line follows by triangular inequality, the second follows by application of the fundamental theorem of calculus since ℓ in is differentiable by Assumption (E). The third uses Cauchy-Schwarz inequality along with a triangular inequality. Finally, the last line follows using that ℓ in is L-smooth in its second argument, locally in ω and uniformly in x and y by Assumption (F). Taking the expectation under P yields:|L in (ω, h)| ≤E P [|ℓ in (ω, h(x), x, y)|] ≤ L + 1 2 ∥h∥ 2 H + E P ∥∂ v ℓ in (ω, 0, x, y)∥ 2 + |ℓ in (ω, 0, x, y)| < +∞,where ∥h∥ H is finite since h ∈ H and expectations under P of ∥∂ v ℓ in (ω, 0, x, y)∥ 2 and |ℓ in (ω, 0, x, y)| are finite by Assumption (D). This shows that L in (ω, h) is well defined on Ω × H.Candidate differential. Fix (ω, h) in Ω × H and consider the following linear form d in in H:d in g := E P g(x) ⊤ ∂ v ℓ in (ω, h(x), x, y) , ∀g ∈ H.We need to show that it is a bounded form. To this end, we will show that d in is a scalar product with some vector D in in H. The following equalities hold:d in g = E P g(x) ⊤ ∂ v ℓ in (ω, h(x), x, y) = E P g(x) ⊤ E P [∂ v ℓ in (ω, h(x), x, y)|x] = E P g(x) ⊤ D in (x) ,where the second line follows by the "tower" property for conditional expectations and where we defineD in (x) := E P [∂ v ℓ in (ω, h(x), x, y)|x] in the last line. D in is a the candidate representation of d in in H.We simply need to check that D in is an element of H. To see this, we use the following upper-bounds:E P ∥D in (x)∥ 2 ≤ E P E P ∥∂ v ℓ in (ω, h(x), x, y)∥ 2 x = E P ∥∂ v ℓ in (ω, h(x), x, y)∥ 2 ≤ 2E P ∥∂ v ℓ in (ω, h(x), x, y) -∂ v ℓ in (ω, 0, x, y)∥ 2 + 2E P ∥∂ v ℓ in (ω, 0, x, y)∥ 2 ≤ 2L 2 E P ∥h(x)∥ 2 + 2E P ∥∂ v ℓ in (ω, 0, x, y)∥ 2 < +∞.The first inequality is an application of Jensen's inequality by convexity of the squared norm. The second line follows by the "tower" property for conditional probability distributions while the third follows by triangular inequality and Jensen's inequality applied to the square function. The last line uses that ℓ in is L-smooth in its second argument, locally in ω and uniformly in x, y by Assumption (F). Since h is square integrable under P by construction and ∥∂ v ℓ in (ω, 0, x, y)∥ is also square integrable by Assumption (D), we deduce from the above upper-bounds that D in (x) must also be square integrable and thus an element of H. Therefore, we have shown that d in is a continuous linear form admitting the following representation:EQUATION)Differentiability of h → L in (ω, h). To prove differentiability, we simply control the first order error E(g) defined as:EQUATION)For a given g ∈ H, the following inequalities hold:E(g) = E P 1 0 dt g(x) ⊤ (∂ v ℓ in (ω, h(x) + tg(x), x, y) -∂ v ℓ in (ω, h(x), x, y)) ≤ E P 1 0 |g(x) ⊤ (∂ v ℓ in (ω, h(x) + tg(x), x, y) -∂ v ℓ in (ω, h(x), x, y)) | dt ≤ L 2 E P ∥g(x)∥ 2 = L 2 ∥g∥ 2 H ,where the first inequality follows by application of the fundamental theorem of calculus since ℓ in is differentiable in its second argument by Assumption (E). The second line follows by Jensen's inequality while the last line uses that v → ∂ℓ in (ω, v, x, y) is L-Lipschitz locally in ω and uniformly in x and y by Assumption (F). Therefore, we have shown that E(g) = o(∥g∥ H ) which precisely means that h → L in (ω, h) is differentiable with differential d in . Moreover, D in is the partial gradient of L in (ω, h) in the second variable:∂ h L in (ω, h) = D in = x → E P [∂ v ℓ in (ω, h(x), x, y)|x] .Lemma D.3 (Differentiability of L out ). Under Assumptions (A), (B) and (J) to (L), L out is jointly differentiable in ω and h. Moreover, its partial derivatives ∂ ω L out (ω, h) and ∂ h L out (ω, h) are elements in Ω and H given by:EQUATIONProof. We follow a similar procedure as in Lemma D.2, where we decompose the proof into three steps: verifying that the objective L out is well-defined, identifying a candidate for the differential and proving that it is the differential of L out .Well-definiteness of the objective. Let (ω, h) be in Ω × H. First, note that by Assumption (B), we have thatEQUATION)The next inequalities control the growth of ℓ out :|ℓ out (ω, h(x), x, y)| ≤ |ℓ out (ω, 0, x, y)| + |ℓ out (ω, h(x), x, y) -ℓ out (ω, 0, x, y)| ≤ |ℓ out (ω, 0, x, y)| + 1 0 dt h(x) ⊤ ∂ v ℓ out (ω, th(x), x, y) ≤ |ℓ out (ω, 0, x, y)| + ∥h(x)∥ 1 0 dt ∥∂ v ℓ out (ω, th(x), x, y)∥ ≤ |ℓ out (ω, 0, x, y)| + C ∥h(x)∥ (1 + ∥h(x)∥ + ∥x∥ + ∥y∥) ≤ |ℓ out (ω, 0, x, y)| + C 1 + 3 ∥h(x)∥ 2 + ∥x∥ 2 + ∥y∥ 2 .The first line is due to the triangular inequality while the second line follows by differentiability of ∂ v ℓ out in its second argument (Assumption (K)). The third line follows by Cauchy-Scwharz inequality wile the fourth line uses that ℓ out has at most a linear growth in its last three arguments by Assumption (L). Using the above inequalities, we get the following upper-bound on L out :EQUATIONIn the above upper-bound, E Q ∥h(x)∥ 2 is finite by Equation (32). Additionally, E Q ∥x∥ 2 + ∥y∥ 2 is finite since Q has finite second moments by Assumption (A) while E Q [|ℓ out (ω, 0, x, y)|] is also finite by Assumption (J). Therefore, L out is well defined over Ω × H.Candidate differential. Fix (ω, h) in Ω × H and define the following linear form:d out (ϵ, g) := ϵ ⊤ E Q [∂ ω ℓ out (ω, h(x), x, y)] + E Q g(x) ⊤ ∂ v ℓ out (ω, h(x), x, y)Define D out = (D ω , D h ) to be:D ω := E Q [∂ ω ℓ out (ω, h(x), x, y)] D h := x → r(x)E Q [∂ v ℓ out (ω, h(x), x, y)|x] .By an argument similar to the one in Lemma D.2, we see thatd out (ϵ, g) = ⟨g, D h ⟩ H + ϵ ⊤ D ω .We now need to show that D ω and D h are well defined elements of Ω and H.Square integrability of D h . We use the following upper-bounds:E P ∥D h (x)∥ 2 ≤ E P r(x) 2 E Q [∥∂ v ℓ out (ω, h(x), x, y)∥|x] 2 ≤ E P r(x) 2 E Q ∥∂ v ℓ out (ω, h(x), x, y)∥ 2 x ≤ M E P r(x)E Q ∥∂ v ℓ out (ω, h(x), x, y)∥ 2 x = M E Q ∥∂ v ℓ out (ω, h(x), x, y)∥ 2 ≤ 4M C 1 + E Q ∥h(x)∥ 2 + E Q ∥x∥ 2 + ∥y∥ 2 .The first inequality is an application of Jensen's inequality by convexity of the norm, while the second one is an application of Cauchy-Schwarz inequality. The third line uses that r(x) is upper-bounded by a constant M by Assumption (B), and the fourth line follows from the "tower" property for conditional probability distributions. Finally, the last line follows by Assumption (L) which ensures that ∂ v ℓ out has at most a linear growth in its last three arguments. By Equation ( 32), we have thatE Q ∥h(x)∥ 2 < +∞.Moreover, since Q has finite second order moment by Assumption (A), we also have that E Q ∥x∥ 2 + ∥y∥ 2 < +∞. We therefore conclude that E P ∥D h (x)∥ 2 is finite which ensure that D h belongs to H.Well-definiteness of D ω . To show that D ω is well defined, we need to prove that (x, y) → ∂ ω ℓ out (ω, h(x), x, y) is integrable under Q. By Assumption (L), we know that ∂ ω ℓ out has at most a quadratic growth in it last three arguments so that the following inequality holds.∥∂ ω ℓ out (ω, h(x), x, y)∥ ≤ C 1 + ∥h(x)∥ 2 + ∥x∥ 2 + ∥y∥ 2 .We can directly conclude by taking the expectation under Q in the above inequality and recalling that E Q ∥h(x)∥ 2 is finite by Equation ( 32), and that Q has finite second-order moments by Assumption (A).Differentiability of L out . Since differentiability is a local notion, we may assume without loss of generality that ∥ϵ∥ 2 +∥g∥ 2 H ≤ 1. Introduce the functions ∆ 1 and ∆ 2 defined over Ω×H, X ×Y×[0, 1] as follows:∆ 1 (ϵ, g, x, y, t) := ∂ v ℓ out (ω + tϵ, h(x) + tg(x), x, y) -∂ v ℓ out (ω + tϵ, h(x), x, y) ∆ ′ 1 (ϵ, g, x, y, t) := ∂ v ℓ out (ω + tϵ, h(x), x, y) -∂ v ℓ out (ω, h(x), x, y) ∆ 2 (ϵ, g, x, y, t) := ∂ ω ℓ out (ω + tϵ, h(x) + tg(x), x, y) -∂ ω ℓ out (ω + tϵ, h(x), x, y) ∆ ′ 2 (ϵ, g, x, y, t) := ∂ ω ℓ out (ω + tϵ, h(x), x, y) -∂ ω ℓ out (ω, h(x), x, y) .We consider the first-order error E(ϵ, g) which admits the following upper-bounds:E(ϵ, g) := |L out (ω + ϵ, h + g) -L out (ω, h) -d out (ϵ, g)| = E Q 1 0 dt g(x) ⊤ (∆ 1 + ∆ ′ 1 )(ϵ, g, x, y, t) + ϵ ⊤ (∆ 2 + ∆ ′ 2 )(ϵ, g, x, y, t) ≤E Q 1 0 dt ∥g(x)∥ (∥∆ 1 (ϵ, g, x, y, t)∥ + ∥∆ ′ 1 (ϵ, g, x, y, t)∥) + E Q 1 0 dt ∥ϵ∥ (∥∆ 2 (ϵ, g, x, y, t)∥ + ∥∆ ′ 2 (ϵ, g, x, y, t)∥) ≤E Q ∥g(x)∥ 2 1 2      E Q 1 0 dt ∥∆ 1 (ϵ, g, x, y, t)∥ 2 1 2 A1(ϵ,g) + E Q 1 0 dt ∥∆ ′ 1 (ϵ, g, x, y, t)∥ 2 1 2 A2(ϵ,g)      + ∥ϵ∥      E Q 1 0 dt ∥∆ 2 (ϵ, g, x, y, t)∥ A3(ϵ,g) + E Q 1 0 dt ∥∆ ′ 2 (ϵ, g, x, y, t)∥ A4(ϵ,g)      ≤M ∥g∥ H (A 1 (ϵ, g) + A 2 (ϵ, g)) + ∥ϵ∥ (A 3 (ϵ, g) + A 4 (ϵ, g)) .The second line uses differentiability of ℓ out (Assumption (K)). The third uses the triangular inequality, while the fourth line uses Cauchy-Schwarz inequality. Finally, the last line uses Equation (32).We simply need to show that each of the terms A 1 , A 2 , A 3 and A 4 converge to 0 as ϵ and g converge to 0. We treat each term separately.Controlling A 1 and A 3 . For ϵ small enough so that Assumption (L) holds, the following upperbounds on A 1 and A 2 hold:A 1 (ϵ, g) ≤CE Q ∥g(x)∥ 2 1 2 ≤CM 1 2 ∥g∥ H A 3 (ϵ, g) ≤CE Q [(1 + ∥h(x) + tg(x)∥ + ∥h(x)∥ + ∥x∥ + ∥y∥) ∥g(x)∥] ≤CE Q ∥g(x)∥ 2 1 2 1 + 2E Q ∥h(x)∥ 2 1 2 + E Q ∥g(x)∥ 2 1 2 + E Q ∥x∥ 2 + ∥y∥ 2 1 2 ≤CM 1 2 ∥g∥ H 1 + M 1 2 + 2E Q ∥h(x)∥ 2 1 2 + E Q ∥x∥ 2 + ∥y∥ 2 1 2.For A 1 , we used that ∂ v ℓ o ut has is Lipschitz continuous in its second argument for any x, y ∈ X × Y and locally in ω by Assumption (L). The second upper-bound on A 1 uses Equation (32). For A 2 , we used the locally Lipschitz property of ∂ ω ℓ out from Assumption (L), followed by Cauchy-Schwarz inequality and Equation (32). For the last line, we also used that ∥g∥ H ≤ 1 by assumption. The above upper-bounds on A 1 and A 3 ensure that these quantities converge to 0 as ϵ and g approach 0.Controlling A 2 and A 4 . To show that A 2 and A 4 converge to 0, we will use the dominated convergence theorem. It is easy to see that ∆ ′ 1 (ϵ, g, x, y, t) and ∆ ′ 2 (ϵ, g, x, y, t) converge point-wise to 0 when ϵ and g converge to 0 since (ω, v) → ∂ v ℓ out (ω, v, x, y) and (ω, v) → ∂ ω ℓ out (ω, v, x, y) are continuous by Assumption (K). It remains to dominate these functions. For ϵ small enough so that Assumption (L) holds, we have that:∆ ′ 1 (ϵ, g, x, y, t) 2 ≤16C 2 1 + ∥h(x)∥ 2 + ∥x∥ 2 + ∥y∥ 2 ∆ ′ 2 (ϵ, g, x, y, t) ≤2C 1 + ∥h(x)∥ 2 + ∥x∥ 2 + ∥y∥ 2 .Both upper-bounds are integrable under Q since E Q ∥h(x)∥ 2 < +∞ by Equation ( 32) and Q has finite second-order moment by Assumption (A). Therefore, by the dominated convergence theorem, we deduce that A 2 and A 4 converge to 0 as ϵ and g approach 0.Finally, we have shown that E(ϵ, g) = o (∥ϵ∥ + ∥g∥ H ) which allows to conclude that L out is differentiable with the partial derivatives given by Equation (31).Lemma D.4 (Differentiability of ∂ h L in ). Under Assumptions (A) and (E) to (I), the differential mapEQUATIONwhere∂ 2 h L in (ω, h) : H → H is a linear symmetric operator representing the partial derivative of ∂ h L in (ω, h) w.r.t h and (∂ ω,h L in (ω, h)) ⋆ is the adjoint of ∂ ω,h L in (ω, h) : H → Ω which represents the partial derivative of ∂ h L in (ω, h) w.r.t ω. Moreover, ∂ 2 h L in (ω, h) and ∂ ω,h L in (ω, h) are given by:EQUATIONProof. Let (ω, h) be in Ω × H. To show that ∂ ω L in is Hadamard differentiable, we proceed in two steps: we first identify a candidate differential and show that it is a bounded operator, then we prove Hadamard differentiability.Candidate differential. For a given (ω, h) ∈ Ω × H, we consider the following linear operators C w,h : H → H and B w,h : H → Ω:C ω,h g = E P ∂ 2 v ℓ in (ω, h(x), x, y) x g(x), B ω,h g = E P [∂ ω,v ℓ in (ω, h(x), x, y)g(x)] ,where the expectations are over y conditionally on x. Next, we show that C ω,h and B ω,h are well-defined and bounded.Well-definiteness of the operator C ω,h . The first step is to show that the image C ω,h g of any element g ∈ H by C ω,h is also an element in H. To this end, we simply need to find a finite upper-bound on ∥C ω,h g∥ H for a given g ∈ H:∥C ω,h g∥2 H = E P E P ∂ 2 v ℓ in (ω, h(x), x, y) x g(x) 2 ≤ E P E P ∂ 2 v ℓ in (ω, h(x), x, y) x 2 op ∥g(x)∥ 2 ≤ E P E P ∂ 2 v ℓ in (ω, h(x), x, y) op x 2 ∥g(x)∥ 2 ≤ E P ∂ 2 v ℓ in (ω, h(x), x, y) 2 op ∥g(x)∥ 2 ≤ L 2 ∥g∥ 2 H .The second line follows using the operator norm inequality, the third line follows by Jensen's inequality applied to the norm, while the fourth uses the "tower" property for conditional distributions. Finally, the last line uses that ∂ 2 v ℓ in is upper-bounded uniformly in x and y by Assumption (F). Therefore, we conclude that C ω,h g belongs to H. Moreover, the inequality ∥C ω,h g∥ H ≤ L ∥g∥ H also establishes the continuity of the operator C ω,h .Well-definiteness of the operator B ω,h . We first show that the image B ω,h is bounded. For a given g in H, we write:∥B ω,h g∥ = ∥E P [∂ ω,v ℓ in (ω, h(x), x, y)g(x)]∥ ≤E P [∥∂ ω,v ℓ in (ω, h(x), x, y)g(x)∥] ≤E P ∥∂ ω,v ℓ in (ω, h(x), x, y)∥ op ∥g(x)∥ ≤ ∥g∥ H E P ∥∂ ω,v ℓ in (ω, h(x), x, y)∥ 2 op 1 2 ≤ ∥g∥ H E P ∥∂ ω,v ℓ in (ω, h(x), x, y) -∂ ω,v ℓ in (ω, 0, x, y)∥ 2 op 1 2 + ∥g∥ H E P ∥∂ ω,v ℓ in (ω, 0, x, y)∥ 2 op 1 2 ≤C ∥g∥ H E P ∥h(x)∥ 2 1 2 + E P (1 + ∥x∥ + ∥y∥) 2 1 2 ≤C ∥g∥ H ∥h∥ H + 2E P 1 + ∥x∥ 2 + ∥y∥ 2 < +∞.In the above expression, the second line is due to Jensen's inequality applied to the norm function, the third line follows from the operator norm inequality, while the fourth follows by Cauchy-Schwarz. The fifth line is due to the triangular inequality. Finally, the sixth line relies on two facts: 1) that v → ∂ ω,v ℓ in (ω, v, x, y) is Lipschitz uniformly in x and y and locally in ω by Assumption (I), and, 2) that ∥∂ ω,v ℓ in (ω, 0, x, y)∥ has at most a linear growth in x and y locally in ω by Assumption (H). Since P has finite second order moments by Assumption (A) and both h and g are square integrable, we conclude that the constant ∥B ω,h ∥ is finite. Moreover, the last inequality establishes that B ω,h is a continuous linear operator from H to Ω. One can then see that the adjoint of B ω,h admits a representation of the form:(B ω,h ) ⋆ ϵ := (∂ ω,h L in (ω, h)) ⋆ ϵ = x → E P (∂ ω,v ℓ in (ω, h(x), x, y)) ⊤ x ϵ.Therefore, we can consider the following candidate operator d 2 in for the differential of ∂ h L in :d 2 in (ϵ, g) := C ω,h g + (B ω,h ) ⋆ ϵ.Differentiablity of ∂ h L in . We will show that ∂ h L in is jointly Hadamard differentiable at (ω, h) with differential operator given by:EQUATION)To this end, we consider a sequence (ϵ k , g k ) k≥1 converging in Ω × H towards an element (ϵ, g) ∈ Ω × H and a non-vanishing real valued sequence t k converging to 0. Define the first-order error E k as follows:E k := 1 t k (∂ h L in (ω + t k ϵ k , h + t k g k ) -∂ h L in (ω, h)) -C ω,h g -(B ω,h ) ⋆ ϵ 2 H .Introduce the functions P 1 , P 2 , ∆ 1 and ∆ 2 defined over N ⋆ , X × Y × [0, 1] as follows:P 1 (k, x, y, s) = ∂ 2 v ℓ in (ω + st k ϵ k , h(x) + st k g k (x), x, y), k ≥ 1 ∂ 2 v ℓ in (ω, h(x), x, y), k = 0 P 2 (k, x, y, s) = (∂ ω,v ℓ in (ω + st k ϵ k , h(x) + st k g k (x), x, y)) ⊤ k ≥ 1 (∂ ω,v ℓ in (ω, h(x), x, y)) ⊤ , k = 0. ∆ 1 (k, x, y, s) = P 1 (k, x, y, s) -P 1 (0, x, y, s), ∆ 2 (k, x, y, s) = P 2 (k, x, y, s) -P 2 (0, x, y, s).By joint differentiability of (ω, v) → ∂ v ℓ in (ω, v, x, y) (Assumption (E)), we use the fundamental theorem of calculus to express E k in terms of ∆ 1 and ∆ 2 :E k =E P E P 1 0 dt (P 1 (k, x, y, s)g k (x)-P 1 (0, x, y, s)g(x) + P 2 (k, x, y, s)ϵ k -P 2 (0, x, y, s)ϵ) x 2 ≤E P E P 1 0 dt ∥P 1 (k, x, y, s)g k (x)-P 1 (0, x, y, s)g(x) + P 2 (k, x, y, s)ϵ k -P 2 (0, x, y, s)ϵ∥ 2 x =E P 1 0 dt ∥P 1 (k, x, y, s)g k (x) -P 1 (0, x, y, s)g(x) + P 2 (k, x, y, s)ϵ k -P 2 (0, x, y, s)ϵ∥ 2 ≤4 E P 1 0 dt ∥∆ 1 (k, x, y, t)∥ 2 op ∥g(x)∥ 2 A (1) k +4 ∥ϵ∥ 2 E P 1 0 dt ∥∆ 2 (k, x, y, t)∥ 2 op B (1) k + 4 E P 1 0 dt ∥P 1 (k, x, y, t)∥ 2 op ∥g(x) -g k (x)∥ 2 A (2) k +4 ∥ϵ -ϵ k ∥ 2 E P 1 0 dt ∥P 2 (k, x, y, t)∥ 2 op B (2) k .The second line uses Jensen's inequality applied to the squared norm, the fourth line results from the "tower" property of conditional distributions. The fifth line uses Jensen's inequality for the square function followed by the operator norm inequality. It remains to show that A(1)EQUATIONk and A(2) k converge to 0 and that B(2)k is bounded. Upper-bound on A (1)k . We will use the dominated convergence theorem. Assumption (F) ensures the existence of a positive constant L and a neighborhoodB of ω so that v → ∂ 2 v ℓ in (ω ′ , v, x, y) op is bounded by L for any ω ′ , x, y ∈ B × X × Y. Since ω + t k ϵ k → ω,then there exists some K 0 so that, for any k ≥ K 0 , we can ensure that ω + t k ϵ k ∈ B. This allows us to deduce that:EQUATION)for any k ≥ K 0 and any x, y ∈ X × Y, with ∥g(x)∥ 2 being integrable under P.Moreover, we also have the following point-wise convergence for P-almost all x ∈ X :EQUATION)Equation ( 39) follows by noting that ω + t k ϵ k → ω and that h(x) + t k g k (x) → h(x) for P-almost all x ∈ X , since t k converges to 0, ϵ k converges to ϵ and g k converges to g in H (a fortiori converges point-wise for P-almost all x ∈ X ). Additionally, the map (ω, v) → ∂ 2 v ℓ in (ω, v, x, y) op is continuous by Assumption (G), which allows to establish Equation (39). From Equations ( 38) and ( 39) we can apply the dominated convergence theorem which allows to deduce that A(1)k → 0. Upper-bound on A (2)k . By a similar argument as for A(1) k and using Assumption (F), we know that there exists K 0 > 0 so that for any k ≥ K 0 :EQUATION)Therefore, we directly get that:EQUATIONwhere we used that g k → g by construction.Upper-bound on B(2)k . We will show that ∥P 2 (k, x, y, t)∥ op is upper-bounded by a square integrable function under P. By Assumptions (H) and (I), there exists a neighborhood B and a positive constantC such that, for all ω ′ , v 1 , v 2 , x, y ∈ B × V × V × X × Y: ∥∂ ω,v1 ℓ in (ω ′ , 0, x, y)∥ ≤ C (1 + ∥x∥ + ∥y∥) (42) ∥∂ ω,v ℓ in (ω ′ , v 1 , x, y) -∂ ω,v ℓ in (ω ′ , v 2 , x, y)∥ ≤ C ∥v 1 -v 2 ∥ (43)By a similar argument as for A(1) k , there exists K 0 so that for any k ≥ K 0 , the above inequalities hold when choosing ω ′ = ω + t k ϵ k . Using this fact, we obtain the following upper-bound on ∥P 2 (k, x, y, t)∥ op for k ≥ K 0 :∥P 2 (k, x, y, t)∥ op ≤ ∥∂ ω,v ℓ in (ω + st k ϵ k , h(x) + st k g k (x), x, y) -∂ ω,v ℓ in (ω + st k ϵ k , 0, x, y)∥ op + ∥∂ ω,v ℓ in (ω + st k ϵ k , 0, x, y)∥ op ≤C (1 + ∥h(x) + st k g k (x)∥ + ∥x∥ + ∥y∥) ≤C (1 + ∥h(x)∥ + t k ∥g k (x)∥ + ∥x∥ + ∥y∥)Therefore, by taking expectations and integrating over t, it follows:B (2) k ≤ C 2 E P (1 + ∥h(x)∥ + t k ∥g k (x)∥ + ∥x∥ + ∥y∥) 2 ≤ 4C 2 E P 1 + ∥h(x)∥ 2 + t 2 k ∥g k (x)∥ 2 + ∥x∥ 2 + ∥y∥ 2 .By construction t 2 k ∥g k (x)∥ 2 → 0 and is therefore a bounded sequence. Moreover, E P ∥h(x)∥ 2 is finite since h belongs to H. Finally, E P ∥x∥ 2 + ∥y∥ 2 < +∞ by Assumption (A). Therefore, we have shown that B(2)k is bounded. Upper-bound on B (1)k . By a similar argument as for B(2) k and using again Assumptions (H) and (I), there exists K 0 so that for any k ≥ K 0 :∥∆ 2 (k, x, y, t)∥ op ≤ ∥∂ ω,v ℓ in (ω + st k ϵ k , h(x) + st k g k (x), x, y) -∂ ω,v ℓ in (ω + st k ϵ k , h(x), x, y)∥ op + ∥∂ ω,v ℓ in (ω + st k ϵ k , h(x), x, y) -∂ ω,v ℓ in (ω, h(x), x, y)∥ op ≤Ct k ∥g k (x)∥ + ∥∂ ω,v ℓ in (ω + st k ϵ k , h(x), x, y) -∂ ω,v ℓ in (ω, h(x), x, y)∥ op ,where we used Equation ( 43) to get an upper-bound on the first terms. By squaring the above inequality and taking the expectation under P we get:B (1) k ≤ 2Ct k ∥g k ∥ 2 H + 2E P   ∥∂ω,vℓin(ω + st k ϵ k , h(x), x, y) -∂ ω,v ℓ in (ω, h(x), x, y)∥ 2 op e k (x,y)    . (44)We only need to show that E P [e k (x, y)] converges to 0 since the first term 2Ct k ∥g k ∥ 2 H already converges to 0 by construction of t k and g k . To achieve this, we will use the dominated convergence theorem. It is easy to see that e k (x, y) converges to 0 point-wise by continuity of ω → ∂ ω,v ℓ in (ω, v, x, y) (Assumption (G)). Therefore, we only need to show that e k (x, y) is dominated by an integrable function. Provided that k ≥ K 0 , we can use Equations ( 42) and ( 43) to get the following upper-bounds:1 4 e k (x, y) ≤ ∥∂ ω,v ℓ in (ω + st k ϵ k , h(x), x, y) -∂ ω,v ℓ in (ω + st k ϵ k , 0, x, y)∥ 2 op + ∥∂ ω,v ℓ in (ω, h(x), x, y) -∂ ω,v ℓ in (ω, 0, x, y)∥ 2 op + ∥∂ ω,v ℓ in (ω, 0, x, y)∥ 2 op + ∥∂ ω,v ℓ in (ω + st k ϵ k , 0, x, y)∥ 2 op ≤2C 2 1 + ∥h(x)∥ 2 + ∥x∥ 2 + ∥y∥ 2 .The l.h.s. of the last line is an integrable function that is independent of k, since h is square integrable by definition and ∥x∥ 2 + ∥y∥ 2 are integrable by Assumption (A). Therefore, by application of the dominated convergence theorem, it follows that E P [e k (x, y)] → 0, we have shown that B(1) k → 0. To conclude, we have shown that the first-order error E k converges to 0 which means that (ω, h) → ∂ h L in (ω, h) is jointly differentiable on Ω × H, with differential given by Equations ( 34) and ( 35).Proof. The strategy is to show that the conditions on L in and L out stated in Proposition 2.2 hold. By Assumption (C), for any ω ∈ Ω, there exists a positive constant µ and a neighborhood B of ω on which the function ℓ in (ω ′ , v, x, y) is µ-strongly convex in v for any (ω ′ , x, y) ∈ B × X × Y. Therefore, by integration, we directly deduce that h → L in (ω ′ , h) is µ strongly convex in h for any ω ′ ∈ B. By Lemmas D.2 and D.4, h → L in (ω, h) is differentiable on H for all ω ∈ Ω and ∂ h L in is Hadamard differentiable on Ω × H. Additionally, L out is jointly differentiable in ω and h by Lemma D.3. Therefore, the conditions on L in and L out for applying Proposition 2.2 hold. Using the notations from Proposition 2.2, we have that the total gradient ∇F(ω) can be expressed as:∇F(ω) = g ω + B ω a ⋆ ω (45)whereg ω = ∂ ω L out (ω, h ⋆ ω ), B ω = ∂ ω,h L in (ω, h ⋆ ω )and where a ⋆ ω is the minimizer of the adjoint objective L adj :L adj (ω, a) := 1 2 a ⊤ C ω a + a ⊤ d ω , with C ω = ∂ 2 h L in (ω, h ⋆ ω ) and d ω = ∂ h L out (ω, h ⋆ ω ).Recalling the expressions of the first and second order differential operators from Lemmas D.2 and D.4, we deduce the expression of the adjoint objective as a sum of two expectations under P and Q given the optimal prediction functionL adj (ω, a) = 1 2 E (x,y)∼P a(x) ⊤ ∂ 2 v ℓ in (ω, h ⋆ ω (x), x, y) a(x) + E (x,y)∼Q a(x) ⊤ ∂ v ℓ out (ω, h ⋆ ω (x), x, y) .Furthermore, the vectors g ω and B ω a ⋆ ω appearing in Equation ( 45) can also be expressed as expectations:g ω = E (x,y)∼Q [∂ ω ℓ out (ω, h ⋆ ω (x), x, y)] B ω a ⋆ ω = E (x,y)∼P [∂ ω,v ℓ in (ω, h ⋆ ω (x), x, y) a ⋆ ω (x)] .We provide a convergence result of Algorithm 1 to stationary points of F. Our analysis uses the framework of biased stochastic gradient descent (Biased SGD) [Demidovich et al., 2024] where the bias arises from suboptimality errors when solving the inner-level and adjoint problems.Gradient estimators. Recall that the total gradient ∇F(ω) admits the following expression under Assumptions (A) to (L):∇F(ω) = E Q [∂ ω ℓ out (ω, h ⋆ ω (x), x, y)] + E P [∂ ω,v ℓ in (ω, h ⋆ ω (x), x, y) a ⋆ ω (x)] ,We denote by ĝ(ω) the gradient estimator, i.e., the mapping ĝ : Ω → Ω computed by Algorithm 4 and which admits the following expression:ĝ(ω) = 1 |Bout| (x,ỹ)∈Bout ∂ ω ℓ out ω, ĥω (x), x, ỹ + 1 |Bin| (x,y)∈Bin ∂ ω,v ℓ in ω, ĥω (x), x, y âω (x),where B in and B out are samples from P and Q independent from ĥω and âω and independent from each other (i.e. B in ⊥ B out ). Here, there are three independent sources of randomness when computing ĝ(ω): estimation of ĥω and âω in Algorithms 2 and 3, as well as random batches B in and B out . We denote by E[•] the expectation with respect to all random variables appearing in the expression of ĝ(ω), and by E[•| ĥω ] and E[•| ĥω , âω ] the conditional expectations knowing ĥω only and both ĥω and âω . ĝ(ω) is a biased estimator of ∇F(ω) (i.e., E[ĝ(ω)] is not equal to ∇F(ω)), as the bias is due to using sub-optimal solutions ĥω and âω instead of h ⋆ ω and a ⋆ ω in the expression of ĝ(ω). Furthermore, we define G(ω) := E ĝ(ω)| ĥω , âω to be the conditional expectation of ĝ(ω) at ω given estimates ĥω and âω of h ⋆ ω and a ⋆ ω . By independence of B in and B out from ĥω and âω , G(ω) admits the following expression:G(ω) =E Q ∂ ω ℓ out ω, ĥω (x), x, y + E P ∂ ω,v ℓ in ω, ĥω (x), x, y âω (x) ,where the expectation is taken w.r.t. (x, y) ∼ Q.We introduce the approximate adjoint objective Ladj (ω, a) where h ⋆ ω is replaced by ĥω :EQUATION)By independence of the estimator ĥ and the samples B used for computing Ladj (ω, a, ĥω , B) it is easy to see that E Ladj (ω, a, ĥω , B) ĥω = Ladj (ω, a). Hence, it is natural to think of âω as an approximation to the minimizer ãω of Ladj (ω, a) in H.Sub-optimality assumption.The following assumption quantifies the sub-optimality errors made by ĥω and âω .(a) For some positive constants ϵ in and ϵ adj and for all ω ∈ Ω, ĥω and âω satisfy:E L in (ω, ĥω ) -L in (ω, h ⋆ ω ) ≤ ϵ in , E Ladj (ω, âω ) -Ladj (ω, ãw ) ≤ ϵ adj .Assumptions on ℓ in .(b) (Strong convexity) ℓ in (ω, v, x, y) is µ-strongly convex in v ∈ V. (c) v → ∂ 2 v ℓ in (ω, v, x, y) is C 1 -Lipschitz on Ω × X × Y. (d) v → ∂ ω,v ℓ in (ω, v, x, y) is differentiable and C 2 -Lipschitz and bounded by B 2 ∈ R. Assumptions on ℓ out . (e) ∥∂ h L out (ω, h ⋆ ω )∥ H is bounded by a positive constant B 3 . (f) v → ∂ v ℓ out (ω, v, x, y) is C 4 -Lipschitz for all (ω, x, y) ∈ Ω × X × Y. (g) v → ∂ ω ℓ out (ω, v, x, y) is differentiable and C 3 -Lipschitz for all (ω, x, y) ∈ Ω × X × Y. (h) (x, y) → ∂ ω ℓ out (ω, h ⋆ ω (x),x, y) has a variance bounded by σ 2 out for all ω ∈ Ω. Smoothness of the total objective.(i) Function F(ω) is L-smooth for all ω ∈ Ω, and bounded from below by F ⋆ ∈ R. Remark E.1. Assumption (a) reflects the generalization errors made when optimizing L in and Ladj using Algorithms 2 and 3. In the case of over-parameterized networks, these errors can be made smaller by increasing network capacity, number of steps and the number of samples [Allen-Zhu et al., 2019 , Du et al., 2019 , Zou et al., 2020] . Remark E.2. Assumptions (b) to (h) are similar in spirit to those used for analyzing bi-level optimization algorithms (ex: [Arbel and Mairal, 2022a , Assumptions 1 to 5]). In particular, Assumption (e) is even weaker than [Arbel and Mairal, 2022a, Assumptions 2] where ∂ h L out (ω, h) needs to be bounded for all ω, h in Ω × H. For instance, Assumption (e) trivially holds when∂ h L out (ω, h) is a linear transformation of ∂ h L in (ω, h), as is the case for min-max problems. There ∂ h L in (ω, h ⋆ ω ) = 0 so that ∂ h L out (ω, h ⋆ ω ) = 0 is bounded, while ∂ h L out (ω, h)might not be bounded in general.The general strategy of the proof is to first show that the conditions for applying the general convergence result for biased SGD of Demidovich et al. [2024, Theorem 3] hold. We start with Proposition E.3 which shows that the biased gradient estimate ĝ(ω) satisfies the conditions of Demidovich et al. [2024, Assumption 9] . Proposition E.3. Let Assumptions (a) to (h) and Assumptions (A) to (L) hold.EQUATIONEQUATIONwhere c 1 and c 2 are constants defined in Equation (50) and σ 0 is a positive constant defined in Equation (51).Proof. We prove each bound separately.Lower-bound on ∇F(ω) ⊤ E[ĝ(ω)]. Fist note that E[ĝ(ω)] = E[G(ω)].Hence, by direct calculation, we have that:∇F(ω) ⊤ E[ĝ(ω)] =∇F(ω) ⊤ E[G(ω)] = 1 2 ∥∇F(ω)∥ 2 + E ∥G(ω)∥ 2 -1 2 E ∥G(ω) -∇F(ω)∥ 2 ≥ 1 2 ∥∇F(ω)∥ 2 -1 2 (c 1 ϵ in + c 2 ϵ out ), where we use Lemma E.4 to get the last lower-bound.Upper-bound on E[∥ĝ(ω)∥ 2 ]. E ∥ĝ(ω)∥ 2 = E ∥ĝ(ω) -G(ω)∥ 2 + E ∥G(ω) -∇F(ω) + ∇F(ω)∥ 2 ≤E ∥ĝ(ω) -G(ω)∥ 2 + 2E ∥G(ω) -∇F(ω)∥ 2 + 2 ∥∇F(ω)∥ 2 ≤σ 2 0 + 2 (c 1 ϵ in + c 2 ϵ adj ) + 2 ∥∇F(ω)∥ 2 ,where the first line uses that E[ĝ(ω)] = E[G(ω)] and the last line uses Lemmas E.4 and E.5.We can now directly use Proposition E.3 and Assumption (i) on F to prove the Theorem 3.1 using the biased SGD convergence result in [Demidovich et al., 2024, Theorem 3] .Proof of Theorem 3.1. The proof is a direct application of [Demidovich et al., 2024, Theorem 3] given that the variance and bias conditions on the estimator ĝ(ω) are satisfied by Proposition E.3 and that F is L-smooth and has a finite lower-bound F ⋆ ∈ R by Assumption (i).E.3 Bias-variance decomposition.Lemma E.4 (Bias control). Let Assumptions (a) to (g) and Assumptions (A) to (L) hold.EQUATIONwhere c 1 and c 2 are non-negative constants defined in Equation (50).Proof.G(ω) -∇F(ω) = E Q ∂ ω ℓ out ω, ĥω (x), x, y -∂ ω ℓ out (ω, h ⋆ ω (x), x, y) A1 + E P ∂ ω,v ℓ in ω, ĥω (x), x, y âω (x) -∂ ω,v ℓ in (ω, h ⋆ ω (x), x, y) a ⋆ ω (x) =A 1 + E P ∂ ω,v ℓ in ω, ĥω (x), x, y -∂ ω,v ℓ in (ω, h ⋆ ω (x), x, y) a ⋆ ω (x) A2 + E P ∂ ω,v ℓ in ω, ĥω (x), x, y (â ω (x) -a ⋆ ω (x)) A3 .We have the following:∥A 1 ∥ ≤M C 3 ĥω -h ⋆ ω H ∥A 2 ∥ ≤C 2 ∥a ⋆ ω ∥ H ĥω -h ⋆ ω H ≤ C 2 B 3 µ -1 ĥω -h ⋆ ω H ∥A 3 ∥ ≤B 2 E P [∥â ω (x) -a ⋆ ω (x)∥] ≤B 2 (E P [∥â ω (x) -ãω (x)∥] + E P [∥a ⋆ ω (x) -ãω (x)∥]) ≤B 2 ∥â ω -ãω ∥ H + ∥ã ω -a ⋆ ω ∥ L1 .The first inequality holds since v → ∂ ω ℓ out (ω, v, x, y) is C 3 -Lipschitz by Assumption (g) and dQ(•, Y) admits a density w.r.t dP(•, Y) bounded by a positive constant M by Assumption (B). The second inequality holds since v → ∂ ω,v ℓ in (ω, v, x, y) is C 2 -Lipschitz by Assumption (d) and ∥a ⋆ ω ∥ H is upper-bounded by µ -1 B 3 as a consequence of Lemma E.6. Finally, the inequality on ∥A 3 ∥ holds since ∂ ω,v ℓ in (ω, v, x, y) is bounded by a constant B 2 by Assumption (d). Therefore, it holds that the difference between G(ω) and ∇F satisfies:∥G(ω) -∇F(ω)∥ 2 ≤ 3 M C 3 + C 2 B 3 µ 2 ĥω -h ⋆ ω 2 H + 3B 2 2 ∥â ω -ãω ∥ 2 H + 3B 2 2 ∥ã ω -a ⋆ ω ∥ 2 L1 .Taking the expectation over ĥω and âω and using the bounds in Lemma E.6 yields:E ∥G(ω) -∇F(ω)∥ 2 ≤6 M C 3 + C 2 B 3 µ -1 2 µ -1 ϵ in + 6B 2 2 µ -1 ϵ adj + 6B 2 2 µ -1 µ -1 C 4 M + µ -2 C 1 B 3 2 ϵ in .Finally, the upper bound on the bias holds with c 1 and c 2 defined as:EQUATION)Lemma E.5 (Variance control). Let Assumptions (a) to (h) and Assumptions (A) to (L) hold, then the variance is upper-bounded as follows:E ∥ĝ(ω) -G(ω)∥ 2 ≤ σ 2 0, where σ 0 is a positive constant given by:σ 2 0 := 2 |B out | 2C 2 3 µ -1 ϵ in + σ 2 out + 4B 2 2 |B in | µ -1 ϵ adj + 2µ -3 C 2 4 M 2 ϵ in + µ -2 B 2 3 (51)Proof. By definition of ĝ(ω) and G(ω) we have that:E ∥ĝ(ω) -G(ω)∥ 2 = 1 |B out | E ∂ ω ℓ out ω, ĥω (x), x, y -E Q ∂ ω ℓ out ω, ĥω (x ′ ), x ′ , y ′ ĥω 2 V1 + 1 |B in | E ∂ ω,v ℓ in ω, ĥω (x), x, y âω (x) 2 V2 - 1 |B in | E E P ∂ ω,v ℓ in ω, ĥω (x), x, y âω (x) ĥω , âω 2 V3 ≤ V 1 |B out | + V 2 |B in | .Where the first line is a direct consequence of the independence of B in and B out . Moreover, we can bound V 1 and V 2 as follows:V 1 ≤2E ∂ ω ℓ out ω, ĥω (x), x, y -∂ ω ℓ out (ω, h ⋆ ω (x), x, y) 2 + 2E Q ∥∂ ω ℓ out (ω, h ⋆ ω (x), x, y) -E Q [∂ ω ℓ out (ω, h ⋆ ω (x ′ ), x ′ , y ′ )]∥ 2 ≤2C 2 3 E ĥω -h ⋆ ω 2 + 2σ 2 out ≤ 4C 2 3 µ -1 ϵ in + 2σ 2 out V 2 ≤B 2 2 E ∥â ω ∥ 2 H ≤ 2B 2 2 E ∥â ω -ãω ∥ 2 H + E ∥ã ω ∥ 2 H ≤4B 2 2 µ -1 ϵ adj + 2µ -3 C 2 4 M 2 ϵ in + µ -2 B 2 3 .For V 1 , we used that v → ∂ ω ℓ out (ω, v, x, y) is C 3 -Lipschitz uniformly in ω, x and y by Assumption (g), and that (x, y) → ∂ ω ℓ out (ω, h ⋆ ω (x), x, y) has a variance uniformly bounded by σ 2 out as a consequence of Assumption (h). For V 2 , we use Assumption (d) where we have that ∂ ω,v ℓ in is uniformly bounded by a constant B 2 and apply the bounds onE ∥ã ω ∥ 2 H and E ∥â ω -ãω ∥ 2 Hfrom Lemma E.6.Lemma E.6. For the optimal prediction and adjoint functions h ⋆ ω , a ⋆ ω defined in Equations ( 6) and (7), as well as their estimated versions ĥω , âω given in Algorithm 1. Under Assumptions (a) to (f) and Assumptions (A) to (L), the following holds:E ĥω -h ⋆ ω 2 H ≤ 2µ -1 ϵ in , E ∥â ω -ãω ∥ 2 H ≤ 2µ -1 ϵ adj , E ∥ã ω -a ⋆ ω ∥ 2 L1 ≤ 2 µ -2 C 1 B 3 + µ -1 C 4 M 2 µ -1 ϵ in , E ∥ã ω ∥ 2 H ≤ 4µ -3 C 2 4 M 2 ϵ in + 2µ -2 B 2 3 , ∥a ⋆ ω ∥ H ≤ µ -1 B 3 .Proof. We show each of the upper bounds separately.Upper-bound on E ĥω -h ⋆ ω 2 H. We use the strong-convexity Assumption (b) and Assumption (a) to show the first bound:EQUATIONUpper-bound on E ∥â ω -ãω ∥ 2 H . The second bound can be proven in the same way as the first, using that, by definition, Ladj is continuous and µ-strongly convex in a together with Assumption (a).E ∥ã w -a ⋆ ω ∥ 2 L1 .We exploit the closed form expressions of ãω and a ⋆ ω :ãω = -     ∂ 2 h L in ω, ĥω Ĥω     -1 ∂ h L out ω, ĥω bω , a ⋆ ω = -   ∂ 2 h L in (ω, h ⋆ ω ) Hω   -1 ∂ h L out (ω, h ⋆ ω ) bω .By standard linear algebra, the difference ãω -a ⋆ ω can be expressed as:ãω -a ⋆ ω = H -1 ω Ĥω -H ω Ĥ-1 ω b ω + Ĥ-1 ω b ω -bω .By taking the L 1 (P) norm of the above and using the upper-bounds from Lemma E.7, we can write:∥ã ω -a ⋆ ω ∥ L1 ≤ H -1 ω Ĥω -H ω Ĥ-1 ω b ω L1 + Ĥ-1 ω b ω -bω L1 ≤ µ -1 Ĥω -H ω Ĥ-1 ω b ω L1 + µ -1 b ω -bω L1 ≤ µ -1 C 1 ĥω -h ⋆ ω H Ĥ-1 ω b ω H + µ -1 b ω -bω H ≤ µ -2 C 1 ĥω -h ⋆ ω H ∥b ω ∥ H + µ -1 C 4 M ĥω -h ⋆ ω H ≤ µ -2 C 1 B 3 + µ -1 C 4 M ĥω -h ⋆ ω H .The first line we used the triangular inequality, the second line from Equation ( 53 52), we get:E ∥ã ω -a ⋆ ω ∥ 2 L1 ≤ µ -2 C 1 B 3 + µ -1 C 4 M 2 E ĥω -h ⋆ ω 2 H ≤ 2 µ -2 C 1 B 3 + µ -1 C 4 M 2 µ -1 ϵ in .Upper-bound on E ∥ã ω ∥ 2 H . We use the closed-form expression of ãω :∥ã ω ∥ H = Ĥ-1 ω bω H ≤ µ -1 bω H ≤ µ -1 bω -b ω H + ∥b ω ∥ H ≤ µ -1 C 4 M ĥω -h ⋆ ω H + B 3 ,where the first line uses that Ĥ-1 ω op ≤ µ -1 by Assumption (b), the second line follows by triangular inequality while the last line uses Equation (55) from Lemma E.7 for the first term and Assumption (e) for the second terms. By squaring the above bound and taking expectation w.r.t. ĥω , we get:E ∥ã ω ∥ 2 H ≤ E µ -1 C 4 M ĥω -h ⋆ ω + µ -1 B 3 2 ≤ 2µ -2 C 2 4 M 2 E ĥω -h ⋆ ω 2 + 2µ -2 B 2 3 ≤ 4µ -3 C 2 4 M 2 ϵ in + 2µ -2 B 2 3 ,where the last line uses Equation ( 52).Upper-bound on ∥a ⋆ ω ∥ H . Using the closed-form expression of a ⋆ , it holds that:∥a ⋆ ω ∥ H = H -1 ω b ω H ≤ µ -1 B 3 ,where we used that H -1 ω op ≤ µ -1 by Assumption (b) and that ∥b ω ∥ H ≤ B 2 by Assumption (d).Lemma E.7. Consider the operators H ω , Ĥω , b ω , bω defined in Lemma E.6. Under Assumptions (b), (c) and (f) and Assumptions (A) to (L), the following holds for any (ω, s) ∈ Ω × H:EQUATIONb ω -bωEQUATION)Proof. We show each of the upper bounds separately.Upper-bound on H -1 ω s L1 and Ĥ-1 ω s. Using strong convexity Assumption (b), we have the following inequality for the Hessian operator H ω acting on a function c ∈ H ⊂ L 1 (P):∥H ω c∥ L1 = E P E P ∂ 2 v ℓ in (ω, h ⋆ ω (x), x, y) x c(x) ≥ µ ∥c∥ L1 .by positive-definiteness of H ω , we take c = H -1 ω s for some s ∈ H:H -1 ω s L1 ≤ µ -1 ∥s∥ L1 .By the same arguments, the above bound applies to Ĥ-1 ω s.Ĥω -H ω c. We can express the operator Ĥω -H ω acting on some s ∈ H as follows:Ĥω -H ω s = E P E P ∂ 2 v ℓ in (ω, ĥω (x), x, y) -∂ 2 v ℓ in (ω, h ⋆ ω (x), x, y) x s(x) .Using Assumption (c) we upper-bound the L 1 norm of the above quantity as follows:Ĥω -H ω s L1 ≤ C 1 E P ĥω (x)) -h ⋆ ω (x) ∥s(x)∥ ≤ C 1 ĥω -h ⋆ ω H ∥s∥ H ,where we used Cauchy-Schwarz inequality to get the last inequality.Upper-bound on b ω -bω. Using Lemma D.3 to get an expression of b ω and bω , we obtain the following upper-bound:b ω -bω 2 H = E P E Q ∂ v ℓ out (ω, h ⋆ ω (x), x, y) -∂ v ℓ out (ω, ĥω (x), x, y) |x r(x) 2 ≤ C 2 4 M 2 ĥω -h ⋆ ω 2 H ,where we used that the density r(x) is upper-bounded by a positive constant M by Assumption (B) and that ∂ v ℓ out is C 4 Lipschitz in its second argument by Assumption (f).F.1 Parametric approximation of the functional bilevel problemIn this section, we approximate the functional problem in Equation (FBO) with a parametric bilevel problem where inner-level functions are parametrized as h(x) = τ (θ)(x) with parameters θ. Here, the inner-level variable is θ instead of the function h. Standard bilevel optimization algorithms like AID can be applied, which involve differentiating twice with respect to the parametric model. However, for models like deep neural networks, the inner objective may not be strongly convex in θ, leading to a non-positive or degenerate Hessian (Proposition F.1). This can cause numerical instabilities and divergence from the gradient in Equation (4) (Proposition F.2), especially when using AID, which relies on solving a quadratic problem defined by the Hessian.If the model has multiple solutions, the Hessian may be degenerate, making the implicit function theorem inapplicable. In contrast, functional implicit differentiation requires solving a positive definite quadratic problem in H to find an adjoint function a ⋆ ω , ensuring a solution even when h ⋆ ω is sub-optimal, due to the strong convexity of L in (ω, h). This stability with sub-optimal solutions is crucial for practical algorithms like the one in Section 3, where the optimal prediction function is approximated within a parametric family, such as neural networks.Formally, to establish a connection with parametric implicit differentiation, let us consider τ : Θ → H to be a map from a finite dimensional set of parameters Θ to the functional Hilbert space H and define a parametric version of the outer and inner objectives in Equation (FBO) restricted to functions inH Θ := {τ (θ) | θ ∈ Θ}: G out (ω, θ) := L out (ω, τ (θ)) G in (ω, θ) := L in (ω, τ (θ)) .(56) The map τ can typically be a neural network parameterization and allows to obtain a "more tractable" approximation to the abstract solution h ⋆ ω in H where the function space H is often too large to perform optimization. This is typically the case when H is an L 2 -space of functions as we discuss in more details in Section 3. When H is a Reproducing Kernel Hilbert Space (RKHS), τ may also correspond to the Nyström approximation [Williams and Seeger, 2000] , which performs the optimization on a finite-dimensional subspace of an RKHS spanned by a few data points.The corresponding parametric version of the problem (FBO) is then formally defined as:min ω∈Ω G tot (ω) := G out (ω, θ ⋆ ω ) s.t. θ ⋆ ω ∈ arg min θ∈Θ G in (ω, θ). (PBO)The resulting bilevel problem in Equation (PBO) often arises in machine learning but is generally ambiguously defined without further assumptions on the map τ as the inner-level problem might admit multiple solutions [Arbel and Mairal, 2022b] . Under the assumption that τ is twice continuously differentiable and the rather strong assumption that the parametric Hessian ∂ 2 θ G in (ω, θ ⋆ ω ) is invertible for a given ω, the expression for the total gradient ∇ ω G tot (ω) follows by direct application of the parametric implicit function theorem [Pedregosa, 2016] :EQUATIONwhere u ⋆ ω is the adjoint vector in Θ. Without further assumptions, the expression of the gradient in Equation ( 57) is generally different from the one obtained in Proposition 2.2 using the functional point of view. Nevertheless, a precise connection between the functional and parametric implicit gradients can be obtained under expressiveness assumptions on the parameterization τ , as discussed in the next two propositions. Proposition F.1. Under the same assumptions as in Proposition 2.2 and assuming that τ is twice continuously differentiable, the following expression holds for any (ω, θ) ∈ Ω × Θ:EQUATIONwhere ∂ 2 θ τ (θ) is a linear operator measuring the distortion induced by the parameterization and acts on functions in H by mapping them to a matrix p × p where p is the dimension of the parameter space Θ. If, in addition, τ is expressive enough so that τ (θ ⋆ ω ) = h ⋆ ω , then the above expression simplifies to:EQUATION)Proposition F.1 follows by direct application of the chain rule, noting that the distortion term on the right of (58) vanishes when θ= θ ⋆ ω since ∂ h L in (ω, τ (θ ⋆ ω )) = ∂ h L in (ω, h ⋆ ω ) = 0 by optimality of h ⋆ ω .A consequence is that, for an optimal parameter θ ⋆ ω , the parametric Hessian is necessarily symmetric positive semi-definite. However, for an arbitrary parameter θ, the distortion does not vanish in general, making the Hessian possibly non-positive. This can result in numerical instability when using algorithms such as AID for which an adjoint vector is obtained by solving a quadratic problem defined by the Hessian matrix ∂ 2 θ G in evaluated on approximate minimizers of the inner-level problem. Moreover, if the model admits multiple solutions θ ⋆ ω , the Hessian is likely to be degenerate making the implicit function theorem inapplicable and the bilevel problem in Equation (PBO) ambiguously defined1 . On the other hand, the functional implicit differentiation requires finding an adjoint function a ⋆ ω by solving a positive definite quadratic problem in H which is always guaranteed to have a solution even when the inner-level prediction function is only approximately optimal. Proposition F.2. Assuming that τ is twice continuously differentiable and that for a fixed ω ∈ Ω we have τ (θ ⋆ ω ) = h ⋆ ω , and J ω := ∂ θ τ (θ ⋆ ω ) has a full rank, then, under the same assumptions as in Proposition 2.2, ∇ ω G tot (ω) is given by:EQUATIONwhere P ω : H → H is a projection operator of rank dim(Θ). If, in addition, the equality τ (θ⋆ ω ′ ) = h ⋆ ω ′ holds for all ω ′ in a neighborhood of ω, then ∇ ω G tot (ω) := ∇F(ω) = g ω + B ω a ⋆ ω .Proposition F.2, which is proven below, shows that, even when the parametric family is expressive enough to recover the optimal prediction function h ⋆ ω at a single value ω, the expression of the total gradient in Equation ( 60) using parametric implicit differentiation might generally differ from the one obtained using its functional counterpart. Indeed the projector P ω , which has a rank equal to dim(Θ), biases the adjoint function by projecting it into a finite dimensional space before applying the cross derivative operator. Only under a much stronger assumption on τ , requiring it to recover the optimal prediction function h ⋆ ω in a neighborhood of the outer-level variable ω, both parametric and functional implicit differentiation recover the same expression for the total gradient. In this case, the projector operator aligns with the cross-derivative operator so that B ω P ω = B ω . Finally, note that the expressiveness assumptions on τ made in Propositions F.1 and F.2 are only used here to discuss the connection with the parametric implicit gradient and are not required by the method we introduce in Section 3.Proof of Proposition F.2. Here we want to show the connection between the parametric gradient of the outer variable ∇ ω G tot (ω) usually used in approximate differentiation methods and the functional gradient of the outer variable ∇F(ω) derived from the functional bilevel problem definition in Equation (FBO). Recall the definition of the parametric inner objective G in (ω, θ) := L in (ω, τ (θ)). According to Proposition F.1, we have the following relation∂ 2 θ G in (ω, θ ⋆ ω ) :=J ω C ω J ⊤ ω with J ω := ∂ θ τ (θ ⋆ ω ).By assumption, J ω has a full rank which matches the dimension of the parameter space Θ. Recall from the assumptions of Theorem 2.1 that the Hessian operator C ω is positive definite by the strong convexity of the inner-objective L in in the second argument. We deduce that ∂ 2 θ G in (ω, θ ⋆ ω ) must be invertible, since, by construction, the dimension of Θ is smaller than that of the Hilbert space H which has possibly infinite dimension. Recall from Theorem 2.1, B ω := ∂ ω,h L in (ω, h ⋆ ω ) and the assumption that τ (θ ⋆ ω ) = h ⋆ ω . We apply the parametric implicit function theorem to get the following expression of the Jacobian ∂ ω θ ⋆ ω :∂ ω θ ⋆ ω := -B ω J ⊤ ω J ω C ω J ⊤ ω -1 .Hence, differentiating the total objective G tot (ω) := G out (ω, θ ⋆ ω ) = L out (ω, τ (θ ⋆ ω )) and applying the chain rule directly results in the following expression:EQUATIONwith previously defined gω := ∂ ω L out (ω, h ⋆ ω ) and d ω := ∂ h L out (ω, h ⋆ ω ).We now introduce the operator Pω := J ⊤ ω J ω C ω J ⊤ ω -1 J ω C ω .The operator P ω is a projector as it satisfies P 2 ω = P ω . Hence, using the fact that the Hessian operator is invertible, and recalling that the adjoint function is given by a ⋆ ω = -C -1 ω d ω , we directly get form Equation (61) that:∇ ω G tot (ω) := g ω + B ω P ω a ⋆ ω .If we further assume that τ (θ ⋆ ω ′ ) = h ⋆ ω ′ holds for all ω ′ in a neighborhood of ω, then differentiating with respect to ω results in the following identity:∂ ω θ ⋆ ω J ω = ∂ ω h ⋆ ω .Using the expression of ∂ ω h ⋆ ω from Equation (3), we have the following identity:-∂ ω θ ⋆ ω J ω C ω = B ω .In other words, B ω is of the form B ω := DJ ω C ω for some finite dimensional matrix D of size dim(Ω)×dim(Θ). Recalling the expression of the total gradient, we can deduce the equality between parametric and functional gradients:∇ ω G tot (ω) = g ω -B ω J ⊤ ω J ω C ω J ⊤ ω -1 J ω d ω = g ω -DJ ω C ω J ⊤ ω J ω C ω J ⊤ ω -1 J ω d ω = g ω -DJ ω d ω = g ω -DJ ω C ω C -1 ω d ω = g ω + B ω a ⋆ ω = ∇F(ω).The first equality follows from the general expression of the total gradient ∇ ω G tot (ω). In the second line we use the expression of B ω which then allows to simplify the expression in the third line. Then, recalling that the Hessian operator C ω is invertible, we get the fourth line. Finally, the result follows by using again the expression of B ω and recalling the definition of the adjoint function a ⋆ ω .The optimization of the prediction function ĥω in the inner-level optimization loop is similar to AID, although the total gradient computation differs significantly. Unlike AID, Algorithm 1 does not require differentiating through the parameters of the prediction model when estimating the total gradient ∇F(ω). This property results in an improved cost in time and memory in most practical cases as shown in Table 1 and Figure 4 . More precisely, AID requires computing Hessian-vector products of size p in , which corresponds to the number of hidden layer weights of the neural network ĥω . While FuncID only requires Hessian-vector products of size d v , i.e. the output dimension of ĥω . In many practical cases, the network's parameter dimension p in is much larger than its output size d v , which results in considerable benefits in terms of memory when using FuncID rather than AID, as shown in Figure 4 (left). Furthermore, unlike AID, the overhead of evaluating Hessian-vector products in FuncID is not affected by the time cost for evaluating the prediction network. When ĥω is a deep network, such an overhead increases significantly with the network size, making AID significantly slower (Figure 4 (right)).γ(T Lin + T h ) βp in + M h FuncID γT Lin + (2 + δ)T a + T h βd v + M aTable 1 : Cost in time and memory for performing a single total gradient estimation using either AID or FuncID and assuming the prediction model is learned. Time cost: T h and T a represent the time cost of evaluating both prediction and adjoint models h and a, while T in is the time cost for evaluating the inner objective once the outputs of h are computed. The factors γ and δ are multiplicative overheads for evaluating hessian-vector products and gradient. Memory cost: M h and M a represent the memory cost of storing the intermediate outputs of h and a, p in and d v are the memory costs of storing the Hessian-vector product for AID and FuncID respectively and β is a multiplicative constant that depends on a particular implementation. The continuous lines are experimental results obtained using a JAX implementation [Bradbury et al., 2018] running on a GPU. The dashed lines correspond to theoretical estimates obtained using the algorithmic costs given in Table 1 with γ = 12, δ = 2 for time, and the constant factors in the memory cost fitted to the data.We closely follow the experimental setting of the state-of-the-art method DFIV [Xu et al., 2021a] . The goal of this experiment is to learn a model f ω approximating the structural function f struct that accurately describes the effect of the treatment t on the outcome o with the help of an instrument x, as illustrated in Figure 5 . We follow the exact same data generation procedure as in Xu et al. [2021a, Appendix E.3 ]. From the dsprites dataset [Matthey et al., 2017] , we generate the treatment t and outcome o as follows:1. Uniformly sample latent parameters scale, rotation, posX, posY from dsprites. Here, function Fig returns the corresponding image to the latent parameters, and λ, ε are noise variables generated from λ ∼ N (0.0, 0.1I) and ε ∼ N (0.0, 0.5). Each element of the matrix A ∈ R 10×4096 is generated from Unif(0.0, 1.0) and fixed throughout the experiment. From the data generation process, we can see that t and o are confounded by posY. We use the instrumental variable x =(scale, rotation, posX) ∈ R 3 , and figures with random noise as treatment variable t. The variable posY is not revealed to the model, and there is no observable confounder. The structural function for this setting isf struct (t) = ∥At∥ 2 2 -5000 1000 .Test data points are generated from grid points of latent variables. The grid consist of 7 evenly spaced values for posX, posY, 3 evenly spaced values for scale, and 4 evenly spaced values for orientation.All results are reported over an average of 20 runs with different seeds on 24GB NVIDIA RTX A5000 GPUs.Feature maps. As in the DFIV setting, we approximate the true structural function f struct with f ω = u ⊤ ψ χ (t) where ψ χ is a feature map of the treatment t, u is a vector in R d2 , and f ω is parameterized by ω = (u, χ). To solve the inner-problem of the bilevel formulation in Section 4.1, the inner prediction function h ω is optimized over functions of the form h(x) = V ϕ(x) where we denote ϕ the feature map of the instrument x and V is a matrix in R d1×d1 . The feature maps ψ χ and ϕ are neural networks (Table 2 ) that are optimized using empirical objectives from Section 3.1 and synthetic dsprites data, the linear weights V and u are fitted exactly at each iteration.Choice of the adjoint function in FuncID. In the dsprites experiment, we call linear FuncID the functional implicit diff. method with a linear choice of the adjoint function. Linear FuncID uses an adjoint function of the form a ⋆ ω (x) = W ϕ(x) with W ∈ R d1×d1 . In other words, to find a ⋆ ω , the features ϕ are fixed and only the optimal linear weight W is computed in closed-form. In the FuncID method, the adjoint function lives in the same function space as h ω . This is achieved by approximating a ⋆ ω with a separate neural network with the same architecture as h ω . Hyper-parameter tuning. As in the setup of DFIV, for training all methods, we use 100 outer iterations (N in Algorithm 1), and 20 inner iterations (M in Algorithm 1) per outer iteration with full-batch. We select the hyper-parameters based on the best validation loss, which we obtain using a validation set with instances of all three variables (t, o, x) [Xu et al., 2021a, Appendix A] . Because of the number of linear solvers, the grid search performed for AID is very large, so we only run it with one seed. For other methods, we run the grid search on 4 different seeds and take the ones with the highest average validation loss. Additionally, for the hyper-parameters that are not tuned, we take the ones reported in Xu et al. [2021a] .• Deep Feature Instrumental Variable Regression: All DFIV hyper-parameters are set based on the best ones reported in Xu et al. [2021a] .• Approximate Implicit Differentiation: We perform a grid search over 5 linear solvers (two variants of gradient descent, two variants of conjugate gradient and an identity heuristic solver), linear solver learning rate 10 -n with n ∈ {3, 4, 5}, linear solver number of iterations {2, 10, 20}, inner optimizer learning rate 10 -n with n ∈ {2, 3, 4}, inner optimizer weight decay 10 -n with n ∈ {1, 2, 3} and outer optimizer learning rate 10 -n with n ∈ {2, 3, 4}.• Iterative Differentiation: We perform a grid search over number of "unrolled" inner iterations {2, 5} (this is chosen because of memory constraints since "unrolling" an iteration is memory-heavy), number of warm-start inner iterations {18, 15}, inner optimizer learning rate 10 -n with n ∈ {2, 3, 4}, inner optimizer weight decay 10 -n with n ∈ {1, 2, 3} and outer optimizer learning rate 10 -n with n ∈ {2, 3, 4}.• Gradient Penalty: The method is based on Eq. 5.1 in Shen and Chen [2023] , for this singlelevel method we perform a grid search on the learning rate 10 -n with n ∈ {2, 3, 4, 5, 6}, weight decay 10 -n with n ∈ {1, 2, 3}, and the penalty weight 10 -n with n ∈ {0, 1, 2, 3}.Since the method has only a single optimization loop, we increase the number of total iterations to 2000 compared to the other methods (100 outer-iterations and 20 inner iterations).• Value Function Penalty: The method is based on Eq. 3.2 in Shen and Chen [2023] , for this method we perform the same grid search as for the Gradient Penalty method. However, since this method has an inner loop, we perform 100 outer iterations and perform a grid search on the number of inner iterations with 10 and 20.• FuncID: We perform a grid search over the number of iterations for learning the adjoint network {10, 20}, adjoint optimizer learning rate 10 -n with n ∈ {2, 3, 4, 5, 6} and adjoint optimizer weight decay 10 -n with n ∈ {1, 2, 3}. The rest of the parameters are the same as for DFIV since the inner and outer models are almost equivalent to the treatment and instrumental networks used in their experiments.We run an additional experiment with 10k training points using the same setting described above to illustrate the effect of the sample size on the methods. Figure 6 shows that a similar conclusion can be drawn when increasing the training sample size from 5k to 10k, thus illustrating the robustness of the obtained results. The key observation here is that the same batches of data are used for both the inner and outer problems, i.e. B in = B out . Therefore, we only need to evaluate the function a on a finite set of points x where (x, y) ∈ B in . Without restricting the solution set of a or adding regularization to Ladj , the optimal solution a ⋆ ω simply matches -∂ v f ( ĥω (x), y) on the set of points x s.t. (x, y) ∈ B in . Our implementation directly exploits this observation and uses the following expression for the total gradient estimation:g out = - (x,y)∈Bin∂ ω,v f ( ĥω (x), r ω (x), s ω (x))∂ v f ( ĥω (x), y).(64)As in the experiments of Nikishin et al. [2022] , we use the CartPole environment with 2 actions, 4-dimensional continuous state space, and optimal returns of 500. For evaluation, we use a separate copy of the environment. The reported return is an average of 10 runs with different seeds.Networks. We us the same neural network architectures that are used in the CartPole experiment of Nikishin et al. [2022, Appendix D] . All networks have two hidden layers and ReLU activations. Both hidden layers in all networks have dimension 32. In the misspecified setting with the limited model class capacity, we set the hidden layer dimension to 3 for the dynamics and reward networks.Hyper-parameters. We perform 200000 environment steps (outer-level steps) and set the number of inner-level iterations to M = 1 for both OMD and FuncID. for MLE, we perform a single update to the state-value function for each update to the model. For training, we use a replay buffer with a batch size of 256, and set the discount factor γ to 0.99. When sampling actions, we use a temperature parameter α = 0.01 as in Nikishin et al. [2022] . The learning rate for outer parameters ω is set to 10 -3 . For the learning rate of the inner neural network and the moving average coefficient τ , we perform a grid search over 10 -4 , 10 -3 , 3 • 10 -3 and {5 • 10 -3 , 10 -2 } as in Nikishin et al. [2022] .Time comparison. Figure 7 shows the average reward on the evaluation environment as a function of training time in seconds. We observe that our model is the fastest to reach best performance both in the well-specified and misspecified settings. MDP model comparison. Figure 8 shows the average prediction error of different methods during training. The differences in average prediction error between the bilevel approaches (OMD, FuncID) and MLE reflect their distinct optimization objectives and trade-offs. OMD and FuncID focus on maximizing performance in the task environment, while MLE emphasizes accurate representation of all aspects of the environment, which can lead to smaller prediction errors but may not necessarily correlate with superior evaluation performance. We also observe that FuncID has a stable prediction error in both settings meanwhile OMD and MLE exhibit some instability. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Our new optimization method does not have an identifiable risk for misuse. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: Our code includes an open access licence in the README file. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.although a generalized version of such a theorem was recently provided under restrictive assumptions[Arbel and Mairal,2022b]. 
 
 Appendix:  


Decision: This paper develops a new functional point of view for bilevel optimization problems in machine learning that provides a new insight into solving bilevel optimization problems with nonconvex lower-level objectives by leveraging their strong convexity in the functional space. This development is noteworthy because, although neural networks are nonconvex, the loss function in model training can be convex or strongly convex. The authors have developed a complete paper with theory, practical insights, and plenty of computational tests. Bilevel optimization is a tough topic, but the authors have also made the paper notably accessible.

Reviews: review: paper_summary: paper_summary: This paper offers a new functional point of view for bilevel optimization problems in machine learning. This functional approach allows the use of an overparameterized neural network as inner prediction function while previous works have used an inner objective that is strongly convex with respect to the parameters of the prediction function. For the inner problem, the prediction function is a function that lies in a Hilbert space of square integrable functions ($L^2$). The authors develop the theory of Functional Implicit Differentiation, which is a flexible class of algorithms to do functional bilevel optimization over $L^2$ spaces. First, they show that strong convexity assumption with respect to the prediction function as opposed to the model parameters ensures the existence and uniqueness of the optimal prediction function $h_\omega^\ast$ for the inner-level objective. Second, they show that differentiability assumptions on the inner-level objective and its Fr\'echet derivative with respect to prediction function $h$ ensure the differentiability of the map $\omega \rightarrow h_\omega^\ast$. Given further assumption about joint differentiability of outer objective, they show that it is possible to compute total objective $\mathcal{F}$ using the adjoint function that minimizes a quadratic objective over a Hilbert space. Finally, assuming that the inner objective and outer objective are defined over the distributions and that a batch of data samples from each distribution can be sampled, an iterative algorithm is proposed that has the following three steps: (1) approximation of inner objective, outer objective and quadratic objective for adjoint function using the batch samples, (2) do a gradient-based update for the parameters of the prediction function and adjoint function assuming they have a fixed parametric form and (3) total gradient approximation and gradient-based update of the functional parameter $\omega$ obtained by solving the outer objective. They show that these class of algorithms converge to a stationary point at $\mathcal{O}(1/N)$ rate. Experiments are performed using Two-stage least squares regression and Model-based reinforcement learning as use cases.
strength_weakness: strength_weakness: The paper is very technical and requires a good understanding of monotone operator theory and theory of Fr\'echet and Hadamard differentiability to understand it fully. Still, the technical details about the existence and uniqueness of the prediction function and the map $\omega \rightarrow h_\omega^\ast$ are deferred to the appendix for interested readers. I would recommend the authors simplify the notation in the main paper a bit to reduce the clutter. For instance, $h_\omega^\ast(x)$ can be represented simply as $h_\omega$. It may help to have a table that shows the functional arguments (e.g. $h_\omega$, $a_\omega$, etc.), functional parameters $\theta$, $\xi$ and $\omega$ and then function variables $x$ and $y$.
questions: questions: 1. On Line 126, there is a typo: "elemen" should be replaced by "element". 
2. On Line 629, capitalize the word "euclidean".
3. What purpose does the variable $y$ and space $\mathcal{Y}$ serve in general? Is it a variable or a parameter? Your prediction function $h_\omega$ is a function of $x$. In the case of 2SLS, you have $y = (t, o)$ and $x$ is the instrumental variable. The function you are interested in is $f_\omega(t)$ and is a function of $t$, which is a subvariable of $y$ and $h_\omega(x) = \mathbb{E}[f_\omega(t) | x]$. Maybe it will help to clarify what $y$ and $\mathcal{Y}$ are in the theoretical section. 
4. What is the rationale for the label "FuncID linear" in Figure 1? Why do you refer to it as "linear"?
5. Why do you feel "FuncID linear" converge faster in terms of outer iterations than "FuncID" in Figure 1?
limitations: limitations: The authors didn't go into a discussion about the limitations nor did they provide a conclusion section and simply concluded with a results section.
paper_summary: paper_summary: This paper introduces a novel functional perspective on bilevel optimization, where the inner objective is defined over a function space. The authors developed functional implicit differentiation and functional adjoint sensitivity, which together facilitate the establishment of a gradient-based algorithm in the functional space. They also analyze the convergence rate of the proposed algorithm and apply it to two-stage least squares regression and model-based reinforcement learning. Experimental results validate the effectiveness of the proposed method.
strength_weakness: strength_weakness: 1. The convergence analysis in this paper is based on the stochastic biased gradient descent framework, making the results explicitly dependent on the sub-optimality constant $\epsilon_{in}$ and $\epsilon_{adj}$. However, it is unclear how these errors relate to the inner loop $M$ and $K$. It might be beneficial to leverage the strong convexity of the inner and adjusted objective functions to clarify these dependencies. See similar techniques used in [1]-[3]. 

[1] K. Ji, J. Yang, and Y. Liang. Bilevel optimization: Convergence analysis and enhanced design. ICML, 2021. 

[2] T. Chen, Y. Sun, and W. Yin. Closing the Gap: Tighter Analysis of Alternating Stochastic Gradient Methods for Bilevel Problems. NeurIPS, 2021. 

[3] M. Dagréou, P. Ablin, S. Vaiter, T. Moreau. A framework for bilevel optimization that enables stochastic and global variance reduction algorithms. NeurIPS, 2022. 


2. Minor issues: experimental baselines. It might be better to also compare with those bilevel methods that are capable to solve nonconvex lower-level problem [4]-[5] as they can also potentially solve the two-stage least squares regression and model-based reinforcement learning problem. 

[4] J. Kwon, D. Kwon, S. Wright, and R. D. Nowak. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. ICLR 2024. 

[5] H. Shen, and T. Chen. On Penalty-based Bilevel Gradient Descent Method. ICML 2023.
questions: questions: Same as weakness.
limitations: limitations: Same as weakness.
paper_summary: paper_summary: The paper proposes a functional approach to bilevel optimization for machine learning, focusing on inner-level problems defined over function spaces rather than traditional parametric settings. This allows the application of the proposed method to machine learning tasks without requiring the strong convexity typically assumed in bilevel optimization.
strength_weakness: strength_weakness: - The paper does not sufficiently compare the proposed methods against a broad spectrum of existing algorithms, particularly the latest advancements in the field of bilevel optimization problem. This lack of comprehensive benchmarking restricts the ability to fully evaluate the performance enhancements or potential drawbacks of the proposed methods relative to the state-of-the-art.

- The paper lacks a thorough analysis of how the proposed methods perform across varied settings and parameter configurations. It is not clear if there are specific scenarios where the methods might underperform or fail to converge. Additional details on the robustness of the methods in diverse operational environments would be beneficial.

- The assumptions necessary for the theoretical framework might not be easily verifiable in practical scenarios, potentially limiting the method's applicability.
questions: questions: 1. Can the proposed functional framework be generalized to other types of bilevel problems that involve constraints on the lower-level problem? How would you address potential non-smoothness or non-convexity in these extended settings?

2. Please explain how does the choice of function space impact the stability and convergence of the FuncID method?
limitations: limitations: Authors addressed their work limitation.
paper_summary: paper_summary: The authors propose a functional view of bilevel optimization in machine learning, in which the inner objective is often strongly convex for many problems of interest. The authors prove the resulting optimization algorithm's convergence and benchmark their approach for regression and model-based RL tasks.
strength_weakness: strength_weakness: 1. The authors seem to do a "bait-and-switch" when going from function space to parameter space. Namely, the authors repeatedly emphasize that in function space, the inner objective is generally strongly convex. However, as the authors point out, this is no longer true when moving to function parameterizations -- and all the experiments concern concrete parameterizations. I'm fairly lost here as to how the authors theoretically handle this jump, as all the theory seems to assume that we are operating in function space.
2. The assumptions in Theorem 3.1 seem very strong. Namely, it is assumed that the inner and adjoint problems are solved close to optimality; however, the inner optimization is nonconvex in the parameters (as discussed above).
3. The authors repeatedly scatter in comparisons to AID throughout the text. I think it would be good to have these all summarized in one table in the appendix (runtime, convergence guarantees, assumptions, etc.)
4. The experimental results are quite weak. Namely, buried in H.3 the authors experiment with boosting the training sample size by a factor of two, and in this setting FuncID seems to underperform DFIV (although FuncID linear is now better). It seems that the improvement over the state of the art is not very robust.
5. FuncID seems to require more outer iterations to converge than baselines (Figure 1b), which I don't see being discussed by the authors.
questions: questions: 1. Can you elaborate on $\epsilon_{in}$ and $\epsilon_{adjoint}$ in Theorem 3.1? I'm not sure what the expectation is taken over in the definition of line 1100.
2. Why is the prediction error of FuncID worse than MLE in Figure 7b?

Notes:
1. In (FBO), there's a space between the colon and the equals (use \coloneqq).
2. In line 94, link to where the assumption is discussed.
3. Font size in the figures is very small and hard to read.
limitations: limitations: I'd like to see an explicit limitations section. The authors' response to the limitations question in the checklist is that they state theorem assumptions -- this is not a comprehensive discussion of limitations, which should include practical experimental considerations as well.

Comments: title: Discussion period authors' comment
Reviewer response?
Discussion period
comment: Thank you to all reviewers for responding to our rebuttal, and to the area chair for encouraging the discussion on our submission. We will ensure that the changes mentioned in the rebuttal and confirmed by the reviewers are integrated into the final version of the paper, if it is accepted.
I thank the authors for their detailed responses. It solves all of my concerns so that I will raise my score.
Reviewer MXrv, could you please review the authors' response and see whether it addresses your questions? Please acknowledge having done so in a comment. Thanks.
Thank you for your responses to my questions and for the clarifications! Good to hear that the typos will be fixed in the final draft and that the missing details will be added to the final draft! I am also glad that you are adding a detailed discussion/conclusion at the end as it is important to leave the reader with some take-away points from your paper. Good luck!
Appreciate authors for their complete response and clarification. I am satisfied with their response and maintain my score.
Thank you to the authors for their clarification. The comparison in points A and B above is really important, and it gets lost in the technical details of the paper. I strongly suggest featuring this prominently in the main paper body, and perhaps cutting from section 3.

I've raised my score.
Dear Reviewers,

Thanks for the contributions in the review period and looking forward to a fruitful discussion on Submission 6898 in the next few days.

Kind regards,

Your AC