# Learning Theory for Kernel Bilevel Optimization

Fares El Khoury, Edouard Pauwels, Samuel Vaiter, Michael Arbel

► To cite this version:

Fares El Khoury, Edouard Pauwels, Samuel Vaiter, Michael Arbel. Learning Theory for Kernel Bilevel Optimization. NeurIPS 2025 - 39th Annual Conference on Neural Information Processing Systems, Dec 2025, San Diego, United States. pp.1-47. hal-04950585v3

HAL Id: hal-04950585

https://inria.hal.science/hal-04950585v3

Submitted on 15 Nov 2025

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

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

<figure><img src="image_1.png" /></figure>
Learning Theory for Kernel Bilevel Optimization

**Fares El Khoury***

Université Grenoble Alpes, Inria,
CNRS, Grenoble INP, LJK,
38000 Grenoble, France

**Edouard Pauwels**

Toulouse School of Economics,
Université Toulouse Capitole,
31080 Toulouse, France

**Samuel Vaiter**

CNRS & Université Côte d'Azur,
Laboratoire J. A. Dieudonné,
06108 Nice, France

**Michael Arbel**

Université Grenoble Alpes, Inria,
CNRS, Grenoble INP, LJK,
38000 Grenoble, France

Abstract

Bilevel optimization has emerged as a technique for addressing a wide range of machine learning problems that involve an outer objective implicitly determined by the minimizer of an inner problem. While prior works have primarily focused on the *parametric* setting, a *learning-theoretic* foundation for bilevel optimization in the *nonparametric* case remains relatively unexplored. In this paper, we take a first step toward bridging this gap by studying *Kernel Bilevel Optimization* (KBO), where the inner objective is optimized over a reproducing kernel Hilbert space. This setting enables rich function approximation while providing a foundation for rigorous theoretical analysis. In this context, we derive novel *finite-sample generalization bounds* for KBO, leveraging tools from empirical process theory. These bounds further allow us to assess the statistical accuracy of gradient-based methods applied to the empirical discretization of KBO. We numerically illustrate our theoretical findings on a synthetic instrumental variable regression task.

# 1 Introduction

Bilevel optimization involves a nested structure where one optimization problem, called *outer-level*, is constrained by the solution of another one, called *inner-level* [19]. This formulation has found applications in a broad spectrum of machine learning fields, including hyperparameter tuning [43, 13, 28], meta-learning [14, 59], inverse problems [34], and reinforcement learning [35, 48], making it a powerful tool in theoretical and practical contexts. Its widespread use naturally raises fundamental questions about the *generalization properties* of models learned through this procedure as the number of data samples increases. Several existing works have studied the generalization and convergence of bilevel algorithms under the assumption that the inner-level problem is *strongly convex* and that its parameters lie in a *finite-dimensional* space. These include analyses of the convergence of stochastic bilevel optimization algorithms [3, 24, 29, 38] and approaches based on algorithmic stability [7, 78]. The strong convexity assumption ensures a *unique* inner-level solution, which is crucial for stability and convergence analysis in bilevel optimization. Moreover, restricting the inner-level parameters to a finite-dimensional space does not cover possibly *richer infinite-dimensional* spaces, as in kernel methods, where the parameter's dimension may grow with the sample size. This sample size dependence in *nonparametric* methods poses additional challenges as solutions at different sample sizes are not directly comparable. In contrast, in the finite-dimensional setting, generalization bounds can be derived by quantifying the convergence of inner-level finite-sample solutions toward the infinite samples *population solution* limit, within a fixed Euclidean space.

<footnote>*Correspondence to: fares.el-khoury@inria.fr.</footnote>
Albeit theoretically convenient, strong convexity and finite-dimensionality limit models expressive-ness, restricting them to linear functions. Moving beyond linear models requires either relaxing the strong convexity assumption to accommodate more expressive models, such as deep neural networks [30], or considering nonparametric bilevel problems, where the inner-level variable lies in an expres- sive infinite-dimensional function space, such as a Reproducing Kernel Hilbert Space (RKHS) [63]. Early works in bilevel optimization for machine learning followed the latter approach, developing methods for hyperparameter selection in kernel-based models [39, 41]. These works leverage the representer theorem [64] to transform the infinite-dimensional problem into a finite-dimensional one, with dimension depending on the sample size. However, they do not address how sample size impacts generalization. Another line of research instead focuses on relaxing the strong convexity assumption, proposing bilevel algorithms in the absence of convexity [4, 42, 66]. Yet, non-convex bilevel optimization is a very hard problem in general [47, 4, 17], and strong generalization guarantees in this setting remain out of reach due to instabilities of the inner-level solutions. Overall, learning theory for bilevel problems beyond the strongly convex parametric setting is essentially lacking.

In the present work, we take an initial step toward developing a learning theory that goes beyond the finite-dimensional setting. Specifically, we propose to study Kernel Bilevel Optimization (KBO) problems, where the inner objective $L_{in} : \mathbb{R}^d \times \mathcal{H} \to \mathbb{R}$ finds an optimal inner solution $h_{\omega}^*$ in an RKHS $\mathcal{H}$ for a given parameter $\omega$ in $\mathbb{R}^d$, while the outer objective $L_{out} : \mathbb{R}^d \times \mathcal{H} \to \mathbb{R}$ optimizes the parameter $\omega$ over a closed subset $C$ of $\mathbb{R}^d$, given the inner solution $h_{\omega}^*$:

$$ \min_{\omega \in C} \mathcal{F}(\omega) := L_{out}(\omega, h_{\omega}^*) \quad \text{s.t.} \quad h_{\omega}^* = \arg \min_{h \in \mathcal{H}} L_{in}(\omega, h). \qquad (\text{KBO}) $$

In particular, we focus on objectives that are expectations of point-wise losses, a common setting in learning theory. RKHS provides a natural framework to study learning-theoretic arguments, and has been instrumental for many fruitful results in pattern recognition and machine learning. They allow to describe very expressive non-linear models with simple and stable algorithms, while enabling a rich statistical analysis and featuring adaptivity to the regularity of the population problem [65, 63, 33]. Our choice is also motivated by the relevance of kernel methods, even in the deep learning era. They remain competitive for some prediction problems, such as those involving physics [26, 44]. Additionally, the mathematics of kernel methods are useful to describe the limiting behavior of deep network training for very large models [37, 11]. In this limit, the problem becomes (strongly) convex in an infinite-dimensional function space, simplifying the difficulties of non-convex model parameterizations, a major bottleneck in the analysis of such models. This point of view was leveraged by Petrulionytė et al. [58] who introduced functional bilevel optimization, and our setting can be seen as a special case for which the underlying function space is an RKHS. From a practical perspective, our setting is amenable to first-order methods using implicit differentiation techniques [31, 6, 15].

**Contributions.** We leverage empirical process theory and its extension to U-processes [67] to derive uniform generalization bounds for the value function of (KBO), quantifying the discrepancy between $\mathcal{F}$ and its plug-in estimator $\hat{\mathcal{F}}$ in terms of both their values and gradients. Classical empirical process results [72] do not directly apply here, as our functional setting involves processes taking values in an infinite-dimensional space rather than real numbers. To address this, we exploit the RKHS structure to represent the resulting error terms as real-valued U-processes, enabling the use of results from Sherman [67]. The control in terms of gradients is crucial to study first-order optimization methods, since $\hat{\mathcal{F}}$ is typically non-convex and iterative methods seek to find approximate critical points where $\|\nabla \hat{\mathcal{F}}\|$ is small. Our result relies on an equivalence we establish between $\nabla \hat{\mathcal{F}}$ and a plug-in statistical estimate of $\nabla \mathcal{F}$ that is more amenable to a statistical analysis. We then use our uniform bounds to provide generalization guarantees for gradient descent and projected gradient descent applied to $\hat{\mathcal{F}}$. Under specific assumptions, we show convergence rates for sub-optimality measures that depend on the sample sizes and the number of algorithmic iterations. This illustrates the practical relevance of our generalization bounds on simple bilevel algorithms. For a large number of steps, gradient algorithms applied to the empirical (KBO) find approximate critical points of the population (KBO) up to a statistical error which we control.

**Organization of the paper.** In Section 2, we describe (KBO), give two application examples, and explain implicit differentiation in an RKHS. In Section 3, we present the empirical (KBO) and state our first main result on the gradient of its value function. Section 4 provides uniform generalization bounds for (KBO), with applications to bilevel gradient methods, as well as a sketch of the proof of our main result. Finally, in Section 5, we illustrate our theoretical findings with experiments on synthetic data for the instrumental variable regression problem.
# 2 Kernel bilevel optimization

## 2.1 Problem formulation

We consider the (KBO) problem with an RKHS $\mathcal{H}$, which is a space of real-valued functions defined on a Borel input space $\mathcal{X} \subset \mathbb{R}^p$ and associated with a reproducing kernel $K: \mathcal{X} \times \mathcal{X} \to \mathbb{R}$. We are interested, in particular, in (regularized) objectives expressed as expectations of point-wise loss functions, a formulation widely adopted in machine learning as it allows the loss functions to represent the average performance over some data distribution. Specifically, given two probability distributions $\mathbb{P}$ and $\mathbb{Q}$ supported on $\mathcal{X} \times \mathcal{Y}$ for some target space $\mathcal{Y} \subset \mathbb{R}^q$, we consider objectives of the form:

$$L_{out}(\omega, h) = \mathbb{E}_{\mathbb{Q}}[\ell_{out}(\omega, h(x), y)], \quad L_{in}(\omega, h) = \mathbb{E}_{\mathbb{P}}[\ell_{in}(\omega, h(x), y)] + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^2,$$

where $\ell_{in}$ and $\ell_{out}: \mathbb{R}^d \times \mathbb{R} \times \mathbb{R}^q \to \mathbb{R}$ represent the inner and outer point-wise loss functions, $\lambda > 0$ is the regularization parameter which is fixed through this work, and $\| \cdot \|_{\mathcal{H}}$ denotes the norm in the RKHS $\mathcal{H}$. The regularization term in $L_{in}$ is often used in practice to prevent overfitting by penalizing overly complex models. In our setting, it ensures strong convexity of $h \mapsto L_{in}(\omega, h)$ under mild assumptions on $\ell_{in}$, which will be critical to leverage functional implicit differentiation.

**Assumptions.** Through the paper, we make the following five assumptions to derive generalization bounds while retaining a simple and modular presentation.

(A) (Measurability of $K$). $K$ is measurable on $\mathcal{X} \times \mathcal{X}$.

(B) (Boundedness of $K$). There exists a constant $\kappa > 0$ such that $K(x, x) \le \kappa$, for any $x \in \mathcal{X}$.

(C) (Compactness of $\mathcal{Y}$). The subset $\mathcal{Y}$ of $\mathbb{R}^q$ is compact.

(D) (Regularity of $\ell_{in}$ and $\ell_{out}$). The functions $\ell_{in}$ and $\ell_{out}$ are of class $C^3$ jointly in their first two arguments $(\omega, v)$, and their derivatives are jointly continuous in $(\omega, v, y)$.

(E) (Convexity of $\ell_{in}$). For any $(\omega, y) \in \mathbb{R}^d \times \mathbb{R}^q$, the map $v \mapsto \ell_{in}(\omega, v, y)$ is convex.

**Assumptions (A) and (B) on $K$** hold for a wide class of kernels, such as the Gaussian, Laplacian, and Matérn [61] kernels. They also hold if $K$ is a Mercer kernel [50], i.e., a continuous, positive-definite kernel on a compact domain $\mathcal{X}$. Specifically, a kernel built using neural network features, e.g., neural tangent kernel [37], satisfies these assumptions on the space of images, which is compact since the pixel values have a bounded range. Assumption (C) on $\mathcal{Y}$ is a mild assumption that holds in most supervised learning applications, such as classification where $\mathcal{Y}$ is finite, or cases where $\mathcal{Y} = [0, 1]^q$, enabling the representation of complex data, like images. Assumption (D) on the point-wise objectives is a mild regularity assumption, which is met for the most commonly used loss functions in practice, including the squared loss, logistic loss, cross-entropy loss, and KL divergence. Finally, Assumption (E) is essential to ensure the existence and uniqueness of a smooth minimizer $h^*$. It is a relatively weak assumption that was recently considered in [58] in the context of functional bilevel optimization, and that holds in many cases of interest, as discussed in Section 2.2.

**Remark 2.1.** Assumptions (B) to (D) can be relaxed at the expense of weaker yet more technical assumptions, such as finite moment assumptions on $\mathbb{P}$ and $\mathbb{Q}$, and suitable polynomial growth of the kernel and some partial derivatives of $\ell_{in}$ and $\ell_{out}$. It is also sufficient to require that Assumption (D) holds on $\mathcal{U} \times \mathbb{R} \times \mathbb{R}^q$ where $\mathcal{U}$ is an open neighborhood of $\mathcal{C}$, and that Assumption (E) holds for any $\omega \in \mathcal{C}$ and $y \in \mathcal{Y}$. We prefer to keep these stronger yet simpler assumptions for clarity.

## 2.2 Examples of (KBO) in machine learning

To illustrate the relevance of (KBO), we consider two examples that highlight its applicability.

**Hyperparameter selection under distribution shift.** In this application, the aim is to select the best hyperparameters for a machine learning model, e.g., regularization parameters, while accounting for distribution shift between the training and test data, i.e., when the training and test data distributions $\mathcal{D}_{train}$ and $\mathcal{D}_{test}$ are different [57, 28]. This can be viewed as an instance of (KBO) when using models in an RKHS. At the inner-level, the model $h$ is trained to minimize the regularized training squared error loss, with the hyperparameter $\omega > 0$ representing the weight for the data fitting term. At the outer-level, the task is to select the hyperparameter $\omega$ that maximizes the model's performance on the
distribution-shifted test data. Both inner and outer objectives can thus be formulated as:

$$L_{out}(\omega, h) = \frac{1}{2} \mathbb{E}_{(x,y) \sim \mathcal{D}_{\text{test}}} [|h(x) - y|^2], \quad L_{in}(\omega, h) = \frac{\omega}{2} \mathbb{E}_{(x,y) \sim \mathcal{D}_{\text{train}}} [|h(x) - y|^2] + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^2.$$

This formulation could be used for domain adaptation [12] or domain generalization [73] to choose hyperparameters that perform well on the distribution-shifted test data.

**Instrumental variable regression.** It is a technique used to address endogeneity in statistical modeling by leveraging instruments to estimate causal relationships [53]. The goal is to estimate a function $t \mapsto f_\omega(t)$ parameterized by a vector $\omega$, that satisfies $y = f_\omega(t) + \epsilon$, where $y \in \mathbb{R}$ is the observed outcome, $t$ is the treatment, and $\epsilon$ is the error term. The key issue is that $t$ is endogenous, which means that it is correlated with $\epsilon$, making direct regression inconsistent. Indeed, such correlation leads to biased estimates of $f_\omega(t)$ as the assumption of exogeneity, i.e., independence of $t$ and $\epsilon$, is violated. To resolve this, one can use an instrumental variable $x$, uncorrelated with $\epsilon$ but correlated with $t$, to recover the relationship between $y$ and $t$ without being directly affected by the bias introduced by $\epsilon$, typically via two-stage least squares regression [68, 51]. As shown in [58], this approach can be naturally expressed as a bilevel problem with inner and outer objectives of the form:

$$L_{out}(\omega, h) = \frac{1}{2} \mathbb{E}_{x,y} [ |h(x) - y|^2 ], \quad L_{in}(\omega, h) = \frac{1}{2} \mathbb{E}_{x,t} [ |h(x) - f_\omega(t)|^2 ] + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^2,$$

where $h$ can be chosen to be in an RKHS to allow flexibility in the estimation while retaining uniqueness of the solution $h_\omega^*$, a key property in bilevel optimization.

## 2.3 Implicit differentiation in an RKHS

A stationarity measure in (KBO) is the gradient $\nabla \mathcal{F}(\omega)$ of the value function $\mathcal{F}$. Computing $\nabla \mathcal{F}(\omega)$, however, is challenging and will be addressed in this section. At a high level, our approach proceeds in two steps. First, we derive an abstract, a priori intractable, expression for $\nabla \mathcal{F}(\omega)$ using implicit differentiation in an RKHS, which is the main source of difficulty. Then, we leverage the structure of our problem to reformulate the gradient in Proposition 2.2 using the solution of a regression problem in the RKHS (the adjoint problem). This more concrete formulation can be approximated with finite samples and will later serve as the foundation of our statistical analysis. Formally, evaluating the gradient requires computing the Jacobian $\partial_\omega h_\omega^*$, which can be viewed as a linear operator from $\mathcal{H}$ to $\mathbb{R}^d$. Indeed, $h_\omega^*$ depends implicitly on $\omega$. A key ingredient for computing $\partial_\omega h_\omega^*$ is the implicit function theorem [36], which guarantees the differentiability of the implicit function $\omega \mapsto h_\omega^*$ and allows characterizing $\partial_\omega h_\omega^*$ as the unique solution of a linear system of the form:

$$\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^*) + \partial_{\omega} h_{\omega}^* \partial_h^2 L_{in}(\omega, h_{\omega}^*) = 0, \quad (1)$$

where $\partial_h^2 L_{in}(\omega, h_\omega^*)$ is an operator from $\mathcal{H}$ to itself representing the Hessian of $L_{in}$ w.r.t. $h$, while $\partial_{\omega,h}^2 L_{in}(\omega, h_\omega^*)$ is an operator from $\mathcal{H}$ to $\mathbb{R}^d$ representing the cross derivatives of $L_{in}$ w.r.t. to $\omega$ and $h$. Applying such result requires $h \mapsto L_{in}(\omega, h)$ to be Fréchet differentiable with invertible Hessian operator and jointly Fréchet differentiable gradient map $(\omega, h) \mapsto \partial_h L_{in}(\omega, h)$. All these properties are satisfied in our setting under Assumptions (A) to (E) as shown in Propositions B.1 to B.3 of Appendix B.1. Furthermore, when $L_{out}$ is Fréchet differentiable, which is our case under Assumptions (A) to (D) (Proposition B.1 of Appendix B.1), then by composition with $\omega \mapsto (\omega, h_\omega^*)$, the map $\omega \mapsto \mathcal{F}(\omega)$ must also be differentiable with gradient obtained using the chain rule:

$$\nabla \mathcal{F}(\omega) = \partial_{\omega} L_{out}(\omega, h_{\omega}^*) + \partial_{\omega} h_{\omega}^* \partial_{h} L_{out}(\omega, h_{\omega}^*).$$

The above expression for the gradient is intractable as it involves abstract operators, namely the derivatives $\partial_h$, $\partial_h^2$, and $\partial_{\omega,h}^2$, the last two of which arise when replacing $\partial_\omega h_\omega^*$ by its expression in Equation (1). In Proposition 2.2 below, we derive an explicit expression for $\nabla \mathcal{F}(\omega)$ which exploits the particular structure of the objectives $L_{in}$ and $L_{out}$ as expectations of point-wise losses.

**Proposition 2.2 (Expression of the total gradient).** Under Assumptions (A) to (E), $\mathcal{F}$ is differentiable on $\mathbb{R}^d$, with gradient $\nabla \mathcal{F}(\omega)$, for any $\omega \in \mathbb{R}^d$, given by:

$$\nabla \mathcal{F}(\omega) = \mathbb{E}_{\mathbb{Q}} [\partial_{\omega} \ell_{out}(\omega, h_{\omega}^*(x), y)] + \mathbb{E}_{\mathbb{P}} [\partial_{\omega,v}^2 \ell_{in}(\omega, h_{\omega}^*(x), y) a_{\omega}^*(x)], \quad (2)$$

where the adjoint function $a_\omega^* \in \mathcal{H}$ is the unique minimizer of a strongly convex quadratic objective $a \mapsto L_{adj}(\omega, a)$ defined on $\mathcal{H}$ as:

$$L_{adj}(\omega, a) := \frac{1}{2} \mathbb{E}_{\mathbb{P}} [\partial_v^2 \ell_{in}(\omega, h_\omega^*(x), y) a^2(x)] + \mathbb{E}_{\mathbb{Q}} [\partial_v \ell_{out}(\omega, h_\omega^*(x), y) a(x)] + \frac{\lambda}{2} \|a\|_{\mathcal{H}}^2, \quad (3)$$
where $\partial_{\omega}l_{out}$ and $\partial_{v}l_{out}$ are the first-order partial derivatives of $l_{out}$ w.r.t. $\omega$ and $v$, while $\partial_{\omega,v}^2l_{in}$ and $\partial_{v}^2l_{in}$ denote the second-order partial derivatives of $l_{in}$ w.r.t. $\omega$ and $v$.

Proposition 2.2 is proved in Appendix B and relies essentially on proving Bochner's integrability [25, Definition 1, Chapter 2] of some suitable operators on $\mathcal{H}$, and then applying Lebesgue's dominated convergence theorem for Bochner's integral [25, Theorem 3, Chapter 2] to interchange derivatives and expectations. The expression in Proposition 2.2 provides a natural way for approximating $\nabla \mathcal{F}(\omega)$ by estimating all expectations using finite-sample averages, as we further discuss in Section 3.

## 3 Finite-sample approximation of (KBO)

In this section, we consider an approximation of (KBO) when only a finite number of i.i.d. samples $(x_i, y_i)_{1\le i\le n}$ and $(\tilde{x}_j, \tilde{y}_j)_{1\le j\le m}$ from $\mathbb{P}$ and $\mathbb{Q}$ are available. This setting is ubiquitous in machine learning as it allows finding tractable approximate solutions to the original problem. As we are interested in approximately solving (KBO) using gradient methods, our focus here is to derive estimators for both the value function $\mathcal{F}(\omega)$ and its gradient $\nabla \mathcal{F}(\omega)$, whose generalization properties will be studied in Section 4.

In Section 3.1, we follow a commonly used approach of first deriving a plug-in estimator $\hat{\mathcal{F}}$ of the value function, then considering its gradient $\nabla \hat{\mathcal{F}}(\omega)$ as an approximation to $\nabla \mathcal{F}(\omega)$. In Section 3.2, we show that this approximation is equivalent to a second estimator, more amenable to a statistical analysis, obtained by directly computing a plug-in estimator of $\nabla \mathcal{F}$ based on its expression in Equation (2). Figure 1 summarizes such equivalence.

<figure><img src="image_2.png" /><figcaption>Figure 1: A commutative diagram illustrating that plug-in statistical estimation and differentiation can be interchanged for $\mathcal{F}$ and $\hat{\mathcal{F}}$ resulting in a single gradient estimator.</figcaption></figure>

### 3.1 Value function: plug-in estimator and its gradient

A natural approach for finding approximate solutions to (KBO) is to consider an approximate problem obtained after replacing the objectives $L_{in}$ and $L_{out}$ by their empirical approximations $\hat{L}_{in}$ and $\hat{L}_{out}$:

$$ \hat{L}_{out}(\omega, h) = \frac{1}{m} \sum_{j=1}^{m} \ell_{out}(\omega, h(\tilde{x}_j), \tilde{y}_j), \quad \hat{L}_{in}(\omega, h) = \frac{1}{n} \sum_{i=1}^{n} \ell_{in}(\omega, h(x_i), y_i) + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^2. $$

A plug-in estimator $\omega \mapsto \hat{\mathcal{F}}(\omega)$ is then obtained by first finding a solution $\hat{h}_{\omega}$ minimizing $h \mapsto \hat{L}_{in}(\omega, h)$, that is meant to approximate the optimal inner solution $h_{\omega}^{*}$, and subsequently plugging it into $\hat{L}_{out}$. This procedure results in the following empirical version of (KBO):

$$ \min_{\omega \in \mathcal{C}} \hat{\mathcal{F}}(\omega) := \hat{L}_{out}(\omega, \hat{h}_{\omega}) \quad \text{s.t.} \quad \hat{h}_{\omega} = \arg \min_{h \in \mathcal{H}} \hat{L}_{in}(\omega, h). $$

The inner problem still requires optimizing over a, potentially infinite-dimensional, RKHS. However, its finite-sum structure allows equivalently expressing it as a finite-dimensional bilevel optimization, by application of the so-called representer theorem [64]:

$$ \begin{array}{ll} \min_{\omega \in \mathcal{C}} \hat{\mathcal{F}}(\omega) := \frac{1}{m} \sum_{j=1}^{m} \ell_{out}(\omega, (\bar{\mathbf{K}} \hat{\gamma}_{\omega})_j, \tilde{y}_j) & (\text{KBO}) \\ \text{s.t.} \quad \hat{\gamma}_{\omega} = \arg \min_{\gamma \in \mathbb{R}^n} \frac{1}{n} \sum_{i=1}^{n} \ell_{in}(\omega, (\mathbf{K} \gamma)_i, y_i) + \frac{\lambda}{2} \gamma^{\top} \mathbf{K} \gamma. & \end{array} $$

Here, $\mathbf{K} \in \mathbb{R}^{n \times n}$ and $\bar{\mathbf{K}} \in \mathbb{R}^{m \times n}$ are matrices containing the pairwise kernel similarities between the data points, i.e., $\mathbf{K}_{ij} := K(x_i, x_j)$ and $\bar{\mathbf{K}}_{ij} := K(\tilde{x}_i, \tilde{y}_j)$, while $\gamma$ is a parameter vector in $\mathbb{R}^n$ representing the inner-level variables. The optimal solution $\hat{\gamma}_{\omega}$ enables recovering the prediction function $\hat{h}_{\omega}$ by linearly combining kernel evaluations at inner-level samples, i.e., $\hat{h}_{\omega} = \sum_{i=1}^{n} (\hat{\gamma}_{\omega})_i K(x_i, \cdot)$. The formulation in (KBO) enables deriving an expression for the gradient $\nabla \hat{\mathcal{F}}(\omega)$ in terms of
the Jacobian $\partial_{\omega}\hat{\gamma}_{\omega}$ by direct application of the chain rule. Unlike $\partial_{\omega}h_{\omega}^{*}$, which requires solving the infinite-dimensional linear system in Equation (1), $\partial_{\omega}\hat{\gamma}_{\omega}$ can be obtained by solving a finite-dimensional linear system using the implicit function theorem (see Proposition C.1 of Appendix C). Hence, (KBØ) falls into a class of optimization problems for which a rich body of literature has proposed practical and scalable algorithms, leveraging the expression of $\nabla\hat{\mathcal{F}}(\omega)$ [38, 3, 24]. Solving (KBØ) thus provides a practical way to approximate the solution of the original population problem (KBO), as proposed by several prior works on bilevel optimization involving kernel methods [39, 41].

**Non-applicability of existing results.** Despite its practical advantages, the above approach yields algorithms that are not directly amenable to a statistical analysis. The key challenge is to be able to control the approximation error between the true gradient $\nabla\mathcal{F}(\omega)$ and its approximation $\nabla\hat{\mathcal{F}}(\omega)$ as the sample sizes $n$ and $m$ increase. Existing statistical analyses for bilevel optimization, such as [7, 78], consider objectives that are expectations or finite sums of point-wise losses, as we do here. While these results can be applied to our setting for each fixed $n$, they do not capture the generalization behavior as $n$ grows. In particular, they require both the inner- and outer-level parameters to lie in spaces of fixed dimensions, that are independent of $n$ and $m$. That is because these parameters are expected to converge to some fixed vectors as $n, m \to +\infty$. In contrast, in our setting, the inner-level parameter $\gamma$ lies in $\mathbb{R}^n$, so its dimension grows with $n$ and is not expected to converge to any well-defined object. Existing non-kernel generalization bounds are discussed in Appendix A.

**Relation to instrumental variable regression.** Our formulation is related to instrumental variable regression, which is a special case, but differs in that we study regularized inner problems with a fixed $\lambda$, independent of the sample sizes. In contrast, the instrumental variable regression literature typically considers un-regularized population problems ($\lambda = 0$), for which a closed-form expression of the inner minimizer is available [32, 68, 76, 45]. Our contribution lies in an orthogonal direction: we handle more general inner objectives, for which even the analysis of regularized problems raises new obstacles that had not been tackled before. Moreover, some prior works have provided convergence rates in the instrumental variable regression setting [1, 2, 22, 23]. Yet, these studies focus on asymptotic results with sieve estimators, meanwhile we leverage the RKHS structure to provide finite-sample bounds. More recently, Meunier et al. [51] established minimax optimal rates under source assumptions by exploiting bounds for vector-valued kernel ridge regression [46] via a spectral filtering technique [52]. However, this approach is not applicable to our case, as it requires the losses to be quadratic in $\omega$, as further discussed in Appendix A.

Next, we provide an equivalent expression for $\nabla\hat{\mathcal{F}}(\omega)$, essential for our analysis in Section 4.

## 3.2 Plug-in estimator of the total gradient

We now consider an a priori different approach for approximating the total gradient $\nabla\mathcal{F}(\omega)$ based on direct plug-in estimation from Equation (2), and show that it recovers the previously introduced estimator $\nabla\hat{\mathcal{F}}(\omega)$. Such approach consists in replacing all expectations in Equation (2) by empirical averages, then replacing $h_{\omega}^{*}$ and $a_{\omega}^{*}$ by their finite-sample estimates $\hat{h}_{\omega}$ and $\hat{a}_{\omega}$. This yields the following estimator of the total gradient:

$$ \widehat{\nabla\mathcal{F}}(\omega) = \frac{1}{m} \sum_{j=1}^{m} \partial_{\omega} \ell_{out}(\omega, \hat{h}_{\omega}(\tilde{x}_{j}), \tilde{y}_{j}) + \frac{1}{n} \sum_{i=1}^{n} \partial_{\omega,v}^{2} \ell_{in}(\omega, \hat{h}_{\omega}(x_{i}), y_{i}) \hat{a}_{\omega}(x_{i}). \quad (4) $$

Just as in Section 3.1, $h_{\omega}^{*}$ can be estimated by $\hat{h}_{\omega}$, the minimizer of the empirical objective $h \mapsto \hat{L}_{in}(\omega, h)$. Similarly, $a_{\omega}^{*}$ can be approximated by $\hat{a}_{\omega}$, the minimizer of $a \mapsto \hat{L}_{adj}(\omega, a)$ defined as:

$$ \hat{L}_{adj}(\omega, a) = \frac{1}{2n} \sum_{i=1}^{n} \partial_{v}^{2} \ell_{in}(\omega, \hat{h}_{\omega}(x_{i}), y_{i}) a^{2}(x_{i}) + \frac{1}{m} \sum_{j=1}^{m} \partial_{v} \ell_{out}(\omega, \hat{h}_{\omega}(\tilde{x}_{j}), \tilde{y}_{j}) a(\tilde{x}_{j}) + \frac{\lambda}{2} \|a\|_{H}^{2}, \quad (5) $$

which serves as the empirical counterpart of the adjoint objective $L_{adj}$ given in Equation (3). Both functions $\hat{h}_{\omega}$ and $\hat{a}_{\omega}$ can be expressed as linear combinations of kernel evaluations with some given parameter vectors whose dimensions increase with the sample size $n$ (see Proposition C.1 for $\hat{h}_{\omega}$ and Lemma C.2 for $\hat{a}_{\omega}$, both in Appendix C). However, these parameters are not required to compute the plug-in estimator $\nabla\hat{\mathcal{F}}(\omega)$ in Equation (4), since only the function values of $\hat{h}_{\omega}$ and $\hat{a}_{\omega}$ are needed.
This property is precisely what makes $\nabla \widehat{\mathcal{F}}(\omega)$ suitable for a statistical analysis. Indeed, its estimation error depends on the approximation errors of $\hat{h}_\omega$ and $\hat{a}_\omega$, which always belong to the same space $\mathcal{H}$ regardless of the sample size, and are expected to approach their population counterparts. This contrasts with $\nabla \widehat{\mathcal{F}}(\omega)$ obtained by implicit differentiation, whose analysis would need controlling the behavior of the vector $\hat{\gamma}_\omega$ that resides in a growing-dimensional space as $n \to +\infty$.

The next proposition establishes a link between practical applications and theoretical analysis by demonstrating that, surprisingly, both estimators $\nabla \widehat{\mathcal{F}}(\omega)$ and $\nabla \widehat{\mathcal{F}}(\omega)$ are precisely equal.

**Proposition 3.1.** Under Assumptions (A) to (E), the gradient $\nabla \widehat{\mathcal{F}}(\omega)$ of the plug-in estimator $\widehat{\mathcal{F}}(\omega)$ of $\mathcal{F}(\omega)$ defined in (KBO) is equal to the plug-in estimator $\nabla \widehat{\mathcal{F}}(\omega)$ of the total gradient $\nabla \mathcal{F}(\omega)$ introduced in Equation (4).

Proposition 3.1 is proved in Appendix C and relies on an application of the representer theorem [64] to provide explicit expressions for both estimators in terms of $\hat{\gamma}_\omega$, kernel matrices $K$ and $\bar{K}$ and partial derivatives of the point-wise objectives $\ell_{in}$ and $\ell_{out}$. Both expressions are then shown to be equal using optimality conditions on the parameters defining $\hat{a}_\omega$. The result in Proposition 3.1 precisely says that the operations of differentiation and plug-in estimation commute in the case of (KBO). Such a commutativity property does not necessarily hold anymore if one considers spaces other than an RKHS, such as $L_2$ [58, Appendix F]. The main difficulty arises from the argmin constraint and the use of implicit differentiation, which may introduce non-linear dependencies between inner- and outer-level variables, making the exchange of differentiation and discretization nontrivial. Next, we leverage the expression of the plug-in estimator $\nabla \widehat{\mathcal{F}}(\omega)$ to provide generalization bounds.

# 4 Generalization bounds for (KBO)

In this section, we present our main result: a maximal inequality that controls how well both $\mathcal{F}$ and $\nabla \mathcal{F}$ are approximated by their empirical counterparts, uniformly over a compact subset $\Omega$ of $\mathbb{R}^d$.

## 4.1 Maximal inequalities for (KBO)

The following theorem provides finite-sample bounds on the uniform approximation errors on the objective and its gradient in expectation over both inner- and outer-level samples.

**Theorem 4.1 (Maximal inequalities).** Fix any compact subset $\Omega$ of $\mathbb{R}^d$. Under Assumptions (A) to (E), the following maximal inequalities hold:

$$ \mathbb{E} \left[ \sup_{\omega \in \Omega} \left| \mathcal{F}(\omega) - \widehat{\mathcal{F}}(\omega) \right| \right] \le C \left( \frac{1}{\sqrt{m}} + \frac{1}{\sqrt{n}} \right), \quad \mathbb{E} \left[ \sup_{\omega \in \Omega} \left\| \nabla \mathcal{F}(\omega) - \nabla \widehat{\mathcal{F}}(\omega) \right\| \right] \le C \left( \frac{1}{\sqrt{m}} + \frac{1}{\sqrt{n}} \right) $$

where the expectation is taken over the finite samples, and $C$ is a constant that depends only on $\Omega$, the dimension $d$, the regularization parameter $\lambda$, $\kappa$, and local upper bounds on $\ell_{in}$, $\ell_{out}$, and their partial derivatives over suitable compact sets.

**Theorem 4.1** states that the estimation error can be decomposed into two contributions each resulting from finite-sample approximation of $L_{in}$ and $L_{out}$ with a parametric rate of $1/\sqrt{m}$ and $1/\sqrt{n}$ up to a constant factor $C$. We provide a detailed expression for the constant in Theorem E.7 of Appendix E. The restriction to a compact subset $\Omega$ instead of the whole space $\mathbb{R}^d$ allows controlling the complexity of some function classes indexed by the parameter $\omega$. Without further assumptions on the objectives, we obtain a constant $C$ that can grow with the diameter of the subset $\Omega$.

**Role of $\lambda$.** The regularization parameter $\lambda$ simultaneously controls the strong convexity of the inner objective (see Proposition B.3 of Appendix B.1), the boundedness and Lipschitz continuity of the inner solutions (see Appendix D.1), the smoothness of the outer objective (see Proposition D.4 of Appendix D.2), the modulus of continuity of $L_{in}$, $\hat{L}_{in}$, $L_{out}$, $\hat{L}_{out}$ and their partial derivatives (see Appendix E.1), and maximal inequalities for certain processes (see Appendix E.2). Larger values of $\lambda$ yield smoother problems that are faster to optimize with larger step sizes, but introduce a larger statistical bias, while smaller values of $\lambda$ reduce bias but make optimization and generalization more delicate, with the error tending to $+\infty$ as $\lambda \to 0$. Selecting $\lambda$ therefore involves a trade-off between bias (regularized vs un-regularized problems), variance (finite sample vs population, as
done in our study), and *optimization efficiency* (step size). Under our assumptions, the dependence of the constants on $\lambda$ cannot be quantified, precluding generalization guarantees as $\lambda \to 0$. Stronger assumptions are required to obtain a quantitative control on $\lambda$.

**Probabilistic and variance bounds.** The most difficult quantity to control is the expectation of the maximal differences, which we have established. Once this expectation is bounded, a high-probability bound on the maximal differences can be derived via Markov's inequality. Moreover, since these differences are bounded (see Proposition E.4), we can also bound their variance. Indeed, let $Z \in [0, z]$ be a random variable representing the maximal difference between $\mathcal{F}$ or $\nabla \mathcal{F}$ and their respective plug-in estimators. Then, $\mathrm{Var}(Z/z) \le \mathbb{E}[(Z/z)^2] \le \mathbb{E}[Z/z]$, which implies that $\mathrm{Var}(Z) \le z\mathbb{E}[Z]$.

We outline the general proof strategy for Theorem 4.1 in the following section, with a full proof provided in Appendix E.

## 4.2 General proof strategy for Theorem 4.1

The main strategy behind the proof of Theorem 4.1 in Appendix E consists of three steps: (step 1) obtaining a point-wise error decomposition of the errors into manageable error terms that holds almost surely for any $\omega \in \Omega$, then applying maximal inequalities to suitable empirical processes (step 2) and some degenerate U-processes (step 3) to control each of these terms. The final error bounds are obtained by combining all these bounds as shown in Appendix E.3.

**Step 1: point-wise error decomposition.** A main challenge in controlling the errors in Theorem 4.1 is the non-linear dependence of both estimators $\hat{\mathcal{F}}(\omega)$ and $\nabla \hat{\mathcal{F}}(\omega)$ on the empirical distributions, as they are obtained via a plug-in procedure. We address this by breaking down the errors into components based on the discrepancies between expected values and their empirical counterparts of individual point-wise losses and their derivatives, all evaluated at the optimal solution $h_{\omega}^{*}$. Specifically, we denote by $\delta_{\omega}^{\text{out}}$ and $\delta_{\omega}^{\text{in}}$ the errors on the objectives defined as $\delta_{\omega}^{\text{out}} := |L_{\text{out}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{out}}(\omega, h_{\omega}^{*})|$ and $\delta_{\omega}^{\text{in}} := |L_{\text{in}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{in}}(\omega, h_{\omega}^{*})|$. Moreover, we quantify the errors between the partial derivatives of these objectives and their empirical counterparts. To simplify our proof outline, we slightly abuse notation by denoting $\partial_h \delta_{\omega}^{\text{out}}$, $\partial_h \delta_{\omega}^{\text{in}}$, $\partial_{\omega} \delta_{\omega}^{\text{out}}$, $\partial_{\omega,h}^2 \delta_{\omega}^{\text{in}}$, and $\partial_h^2 \delta_{\omega}^{\text{in}}$ to refer to these errors in terms of partial derivatives. For instance, $\partial_h \delta_{\omega}^{\text{out}}$ is defined as $\|\partial_h L_{\text{out}}(\omega, h_{\omega}^{*}) - \partial_h \hat{L}_{\text{out}}(\omega, h_{\omega}^{*})\|_{\mathcal{H}}$, with similar definitions for the other terms (see Appendix E.1). We get $|\mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega)| \le C(\delta_{\omega}^{\text{out}} + \partial_h \delta_{\omega}^{\text{in}})$ and $\|\nabla \mathcal{F}(\omega) - \nabla \hat{\mathcal{F}}(\omega)\| \le C(\partial_{\omega} \delta_{\omega}^{\text{out}} + \partial_h \delta_{\omega}^{\text{out}} + \partial_h^2 \delta_{\omega}^{\text{in}} + \partial_{\omega,h}^2 \delta_{\omega}^{\text{in}})$. Proposition E.4 formalizes this step and includes the exact constants. The error terms in both decompositions are amenable to a statistical analysis using empirical process theory as we discuss next.

**Step 2: maximal inequalities for empirical processes.** Some of the error terms, namely $\delta_{\omega}^{\text{out}}$ and $\partial_{\omega} \delta_{\omega}^{\text{out}}$, can be controlled directly using empirical process theory. For example, $\delta_{\omega}^{\text{out}}$ is associated to the family of random functions $\sqrt{m}(L_{\text{out}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{out}}(\omega, h_{\omega}^{*}))_{\omega \in \Omega}$, which defines an empirical process, a scaled and centered empirical average of real-valued functions indexed by the parameter $\omega$. Thus, provided that suitable estimates of the class complexity are available (as measured by its packing number in Proposition F.1 of Appendix F), which are easy to obtain in our setting, we show in Proposition E.5 of Appendix E.2 that a maximal inequality of the following form follows from classical results on empirical processes: $\mathbb{E}_{\mathcal{Q}}[\sup_{\omega \in \Omega} \delta_{\omega}^{\text{out}}] \le C/\sqrt{m}$.

**Step 3: maximal inequalities for degenerate U-processes.** Step 2 cannot be readily applied to the remaining terms involving partial derivatives w.r.t. $h(\partial_h \delta_{\omega}^{\text{out}}, \partial_h \delta_{\omega}^{\text{in}}, \partial_{\omega,h}^2 \delta_{\omega}^{\text{in}}, \partial_h^2 \delta_{\omega}^{\text{in}}) := D_h$. These are associated to processes that are not real-valued anymore, but take values in an infinite-dimensional space. In fact, one could apply step 2 to get an error per dimension, but then summing the errors yields a divergent sum. While the recent work in [56] develops an empirical process theory for functions taking values in a vector space, the provided complexity estimates would result in an unfavorable dependence on the sample size. Instead, we leverage the structure of the RKHS to control these errors using maximal inequalities for suitable degenerate U-processes of order 2 indexed by the parameter $\omega$ and for which such inequalities were provided in the seminal works of Sherman [67], Nolan and Pollard [54]. U-processes of order 2 are generalization of empirical processes and involve empirical averages of real-valued functions which depend on pairs of samples, instead of a single one as in empirical processes. In our case, these functions arise when taking the square of any term in $D_h$ and exploiting the reproducing property of the RKHS. This approach, presented in Proposition E.6 of Appendix E.2, allows us to obtain maximal inequalities for the terms in $D_h$. For example, it is of the
following form for $\partial_h \delta_{\omega}^{out}: \mathbb{E}_{\mathbb{Q}}[\sup_{\omega \in \Omega} \partial_h \delta_{\omega}^{out}] \leq C/\sqrt{m}$. Combining the maximal inequalities from steps 2 and 3 with the error decomposition from step 1 allows to obtain the result of Theorem 4.1.

**Discussion.** Alternative approaches to $U$-processes could be used to derive generalization bounds, although these would result in a degraded sample dependence. Specifically, one could employ a variational formulation of the RKHS norm appearing in some of the error terms, such as $\partial_h \delta_{\omega}^{out}$, to express them as the error of some real-valued empirical process to which standard results could be applied. However, this comes at the cost of considering processes indexed not only by the finite-dimensional parameter $\omega$, but also by functions in the unit RKHS ball. As a result, these families have much larger complexities as measured by their covering/packing numbers [77, Lemma D.2], which directly impacts the generalization rate. In contrast, our proposed approach bypasses this challenge by using real-valued $U$-processes indexed by finite-dimensional parameters, at the expense of employing a more general empirical process theory for degenerate $U$-processes [67].

To illustrate the implications of Theorem 4.1, we next provide convergence results for bilevel gradient methods.

## 4.3 Applications to empirical bilevel gradient methods

A typical strategy to solve (KBO) is to obtain empirical samples and apply a bilevel optimization algorithm to (KBO), for which our results offer statistical guarantees. Below, we present the generalization error for bilevel gradient descent. Generalization results for the projected bilevel gradient descent are directed to Appendix E.4.

**Bilevel gradient descent.** It is the simplest gradient-based method for solving the unconstrained (KBO) problem, i.e., when $C = \mathbb{R}^d$. It performs the update $\omega_{t+1} = \omega_t - \eta \nabla \hat{\mathcal{F}}(\omega_t)$ for all $t \ge 0$, where $\eta > 0$ is the step size. The algorithm requires access to the strongly convex inner-level solution and its derivative, which can be obtained using implicit differentiation.

**Corollary 4.2 (Generalization for bilevel gradient descent).** Consider Assumptions (A) to (E) and fix $\lambda > 0$. Assume further that $K$ in (KBO) is almost surely definite, and that there exists $c > 0$ such that $\inf_{\omega,v,y} \ell_{out}(\omega, v, y) - c\|\omega\|^2 > -\infty$. Fix $\omega_0 \in \mathbb{R}^d$ and let $\omega_{t+1} = \omega_t - \eta \nabla \hat{\mathcal{F}}(\omega_t)$ for all $t \ge 0$, where $\eta > 0$ is the step size. Then, there exist constants $\bar{\eta} > 0$ and $\bar{c} > 0$ such that for any $0 < \eta < \bar{\eta}$ any $t > 0$, the following holds:

$$ \mathbb{E} \left[ \min_{i=0,\dots,t} \| \nabla \mathcal{F}(\omega_i) \| \right] \le \bar{c} \left( \frac{1}{\sqrt{m}} + \frac{1}{\sqrt{n}} + \frac{1}{\sqrt{t+1}} \right), \mathbb{E} \left[ \limsup_{i \to \infty} \| \nabla \mathcal{F}(\omega_i) \| \right] \le \bar{c} \left( \frac{1}{\sqrt{m}} + \frac{1}{\sqrt{n}} \right) $$

The additional assumption on $\ell_{out}$ serves as a device to ensure the almost sure boundedness of the sequence a priori. It is rather mild, as it can be enforced by a small perturbation of the form $(\omega, v, y) \mapsto \ell_{out}(\omega, v, y) + c\|\omega\|^2$, assuming $\ell_{out} \ge 0$, which is typical in applications. Any other device that ensures a priori boundedness could also be considered. The assumption on $K$ is satisfied almost surely for most commonly used kernels. The proof of Corollary 4.2 follows from Theorem 4.1 and can be found in Appendix E.4. Corollary 4.2 also highlights a key algorithmic insight: the convergence of the bilevel method requires striking a balance between data availability (sample sizes $n$ and $m$) and computational budget (number of gradient steps). This trade-off arises because the total convergence error combines a statistical component, due to approximating the population gradient from finite samples, and an optimization component, due to performing only a limited number of (projected) gradient steps. Ensuring the right balance when designing practical algorithms prevents either insufficient data or limited computation from dominating the overall error.

## 5 Numerical experiments

**Setup.** To empirically validate our theoretical results, we consider the instrumental variable regression problem discussed in Section 2.2, in which we assume a linear dependence on $\omega$ for the function $f_{\omega}$ of the form $f_{\omega}(t) = \omega^{\top}\phi(t)$, where $\phi(t) \in \mathbb{R}^d$ denotes the feature map. We chose this particular problem because it allows us to derive closed-form expressions of the exact value function and its gradient, as well as their plug-in estimators, which are detailed in Appendix I.2. We use the Gaussian kernel and follow the experimental setup of Singh et al. [68], generating synthetic data that remain fixed across all runs. We vary $n$ between 100 and 5,000, setting $m = n$. The case
$n \neq m$ is analyzed separately in Appendix I.5. For the instrumental variable $x$, we consider two distributions: a $p$-dimensional standard Gaussian and a $p$-dimensional Student's $t$-distribution with degrees of freedom $\nu \in \{2.1, 2.5, 2.9\}$. Further details on the experimental setup are provided in Appendix I.4. We optimize the outer loss in (KBO) using gradient descent, where the step size is selected using backtracking line search and $\omega_0$ is randomly drawn from $U(0, 1)^d$. The stopping criterion is when $\|\nabla \hat{F}(\omega_i)\| \le 10^{-5}$, where $i$ is the iteration index. Our code is available at https://github.com/fareselkhoury/KBO.

**Scalable approximations for $\mathcal{F}$ and $\nabla \mathcal{F}$**. Since the expressions of $\mathcal{F}$ and $\nabla \mathcal{F}$ involve expectations, they are intractable to compute exactly. We approximate them accurately using their plug-in estimators $\hat{\mathcal{F}}$ and $\nabla \hat{\mathcal{F}}$, derived in Appendix I.2, and evaluate them using a large number of samples. To make this computation scalable, we approximate the kernel using a *random Fourier features* approximation [60], as detailed in Appendix I.3. Specifically, we use 1,000,000 samples and 26,000 random features, and handle memory constraints via a block decomposition strategy with a block size of 1,000.

**Results.** The plots in Figure 2 show the generalization behavior as a function of the number of inner samples $n$. (a) and (b) display the generalization error at initialization for the value function and the gradient, respectively. (c) presents the generalization bound for the gradient norm at the final iteration, while (d) shows the bound for the minimum gradient norm across all iterations. These results align with our theoretical findings, as all curves closely follow the expected theoretical slope. Additionally, Figure 3 of Appendix I.5 shows that balanced sample sizes lead to improved optimization behavior.

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of gradient descent on (KBO) for the instrumental variable regression task using synthetic data. The plots are averaged over 50 runs and displayed on a log-log scale. The line represents the mean across all runs, and the shaded region indicates the 95% confidence interval.</figcaption></figure>

# 6 Conclusion and perspectives

**Summary.** In this work, we established the first generalization bounds for (KBO). These results are crucial for understanding the generalization properties of algorithms for solving (KBO). They offer rigorous guarantees on the algorithm's performance on unseen data—a fundamental criterion for any algorithmic design—and help control overfitting. Given that our bounds are of order $O(1/\sqrt{m} + 1/\sqrt{n})$, this highlights the equal importance of both outer- and inner-level sample sizes to the overall generalization error. Our findings can impact current practices, particularly in hyperparameter optimization, where the validation dataset is typically much smaller than the training set.

**Limitations and future work.** This paper takes a first step toward providing generalization results for bilevel gradient-based methods in a nonparametric setting. While our theoretical analysis focused on a full-batch bilevel setting with *exact gradients*, extending this framework to *stochastic* variants, such as those in [3, 29, 21, 24], remains an open challenge. A promising direction would be to consider *approximate kernel representations*, such as random Fourier features or neural tangent kernels, which enable scalable learning using kernel methods while preserving useful theoretical properties. Furthermore, the constants in our bounds are likely *conservative*; we did not investigate their tightness or potential for improvement. A deeper analysis of their optimality could provide valuable insights and constitutes an avenue worth exploring. Additionally, providing generalization guarantees for the un-regularized problem, possibly in the form of minimax optimal rates for (KBO), is a worthwhile future direction. This requires controlling the constants as $\lambda \to 0$, provided additional source assumptions are made, as discussed in [69, 70]. Finally, broadening our framework to cover *non-smooth losses*, such as the hinge loss in SVMs, is an interesting direction for future work.
## Acknowledgments and Disclosure of Funding

This work was supported by the ANR project BONSAI (grant ANR-23-CE23-0012-01). EP acknowledges funding from the AI Interdisciplinary Institute ANITI through the French Investments for the Future - PIA3 program under grant agreement ANR-19-PI3A-0004; the Air Force Office of Scientific Research, Air Force Materiel Command, USAF under grant number FA8655-22-1-7012; TSE-P; the ANR project Chess (grant ANR-17-EURE-0010); the ANR project Regulia; and the Institut Universitaire de France (IUF). EP and SV acknowledge support from the ANR project MAD. SV also thanks the PEPR PDE-AI program (ANR-23-PEIA-0004) and the 3IA Côte d'Azur Chair BOGL.

## References

[1] C. Ai and X. Chen. Efficient estimation of models with conditional moment restrictions containing unknown functions. *Econometrica*, 71(6):1795-1843, 2003.

[2] C. Ai and X. Chen. Estimation of possibly misspecified semiparametric conditional moment restriction models with different conditioning variables. *Journal of Econometrics*, 141(1):5-43, 2007.

[3] M. Arbel and J. Mairal. Amortized implicit differentiation for stochastic bilevel optimization. In *International Conference on Learning Representations (ICLR)*, 2022.

[4] M. Arbel and J. Mairal. Non-convex bilevel games with critical point selection maps. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2022.

[5] S. Arora, S. Du, S. Kakade, Y. Luo, and N. Saunshi. Provable representation learning for imitation learning via bi-level optimization. In *International Conference on Machine Learning (ICML)*, 2020.

[6] S. Bai, J. Z. Kolter, and V. Koltun. Deep equilibrium models. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2019.

[7] F. Bao, G. Wu, C. Li, J. Zhu, and B. Zhang. Stability and generalization of bilevel programming in hyperparameter optimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2021.

[8] F. Bauer, S. V. Pereverzev, and L. Rosasco. On regularization algorithms in learning theory. *Journal of Complexity*, 23(1):52-72, 2007.

[9] A. Beck. *Introduction to nonlinear optimization: Theory, algorithms, and applications with MATLAB*. SIAM, 2014.

[10] A. Beck. *First-order methods in optimization*. SIAM, 2017.

[11] M. Belkin, S. Ma, and S. Mandal. To understand deep learning we need to understand kernel learning. In *International Conference on Machine Learning (ICML)*, 2018.

[12] S. Ben-David, J. Blitzer, K. Crammer, and F. Pereira. Analysis of representations for domain adaptation. In *Advances in Neural Information Processing Systems (NIPS)*, 2006.

[13] Y. Bengio. Gradient-based optimization of hyperparameters. *Neural Computation*, 12(8): 1889-1900, 2000. ISSN 0899-7667, 1530-888X. doi: 10.1162/089976600300015187.

[14] L. Bertinetto, J. F. Henriques, P. Torr, and A. Vedaldi. Meta-learning with differentiable closed-form solvers. In *International Conference on Learning Representations (ICLR)*, 2019.

[15] M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert. Efficient and modular implicit differentiation. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2022.

[16] S. Bochner. *Lectures on Fourier integrals*, volume 42. Princeton University Press, 1959.

[17] J. Bolte, Q.-T. Lê, E. Pauwels, and S. Vaiter. Geometric and computational hardness of bilevel programming. *Mathematical Programming*, pages 1-36, 2025.
[18] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http://github.com/jax-ml/jax.

[19] W. Candler and R. Norton. *Multi-Level Programming and Development Policy*, volume 1. World Bank, 1977.

[20] A. Caponnetto and E. De Vito. Optimal rates for the regularized least-squares algorithm. *Foundations of Computational Mathematics*, 7(3):331–368, 2007. doi: 10.1007/s10208-006-0196-8.

[21] T. Chen, Y. Sun, and W. Yin. Closing the gap: Tighter analysis of alternating stochastic gradient methods for bilevel problems. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2021.

[22] X. Chen and D. Pouzo. Estimation of nonparametric conditional moment models with possibly nonsmooth generalized residuals. *Econometrica*, 80(1):277–321, 2012.

[23] X. Chen and D. Pouzo. Sieve Wald and QLR inferences on semi/nonparametric conditional moment models. *Econometrica*, 83(3):1013–1079, 2015.

[24] M. Dagréou, P. Ablin, S. Vaiter, and T. Moreau. A framework for bilevel optimization that enables stochastic and global variance reduction algorithms. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2022.

[25] J. Diestel and J. J. Uhl. *Vector Measures*, volume 15 of *Mathematical Surveys and Monographs*. American Mathematical Society, Providence, RI, 1977. ISBN 978-0-8218-1515-1.

[26] N. Doumèche, F. Bach, G. Biau, and C. Boyer. Physics-informed kernel learning. *Journal of Machine Learning Research (JMLR)*, 2025.

[27] H. W. Engl, M. Hanke, and A. Neubauer. *Regularization of Inverse Problems*. Springer, 1996.

[28] L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil. Bilevel programming for hyperparameter optimization and meta-learning. In *International Conference on Machine Learning (ICML)*, 2018.

[29] S. Ghadimi and M. Wang. Approximation methods for bilevel programming. *arXiv preprint arXiv:1802.02246*, 2018.

[30] I. Goodfellow, Y. Bengio, and A. Courville. *Deep Learning*. The MIT Press, 2016. http://www.deeplearningbook.org.

[31] A. Griewank and C. Faure. Piggyback differentiation and optimization. In *Large-scale PDE-constrained optimization*, pages 148–164. Springer, 2003.

[32] J. Hartford, G. Lewis, K. Leyton-Brown, and M. Taddy. Deep IV: A flexible approach for counterfactual prediction. In *International Conference on Machine Learning (ICML)*, 2017.

[33] T. Hofmann, B. Schölkopf, and A. Smola. Kernel methods in machine learning. *The Annals of Statistics*, 36(3):1171–1220, 2008.

[34] G. Holler, K. Kunisch, and R. C. Barnard. A bilevel approach for parameter learning in inverse problems. *Inverse Problems*, 34(11):115012, 2018. doi: 10.1088/1361-6420/aae473.

[35] M. Hong, H.-T. Wai, Z. Wang, and Z. Yang. A two-timescale stochastic algorithm framework for bilevel optimization: Complexity analysis and application to actor-critic. *SIAM Journal on Optimization*, 33(1):147–180, 2023.

[36] A. D. Ioffe and V. M. Tihomirov. *Theory of Extremal Problems*. Series: Studies in Mathematics and its Applications 6. Elsevier, 1979.

[37] A. Jacot, F. Gabriel, and C. Hongler. Neural Tangent Kernel: Convergence and Generalization in Neural Networks. In *Advances in Neural Information Processing Systems (NIPS)*, 2018.
[38] K. Ji, J. Yang, and Y. Liang. Bilevel optimization: Convergence analysis and enhanced design. In *International Conference on Machine Learning (ICML)*, 2021.

[39] S. Keerthi, V. Sindhwani, and O. Chapelle. An Efficient Method for Gradient-Based Adaptation of Hyperparameters in SVM Models. In *Advances in Neural Information Processing Systems (NIPS)*, 2006.

[40] M. R. Kosorok. *Introduction to Empirical Processes and Semiparametric Inference*. Springer Series in Statistics. Springer, New York, 2008. URL https://doi.org/10.1007/978-0-387-74978-5.

[41] G. Kunapuli, K. P. Bennett, J. Hu, and J.-S. Pang. Classification model selection via bilevel programming. *Optimization Methods & Software*, 23(4):475-489, 2008.

[42] J. Kwon, D. Kwon, S. Wright, and R. D. Nowak. On penalty methods for nonconvex bilevel optimization and first-order stochastic approximation. In *International Conference on Learning Representations (ICLR)*, 2024.

[43] J. Larsen, L. Hansen, C. Svarer, and M. Ohlsson. Design and regularization of neural networks: The optimal use of a validation set. In *Neural Networks for Signal Processing VI. Proceedings of the 1996 IEEE Signal Processing Society Workshop*, pages 62-71, Kyoto, Japan, 1996. IEEE. doi: 10.1109/NNSP.1996.548300.

[44] M. Letizia, G. Losapio, M. Rando, G. Grosso, A. Wulzer, M. Pierini, M. Zanetti, and L. Rosasco. Learning new physics efficiently with nonparametric methods. *The European Physical Journal C*, 82(10):879, 2022.

[45] Z. Li, H. Lan, V. Syrgkanis, M. Wang, and M. Uehara. Regularized DeepIV with Model Selection. *arXiv preprint arXiv:2403.04236*, 2024.

[46] Z. Li, D. Meunier, M. Mollenhauer, and A. Gretton. Towards Optimal Sobolev Norm Rates for the Vector-Valued Regularized Least-Squares Algorithm. *Journal of Machine Learning Research (JMLR)*, 25(181):1-51, 2024. URL https://jmlr.org/papers/v25/23-1663.html.

[47] R. Liu, Y. Liu, S. Zeng, and J. Zhang. Towards gradient-based bilevel optimization with non-convex followers and beyond. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2021.

[48] R. Liu, X. Liu, S. Zeng, J. Zhang, and Y. Zhang. Value-function-based sequential minimization for bi-level optimization. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 2023.

[49] L. Lo Gerfo, L. Rosasco, F. Odone, E. De Vito, and A. Verri. Spectral Algorithms for Supervised Learning. *Neural Computation*, 20(8):1873-1897, 2008. doi: 10.1162/neco.2008.05-07-517.

[50] J. Mercer. Functions of positive and negative type, and their connection with the theory of integral equations. *Philosophical Transactions of the Royal Society A*, 209:415-446, 1909.

[51] D. Meunier, Z. Li, T. Christensen, and A. Gretton. Nonparametric Instrumental Regression via Kernel Methods is Minimax Optimal. *arXiv preprint arXiv:2411.19653*, 2024.

[52] D. Meunier, Z. Shen, M. Mollenhauer, A. Gretton, and Z. Li. Optimal Rates for Vector-Valued Spectral Regularization Learning Algorithms. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2024.

[53] W. K. Newey and J. L. Powell. Instrumental variable estimation of nonparametric models. *Econometrica*, 71(5):1565-1578, 2003.

[54] D. Nolan and D. Pollard. U-processes: rates of convergence. *The Annals of Statistics*, pages 780-799, 1987.

[55] S. Oymak, M. Li, and M. Soltanolkotabi. Generalization guarantees for neural architecture search with train-validation split. In *International Conference on Machine Learning (ICML)*, 2021.
[56] J. Park and K. Muandet. Towards empirical process theory for vector-valued functions: Metric entropy of smooth function classes. In *International Conference on Algorithmic Learning Theory*, pages 1216–1260. PMLR, 2023.

[57] F. Pedregosa. Hyperparameter optimization with approximate gradient. In *International Conference on Machine Learning (ICML)*, 2016.

[58] I. Petrulionytė, J. Mairal, and M. Arbel. Functional Bilevel Optimization for Machine Learning. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2024.

[59] Q. Pham, C. Liu, D. Sahoo, and H. Steven. Contextual transformation networks for online continual learning. In *International Conference on Learning Representations (ICLR)*, 2021.

[60] A. Rahimi and B. Recht. Random Features for Large-Scale Kernel Machines. In *Advances in Neural Information Processing Systems (NIPS)*, 2007.

[61] C. E. Rasmussen and C. K. I. Williams. *Gaussian Processes for Machine Learning*. The MIT Press, 2006.

[62] W. Rudin. *Real and Complex Analysis*. McGraw-Hill, New York, 3rd edition, 1987. ISBN 978-0070542341.

[63] B. Schölkopf and A. J. Smola. *Learning with kernels: support vector machines, regularization, optimization, and beyond*. The MIT Press, 2002.

[64] B. Schölkopf, R. Herbrich, and A. J. Smola. A Generalized Representer Theorem. In D. Helmbold and B. Williamson, editors, *Computational Learning Theory*, pages 416–426, Berlin, Heidelberg, 2001. Springer Berlin Heidelberg. ISBN 978-3-540-44581-4.

[65] J. Shawe-Taylor and N. Cristianini. *Kernel Methods for Pattern Analysis*. Cambridge University Press, Cambridge, UK, 2004.

[66] H. Shen and T. Chen. On penalty-based bilevel gradient descent method. In *International Conference on Machine Learning (ICML)*, 2023.

[67] R. P. Sherman. Maximal Inequalities for Degenerate *U*-Processes with Applications to Optimization Estimators. *The Annals of Statistics*, 22(1):439 – 459, 1994. URL https://doi.org/10.1214/aos/1176325377.

[68] R. Singh, M. Sahani, and A. Gretton. Kernel Instrumental Variable Regression. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2019.

[69] S. Smale and D.-X. Zhou. Learning theory estimates via integral operators and their approximations. *Constructive Approximation*, 26(2):153–172, 2007.

[70] B. K. Sriperumbudur, K. Fukumizu, A. Gretton, A. Hyvärinen, and R. Kumar. Density estimation in infinite dimensional exponential families. *Journal of Machine Learning Research (JMLR)*, 18(57):1–59, 2017.

[71] A. W. van der Vaart. *Asymptotic Statistics*, volume 3. Cambridge University Press, 2000.

[72] A. W. van der Vaart and J. A. Wellner. *Weak Convergence and Empirical Processes: With Applications to Statistics*. Springer Series in Statistics. Springer-Verlag, New York, 1996. doi: 10.1007/978-1-4757-2545-2.

[73] J. Wang, C. Lan, C. Liu, Y. Ouyang, T. Qin, W. Lu, Y. Chen, W. Zeng, and S. Y. Philip. Generalizing to unseen domains: A survey on domain generalization. *IEEE transactions on knowledge and data engineering*, 35(8):8052–8072, 2022.

[74] R. Wang, C. Zheng, G. Wu, X. Min, X. Zhang, J. Zhou, and C. Li. Lower bounds of uniform stability in gradient-based bilevel algorithms for hyperparameter optimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2024.

[75] C. Williams and M. Seeger. Using the Nyström Method to Speed Up Kernel Machines. In *Advances in Neural Information Processing Systems (NIPS)*, 2000.
[76] L. Xu, Y. Chen, S. Srinivasan, N. de Freitas, A. Doucet, and A. Gretton. Learning Deep Features in Instrumental Variable Regression. In *International Conference on Learning Representations (ICLR)*, 2021.

[77] Z. Yang, C. Jin, Z. Wang, M. Wang, and M. I. Jordan. On Function Approximation in Reinforcement Learning: Optimism in the Face of Large State Spaces. In *Advances in Neural Information Processing Systems (NeurIPS)*, 2020.

[78] X. Zhang, H. Chen, B. Gu, T. Gong, and F. Zheng. Fine-grained analysis of stability and generalization for stochastic bilevel optimization. In *Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence*, pages 5508–5516, 2024.
# Appendices

**Roadmap.** In Appendix A, we review existing non-kernel generalization bounds in the bilevel optimization literature and explain why minimax optimal rates cannot be obtained in our setting. We begin the theoretical appendices by presenting and establishing regularity properties of the objective functions in Appendix B. In Appendix C, we introduce the gradient estimators. Appendix D is dedicated to proving the boundedness and Lipschitz continuity of $h_{\omega}^{*}$ and $\hat{h}_{\omega}$, along with local boundedness and Lipschitz properties of $\ell_{in}$, $\ell_{out}$, and their derivatives. The generalization results are provided in Appendix E. In Appendix F, we establish maximal inequalities for bounded and Lipschitz families of functions. Differentiability properties of the objectives are studied in Appendix G. Appendix H contains auxiliary technical lemmas used throughout the proofs. Finally, further details on the experiments and additional numerical results are provided in Appendix I.

**Notations.** $\| \cdot \|$ denotes the Euclidean norm in $\mathbb{R}^d$, $\| \cdot \|_{\mathcal{H}}$ denotes the norm in the RKHS $\mathcal{H}$, $\| \cdot \|_{op}$ denotes the operator norm, and $\| \cdot \|_{HS}$ denotes the Hilbert-Schmidt norm. $\langle \cdot, \cdot \rangle_{\mathcal{H}}$ denotes the inner product on $\mathcal{H}$, and $\langle \cdot, \cdot \rangle_{HS}$ denotes the Hilbert-Schmidt inner product. $K(x, \cdot)$ denotes the feature map, for any $x \in \mathcal{X}$. For any two normed spaces $E$ and $F$, $\mathcal{L}(E, F)$ denotes the space of continuous linear operators from $E$ to $F$. For any two probability distributions $\mathcal{P}$ and $\mathcal{Q}$, $\mathcal{P} \otimes \mathcal{Q}$ denotes the product measure of $\mathcal{P}$ and $\mathcal{Q}$. Given two Hilbert spaces $(H_1, \langle \cdot, \cdot \rangle_{H_1})$ and $(H_2, \langle \cdot, \cdot \rangle_{H_2})$, the tensor product of $u \in H_1$ and $v \in H_2$, denoted by $u \otimes v$, is an operator from $H_2$ to $H_1$ defined, for any $e \in H_2$, as $(u \otimes v)e = u \langle v, e \rangle_{H_2}$. For any $v_1, \dots, v_n \in \mathbb{R}$, $\text{diag}(v_1, \dots, v_n) \in \mathbb{R}^{n \times n}$ denotes a diagonal matrix of size $n \times n$, where the diagonal entries are $v_1, \dots, v_n$ and all the off-diagonal entries are 0. $\mathbb{1}_m$ denotes a vector of size $m$ where all entries are 1. $\mathbb{1}_{n \times n}$ denotes the identity matrix of size $n$. For any vector space $V$ over $\mathbb{R}$, $\text{Id}_V$ denotes the identity operator on $V$. Given a compact set $\mathcal{K}$, $\text{diam}(\mathcal{K})$ denotes its diameter. $v^{\top}$ denotes the transpose of either a vector or a matrix, depending on the context.

# Contents

**A Further Discussion** 17

**B Regularity and Differentiability Results** 17

B.1 Regularity of the objectives 17

B.2 Differentiability of the value function 18

**C Gradient Estimators** 19

**D Preliminary Results** 21

D.1 Boundedness and Lipschitz continuity of $h_{\omega}^{*}$ and $\hat{h}_{\omega}$ 22

D.2 Local boundedness and Lipschitz properties of $\ell_{in}$, $\ell_{out}$, and their derivatives 23

**E Generalization Properties** 24

E.1 Point-wise estimates 24

E.2 Maximal inequalities 29

E.3 Proof of Theorem 4.1 33

E.4 Generalization for bilevel gradient methods 33

**F Maximal Inequalities for Bounded and Lipschitz Family of Functions** 34

**G Differentiability Results** 38
Similarly, the empirical estimates $\hat{L}_{in}$ and $\hat{L}_{out}$ admit finite values, and are differentiable with gradients admitting similar expressions as above with $\mathbb{P}$ and $\mathbb{Q}$ replaced by their empirical estimates $\hat{\mathbb{P}}_n$ and $\hat{\mathbb{Q}}_m$.

**Proposition B.2 (Differentiability of $\partial_h L_{in}$).** Under Assumptions (A) to (D), for any $(\omega, h) \in \mathbb{R}^d \times \mathcal{H}$, the function $(\omega, h) \mapsto \partial_h L_{in}(\omega, h)$ is differentiable with partial derivatives given by:

$$ \partial_{\omega,h}^2 L_{in}(\omega, h) = \mathbb{E}_{\mathbb{P}} \left[ \partial_{\omega,v}^2 \ell_{in}(\omega, h(x), y) K(x, \cdot) \right] \in \mathcal{L}(\mathcal{H}, \mathbb{R}^d), $$

$$ \partial_h^2 L_{in}(\omega, h) = \mathbb{E}_{\mathbb{P}} \left[ \partial_v^2 \ell_{in}(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot) \right] + \lambda \mathrm{Id}_{\mathcal{H}} \in \mathcal{L}(\mathcal{H}, \mathcal{H}). $$

Moreover, for any $\omega \in \mathbb{R}^d$ and $h \in \mathcal{H}$, the operators $\partial_{\omega,h}^2 L_{in}(\omega, h)$ and $\partial_h^2 L_{in}(\omega, h) - \lambda \mathrm{Id}_{\mathcal{H}}$ are Hilbert-Schmidt, i.e., bounded operators with finite Hilbert-Schmidt norm. The same conclusions hold for the empirical estimate $(\omega, h) \mapsto \partial_h \hat{L}_{in}(\omega, h)$ with partial derivatives admitting similar expressions as above with $\mathbb{P}$ replaced by its empirical estimate $\hat{\mathbb{P}}_n$.

**Proposition B.3 (Strong convexity of the inner objective in its second variable and invertibility of the Hessians).** Under Assumptions (A) to (E), $h \mapsto L_{in}(\omega, h)$ and $h \mapsto \hat{L}_{in}(\omega, h)$ are $\lambda$-strongly convex for any $\omega \in \mathbb{R}^d$. Moreover, for any $\omega \in \mathbb{R}^d$ and $h \in \mathcal{H}$, the Hessian operators $\partial_h^2 L_{in}(\omega, h)$ and $\partial_h^2 \hat{L}_{in}(\omega, h)$ are invertible with their operator norm bounded by $\frac{1}{\lambda}$.

*Proof.* By Assumption (E), we know that $v \mapsto \ell_{in}(\omega, v, y)$ is convex for any $\omega \in \mathbb{R}^d$ and $y \in \mathcal{Y}$. Moreover, by Proposition B.1, $(x, y) \mapsto \ell_{in}(\omega, h(x), y)$ is integrable for any $\omega \in \mathbb{R}^d$ and $h \in \mathcal{H}$. Consequently, by integration, we directly deduce that $h \mapsto \mathbb{E}_{\mathbb{P}}[\ell_{in}(\omega, h(x), y)]$ is convex for any $\omega \in \mathbb{R}^d$. Finally, $h \mapsto L_{in}(\omega, h) := \mathbb{E}_{\mathbb{P}}[\ell_{in}(\omega, h(x), y)] + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^2$ must be $\lambda$-strongly convex, for any $\omega \in \mathbb{R}^d$, as a sum of a convex function and a $\lambda$-strongly convex function. Similarly, we deduce that $h \mapsto \hat{L}_{in}(\omega, h)$ is $\lambda$-strongly convex, for any $\omega \in \mathbb{R}^d$. Invertibility follows from the expression of the Hessian operator in Proposition B.2 $\square$

## B.2 Differentiability of the value function

**Proposition B.4 (Total functional gradient $\nabla \mathcal{F}$).** Assume Assumptions (A), (B), (D) and (E) hold. For any $\omega \in \mathbb{R}^d$, the total functional gradient $\nabla \mathcal{F}(\omega)$ satisfies:

$$ \nabla \mathcal{F}(\omega) = \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) + \partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*}) a_{\omega}^{*} \in \mathbb{R}^d, \quad (6) $$

where $a_{\omega}^{*}$ is the unique minimizer of the following quadratic objective:

$$ L_{adj}(\omega, a) := \frac{1}{2} \langle a, H_{\omega}a \rangle_{\mathcal{H}} + \langle a, d_{\omega} \rangle_{\mathcal{H}}, \quad \text{for any } a \in \mathcal{H}, \quad (7) $$

with $H_{\omega} := \partial_h^2 L_{in}(\omega, h_{\omega}^{*}) : \mathcal{H} \to \mathcal{H}$ being the Hessian operator and $d_{\omega} := \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) \in \mathcal{H}$.

*Proof.* By applying Propositions B.1 and B.3, we know that $h \mapsto L_{in}(\omega, h)$ has finite values, is $\lambda$-strongly convex and Fréchet differentiable. Moreover, by Proposition B.2, $\partial_h L_{in}$ is Fréchet differentiable on $\mathbb{R}^d \times \mathcal{H}$, and, a fortiori, Hadamard differentiable. Therefore, by the functional implicit differentiation theorem [36, 58, Theorem 2.1], we deduce that the map $\omega \mapsto h_{\omega}^{*}$ is uniquely defined and is Fréchet differentiable with Jacobian $\partial_{\omega} h_{\omega}^{*}$ solving the following linear system for any $\omega \in \mathbb{R}^d$:

$$ \partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*}) + \partial_{\omega} h_{\omega}^{*} \partial_h^2 L_{in}(\omega, h_{\omega}^{*}) = 0. $$

Using that $\partial_h^2 L_{in}(\omega, h_{\omega}^{*})$ is invertible by Proposition B.3, we can express $\partial_{\omega} h_{\omega}^{*}$ as:

$$ \partial_{\omega} h_{\omega}^{*} = -\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*}) (\partial_h^2 L_{in}(\omega, h_{\omega}^{*}))^{-1}. $$

Furthermore, $L_{out}$ is jointly Fréchet differentiable by application of Proposition B.1, so that $\omega \mapsto \mathcal{F}(\omega)$ is also differentiable by composition of the functions $(\omega, h) \mapsto L_{out}(\omega, h)$ and $\omega \mapsto (\omega, h_{\omega}^{*})$. For a given $\omega \in \mathbb{R}^d$, the gradient of $\mathcal{F}$ is then given by the chain rule:

$$ \nabla \mathcal{F}(\omega) = \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) + \partial_{\omega} h_{\omega}^{*} \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}). \quad (8) $$

Substituting the expression of $\partial_{\omega} h_{\omega}^{*}$ into Equation (8) yields:

$$ \nabla \mathcal{F}(\omega) = \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*}) (\partial_h^2 L_{in}(\omega, h_{\omega}^{*}))^{-1} \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}). $$
To conclude, it suffices to notice that the function $a_{\omega}^{*}$ appearing in Equation (6) must be equal to $-H_{\omega}^{-1}d_{\omega}$. Indeed, $a_{\omega}^{*}$ is defined as the minimizer of the quadratic objective $L_{adj}(\omega, a)$ in Equation (7) which is strongly convex since the Hessian operator is lower-bounded by $\lambda \mathrm{Id}_{\mathcal{H}}$. Consequently, the minimizer $a_{\omega}^{*}$ exists and is uniquely characterized by the optimality condition:

$$H_{\omega}a_{\omega}^{*} + d_{\omega} = 0.$$

The above equation is a linear system in $\mathcal{H}$ whose solution is given by $a_{\omega}^{*} := -H_{\omega}^{-1}d_{\omega}$. $\square$

## C Gradient Estimators

**Proposition C.1 (Expression of $\nabla \hat{\mathcal{F}}(\omega)$ by implicit differentiation).** Under Assumptions (B) to (E), for any $\omega \in \mathbb{R}^d$, the gradient $\nabla \hat{\mathcal{F}}(\omega)$ of the discretized kernel bilevel optimization problem ($\overline{\text{KBO}}$) is given by:

$$\nabla \hat{\mathcal{F}}(\omega) = \frac{1}{m} \mathbf{D}_{\omega}^{\text{out}} \mathbf{1}_{m} - \frac{1}{m} \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} \mathbf{M}^{-1} \mathbf{u} \in \mathbb{R}^{d},$$

where $\mathbf{K}$ and $\overline{\mathbf{K}}$ are the Gram matrices in $\mathbb{R}^{n \times n}$ and $\mathbb{R}^{m \times n}$ with entries given by $\mathbf{K}_{ij} := K(x_i, x_j)$ and $\overline{\mathbf{K}}_{ij} := K(\tilde{x}_i, x_j)$, and $\mathbf{M} \in \mathbb{R}^{n \times n}$, $\mathbf{u} \in \mathbb{R}^n$, $\mathbf{D}_{\mathbf{v}}^{\text{out}} \in \mathbb{R}^m$, $\mathbf{D}_{\omega}^{\text{out}} \in \mathbb{R}^{d \times m}$, $\mathbf{D}_{\mathbf{v}, \mathbf{v}}^{\text{in}} \in \mathbb{R}^{n \times n}$, and $\mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} \in \mathbb{R}^{d \times n}$ are defined as:

$$\begin{align*}
\mathbf{M} &:= \mathbf{K} \mathbf{D}_{\mathbf{v}, \mathbf{v}}^{\text{in}} + n\lambda \mathbf{1}_{n \times n}, & \mathbf{u} &:= \overline{\mathbf{K}}^{\top} \mathbf{D}_{\mathbf{v}}^{\text{out}}, \\
\mathbf{D}_{\mathbf{v}}^{\text{out}} &:= \left( \partial_{\mathbf{v}} \ell_{\text{out}} (\omega, \hat{h}_{\omega}(\tilde{x}_{\mathbf{v}}), \tilde{y}_{\mathbf{v}}) \right)_{1 \le j \le m}, & \mathbf{D}_{\omega}^{\text{out}} &:= \left( \partial_{\omega} \ell_{\text{out}} (\omega, \hat{h}_{\omega}(\tilde{x}_{\mathbf{v}}), \tilde{y}_{\mathbf{v}}) \right)_{1 \le j \le m}, \\
\mathbf{D}_{\mathbf{v}, \mathbf{v}}^{\text{in}} &:= \text{diag} \left( \left( \partial_{\mathbf{v}}^2 \ell_{\text{in}} (\omega, \hat{h}_{\omega}(x_i), y_i) \right)_{1 \le i \le n} \right), & \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} &:= \left( \partial_{\omega, \mathbf{v}}^2 \ell_{\text{in}} (\omega, \hat{h}_{\omega}(x_i), y_i) \right)_{1 \le i \le n}.
\end{align*}$$

*Proof.* Let $\omega \in \mathbb{R}^d$. Recall the expression of $\hat{\mathcal{F}}(\omega)$:

$$\begin{align*}
\hat{\mathcal{F}}(\omega) &:= \frac{1}{m} \sum_{j=1}^{m} \ell_{\text{out}} (\omega, \hat{h}_{\omega}(\tilde{x}_{j}), \tilde{y}_{j}) \\
\text{s.t.} \quad \hat{h}_{\omega} &= \arg \min_{\hat{h}_{\omega} \in \mathcal{H}} \hat{L}_{\text{in}}(\omega, h) := \frac{1}{n} \sum_{i=1}^{n} \ell_{\text{in}} (\omega, h(x_i), y_i) + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^{2}.
\end{align*}$$

By the representer theorem, it is easy to see that $\hat{h}_{\omega}$ must be a linear combination of $K(x_1, \cdot), \dots, K(x_n, \cdot)$:

$$\hat{h}_{\omega} = \sum_{i=1}^{n} (\hat{\gamma}_{\omega})_{i} K(x_{i}, \cdot). \qquad (9)$$

Hence, finding $\hat{h}_{\omega}$ amounts to minimizing $\hat{L}_{\text{in}}(\omega, h)$ over the span of $(K(x_1, \cdot), \dots, K(x_n, \cdot))$, i.e., over functions $h^{\gamma}$ of the form $h^{\gamma} = \sum_{i=1}^{n} (\gamma)_{i} K(x_i, \cdot)$ for $\gamma \in \mathbb{R}^n$. Restricting the objective to such functions results in the following inner optimization problem which is finite-dimensional:

$$\hat{\gamma}_{\omega} := \arg \min_{\gamma \in \mathbb{R}^n} \frac{1}{n} \sum_{i=1}^{n} \ell_{\text{in}} (\omega, (\mathbf{K}\gamma)_i, y_i) + \frac{\lambda}{2} \gamma^{\top} \mathbf{K} \gamma,$$

where we used that $(h^{\gamma}(x_i))_{1 \le i \le n} = \mathbf{K}\gamma$ and $\|h^{\gamma}\|_{\mathcal{H}}^2 = \gamma^{\top}\mathbf{K}\gamma$. Similarly, using that $(h^{\gamma}(\tilde{x}_j))_{1 \le j \le m} = \overline{\mathbf{K}}\gamma$, we can express $\hat{\mathcal{F}}(\omega)$ as follows:

$$\hat{\mathcal{F}}(\omega) = \frac{1}{m} \sum_{j=1}^{m} \ell_{\text{out}} (\omega, (\overline{\mathbf{K}} \hat{\gamma}_{\omega})_j, \tilde{y}_j).$$

Differentiating the above expression w.r.t. $\omega$ and applying the chain rule result in:

$$\nabla \hat{\mathcal{F}}(\omega) = \frac{1}{m} \sum_{j=1}^{m} \partial_{\omega} \ell_{\text{out}} (\omega, (\overline{\mathbf{K}} \hat{\gamma}_{\omega})_j, \tilde{y}_j) + \frac{1}{m} \sum_{j=1}^{m} (\partial_{\omega} \hat{\gamma}_{\omega} \overline{\mathbf{K}}^{\top})_j \partial_{\mathbf{v}} \ell_{\text{out}} (\omega, (\overline{\mathbf{K}} \hat{\gamma}_{\omega})_j, \tilde{y}_j),$$
where $\partial_{\omega}\hat{\gamma}_{\omega}$ denotes the Jacobian of $\hat{\gamma}_{\omega}$. We can further express the above equation in matrix form to get:

$$ \nabla \hat{\mathcal{F}}(\omega) = \frac{1}{m} \mathbf{D}_{\omega}^{\text{out}} \mathbf{1}_{m} + \frac{1}{m} \partial_{\omega} \hat{\gamma}_{\omega} \bar{\mathbf{K}}^{\top} \mathbf{D}_{\omega}^{\text{out}}. \quad (10) $$

Moreover, an application of the implicit function theorem<fnref n="2" /> allows to directly express the Jacobian $\partial_{\omega}\hat{\gamma}_{\omega}$ as a solution of the following linear system obtained by differentiating the optimality condition for $\hat{\gamma}_{\omega}$ w.r.t. $\omega$:

$$ \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} \mathbf{K} + (\partial_{\omega} \hat{\gamma}_{\omega}) \underbrace{(\mathbf{K} \mathbf{D}_{\mathbf{v}, \mathbf{v}}^{\text{in}} + n\lambda \mathbf{1}_{n \times n})}_{\mathbf{M}} \mathbf{K} = 0. $$

A solution of the form $\partial_{\omega}\hat{\gamma}_{\omega} = -\mathbf{D}_{\omega, \mathbf{v}}^{\text{in}}\mathbf{M}^{-1}$ always exists by invertibility of the matrix $\mathbf{M}$. The result follows after replacing $\partial_{\omega}\hat{\gamma}_{\omega}$ by $-\mathbf{D}_{\omega, \mathbf{v}}^{\text{in}}\mathbf{M}^{-1}$ in Equation (10). □

**Lemma C.2 (Estimator of the total functional gradient).** Let $\omega \in \mathbb{R}^d$. Consider the following functional estimator:

$$ \widehat{\nabla \mathcal{F}}(\omega) = \frac{1}{m} \sum_{j=1}^{m} \partial_{\omega} \ell_{\text{out}}(\omega, \hat{h}_{\omega}(\tilde{x}_j), \tilde{y}_j) + \frac{1}{n} \sum_{i=1}^{n} \partial_{\omega, \mathbf{v}}^2 \ell_{\text{in}}(\omega, \hat{h}_{\omega}(x_i), y_i) \hat{a}_{\omega}(x_i). $$

Then, under Assumptions (A) to (E), $\widehat{\nabla \mathcal{F}}(\omega)$ admits the following expression:

$$ \widehat{\nabla \mathcal{F}}(\omega) = \frac{1}{m} \mathbf{D}_{\omega}^{\text{out}} \mathbf{1}_{m} + \frac{1}{n} \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} [\mathbf{K} \quad \mathbf{u}] \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix} \in \mathbb{R}^d, $$

where $\mathbf{D}_{\omega}^{\text{out}}$, $\mathbf{D}_{\omega, \mathbf{v}}^{\text{in}}$, $\mathbf{D}_{\omega}^{\text{out}}$, and $\mathbf{D}_{\omega, \mathbf{v}}^{\text{in}}$ are the same matrices given in Proposition C.1, while $\hat{\alpha}_{\omega} \in \mathbb{R}^n$ and $\hat{\beta}_{\omega} \in \mathbb{R}$ are solutions to the linear system:

$$ \begin{bmatrix} \mathbf{M} \mathbf{K} & \mathbf{M} \mathbf{u} \\ \mathbf{u}^{\top} \mathbf{M} & p \end{bmatrix} \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix} = -\frac{n}{m} \begin{bmatrix} \mathbf{u} \\ v \end{bmatrix}, \quad (11) $$

where the vector $\mathbf{u}$ and matrix $\mathbf{M}$ are the same as in Proposition C.1, while $p$ and $v$ are non-negative scalars.

*Proof.* Let $\omega \in \mathbb{R}^d$. We start by providing an expression of $\hat{a}_{\omega}$ as a linear combination of the kernel evaluated at the inner training points $x_i$, i.e., $K(x_i, \cdot)$, and some element $\xi \in \mathcal{H}$ that we will characterize shortly. From it, we will obtain the expression of $\widehat{\nabla \mathcal{F}}(\omega)$.

**Expression of $\hat{a}_{\omega}$**. Recall that $\hat{a}_{\omega}$ is the unique minimizer of $\hat{L}_{adj}$ in Equation (5), which admits, for any $a \in \mathcal{H}$, the following simple expression by the reproducing property:

$$ \hat{L}_{adj}(\omega, a) = \frac{1}{2n} \sum_{i=1}^{n} \partial_{\omega}^2 \ell_{\text{in}}(\omega, \hat{h}_{\omega}(x_i), y_i) a^2(x_i) \\ + \frac{1}{m} \left\langle a, \underbrace{\sum_{j=1}^{m} \partial_{\omega} \ell_{\text{out}}(\omega, \hat{h}_{\omega}(\tilde{x}_j), \tilde{y}_j) K(\tilde{x}_j, \cdot)}_{\mathcal{H}} \right\rangle + \frac{\lambda}{2} \|a\|_{\mathcal{H}}^2. $$

Hence, by application of the representer theorem, it follows that $\hat{a}_{\omega}$ admits an expression of the form:

$$ \hat{a}_{\omega} = \sum_{i=1}^{n} (\hat{\alpha}_{\omega})_i K(x_i, \cdot) + \hat{\beta}_{\omega} \xi. \quad (12) $$

Therefore, it is possible to recover $\hat{a}_{\omega}$ by minimizing $a \mapsto L_{adj}(\omega, a)$ over the span of $(\xi, K(x_1, \cdot), \dots, K(x_n, \cdot))$. Hence, to find the optimal coefficients $\hat{\alpha}_{\omega} := ((\hat{\alpha}_{\omega})_i)_{1 \le i \le n}$ and $\hat{\beta}_{\omega}$,

<footnote n="2">²In the case where the matrix $\mathbf{K}$ is non-invertible, one needs to restrict $\gamma$ to the orthogonal complement of the null space of $\mathbf{K}$. Such a restriction is valid since the resulting solution $\hat{h}_{\omega}$ will not depend on the component belonging to the null space of $\mathbf{K}$.</footnote>
we first need to express the objective $L_{adj}$ in terms of the coefficients $\alpha \in \mathbb{R}^n$ and $\beta \in \mathbb{R}$ for a given $a^{\alpha, \beta} \in \mathcal{H}$ of the form $a^{\alpha, \beta} = \sum_{i=1}^n (\alpha)_i K(x_i, \cdot) + \beta\xi$. To this end, note that the vector $(\xi(x_1), \dots, \xi(x_n))$ is exactly equal to $\mathbf{u} = \bar{\mathbf{K}}^{\top} \mathbf{D}_{\mathbf{v}}^{\text{out}}$ as defined in Proposition C.1. Moreover, using the reproducing property, we directly have:

$$ (a^{\alpha, \beta}(x_i))_{1 \le i \le n} = [\mathbf{K} \ \mathbf{u}] \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \langle a^{\alpha, \beta}, \xi \rangle_{\mathcal{H}} = [\mathbf{u}^{\top} \ \|\xi\|_{\mathcal{H}}^2] \begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \\ \|a^{\alpha, \beta}\|_{\mathcal{H}}^2 = [\alpha^{\top} \ \beta] \begin{bmatrix} \mathbf{K} & \mathbf{u} \\ \mathbf{u}^{\top} & \|\xi\|_{\mathcal{H}}^2 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix}. $$

We can therefore express the objective $\hat{L}_{adj}$ as follows:

$$ \hat{L}_{adj}(\omega, a^{\alpha, \beta}) = \frac{1}{2n} [\alpha^{\top} \ \beta] \begin{bmatrix} \mathbf{K} \\ \mathbf{u}^{\top} \end{bmatrix} \mathbf{D}_{\mathbf{v}, \mathbf{v}}^{\text{in}} [\mathbf{K} \ \mathbf{u}] \begin{bmatrix} \alpha \\ \beta \end{bmatrix} \\ + \frac{1}{m} [\mathbf{u}^{\top} \ \|\xi\|_{\mathcal{H}}^2] \begin{bmatrix} \alpha \\ \beta \end{bmatrix} + \frac{\lambda}{2} [\alpha^{\top} \ \beta] \begin{bmatrix} \mathbf{K} & \mathbf{u} \\ \mathbf{u}^{\top} & \|\xi\|_{\mathcal{H}}^2 \end{bmatrix} \begin{bmatrix} \alpha \\ \beta \end{bmatrix}. $$

Hence, the optimal coefficients $\hat{\alpha}_{\omega} := ((\hat{\alpha}_{\omega})_i)_{1 \le i \le n}$ and $\hat{\beta}_{\omega}$ are those minimizing the above quadratic form and are characterized by the following optimality condition:

$$ \left[ \underbrace{\overbrace{\left( \mathbf{K} \mathbf{D}_{\mathbf{v},\mathbf{v}}^{\text{in}} + n\lambda \mathbf{1}_{n \times n} \right)}^{\mathbf{M}} \mathbf{K}}_{\mathbf{M}} \underbrace{\overbrace{\left( \mathbf{K} \mathbf{D}_{\mathbf{v},\mathbf{v}}^{\text{in}} + n\lambda \mathbf{1}_{n \times n} \right)}^{\mathbf{M}} \mathbf{u}}_{p \ge 0} \mathbf{u}^{\top} \mathbf{D}_{\mathbf{v},\mathbf{v}}^{\text{in}} \mathbf{u} + n\lambda \|\xi\|_{\mathcal{H}}^2 \right] \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix} = -\frac{n}{m} \left[ \underbrace{\mathbf{u}}_{v \ge 0} \begin{bmatrix} \|\xi\|_{\mathcal{H}}^2 \\ \|\xi\|_{\mathcal{H}}^2 \end{bmatrix} \right]. $$

**Expression of $\widehat{\nabla \mathcal{F}}(\omega)$**. The result follows directly after expressing $\widehat{\nabla \mathcal{F}}(\omega)$ in vector form using the notations $\mathbf{D}_{\omega}^{\text{out}}$ and $\mathbf{D}_{\omega, \mathbf{v}}^{\text{in}}$ from Proposition C.1 and recalling that $(\hat{a}_{\omega}(x_i))_{1 \le i \le n} = [\mathbf{K} \ \mathbf{u}] \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix}$. $\square$

**Proof of Proposition 3.1.** Let $\omega \in \mathbb{R}^d$. Define

$$ \widehat{\nabla \mathcal{F}}(\omega) = \frac{1}{m} \sum_{j=1}^{m} \partial_{\omega} \ell_{\text{out}}(\omega, \hat{h}_{\omega}(\tilde{x}_j), \tilde{y}_j) + \frac{1}{n} \sum_{i=1}^{n} \partial_{\omega, \mathbf{v}}^2 \ell_{\text{in}}(\omega, \hat{h}_{\omega}(x_i), y_i) \hat{a}_{\omega}(x_i), $$

where $\hat{h}_{\omega}$ and $\hat{a}_{\omega}$ are given by Equations (9) and (12). We will show that $\widehat{\nabla \mathcal{F}}(\omega) = \nabla \widehat{\mathcal{F}}(\omega)$. By Proposition C.1 and Lemma C.2, we know that $\nabla \widehat{\mathcal{F}}(\omega)$ and $\widehat{\nabla \mathcal{F}}(\omega)$ admit the following expressions:

$$ \nabla \widehat{\mathcal{F}}(\omega) = \frac{1}{m} \mathbf{D}_{\omega}^{\text{out}} \mathbf{1}_m - \frac{1}{m} \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} \mathbf{M}^{-1} \mathbf{u} \\ \widehat{\nabla \mathcal{F}}(\omega) = \frac{1}{m} \mathbf{D}_{\omega}^{\text{out}} \mathbf{1}_m + \frac{1}{n} \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} [\mathbf{K} \ \mathbf{u}] \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix}. $$

Taking the difference of the two estimators yields:

$$ \widehat{\nabla \mathcal{F}}(\omega) - \nabla \widehat{\mathcal{F}}(\omega) = \frac{1}{m} \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} \left( \mathbf{M}^{-1} \mathbf{u} + \frac{m}{n} [\mathbf{K} \ \mathbf{u}] \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix} \right) \\ = \frac{1}{m} \mathbf{D}_{\omega, \mathbf{v}}^{\text{in}} \mathbf{M}^{-1} \underbrace{\left( \mathbf{u} + \frac{m}{n} \mathbf{M} [\mathbf{K} \ \mathbf{u}] \begin{bmatrix} \hat{\alpha}_{\omega} \\ \hat{\beta}_{\omega} \end{bmatrix} \right)}_{=0}, $$

where the term $\mathbf{u} + \frac{m}{n}(\mathbf{M} \mathbf{K} \hat{\alpha}_{\omega} + \hat{\beta}_{\omega} \mathbf{M} \mathbf{u})$ is equal to 0 by definition of $\hat{\alpha}_{\omega}$ and $\hat{\beta}_{\omega}$ as solutions of the linear system (11) of Lemma C.2. $\square$

# D Preliminary Results

In this section, $\Omega$ is an arbitrary compact subset of $\mathbb{R}^d$ with hull($\Omega$) denoting its convex hull, which is also compact. We also consider an arbitrary fixed positive value $\Lambda$ such that $\lambda \le \Lambda$ as this would allow us to simplify the dependence on $\lambda$ of the boundedness and Lipschitz constants.
### D.1 Boundedness and Lipschitz continuity of $h_{\omega}^*$ and $\hat{h}_{\omega}$

**Proposition D.1 (Boundedness of $h_{\omega}^*$ and $\hat{h}_{\omega}$).** Under Assumptions (A) to (E), the functions $\omega \mapsto \|h_{\omega}^*\|_{\mathcal{H}}$ and $\omega \mapsto \|\hat{h}_{\omega}\|_{\mathcal{H}}$ are bounded over hull($\Omega$) by $\frac{B\sqrt{\kappa}}{\lambda}$, where $B := \sup_{\omega \in \text{hull}(\Omega), y \in \mathcal{Y}} |\partial_v l_{in}(\omega, 0, y)| > 0$. Moreover, for all $\omega \in \text{hull}(\Omega)$ and $x \in \mathcal{X}$, $h_{\omega}^*(x)$ and $\hat{h}_{\omega}(x)$ take value in the compact interval $\mathcal{V} := [-\frac{B\kappa}{\lambda}, \frac{B\kappa}{\lambda}] \subset \mathbb{R}$.

*Proof. Boundedness of $\|h_{\omega}^*\|_{\mathcal{H}}$ and $\|\hat{h}_{\omega}\|_{\mathcal{H}}$.* Let $\omega \in \text{hull}(\Omega)$. Using Lemma H.2, we know, for any $h \in \mathcal{H}$, that:

$$ \|h - h_{\omega}^*\|_{\mathcal{H}} \le \frac{1}{\lambda} \|\partial_h L_{in}(\omega, h)\|_{\mathcal{H}}. $$

This is particularly valid for $h = 0$. Thus,

$$ \|h_{\omega}^*\|_{\mathcal{H}} \le \frac{1}{\lambda} \|\partial_h L_{in}(\omega, 0)\|_{\mathcal{H}}. $$

Using the expression of the partial derivative $\partial_h L_{in}$ established in Proposition B.1, we obtain:

$$ \|h_{\omega}^*\|_{\mathcal{H}} \le \frac{1}{\lambda} \|\mathbb{E}_{\mathbb{P}}[\partial_v l_{in}(\omega, 0, y) K(x, \cdot)]\|_{\mathcal{H}}. $$

By Assumption (B), $K$ is bounded by $\kappa$. Hence, Jensen's inequality yields:

$$ \|h_{\omega}^*\|_{\mathcal{H}} \le \frac{1}{\lambda} \mathbb{E}_{\mathbb{P}} \left[ |\partial_v l_{in}(\omega, 0, y)| \|K(x, \cdot)\|_{\mathcal{H}} \right] \le \frac{\sqrt{\kappa}}{\lambda} \mathbb{E}_{\mathbb{P}} \left[ |\partial_v l_{in}(\omega, 0, y)| \right]. $$

By Assumption (C), $\mathcal{Y}$ is compact, which implies that hull($\Omega$) $\times$ $\mathcal{Y}$ is compact. From Assumption (D), we know that the function $(\omega, y) \mapsto \partial_v l_{in}(\omega, 0, y)$ is continuous. Given that every continuous function on a compact space is bounded, we obtain:

$$ \|h_{\omega}^*\|_{\mathcal{H}} \le \frac{B\sqrt{\kappa}}{\lambda} < +\infty, \quad \text{where } B := \sup_{\omega \in \text{hull}(\Omega), y \in \mathcal{Y}} |\partial_v l_{in}(\omega, 0, y)| > 0. $$

To prove that $\|\hat{h}_{\omega}\|_{\mathcal{H}} \le \frac{B\sqrt{\kappa}}{\lambda}$, we follow a similar approach to that of $\|h_{\omega}^*\|_{\mathcal{H}} \le \frac{B\sqrt{\kappa}}{\lambda}$. More precisely, we investigate the case where the expectation is with respect to the empirical estimate $\hat{\mathbb{P}}_n$ of $\mathbb{P}$.

$h_{\omega}^*(x)$ and $\hat{h}_{\omega}(x)$ belong to $\mathcal{V}$. Let $\omega \in \text{hull}(\Omega)$ and $x \in \mathcal{X}$. By the reproducing property, the Cauchy-Schwarz inequality, and Assumption (B), we have:

$$ |h_{\omega}^*(x)| \le \sqrt{\kappa} \|h_{\omega}^*\| \quad \text{and} \quad |\hat{h}_{\omega}(x)| \le \sqrt{\kappa} \|\hat{h}_{\omega}\|. $$

Using the bound on $\|h_{\omega}^*\|_{\mathcal{H}}$ and $\|\hat{h}_{\omega}\|_{\mathcal{H}}$ already proved in the first part of this proof, we get:

$$ |h_{\omega}^*(x)| \le \frac{B\kappa}{\lambda} \quad \text{and} \quad |\hat{h}_{\omega}(x)| \le \frac{B\kappa}{\lambda}. $$

This concludes the proof. $\square$

**Proposition D.2 (Lipschitz continuity of $\omega \mapsto h_{\omega}^*$).** Under Assumptions (A) to (E), the function $\omega \mapsto h_{\omega}^*$ is $\frac{L\sqrt{\kappa}}{\lambda}$-Lipschitz continuous on hull($\Omega$), where $L := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \|\partial_{\omega,v}^2 l_{in}(\omega, v, y)\| > 0$, and $\mathcal{V}$ is the compact interval introduced in Proposition D.1.

*Proof.* To prove this proposition, we adopt the strategy of finding an upper bound for the Jacobian, which serves as the Lipschitz constant.

Let $\omega \in \text{hull}(\Omega)$. Using Propositions B.1 and B.3, we know that $h \mapsto L_{in}(\omega, h)$ is $\lambda$-strongly convex and Fréchet differentiable. Also, by Proposition B.2, $\partial_h L_{in}$ is Fréchet differentiable on $\mathbb{R}^d \times \mathcal{H}$, and,
a fortiori, Hadamard differentiable. Then, by the functional implicit differentiation theorem [36, 58, Theorem 2.1], the Jacobian $\partial_{\omega}h_{\omega}^{*} : \mathcal{H} \to \mathbb{R}^d$ can be expressed as:

$$ \partial_{\omega}h_{\omega}^{*} = -\partial_{\omega,h}^{2}L_{in}(\omega, h_{\omega}^{*}) (\partial_{h}^{2}L_{in}(\omega, h_{\omega}^{*}))^{-1}. $$

We have:

$$ \begin{aligned} \|\partial_{\omega}h_{\omega}^{*}\|_{\text{op}} &\le \|\partial_{\omega,h}^{2}L_{in}(\omega, h_{\omega}^{*})\|_{\text{op}} \|(\partial_{h}^{2}L_{in}(\omega, h_{\omega}^{*}))^{-1}\|_{\text{op}} \\ &\le \frac{\|\partial_{\omega,h}^{2}L_{in}(\omega, h_{\omega}^{*})\|_{\text{op}}}{\lambda} \\ &= \frac{\|\mathbb{E}_{\mathbb{P}}[\partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)K(x, \cdot)]\|_{\text{op}}}{\lambda} \\ &\le \frac{\|\mathbb{E}_{\mathbb{P}}[\|\partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)\| \|K(x, \cdot)\|_{\mathcal{H}}]}{\lambda} \\ &\le \frac{\sqrt{\kappa}\mathbb{E}_{\mathbb{P}}[\|\partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)\|]}{\lambda}, \end{aligned} \quad (13) $$

where the first line uses the sub-multiplicative property of the operator norm $\| \cdot \|_{\text{op}}$, the second line stems from the fact that $h \mapsto L_{in}(\omega, h)$ is $\lambda$-strongly convex, for any $\omega \in \mathbb{R}^d$, as proved in Proposition B.3, the third line follows from Proposition B.2, the fourth line uses Jensen's inequality, and the last line is a direct consequence of the boundedness of $K$ by $\kappa$ (Assumption (B)). According to Proposition D.1, $h_{\omega}^{*}(x) \in \mathcal{V} := [-\frac{B\kappa}{\lambda}, \frac{B\kappa}{\lambda}]$, which is a compact interval of $\mathbb{R}$, where $B := \sup_{\omega \in \text{hull}(\Omega), y \in \mathcal{Y}} |\partial_{v}\ell_{in}(\omega, 0, y)| > 0$. By Assumption (C), $\mathcal{Y}$ is a compact set, hence $\text{hull}(\Omega) \times \mathcal{V} \times \mathcal{Y}$ is compact. Besides, by Assumption (D), $(\omega, v, y) \mapsto \partial_{v}\ell_{in}(\omega, v, y)$ is continuous over the domain $\text{hull}(\Omega) \times \mathcal{V} \times \mathcal{Y}$. Since every continuous function on a compact set is bounded, this leads to:

$$ \mathbb{E}_{\mathbb{P}}[\|\partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)\|] \le L := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \|\partial_{\omega,v}^{2}\ell_{in}(\omega, v, y)\| < +\infty. $$

Substituting this bound into Equation (13) means that $\frac{L\sqrt{\kappa}}{\lambda}$ is an upper bound on $\|\partial_{\omega}h_{\omega}^{*}\|_{\text{op}}$. Thus, the result follows as desired. $\square$

## D.2 Local boundedness and Lipschitz properties of $\ell_{in}$, $\ell_{out}$, and their derivatives

**Proposition D.3 (Local boundedness).** Under Assumptions (A) to (E), the functions $(\omega, x, y) \mapsto \ell_{out}(\omega, h_{\omega}^{*}(x), y)$, $(\omega, x, y) \mapsto \partial_{\omega}\ell_{out}(\omega, h_{\omega}^{*}(x), y)$, and $(\omega, x, y) \mapsto \partial_{v}\ell_{out}(\omega, h_{\omega}^{*}(x), y)$ are bounded over $\text{hull}(\Omega) \times \mathcal{X} \times \mathcal{Y}$ by some positive constant $M_{out}$. Similarly, the functions $(\omega, x, y) \mapsto \partial_{v}\ell_{in}(\omega, h_{\omega}^{*}(x), y)$, $(\omega, x, y) \mapsto \partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)$, and $(\omega, x, y) \mapsto \partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)$ are bounded over $\text{hull}(\Omega) \times \mathcal{X} \times \mathcal{Y}$ by some positive constant $M_{in}$. The constants $M_{out}$ and $M_{in}$ are defined as:

$$ M_{out} := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \max (|\ell_{out}(\omega, v, y)|, \|\partial_{\omega}\ell_{out}(\omega, v, y)\|, |\partial_{v}\ell_{out}(\omega, v, y)|) > 0, $$

$$ M_{in} := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \max (|\partial_v \ell_{in}(\omega, v, y)|, |\partial_v^2 \ell_{in}(\omega, v, y)|, \|\partial_{\omega,v}^2 \ell_{in}(\omega, v, y)\|) > 0, $$

where $\mathcal{V} \subset \mathbb{R}$ is the compact interval defined in Proposition D.1.

*Proof.* By Proposition D.1, we have that $h_{\omega}^{*}(x) \in \mathcal{V} := [-\frac{B\kappa}{\lambda}, \frac{B\kappa}{\lambda}] \subset \mathbb{R}$, for any $x \in \mathcal{X}$. From Assumption (D), we know that $\ell_{in}$, $\ell_{out}$, and their partial derivatives are all continuous on $\text{hull}(\Omega) \times \mathcal{V} \times \mathcal{Y}$. Also, $\mathcal{Y}$ is compact by Assumption (C). Thus, $\text{hull}(\Omega) \times \mathcal{V} \times \mathcal{Y}$ is compact. As every continuous function defined over a compact space is bounded, we obtain that:

$$ \sup_{\omega \in \text{hull}(\Omega), x \in \mathcal{X}, y \in \mathcal{Y}} |\ell_{out}(\omega, h_{\omega}^{*}(x), y)| \le \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} |\ell_{out}(\omega, v, y)| < +\infty, $$

$$ \sup_{\omega \in \text{hull}(\Omega), x \in \mathcal{X}, y \in \mathcal{Y}} \|\partial_{\bullet}\ell_{\circ}(\omega, h_{\omega}^{*}(x), y)\| \le \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \|\partial_{\bullet}\ell_{\circ}(\omega, v, y)\| < +\infty, $$

where $\bullet \in \{\{v\}, \{w\}, \{w, v\}\}$ and $\circ \in \{in, out\}$. This implies the desired result. $\square$
**Proposition D.4 (Local Lipschitz continuity).** Under Assumptions (A) to (E), there exists a positive constant $Lip_{out}$ so that for any $(x, y)$ in $\mathcal{X} \times \mathcal{Y}$, the functions $\omega \mapsto \ell_{out}(\omega, h_{\omega}^{*}(x), y)$, $\omega \mapsto \partial_{\omega}\ell_{out}(\omega, h_{\omega}^{*}(x), y)$, and $\omega \mapsto \partial_{v}\ell_{out}(\omega, h_{\omega}^{*}(x), y)$ are locally $\frac{\text{Lip}_{out}}{\lambda}$-Lipschitz continuous over hull($\Omega$). Similarly, there exists a positive constant $Lip_{in}$ so that for any $(x, y)$ in $\mathcal{X} \times \mathcal{Y}$, the functions $\omega \mapsto \partial_{v}\ell_{in}(\omega, h_{\omega}^{*}(x), y)$, $\omega \mapsto \partial_{v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)$, and $\omega \mapsto \partial_{\omega,v}^{2}\ell_{in}(\omega, h_{\omega}^{*}(x), y)$ are locally $\frac{\text{Lip}_{in}}{\lambda}$-Lipschitz continuous over hull($\Omega$). The constants $Lip_{out}$ and $Lip_{in}$ are defined, for any $0 < \lambda \le \Lambda$, as:

$$ \text{Lip}_{\text{out}} := (\Lambda + M_{\text{in}}\kappa) \max (M_{\text{out}}, \bar{M}_{\text{out}}) > 0 $$

$$ \text{Lip}_{\text{in}} := (\Lambda + M_{\text{in}}\kappa) \max (M_{\text{in}}, \bar{M}_{\text{in}}) > 0, $$

where:

$$ \bar{M}_{\text{out}} := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \max \left( \left\| \partial_{\omega}^{2} \ell_{\text{out}}(\omega, v, y) \right\|_{\text{op}}, \left\| \partial_{\omega,v}^{2} \ell_{\text{out}}(\omega, v, y) \right\|, \left| \partial_{v}^{2} \ell_{\text{out}}(\omega, v, y) \right| \right) > 0, $$

$$ \bar{M}_{\text{in}} := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \max \left( \left\| \partial_{\omega} \partial_{v}^{2} \ell_{\text{in}}(\omega, v, y) \right\|, \left| \partial_{v}^{3} \ell_{\text{in}}(\omega, v, y) \right|, \left\| \partial_{\omega} \partial_{\omega,v}^{2} \ell_{\text{in}}(\omega, v, y) \right\| \right) > 0, $$

with $M_{\text{out}}$ and $M_{\text{in}}$ being the positive constants defined in Proposition D.3, and $\mathcal{V} \subset \mathbb{R}$ is the compact interval defined in Proposition D.1.

**Proof.** For any $(\omega, x, y) \in \text{hull}(\Omega) \times \mathcal{X} \times \mathcal{Y}$, we have:

$$ \begin{aligned} \| \nabla_{\omega} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \| &= \| \partial_{\omega} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) + \partial_{v} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \| \| \partial_{\omega} h_{\omega}^{*}(x) \| \\ &\le \| \partial_{\omega} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \| + | \partial_{v} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) | \| \partial_{\omega} h_{\omega}^{*} \|_{\text{op}} \| K(x, \cdot) \|_{\mathcal{H}} \\ &\le M_{\text{out}} \left( 1 + \frac{M_{\text{in}}\kappa}{\lambda} \right) \\ &\le \frac{M_{\text{out}} (\Lambda + M_{\text{in}}\kappa)}{\lambda}, \end{aligned} $$

where the first line uses the chain rule, the second line applies the triangle inequality and the reproducing property of the RKHS $\mathcal{H}$, the third line follows from Proposition D.3 to bound the derivatives of $\ell_{out}$, from Proposition D.2, which states that the function $\omega \mapsto h_{\omega}^{*}$ is $\frac{L\sqrt{\kappa}}{\lambda}$-Lipschitz continuous with $L := \sup_{\omega \in \text{hull}(\Omega), v \in \mathcal{V}, y \in \mathcal{Y}} \| \partial_{\omega,v}^{2} \ell_{\text{in}}(\omega, v, y) \| < M_{\text{in}}$, to bound $\| \partial_{\omega} h_{\omega}^{*} \|_{\text{op}}$, and from Assumption (B) to bound $\| K(x, \cdot) \|_{\mathcal{H}}$, and the last line is a direct consequence of $0 < \lambda \le \Lambda$. In a similar way, we obtain:

$$ \begin{aligned} \| \nabla_{\omega} \partial_{\omega} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \|_{\text{op}} &\le \frac{\bar{M}_{\text{out}} (\Lambda + M_{\text{in}}\kappa)}{\lambda}, \quad \| \nabla_{\omega} \partial_{v} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \| \le \frac{\bar{M}_{\text{out}} (\Lambda + M_{\text{in}}\kappa)}{\lambda}, \\ \| \nabla_{\omega} \partial_{v} \ell_{\text{in}}(\omega, h_{\omega}^{*}(x), y) \| &\le \frac{M_{\text{in}} (\Lambda + M_{\text{in}}\kappa)}{\lambda}, \quad \| \nabla_{\omega} \partial_{v}^{2} \ell_{\text{in}}(\omega, h_{\omega}^{*}(x), y) \| \le \frac{\bar{M}_{\text{in}} (\Lambda + M_{\text{in}}\kappa)}{\lambda}, \\ \| \nabla_{\omega} \partial_{\omega,v}^{2} \ell_{\text{in}}(\omega, h_{\omega}^{*}(x), y) \|_{\text{op}} &\le \frac{\bar{M}_{\text{in}} (\Lambda + M_{\text{in}}\kappa)}{\lambda}. \end{aligned} $$

Combining all these bounds concludes the proof. $\square$

# E Generalization Properties

As before, let $\Omega$ be an arbitrary compact subset of $\mathbb{R}^d$.

## E.1 Point-wise estimates

We present a point-wise upper bound on the value error $|\mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega)|$ and gradient error $\|\nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega)\|$. To this end, we introduce the following notation for the error between the inner and outer objectives and their empirical approximations evaluated at the optimal inner solution $h_{\omega}^{*}$:

$$ \delta_{\omega}^{\text{out}} := |L_{\text{out}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{out}}(\omega, h_{\omega}^{*})|, \quad \delta_{\omega}^{\text{in}} := |L_{\text{in}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{in}}(\omega, h_{\omega}^{*})|. $$
By abuse of notation, we introduce the following errors between partial derivatives of $L_{in}$ and $\hat{L}_{in}$ (resp. $L_{out}$ and $\hat{L}_{out}$), evaluated at $(\omega, h_{\omega}^{*})$, i.e.,

$$
\begin{align*}
\partial_h \delta_{\omega}^{\text{out}} &:= \left\| \partial_h L_{\text{out}}(\omega, h_{\omega}^{*}) - \partial_h \hat{L}_{\text{out}}(\omega, h_{\omega}^{*}) \right\|_{\mathcal{H}}, & \partial_{\omega} \delta_{\omega}^{\text{out}} &:= \left\| \partial_{\omega} L_{\text{out}}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{\text{out}}(\omega, h_{\omega}^{*}) \right\|, \\
\partial_h \delta_{\omega}^{\text{in}} &:= \left\| \partial_h L_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_h \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) \right\|_{\mathcal{H}}, & \partial_{\omega} \delta_{\omega}^{\text{in}} &:= \left\| \partial_{\omega} L_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) \right\|, \\
\partial_h^2 \delta_{\omega}^{\text{in}} &:= \left\| \partial_h^2 L_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_h^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) \right\|_{\text{op}}, & \partial_{\omega,h}^2 \delta_{\omega}^{\text{in}} &:= \left\| \partial_{\omega,h}^2 L_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) \right\|_{\text{op}}.
\end{align*}
$$

**Proposition E.1.** Under Assumptions (A) to (E), the following holds for any $\omega \in \Omega$:

$$ \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}} \le \frac{1}{\lambda} \left\| \partial_{h} \hat{L}_{in}(\omega, h_{\omega}^{*}) \right\|_{\mathcal{H}} = \frac{1}{\lambda} \partial_{h} \delta_{\omega}^{in}. $$

*Proof.* Let $\omega \in \Omega$. The function $h \mapsto \hat{L}_{in}(\omega, h)$ is $\lambda$-strongly convex and Fréchet differentiable by Propositions B.1 and B.3. Moreover, $\hat{h}_{\omega}$ is the minimizer of $h \mapsto \hat{L}_{in}(\omega, h)$ by definition. Therefore, using Lemma H.2, we obtain a control on the distance in $\mathcal{H}$ to the optimum $\hat{h}_{\omega}$ of $h \mapsto \hat{L}_{in}(\omega, h)$ in terms of the gradient $\partial_h \hat{L}_{in}(\omega, h)$:

$$ \left\| h - \hat{h}_{\omega} \right\|_{\mathcal{H}} \le \frac{1}{\lambda} \left\| \partial_h \hat{L}_{in}(\omega, h) \right\|_{\mathcal{H}}, \quad \forall h \in \mathcal{H}. $$

In particular, choosing $h = h_{\omega}^{*}$ yields the inequality. The fact that $\left\| \partial_h \hat{L}_{in}(\omega, h_{\omega}^{*}) \right\|_{\mathcal{H}} = \partial_h \delta_{\omega}^{in}$ follows from the optimality of $h_{\omega}^{*}$ which implies that $\partial_h L_{in}(\omega, h_{\omega}^{*}) = 0$. $\square$

**Proposition E.2.** Under Assumptions (A) to (E), the following inequalities hold for any $\omega \in \Omega$:

$$
\begin{align*}
E_{\omega}^{\text{out}} &:= \left| \hat{L}_{\text{out}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{out}}(\omega, \hat{h}_{\omega}) \right| \le C_{\text{out}} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \\
\partial_h E_{\omega}^{\text{out}} &:= \left\| \partial_h \hat{L}_{\text{out}}(\omega, h_{\omega}^{*}) - \partial_h \hat{L}_{\text{out}}(\omega, \hat{h}_{\omega}) \right\|_{\mathcal{H}} \le C_{\text{out}} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \\
\partial_{\omega} E_{\omega}^{\text{out}} &:= \left\| \partial_{\omega} \hat{L}_{\text{out}}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{\text{out}}(\omega, \hat{h}_{\omega}) \right\|_{\mathcal{H}} \le C_{\text{out}} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \\
\partial_h^2 E_{\omega}^{\text{in}} &:= \left\| \partial_h^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_h^2 \hat{L}_{\text{in}}(\omega, \hat{h}_{\omega}) \right\|_{\text{op}} \le C_{\text{in}} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \\
\partial_{\omega,h}^2 E_{\omega}^{\text{in}} &:= \left\| \partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, \hat{h}_{\omega}) \right\|_{\text{op}} \le C_{\text{in}} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}.
\end{align*}
$$

The positive constants $C_{out}$ and $C_{in}$ are defined as:

$$
\begin{align*}
C_{\text{out}} &:= \max \left( M_{\text{out}} \sqrt{\kappa}, \bar{M}_{\text{out}} \kappa, \bar{M}_{\text{out}} \sqrt{\kappa} \right) > 0, \\
C_{\text{in}} &:= \max \left( \bar{M}_{\text{in}} \kappa \sqrt{\kappa}, \bar{M}_{\text{in}} \kappa, M_{\text{in}} \sqrt{d\kappa} \right) > 0,
\end{align*}
$$

where $M_{out}$, $\bar{M}_{out}$, and $\bar{M}_{in}$ are the positive constants defined in Propositions D.3 and D.4.

*Proof.* Lipschitz continuity of some functions of interest. Let $\omega \in \Omega$. According to Proposition D.1, both $h_{\omega}^{*}(x)$ and $\hat{h}_{\omega}(x)$ lie in the compact interval $\mathcal{V} := [-\frac{B\kappa}{\lambda}, \frac{B\kappa}{\lambda}] \subset \mathbb{R}$, for any $x \in \mathcal{X}$, where $B := \sup_{\omega \in \text{hull}(\Omega), y \in \mathcal{Y}} |\partial_v \ell_{in}(\omega, 0, y)| > 0$. By Assumption (C), $\mathcal{Y}$ is a compact set. Hence, $\Omega \times \mathcal{V} \times \mathcal{Y}$ is a compact set as well. Furthermore, by Assumption (D), $(\omega, v, y) \mapsto \ell_{in}(\omega, v, y)$, $(\omega, v, y) \mapsto \ell_{out}(\omega, v, y)$, and their derivatives are all continuous over the compact domain $\Omega \times \mathcal{V} \times \mathcal{Y}$. Therefore, these functions and their derivatives are bounded on this domain. In particular, this also holds when $v$ takes the specific values $h_{\omega}^{*}(x)$ or $\hat{h}_{\omega}(x)$. Let $\bar{v}$ be either $h_{\omega}^{*}(x)$ or $\hat{h}_{\omega}(x)$, for any
$x \in \mathcal{X}$. For any $\omega \in \Omega$ and $y \in \mathcal{Y}$, we have:

$$
\begin{align*}
|\partial_v \ell_{out}(\omega, \bar{v}, y)| &\le \sup_{\omega \in \Omega, v \in \mathcal{V}, y \in \mathcal{Y}} |\partial_v \ell_{out}(\omega, v, y)| \le M_{out} < +\infty, \\
|\partial_v^2 \ell_{out}(\omega, \bar{v}, y)| &\le \sup_{\omega \in \Omega, v \in \mathcal{V}, y \in \mathcal{Y}} |\partial_v^2 \ell_{out}(\omega, v, y)| \le \bar{M}_{out} < +\infty, \\
\|\partial_{\omega,v}^2 \ell_{out}(\omega, \bar{v}, y)\| &\le \sup_{\omega \in \Omega, v \in \mathcal{V}, y \in \mathcal{Y}} \|\partial_{\omega,v}^2 \ell_{out}(\omega, v, y)\| \le \bar{M}_{out} < +\infty, \\
|\partial_v^3 \ell_{in}(\omega, \bar{v}, y)| &\le \sup_{\omega \in \Omega, v \in \mathcal{V}, y \in \mathcal{Y}} |\partial_v^3 \ell_{in}(\omega, v, y)| \le \bar{M}_{in} < +\infty, \\
\|\partial_v \partial_{\omega,v}^2 \ell_{in}(\omega, \bar{v}, y)\|_{\text{op}} &\le \sup_{\omega \in \Omega, v \in \mathcal{V}, y \in \mathcal{Y}} \|\partial_\omega \partial_v^2 \ell_{in}(\omega, v, y)\| \le \bar{M}_{in} < +\infty.
\end{align*}
$$

This means that $v \in \mathcal{V} \mapsto \ell_{out}(\omega, v, y)$, $v \in \mathcal{V} \mapsto \partial_v \ell_{out}(\omega, v, y)$, $v \in \mathcal{V} \mapsto \partial_\omega \ell_{out}(\omega, v, y)$, $v \in \mathcal{V} \mapsto \partial_v^2 \ell_{in}(\omega, v, y)$, and $v \in \mathcal{V} \mapsto \partial_\omega v \ell_{in}(\omega, v, y)$ are Lipschitz continuous, with Lipschitz constants $M_{out}$, $\bar{M}_{out}$, $\bar{M}_{in}$, and $\bar{M}_{in}$, respectively, for any $\omega \in \Omega$ and $y \in \mathcal{Y}$.

**Upper bounds.** We have:

$$
\begin{align*}
E_{\omega}^{\text{out}} := |\hat{L}_{\text{out}}(\omega, h_{\omega}^{*}) - \hat{L}_{\text{out}}(\omega, \hat{h}_{\omega})| &= \left| \frac{1}{m} \sum_{j=1}^{m} \ell_{\text{out}}(\omega, h_{\omega}^{*}(\tilde{x}_j), \tilde{y}_j) - \frac{1}{m} \sum_{j=1}^{m} \ell_{\text{out}}(\omega, \hat{h}_{\omega}(\tilde{x}_j), \tilde{y}_j) \right| \\
&\le \frac{1}{m} \sum_{j=1}^{m} \left| \ell_{\text{out}}(\omega, h_{\omega}^{*}(\tilde{x}_j), \tilde{y}_j) - \ell_{\text{out}}(\omega, \hat{h}_{\omega}(\tilde{x}_j), \tilde{y}_j) \right| \\
&\le \frac{M_{\text{out}}}{m} \sum_{j=1}^{m} \left| h_{\omega}^{*}(\tilde{x}_j) - \hat{h}_{\omega}(\tilde{x}_j) \right| \\
&\le M_{\text{out}} \sqrt{\kappa} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}},
\end{align*}
$$

where the first line uses the definition of $(\omega, h) \mapsto \hat{L}_{out}(\omega, h)$, the second line applies the triangle inequality, the third line leverages the fact that $v \mapsto \ell_{out}(\omega, v, y)$ is $M_{out}$-Lipschitz continuous, for any $\omega \in \Omega$ and $y \in \mathcal{Y}$, and the last line follows from the reproducing property of the RKHS $\mathcal{H}$, Cauchy-Schwarz's inequality, and Assumption (B) to bound $\|K(x, \cdot)\|_{\mathcal{H}}$ by $\sqrt{\kappa}$. Similarly, we obtain:

$$
\begin{align*}
\partial_h E_{\omega}^{\text{out}} \le \bar{M}_{\text{out}} \kappa \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \quad \partial_{\omega} E_{\omega}^{\text{out}} \le \bar{M}_{\text{out}} \sqrt{\kappa} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \\
\partial_h^2 E_{\omega}^{\text{in}} \le \bar{M}_{\text{in}} \kappa \sqrt{\kappa} \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}, \quad \partial_{\omega,h}^2 E_{\omega}^{\text{in}} \le \bar{M}_{\text{in}} \kappa \left\| h_{\omega}^{*} - \hat{h}_{\omega} \right\|_{\mathcal{H}}.
\end{align*}
$$

Combining all the bounds finishes the proof. $\square$

**Proposition E.3.** Under Assumptions (A) to (E), the following inequalities hold for any $\omega \in \Omega$:

$$
\|\partial_h L_{out}(\omega, h_{\omega}^{*})\|_{\mathcal{H}} \le C_{out}, \quad \|\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})\|_{\text{op}} \le C_{in}, \quad \|\partial_{\omega,h}^2 \hat{L}_{in}(\omega, \hat{h}_{\omega})\|_{\text{op}} \le C_{in},
$$

where $C_{out}$ and $C_{in}$ are the positive constants defined in Proposition E.2.

**Proof.** Let $\omega \in \Omega$.

**Upper bound on** $\|\partial_h L_{out}(\omega, h_{\omega}^*)\|_{\mathcal{H}}$. We have:

$$
\begin{align*}
\|\partial_h L_{\text{out}}(\omega, h_{\omega}^*)\|_{\mathcal{H}} &= \left\| \mathbb{E}_{\mathbb{Q}} [\partial_v \ell_{\text{out}}(\omega, h_{\omega}^*(x), y) K(x, \cdot)] \right\|_{\mathcal{H}} \\
&\le \mathbb{E}_{\mathbb{Q}} \left[ |\partial_v \ell_{\text{out}}(\omega, h_{\omega}^*(x), y)| \|K(x, \cdot)\|_{\mathcal{H}} \right] \\
&\le \sqrt{\kappa} \mathbb{E}_{\mathbb{Q}} \left[ |\partial_v \ell_{\text{out}}(\omega, h_{\omega}^*(x), y)| \right],
\end{align*}
$$

where the first line follows from Proposition B.1, the second line results from the triangle inequality, and the last line uses Assumption (B) to bound $\|K(x, \cdot)\|_{\mathcal{H}}$ by $\sqrt{\kappa}$. Furthermore, we know by
Proposition D.1 that $(\omega, h_{\omega}^{*}(x), y)$ belongs to the compact subset $\Omega \times \mathcal{V} \times \mathcal{Y}$, and by Proposition D.3 that $\partial_{\nu}l_{out}(\omega, h_{\omega}^{*}(x), y)$ is bounded by a constant $M_{out}$ on hull($\Omega$) $\times \mathcal{V} \times \mathcal{Y}$. Hence, it follows that:

$$ \|\partial_h L_{out}(\omega, h_{\omega}^{*})\|_{\mathcal{H}} \le \sqrt{\kappa} M_{out} \le C_{out}, $$

where $C_{out}$ is defined in Proposition E.2.

**Upper bound on** $\left\|\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})\right\|_{op}$. According to Proposition B.2, $\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})$ is a Hilbert-Schmidt operator, which points to:

$$ \|\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})\|_{op} \le \|\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})\|_{HS} = \sqrt{\sum_{l=1}^{d} \|\partial_{\omega_l,h}^2 L_{in}(\omega, h_{\omega}^{*})\|_{\mathcal{H}}^2}. \quad (14) $$

This means that to find an upper bound on $\left\|\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})\right\|_{op}$, it suffices to establish an upper bound on $\left\|\partial_{\omega_l,h}^2 L_{in}(\omega, h_{\omega}^{*})\right\|_{\mathcal{H}}^2$ for any $l \in \{1, \dots, d\}$. For a fixed $l \in \{1, \dots, d\}$, we have:

$$ \begin{aligned} \|\partial_{\omega_l,h}^2 L_{in}(\omega, h_{\omega}^{*})\|_{\mathcal{H}}^2 &= \|\mathbb{E}_{\mathbb{P}} [\partial_{\omega_l,v}^2 l_{in}(\omega, h_{\omega}^{*}(x), y) K(x, \cdot)]\|_{\mathcal{H}}^2 \\ &\le \mathbb{E}_{\mathbb{P}} [\|\partial_{\omega_l,v}^2 l_{in}(\omega, h_{\omega}^{*}(x), y)\|^2 \|K(x, \cdot)\|_{\mathcal{H}}^2] \\ &\le \mathbb{E}_{\mathbb{P}} [\|\partial_{\omega,v}^2 l_{in}(\omega, h_{\omega}^{*}(x), y)\|^2] \kappa, \end{aligned} $$

where the first line follows from Proposition B.2, the second line is a consequence of Jensen's inequality applied on the convex function $\|\cdot\|^2$, and the last line applies Assumption (B) to bound $\|K(x, \cdot)\|_{\mathcal{H}}^2$ by $\kappa$. Incorporating this upper bound into Equation (14) yields:

$$ \|\partial_{\omega,h}^2 L_{in}(\omega, h_{\omega}^{*})\|_{op} \le \sqrt{\mathbb{E}_{\mathbb{P}} [\|\partial_{\omega,v}^2 l_{in}(\omega, h_{\omega}^{*}(x), y)\|^2]} d\kappa \le M_{in} \sqrt{d\kappa} \le C_{in}, $$

where we used Proposition D.3 to bound $\partial_{\omega,v}^2 l_{in}(\omega, h_{\omega}^{*}(x), y)$ by the constant $M_{in}$.

**Upper bound on** $\left\|\partial_{\omega,h}^2 \hat{L}_{in}(\omega, \hat{h}_{\omega})\right\|_{op}$. The derivation of this upper bound follows the same steps as the previous one, with the only differences being the use of $\hat{L}_{in}$ instead of $L_{in}$, and $\hat{h}_{\omega}$ instead of $h_{\omega}^{*}$.

Note that in the last step of each of the three upper bounds, we used the fact that the functions we are dealing with are continuous by Assumption (D) on $\Omega \times \mathcal{V} \times \mathcal{Y}$, which is compact because $\Omega$ is compact, $\mathcal{Y}$ is compact by Assumption (C), and $\mathcal{V}$ is a compact interval of $\mathbb{R}$ defined in Proposition D.1. Hence, those functions are bounded. $\square$

**Proposition E.4 (Approximation bounds).** Under Assumptions (A) to (E), the following holds for any $\omega \in \Omega$:

$$ \begin{aligned} & \left|\mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega)\right| \le \delta_{\omega}^{out} + \frac{C_{out}}{\lambda} \partial_h \delta_{\omega}^{in}, \\ & \left\|\nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega)\right\| \le \partial_{\omega} \delta_{\omega}^{out} + \frac{C_{in}}{\lambda} \partial_h \delta_{\omega}^{out} + \frac{C_{out} C_{in}}{\lambda^2} \partial_h^2 \delta_{\omega}^{in} \\ & \qquad + \frac{C_{out}}{\lambda} \partial_{\omega,h}^2 \delta_{\omega}^{in} + \frac{C_{out}}{\lambda} \left(1 + 2 \frac{C_{in}}{\lambda} + \frac{C_{in}^2}{\lambda^2}\right) \partial_h \delta_{\omega}^{in}, \end{aligned} $$

where the constants $C_{in}$ and $C_{out}$ are given in Proposition E.2.

*Proof.* In all what follows, we fix a value for $\omega$ in $\Omega$. We start by controlling the value function, then its gradient.

**Control on the value function.** By the triangle inequality, we have:

$$ \left|\mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega)\right| \le \underbrace{\left|L_{out}(\omega, h_{\omega}^{*}) - \hat{L}_{out}(\omega, h_{\omega}^{*})\right|}_{\delta_{\omega}^{out}} + \underbrace{\left|\hat{L}_{out}(\omega, h_{\omega}^{*}) - \hat{L}_{out}(\omega, \hat{h}_{\omega})\right|}_{E_{\omega}^{out}}. \quad (15) $$
According to Proposition E.2, the error term $E_{\omega}^{out}$ is controlled by the norm of the difference $h_{\omega}^{*} - \hat{h}_{\omega}$, i.e., $E_{\omega}^{out} \le C_{out} \|h_{\omega}^{*} - \hat{h}_{\omega}\|_{\mathcal{H}}$. Moreover, by Proposition E.1, we know that $\|h_{\omega}^{*} - \hat{h}_{\omega}\|_{\mathcal{H}} \le \frac{1}{\lambda} \partial_{h} \delta_{\omega}^{in}$. Therefore, combining both bounds yields $E_{\omega}^{out} \le \frac{C_{out}}{\lambda} \partial_{h} \delta_{\omega}^{in}$. The upper bound on the value function follows by substituting the previous inequality into Equation (15).

**Control on the gradient.** By Proposition B.4, we have the following expression for the total gradient $\nabla \mathcal{F}$:

$$\nabla \mathcal{F}(\omega) = \partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^{2} L_{in}(\omega, h_{\omega}^{*}) (\partial_{h}^{2} L_{in}(\omega, h_{\omega}^{*}))^{-1} \partial_{h} L_{out}(\omega, h_{\omega}^{*}).$$

Similarly, the gradient estimator $\widehat{\nabla \mathcal{F}}$ is defined by replacing $L_{out}$ and $L_{in}$ by their empirical versions $\hat{L}_{out}$ and $\hat{L}_{in}$, and $h_{\omega}^{*}$ by $\hat{h}_{\omega} := \arg \min_{h \in \mathcal{H}} \hat{L}_{in}(\omega, h)$ in the above expression, i.e.,

$$\widehat{\nabla \mathcal{F}}(\omega) = \partial_{\omega} \hat{L}_{out}(\omega, \hat{h}_{\omega}) - \partial_{\omega,h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega}) (\partial_{h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega}))^{-1} \partial_{h} \hat{L}_{out}(\omega, \hat{h}_{\omega}).$$

To simplify notations, for any $h \in \mathcal{H}$, we introduce the following operators $R(h), \hat{R}(h) : \mathcal{H} \to \Omega$:

$$R(h) = \partial_{\omega,h}^{2} L_{in}(\omega, h) (\partial_{h}^{2} L_{in}(\omega, h))^{-1} \quad \text{and} \quad \hat{R}(h) = \partial_{\omega,h}^{2} \hat{L}_{in}(\omega, h) (\partial_{h}^{2} \hat{L}_{in}(\omega, h))^{-1}.$$

The difference $\nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega)$ can be decomposed as:

$$\begin{aligned} \nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega) & = (\partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{out}(\omega, h_{\omega}^{*})) + (\partial_{\omega} \hat{L}_{out}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{out}(\omega, \hat{h}_{\omega})) \\ & \quad - \hat{R}(\hat{h}_{\omega}) ( (\partial_{h} L_{out}(\omega, h_{\omega}^{*}) - \partial_{h} \hat{L}_{out}(\omega, h_{\omega}^{*})) + (\partial_{h} \hat{L}_{out}(\omega, h_{\omega}^{*}) - \partial_{h} \hat{L}_{out}(\omega, \hat{h}_{\omega})) ) \\ & \quad - (R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})) \partial_{h} L_{out}(\omega, h_{\omega}^{*}). \end{aligned}$$

By taking the norm of the above equality and using the triangle inequality, we obtain the following upper bound:

$$\begin{aligned} & \|\nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega)\| \\ \le & \underbrace{\|\partial_{\omega} L_{out}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{out}(\omega, h_{\omega}^{*})\|}_{\partial_{\omega} \delta_{\omega}^{out}} + \underbrace{\|\partial_{\omega} \hat{L}_{out}(\omega, h_{\omega}^{*}) - \partial_{\omega} \hat{L}_{out}(\omega, \hat{h}_{\omega})\|}_{\partial_{\omega} E_{\omega}^{out}} \\ & + \|\hat{R}(\hat{h}_{\omega})\|_{\text{op}} \left( \underbrace{\|\partial_{h} L_{out}(\omega, h_{\omega}^{*}) - \partial_{h} \hat{L}_{out}(\omega, h_{\omega}^{*})\|_{\mathcal{H}}}_{\partial_{h} \delta_{\omega}^{out}} + \underbrace{\|\partial_{h} \hat{L}_{out}(\omega, h_{\omega}^{*}) - \partial_{h} \hat{L}_{out}(\omega, \hat{h}_{\omega})\|_{\mathcal{H}}}_{\partial_{h} E_{\omega}^{out}} \right) \\ & + \|R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})\|_{\text{op}} \|\partial_{h} L_{out}(\omega, h_{\omega}^{*})\|_{\mathcal{H}}. \end{aligned} \quad (16)$$

Next, we provide upper bounds on $\|R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})\|_{\text{op}}$ and $\|\hat{R}(\hat{h}_{\omega})\|_{\text{op}}$ in terms of derivatives of $L_{in}$ and $\hat{L}_{in}$.

**Upper bounds on** $\|R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})\|_{\text{op}}$ **and** $\|\hat{R}(\hat{h}_{\omega})\|_{\text{op}}$. By application of Propositions B.2 and B.3, we deduce that $\partial_{\omega,h}^{2} L_{in}(\omega, h_{\omega}^{*})$, $\partial_{h}^{2} L_{in}(\omega, h_{\omega}^{*})$, $\partial_{\omega,h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega})$, and $\partial_{h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega})$ are all bounded operators. Moreover, since $L_{in}$ and $\hat{L}_{in}$ are $\lambda$-strongly convex in their second argument by Proposition B.3, it follows that $\partial_{h}^{2} L_{in}(\omega, h_{\omega}^{*}) \ge \lambda \text{Id}_{\mathcal{H}}$ and $\partial_{h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega}) \ge \lambda \text{Id}_{\mathcal{H}}$. We can therefore apply Lemma H.1 which yields the following inequalities:

$$\begin{aligned} \|R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})\|_{\text{op}} & \le \frac{1}{\lambda^2} \|\partial_{\omega,h}^{2} L_{in}(\omega, h_{\omega}^{*})\|_{\text{op}} \|\partial_{h}^{2} L_{in}(\omega, h_{\omega}^{*}) - \partial_{h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega})\|_{\text{op}} \\ & \quad + \frac{1}{\lambda} \|\partial_{\omega,h}^{2} L_{in}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega})\|_{\text{op}}, \\ \|\hat{R}(\hat{h}_{\omega})\|_{\text{op}} & \le \frac{1}{\lambda} \|\partial_{\omega,h}^{2} \hat{L}_{in}(\omega, \hat{h}_{\omega})\|_{\text{op}}. \end{aligned}$$
By applying the triangle inequality to both terms of the first inequality above, we obtain:

$$
\begin{aligned}
& \|R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})\|_{\text{op}} \\
\le & \frac{1}{\lambda^2} \|\partial_{\omega,h}^2 L_{\text{in}}(\omega, h_{\omega}^{*})\|_{\text{op}} \left( \underbrace{\left\|\partial_h^2 L_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_h^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*})\right\|_{\text{op}}}_{\partial_h^2 \delta_{\omega}^{\text{in}}} + \underbrace{\left\|\partial_h^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_h^2 L_{\text{in}}(\omega, \hat{h}_{\omega})\right\|_{\text{op}}}_{\partial_h^2 E_{\omega}^{\text{in}}} \right) \\
& + \frac{1}{\lambda} \left( \underbrace{\left\|\partial_{\omega,h}^2 L_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*})\right\|_{\text{op}}}_{\partial_{\omega,h}^2 \delta_{\omega}^{\text{in}}} + \underbrace{\left\|\partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, h_{\omega}^{*}) - \partial_{\omega,h}^2 L_{\text{in}}(\omega, \hat{h}_{\omega})\right\|_{\text{op}}}_{\partial_{\omega,h}^2 E_{\omega}^{\text{in}}} \right).
\end{aligned}
$$

**Final bound.** We can now substitute the above bounds on $\|R(h_{\omega}^{*}) - \hat{R}(\hat{h}_{\omega})\|_{\text{op}}$ and $\|\hat{R}(\hat{h}_{\omega})\|_{\text{op}}$ into Equation (16) to obtain the following upper bound on the gradient error:

$$
\begin{aligned}
& \|\nabla \mathcal{F}(\omega) - \overline{\nabla \mathcal{F}}(\omega)\| \\
\le & \partial_{\omega} \delta_{\omega}^{\text{out}} + \partial_{\omega} E_{\omega}^{\text{out}} + \frac{1}{\lambda} \left\| \partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, \hat{h}_{\omega}) \right\|_{\text{op}} (\partial_h \delta_{\omega}^{\text{out}} + \partial_h E_{\omega}^{\text{out}}) \\
& + \|\partial_h L_{\text{out}}(\omega, h_{\omega}^{*})\|_{\mathcal{H}} \left( \frac{1}{\lambda^2} \|\partial_{\omega,h}^2 L_{\text{in}}(\omega, h_{\omega}^{*})\|_{\text{op}} (\partial_h^2 \delta_{\omega}^{\text{in}} + \partial_h^2 E_{\omega}^{\text{in}}) + \frac{1}{\lambda} (\partial_{\omega,h}^2 \delta_{\omega}^{\text{in}} + \partial_{\omega,h}^2 E_{\omega}^{\text{in}}) \right).
\end{aligned}
\tag{17}
$$

Furthermore, by Proposition E.3, we have the following upper bounds on the derivatives of $L_{\text{in}}$ and $L_{\text{out}}$:

$$
\|\partial_h L_{\text{out}}(\omega, h_{\omega}^*)\|_{\mathcal{H}} \le C_{\text{out}}, \quad \|\partial_{\omega,h}^2 L_{\text{in}}(\omega, h_{\omega}^*)\|_{\text{op}} \le C_{\text{in}}, \quad \|\partial_{\omega,h}^2 \hat{L}_{\text{in}}(\omega, \hat{h}_{\omega})\|_{\text{op}} \le C_{\text{in}}.
$$

Incorporating the above bounds into Equation (17), we further get:

$$
\begin{aligned}
\|\nabla \mathcal{F}(\omega) - \overline{\nabla \mathcal{F}}(\omega)\| &\le \partial_{\omega} \delta_{\omega}^{\text{out}} + \partial_{\omega} E_{\omega}^{\text{out}} + \frac{C_{\text{in}}}{\lambda} (\partial_{\omega} \delta_{\omega}^{\text{out}} + \partial_{\omega} E_{\omega}^{\text{out}}) \\
&\quad + C_{\text{out}} \left( \frac{C_{\text{in}}}{\lambda^2} (\partial_{\omega}^2 \delta_{\omega}^{\text{in}} + \partial_{\omega}^2 E_{\omega}^{\text{in}}) + \frac{1}{\lambda} (\partial_{\omega,h}^2 \delta_{\omega}^{\text{in}} + \partial_{\omega,h}^2 E_{\omega}^{\text{in}}) \right).
\end{aligned}
$$

By Proposition E.2, we can upper-bound the error terms $\partial_{\omega} E_{\omega}^{\text{out}}$ and $\partial_{\omega} E_{\omega}^{\text{out}}$ by $C_{\text{out}} \|h_{\omega}^{*} - \hat{h}_{\omega}\|_{\mathcal{H}}$, and $\partial_{\omega}^2 E_{\omega}^{\text{in}}$ and $\partial_{\omega,h}^2 E_{\omega}^{\text{in}}$ by $C_{\text{in}} \|h_{\omega}^{*} - \hat{h}_{\omega}\|_{\mathcal{H}}$. Furthermore, since $\|h_{\omega}^{*} - \hat{h}_{\omega}\|_{\mathcal{H}} \le \frac{1}{\lambda} \partial_{\omega} \delta_{\omega}^{\text{in}}$ by Proposition E.1, we can further show that the gradient error satisfies the desired bound:

$$
\begin{aligned}
\|\nabla \mathcal{F}(\omega) - \overline{\nabla \mathcal{F}}(\omega)\| &\le \partial_{\omega} \delta_{\omega}^{\text{out}} + \frac{C_{\text{in}}}{\lambda} \partial_{\omega} \delta_{\omega}^{\text{out}} + C_{\text{out}} \left( \frac{C_{\text{in}}}{\lambda^2} \partial_{\omega}^2 \delta_{\omega}^{\text{in}} + \frac{1}{\lambda} \partial_{\omega,h}^2 \delta_{\omega}^{\text{in}} \right) \\
&\quad + \frac{C_{\text{out}}}{\lambda} \left( 1 + 2 \frac{C_{\text{in}}}{\lambda} + \frac{C_{\text{in}}^2}{\lambda^2} \right) \partial_{\omega} \delta_{\omega}^{\text{in}}.
\end{aligned}
$$

□

## E.2 Maximal inequalities

**Proposition E.5 (Maximal inequalities for empirical processes).** Let $\Lambda$ be a positive constant. Under Assumptions (A) to (E), the following maximal inequalities hold for any $0 < \lambda \le \Lambda$:

$$
\mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \delta_{\omega}^{\text{out}} \right] \le \sqrt{\frac{1}{\lambda^2 m}} c(\Omega) \max(M_{\text{out}} \text{Lip}_{\text{out}} \text{diam}(\Omega), \Lambda M_{\text{out}}^2),
$$

$$
\mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_{\omega} \delta_{\omega}^{\text{out}} \right] \le \sqrt{\frac{d}{\lambda^2 m}} c(\Omega) \max(M_{\text{out}} \text{Lip}_{\text{out}} \text{diam}(\Omega), \Lambda M_{\text{out}}^2),
$$

where $c(\Omega)$ is a positive constant greater than 1 that depends only on $\Omega$ and $d$, while $\text{Lip}_{\text{out}}$ and $M_{\text{out}}$ are positive constants defined in Propositions D.3 and D.4.
**Proof.** We will apply the result of Proposition F.3 which provides maximal inequalities for real-valued empirical processes that are uniformly bounded and Lipschitz in their parameter. To this end, consider the parametric families:

$$ \mathcal{T}_l^{\text{out}} := \{\mathcal{X} \times \mathcal{Y} \ni (x, y) \mapsto \partial_{w_l} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \mid \omega \in \Omega\}, \quad 1 \le l \le d $$

$$ \mathcal{T}_0^{\text{out}} := \{\mathcal{X} \times \mathcal{Y} \ni (x, y) \mapsto \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y) \mid \omega \in \Omega\}. $$

For any $0 \le l \le d$, these real-valued functions are uniformly bounded by a positive constant $M_{\text{out}}$, thanks to Proposition D.3. Moreover, by Proposition D.4, the functions $\omega \mapsto \partial_{w_l} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y)$ and $\omega \mapsto \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y)$ are all $\lambda^{-1} \text{ Lip}_{\text{out}}$-Lipschitz for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$. Hence, Proposition F.3 is applicable to each of these families, with $\mathbb{D}$ set to $\mathbb{Q}$ and $\mathbb{Z}$ set to $\mathcal{X} \times \mathcal{Y}$. We treat both $\delta_{\omega}^{\text{out}}$ and $\partial_{\omega} \delta_{\omega}^{\text{out}}$ separately.

**A maximal inequality for $\delta_{\omega}^{\text{out}}$.** For $l = 0$, we readily apply Proposition F.3 with $p = 1$ to get the following maximal inequality for $\delta_{\omega}^{\text{out}}$:

$$ \begin{aligned} \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \delta_{\omega}^{\text{out}} \right] &:= \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \left| \mathbb{E}_{(x,y) \sim \mathbb{Q}} [\ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y)] - \frac{1}{m} \sum_{j=1}^{m} \ell_{\text{out}}(\omega, h_{\omega}^{*}(\tilde{x}_j), \tilde{y}_j) \right| \right] \\ &\le \sqrt{\frac{1}{\lambda^2 m}} c(\Omega) \max(M_{\text{out}} \text{Lip}_{\text{out}} \text{diam}(\Omega), \Lambda M_{\text{out}}^2). \end{aligned} $$

**A maximal inequality for $\partial_{\omega} \delta_{\omega}^{\text{out}}$.** We now turn to $\partial_{\omega} \delta_{\omega}^{\text{out}}$, which involves vector-valued processes (as an error between the gradient and its estimate). While the maximal inequalities in Proposition F.3 hold for real-valued processes, we will first obtain maximal inequalities for each component appearing in $\partial_{\omega} \delta_{\omega}^{\text{out}}$ and then sum these to control $\partial_{\omega} \delta_{\omega}^{\text{out}}$. To this end, we first use the Cauchy-Schwarz inequality which implies that $\mathbb{E}_{\mathbb{Q}}[\sup_{\omega \in \Omega} \partial_{\omega} \delta_{\omega}^{\text{out}}] \le \mathbb{E}_{\mathbb{Q}}[\sup_{\omega \in \Omega} (\partial_{\omega} \delta_{\omega}^{\text{out}})^2]^{\frac{1}{2}}$. Thus we only need to control $\mathbb{E}_{\mathbb{Q}}[\sup_{\omega \in \Omega} (\partial_{\omega} \delta_{\omega}^{\text{out}})^2]$. Simple calculations show that:

$$ \begin{aligned} & \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_{\omega} \delta_{\omega}^{\text{out}} \right]^2 \\ & \le \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} (\partial_{\omega} \delta_{\omega}^{\text{out}})^2 \right] \\ & \le \sum_{l=1}^{d} \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \left| \mathbb{E}_{(x,y) \sim \mathbb{Q}} [\partial_{w_l} \ell_{\text{out}}(\omega, h_{\omega}^{*}(x), y)] - \frac{1}{m} \sum_{j=1}^{m} \partial_{w_l} \ell_{\text{out}}(\omega, h_{\omega}^{*}(\tilde{x}_j), \tilde{y}_j) \right| \right]^2 \\ & \le \left( \sqrt{\frac{d}{\lambda^2 m}} c(\Omega) \max(M_{\text{out}} \text{Lip}_{\text{out}} \text{diam}(\Omega), \Lambda M_{\text{out}}^2) \right)^2, \end{aligned} $$

where the last inequality follows by application of Proposition F.3 with $p = 2$ to each term in the right-hand side of the first inequality for $1 \le l \le d$. We get the desired bound on $\mathbb{E}_{\mathbb{Q}}[\sup_{\omega \in \Omega} \partial_{\omega} \delta_{\omega}^{\text{out}}]$ by taking the square root of the above inequality. $\square$

**Proposition E.6 (Maximal inequalities for RKHS-valued empirical processes).** Let $\Lambda$ be a positive constant. Under Assumptions (A) to (E), the following maximal inequalities hold for any $0 < \lambda \le \Lambda$:

$$ \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_h \delta_{\omega}^{\text{out}} \right] \le \lambda^{-\frac{1}{4}} m^{-\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{\text{out},1} \tilde{L}_{\text{out},1} \text{diam}(\Omega), \Lambda \tilde{M}_{\text{out},1}^2 \right) \right)^{\frac{1}{4}}, $$

$$ \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_h \delta_{\omega}^{\text{in}} \right] \le \lambda^{-\frac{1}{4}} n^{-\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{\text{in},1} \tilde{L}_{\text{in},1} \text{diam}(\Omega), \Lambda \tilde{M}_{\text{in},1}^2 \right) \right)^{\frac{1}{4}}, $$

$$ \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{\omega,h}^2 \delta_{\omega}^{\text{in}} \right] \le \lambda^{-\frac{1}{4}} n^{-\frac{1}{2}} d^{\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{\text{in},1} \tilde{L}_{\text{in},1} \text{diam}(\Omega), \Lambda \tilde{M}_{\text{in},1}^2 \right) \right)^{\frac{1}{4}}, $$

$$ \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_h^2 \delta_{\omega}^{\text{in}} \right] \le \lambda^{-\frac{1}{4}} n^{-\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{\text{in},2}^2 \tilde{L}_{\text{in},2} \text{diam}(\Omega), \Lambda \tilde{M}_{\text{in},2}^2 \right) \right)^{\frac{1}{4}}, $$
where $c(\Omega)$ is a positive constant greater than 1 that depends only on $\Omega$ and $d$, $\tilde{L}_{out,1}$, $\tilde{L}_{in,1}$, $\tilde{L}_{in,2}$, $\tilde{M}_{out,1}$, $\tilde{M}_{in,1}$, and $\tilde{M}_{in,2}$ are positive constants defined as:

$$ \tilde{L}_{out,1} := 2 \mathrm{Lip}_{out} M_{out} \kappa, \quad \tilde{L}_{in,1} := 2 \mathrm{Lip}_{in} M_{in} \kappa, \quad \tilde{L}_{in,2} := 2 \mathrm{Lip}_{in} M_{in} \kappa^2, \\ \tilde{M}_{out,1} := M_{out}^2 \kappa, \quad \tilde{M}_{in,1} := M_{in}^2 \kappa, \quad \tilde{M}_{in,2} := M_{in}^2 \kappa^2, $$

and $\mathrm{Lip}_{out}$, $\mathrm{Lip}_{in}$, $M_{out}$, and $M_{in}$ are positive constants given in Propositions D.3 and D.4.

**Proof.** Consider parametric families of real-valued functions indexed by $\Omega$ of the form:

$$ \mathcal{T}_{s,a} := \{t_{\omega} : ((x, y), (x', y')) \mapsto f_s(\omega, x, y) f_s(\omega, x', y') K^a(x, x') \mid \omega \in \Omega\}, $$

where $a \in \{1, 2\}$, $s$ is an integer satisfying $0 \le s \le d+2$, and $f_s(\omega, x, y)$ are real-valued functions given by:

$$ \begin{aligned} f_0 : (\omega, x, y) &\mapsto \partial_v \ell_{out}(\omega, h_{\omega}^{*}(x), y), & f_1 : (\omega, x, y) &\mapsto \partial_v \ell_{in}(\omega, h_{\omega}^{*}(x), y), \\ f_2 : (\omega, x, y) &\mapsto \partial_v^2 \ell_{in}(\omega, h_{\omega}^{*}(x), y), & f_{2+l} : (\omega, x, y) &\mapsto \partial_{\omega_l, v}^2 \ell_{in}(\omega, h_{\omega}^{*}(x), y), & 1 \le l \le d. \end{aligned} $$

For any $1 \le s \le d+2$, the real-valued functions $f_s$ are uniformly bounded by a positive constant $M_{in}$ thanks to Proposition D.3. Moreover, since the kernel $K$ is bounded by $\kappa$ due to Assumption (B), it follows that all elements $t_{\omega}$ of $\mathcal{T}_{s,a}$ are uniformly bounded by $\tilde{M}_{in,a} := M_{in}^2 \kappa^a$. Moreover, for $1 \le s \le d+2$, the functions $\omega \mapsto f_s(\omega, x, y)$ are $\lambda^{-1} \mathrm{Lip}_{in}$-Lipschitz for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$ by Proposition D.4. Hence, it follows that the map $\omega \mapsto t_{\omega}((x, y), (x', y'))$ is $\lambda^{-1} \tilde{L}_{in,a}$-Lipschitz with $\tilde{L}_{in,a} := 2 \mathrm{Lip}_{in} M_{in} \kappa^a$ for any $(x, y)$ and $(x', y')$ in $\mathcal{X} \times \mathcal{Y}$. Similarly, for $s=0$, we get the same properties, albeit, with different constants, i.e., the family $\mathcal{T}_{0,a}$ is uniformly bounded by a constant $\tilde{M}_{out,a} := M_{out}^2 \kappa^a$ with $M_{out}$ introduced in Proposition D.3, and is $\lambda^{-1} \tilde{L}_{out,a}$-Lipschitz in its parameter with $\tilde{L}_{out,a} := 2 \mathrm{Lip}_{out} M_{out} \kappa^a$ where $\mathrm{Lip}_{out}$ is given in Proposition D.4. Hence, the maximal inequality in Proposition F.4 is applicable to each of these families with $\mathcal{Z}$ set to $\mathcal{X} \times \mathcal{Y}$, and $\mathcal{D}$ set either to $\mathbb{P}$ for $1 \le s \le d+2$, or to $\mathbb{Q}$ for $s=0$. For conciseness, in all what follows, we will write $z = (x, y)$ and $z_i = (x_i, y_i)$ and $\tilde{z}_j = (\tilde{x}_j, \tilde{y}_j)$ for $1 \le i \le n$ and $1 \le j \le m$.

**Maximal inequalities for $\partial_h \delta_{\omega}^{out}$ and $\partial_h \delta_{\omega}^{in}$**. We control $\partial_h \delta_{\omega}^{out}$ first as $\partial_h \delta_{\omega}^{in}$ will be dealt with similarly. Using Cauchy-Schwarz inequality and standard calculus, we have that:

$$ \begin{aligned} & \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_h \delta_{\omega}^{out} \right]^2 \\ & \le \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} (\partial_h \delta_{\omega}^{out})^2 \right] \\ & := \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \left\| \mathbb{E}_{(x,y) \sim \mathbb{Q}} [\partial_v \ell_{out}(\omega, h_{\omega}^{*}(x), y) K(x, \cdot)] - \frac{1}{m} \sum_{j=1}^{m} \partial_v \ell_{out}(\omega, h_{\omega}^{*}(\tilde{x}_j), \tilde{y}_j) K(\tilde{x}_j, \cdot) \right\|_{\mathcal{H}}^2 \right] \\ & = \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \mathbb{E}_{z,z' \sim \mathbb{Q} \otimes \mathbb{Q}} [t_{\omega}(z, z')] + \frac{1}{m^2} \sum_{i,j=1}^{m} t_{\omega}(z_i, z_j) - \frac{2}{m} \sum_{j=1}^{m} \mathbb{E}_{z \sim \mathbb{Q}} [t_{\omega}(z, \tilde{z}_j)] \right], \end{aligned} $$

where $t_{\omega}(z, z') := \partial_v \ell_{out}(\omega, h_{\omega}^{*}(x), y) \partial_v \ell_{out}(\omega, h_{\omega}^{*}(x'), y') K(x, x') \in \mathcal{T}_{0,1}$. The last term is precisely what Proposition F.4 controls when applying it to the family $\mathcal{T}_{0,1}$ and choosing $\mathcal{D}$ to be $\mathbb{Q}$. Therefore, the following maximal inequality holds by application of Proposition F.4:

$$ \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_h \delta_{\omega}^{out} \right] \le \lambda^{-\frac{1}{4}} m^{-\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{out,1} \tilde{L}_{out,1} \operatorname{diam}(\Omega), \Lambda \tilde{M}_{out,1}^2 \right) \right)^{\frac{1}{4}}, $$

where $c(\Omega)$ is a positive constant greater than 1 that depends only on $\Omega$ and $d$. We obtain a similar inequality for $\partial_h \delta_{\omega}^{in}$ by carrying out similar calculations, then applying Proposition F.4 to the family $\mathcal{T}_{1,1}$ and choosing $\mathbb{P}$ for the probability distribution $\mathcal{D}$. The resulting bound is then of the form:

$$ \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_h \delta_{\omega}^{in} \right] \le \lambda^{-\frac{1}{4}} n^{-\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{in,1} \tilde{L}_{in,1} \operatorname{diam}(\Omega), \Lambda \tilde{M}_{in,1}^2 \right) \right)^{\frac{1}{4}}. $$
A maximal inequality for $\partial_{\omega,h}^2 \delta_{\omega}^{in}$. We have:

$$
\begin{align*}
& \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{\omega,h}^2 \delta_{\omega}^{in} \right]^2 \\
\le {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} (\partial_{\omega,h}^2 \delta_{\omega}^{in})^2 \right] \\
:= {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \left\| \mathbb{E}_{(x,y) \sim \mathbb{P}} [\partial_{\omega,v}^2 \ell_{in}(\omega, h_{\omega}^*(x), y) K(x, \cdot)] - \frac{1}{n} \sum_{i=1}^n \partial_{\omega,v}^2 \ell_{in}(\omega, h_{\omega}^*(x_i), y_i) K(x_i, \cdot) \right\|_{\text{op}}^2 \right] \\
\le {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \left\| \mathbb{E}_{(x,y) \sim \mathbb{P}} [\partial_{\omega,v}^2 \ell_{in}(\omega, h_{\omega}^*(x), y) K(x, \cdot)] - \frac{1}{n} \sum_{i=1}^n \partial_{\omega,v}^2 \ell_{in}(\omega, h_{\omega}^*(x_i), y_i) K(x_i, \cdot) \right\|_{\text{HS}}^2 \right] \\
:= {} & \sum_{l=1}^d \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \left\| \mathbb{E}_{(x,y) \sim \mathbb{P}} [\partial_{\omega,l,v}^2 \ell_{in}(\omega, h_{\omega}^*(x), y) K(x, \cdot)] - \frac{1}{n} \sum_{i=1}^n \partial_{\omega,l,v}^2 \ell_{in}(\omega, h_{\omega}^*(x_i), y_i) K(x_i, \cdot) \right\|_{\mathcal{H}}^2 \right] \\
:= {} & \sum_{l=1}^d \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \mathbb{E}_{z,z' \sim \mathbb{P} \otimes \mathbb{P}} [t_{\omega,l}(z, z')] + \frac{1}{n^2} \sum_{i,j=1}^n t_{\omega,l}(z_i, z_j) - \frac{2}{n} \sum_{i=1}^n \mathbb{E}_{z \sim \mathbb{P}} [t_{\omega,l}(z, z_i)] \right],
\end{align*}
$$

where we introduced $t_{\omega,l}(z, z') := \partial_{\omega,l,v}^2 \ell_{in}(\omega, h_{\omega}^*(x), y) \partial_{\omega,l,v}^2 \ell_{in}(\omega, h_{\omega}^*(x'), y') K(x, x') \in \mathcal{T}_{2+l,1}$.
Here, (a) follows from the Cauchy-Schwarz inequality, (b) is obtained by definition of $\partial_{\omega,h}^2 \delta_{\omega}^{in}$, while
(c) uses the general fact that the operator norm of an operator is upper-bounded by its Hilbert-Schmidt
norm which is finite in our case by application of Proposition B.2. Moreover, (d) further uses the
Hilbert-Schmidt norm of an operator in terms of the norm of its rows, while (e) simply expands the
squared RKHS norm and uses the reproducing property in the RKHS $\mathcal{H}$. Each term in the last item
(e) is precisely what Proposition F.4 controls when applying it to the families $\mathcal{T}_{2+l,1}$ for $1 \le l \le d$
and choosing $\mathbb{D}$ to be $\mathbb{P}$. Therefore, the following maximal inequality holds by a direct application of
Proposition F.4:

$$
\mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{\omega,h}^2 \delta_{\omega}^{in} \right] \le \lambda^{-\frac{1}{4}} n^{-\frac{1}{2}} d^{\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{in,1}, \tilde{L}_{in,1} \operatorname{diam}(\Omega), \Lambda \tilde{M}_{in,1}^2 \right) \right)^{\frac{1}{4}},
$$

where $c(\Omega)$ is a positive constant greater than 1 that depends only on $\Omega$ and $d$.

A maximal inequality for $\partial_h^2 \delta_\omega^{in}$. We will use a similar approach as for $\partial_{\omega,h}^2 \delta_\omega^{in}$. We have:

$$
\begin{align*}
& \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_h^2 \delta_{\omega}^{in} \right]^2 \\
\le {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} (\partial_h^2 \delta_{\omega}^{in})^2 \right] \\
:= {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \left\| \mathbb{E}_{(x,y) \sim \mathbb{P}} [\partial_v^2 \ell_{in}(\omega, h_{\omega}^*(x), y) K(x, \cdot) \otimes K(x, \cdot)] \right. \right. \\
& \qquad \left. - \frac{1}{n} \sum_{i=1}^n \partial_v^2 \ell_{in}(\omega, h_{\omega}^*(x_i), y_i) K(x_i, \cdot) \otimes K(x_i, \cdot) \right\|_{\text{op}}^2 \\
\le {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \left\| \mathbb{E}_{(x,y) \sim \mathbb{P}} [\partial_v^2 \ell_{in}(\omega, h_{\omega}^*(x), y) K(x, \cdot) \otimes K(x, \cdot)] \right. \right. \\
& \qquad \left. - \frac{1}{n} \sum_{i=1}^n \partial_v^2 \ell_{in}(\omega, h_{\omega}^*(x_i), y_i) K(x_i, \cdot) \otimes K(x_i, \cdot) \right\|_{\text{HS}}^2 \\
:= {} & \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \mathbb{E}_{z,z' \sim \mathbb{P} \otimes \mathbb{P}} [t_{\omega}(z, z')] + \frac{1}{n^2} \sum_{i,j=1}^n t_{\omega}(z_i, z_j) - \frac{2}{n} \sum_{i=1}^n \mathbb{E}_{z \sim \mathbb{P}} [t_{\omega}(z, z_i)] \right],
\end{align*}
$$
where we introduced $t_{\omega}(z, z') := \partial_{v}^{2}\ell_{in}(\omega, x, y)\partial_{v}^{2}\ell_{in}(\omega, x', y')K^{2}(x, x') \in \mathcal{T}_{2,2}$. Here, (a) follows from the Cauchy-Schwarz inequality, (b) is obtained by definition of $\partial_{h}^{2}\delta_{\omega}^{in}$, while (c) uses the general fact that the operator norm of an operator is upper-bounded by its Hilbert-Schmidt norm which is finite in our case by application of Proposition B.2. Moreover, (d) further uses the identity in Lemma H.3 for computing the Hilbert-Schmidt norm of sum/expectation of tensor-product operators. The last item (d) is precisely what Proposition F.4 controls when applying it to the family $\mathcal{T}_{2,2}$ and choosing $\mathbb{D}$ to be $\mathbb{P}$. Therefore, the following maximal inequality holds by direct application of Proposition F.4:

$$ \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_h^2 \delta_{\omega}^{in} \right] \leq \lambda^{-\frac{1}{4}} n^{-\frac{1}{2}} \left( c(\Omega) \max \left( \tilde{M}_{in,2} \tilde{L}_{in,2} \operatorname{diam}(\Omega), \Lambda \tilde{M}_{in,2}^2 \right) \right)^{\frac{1}{4}}, $$

where $c(\Omega)$ is a positive constant greater than 1 that depends only on $\Omega$ and $d$. $\square$

## E.3 Proof of Theorem 4.1

**Theorem E.7 (Generalization bounds).** The following holds under Assumptions (A) to (E):

$$ \begin{aligned} \mathbb{E} \left[ \sup_{\omega \in \Omega} \left| \mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega) \right| \right] &\le \frac{1}{\lambda m^{\frac{1}{2}}} + \frac{C_{out}}{\lambda^{\frac{5}{4}} n^{\frac{1}{2}}}, \\ \mathbb{E} \left[ \sup_{\omega \in \Omega} \left\| \nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega) \right\| \right] &\le \frac{1}{\lambda} \left( d^{\frac{1}{2}} + \frac{C_{in}}{\lambda^{\frac{1}{4}}} \right) \frac{1}{m^{\frac{1}{2}}} + \frac{C_{out}}{\lambda^{\frac{5}{4}}} \left( 2 + 3 \frac{C_{in}}{\lambda} + \frac{C_{in}^2}{\lambda^2} \right) \frac{1}{n^{\frac{1}{2}}}, \end{aligned} $$

where the constants $C_{in}$ and $C_{out}$ are given in Proposition E.2.

*Proof.* Using the point-wise estimates in Proposition E.4 and taking their supremum over $\Omega$ followed by the expectations over data, the following error bounds hold:

$$ \begin{aligned} \mathbb{E} \left[ \sup_{\omega \in \Omega} \left| \mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega) \right| \right] &\le \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \delta_{\omega}^{\mathrm{out}} \right] + \frac{C_{\mathrm{out}}}{\lambda} \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{h} \delta_{\omega}^{\mathrm{in}} \right], \\ \mathbb{E} \left[ \sup_{\omega \in \Omega} \left\| \nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega) \right\| \right] &\le \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_{\omega} \delta_{\omega}^{\mathrm{out}} \right] + \frac{C_{\mathrm{in}}}{\lambda} \mathbb{E}_{\mathbb{Q}} \left[ \sup_{\omega \in \Omega} \partial_{h} \delta_{\omega}^{\mathrm{out}} \right] \\ &\quad + \frac{C_{\mathrm{out}}}{\lambda} \left( 1 + 2 \frac{C_{\mathrm{in}}}{\lambda} + \frac{C_{\mathrm{in}}^2}{\lambda^2} \right) \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{h} \delta_{\omega}^{\mathrm{in}} \right] \\ &\quad + \frac{C_{\mathrm{out}} C_{\mathrm{in}}}{\lambda^2} \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{h}^2 \delta_{\omega}^{\mathrm{in}} \right] + \frac{C_{\mathrm{out}}}{\lambda} \mathbb{E}_{\mathbb{P}} \left[ \sup_{\omega \in \Omega} \partial_{\omega,h}^2 \delta_{\omega}^{\mathrm{in}} \right]. \end{aligned} $$

Furthermore, we can use the maximal inequalities in Propositions E.5 and E.6 to control each term appearing in the right-hand side of the above inequalities:

$$ \begin{aligned} \mathbb{E} \left[ \sup_{\omega \in \Omega} \left| \mathcal{F}(\omega) - \hat{\mathcal{F}}(\omega) \right| \right] &\le R \left( m^{-\frac{1}{2}} \lambda^{-1} + C_{\mathrm{out}} n^{-\frac{1}{2}} \lambda^{-(1+\frac{1}{4})} \right), \\ \mathbb{E} \left[ \sup_{\omega \in \Omega} \left\| \nabla \mathcal{F}(\omega) - \widehat{\nabla \mathcal{F}}(\omega) \right\| \right] &\le R \left( m^{-\frac{1}{2}} \lambda^{-1} d^{\frac{1}{2}} + C_{\mathrm{in}} m^{-\frac{1}{2}} \lambda^{-(1+\frac{1}{4})} \right. \\ &\qquad \left. + C_{\mathrm{out}} n^{-\frac{1}{2}} \lambda^{-(1+\frac{1}{4})} \left( 1 + 2 \frac{C_{\mathrm{in}}}{\lambda} + \frac{C_{\mathrm{in}}^2}{\lambda^2} \right) \right. \\ &\qquad \left. + C_{\mathrm{out}} C_{\mathrm{in}} n^{-\frac{1}{2}} \lambda^{-(2+\frac{1}{4})} + C_{\mathrm{out}} n^{-\frac{1}{2}} \lambda^{-(1+\frac{1}{4})} \right), \end{aligned} $$

where the constant $R$ depends only on the Lipschitz constants $\mathrm{Lip}_{in}$ and $\mathrm{Lip}_{out}$, the upper bounds $M_{in}$ and $M_{out}$, the bound $\kappa$ on the kernel, the set $\Omega$, and the dimension $d$. Rearranging the obtained upper bounds concludes the proof. $\square$

## E.4 Generalization for bilevel gradient methods

*Proof of Corollary 4.2.* Consider that $\inf_{\omega,v,y} \ell_{out}(\omega, v, y) - c\|\omega\|^2 \ge 0$, which entails $\ell_{out}(\omega, v, y) \ge c\|\omega\|^2$ for all $v, y$. Using Proposition D.1 and setting $B = \sup_{y \in Y} |\partial_v \ell_{in}(\omega_0, 0, y)|$,
we have almost surely that:

$$ \hat{\mathcal{F}}(\omega_0) \leq \max_{|v| \leq \frac{B\kappa}{\lambda}, y \in \mathcal{Y}} \ell_{out}(\omega_0, v, y) =: \bar{\ell}. $$

Therefore, for any $\omega$ such that $\hat{\mathcal{F}}(\omega) \leq \hat{\mathcal{F}}(\omega_0)$, we have $\|\omega\|^2 \leq \bar{\ell}/c$. Define $\Omega$ as the ball of radius $\sqrt{\bar{\ell}/c}$ centered at 0. Using the fact that $\nabla \hat{\mathcal{F}} = \nabla \mathcal{F}$ in Proposition 3.1 and the representation in Equation (4), it is clear from Proposition D.2 and Assumption (D) that $\nabla \hat{\mathcal{F}}$ is Lipschitz on $\Omega$ with a deterministic constant $L$. It follows from standard results on gradient descent for nonconvex $\hat{\mathcal{F}}$ with Lipschitz gradient (see, e.g., [9, Theorems 4.25, 4.26]) that if we take $\bar{\eta} = 1/L$, then almost surely:

* $\hat{\mathcal{F}}(\omega_t) \leq \hat{\mathcal{F}}(\omega_0)$ and $\omega_t \in \Omega$ for all $t \geq 0$.

* $\nabla \hat{\mathcal{F}}(\omega_t) \to 0$ as $t \to \infty$.

* $\min_{i=0,...,t} \|\nabla \hat{\mathcal{F}}(\omega_i)\| \leq \bar{c}/\sqrt{t+1}$ for all $t \geq 0$, where $\bar{c}$ is a deterministic constant.

The corollary then follows by combining Proposition 3.1 and the uniform bound in Theorem 4.1. $\square$

**Bilevel projected gradient descent.** Considering the constrained (KBO) problem and assuming that $C$ is convex and compact, the projected gradient descent initialized at $\omega_0 \in C$ iterates the following recursion $\omega_{t+1} = \Pi_C(\omega_t - \eta \nabla \hat{\mathcal{F}}(\omega_t))$ for all $t \geq 0$, where $\Pi_C$ denotes the orthogonal projection onto $C$ and $\eta > 0$ is the step size. The algorithmic requirements are the same as the gradient descent algorithm, with the additional cost of computing the projection, which is typically cheap for basic sets such as balls. In the constrained setting, the optimality condition should take the constraints into account. To this end, we consider the gradient mappings $\hat{G}_\eta: \omega \mapsto \frac{1}{\eta}(\omega - \Pi_C(\omega - \eta \nabla \hat{\mathcal{F}}(\omega)))$ and $G_\eta: \omega \mapsto \frac{1}{\eta}(\omega - \Pi_C(\omega - \eta \nabla \mathcal{F}(\omega)))$ [10, Section 10.3]. This captures the stationarity of the recursion, and any local minimum of $\mathcal{F}$ on $C$ satisfies $G_\eta = 0$ for all $\eta > 0$.

**Corollary E.8 (Generalization for bilevel projected gradient descent).** Consider Assumptions (A) to (E) and fix $\lambda > 0$. Assume further that $K$ in (KBO) is almost surely definite, and that $C$ is convex and compact. Fix $\omega_0 \in C$ and let $\omega_{t+1} = \Pi_C(\omega_t - \eta \nabla \hat{\mathcal{F}}(\omega_t))$, where $\eta > 0$ is the step size and $t \geq 0$ is the iteration index. Then, there exist constants $\bar{\eta} > 0$ and $\bar{c} > 0$ such that for any $0 < \eta < \bar{\eta}$ and $t > 0$, the following holds:

$$ \mathbb{E}\left[\min_{i=0,...,t} \|G_{\eta}(\omega_i)\|\right] \leq \bar{c}\left(\frac{1}{\sqrt{m}} + \frac{1}{\sqrt{n}} + \frac{1}{\sqrt{t+1}}\right), \quad \mathbb{E}\left[\limsup_{i \to \infty} \|G_{\eta}(\omega_i)\|\right] \leq \bar{c}\left(\frac{1}{\sqrt{m}} + \frac{1}{\sqrt{n}}\right). $$

*Proof.* We choose $\Omega = C$. All iterates obviously remain in $\Omega$. Similarly as in the proof of Corollary 4.2, we know that $\nabla \hat{\mathcal{F}}$ is Lipschitz and that $\mathcal{F}$ is bounded on $\Omega$ with deterministic constants. It then follows from classical analysis on the nonconvex projected gradient algorithm (see, e.g., [10, Theorem 10.15]), that for sufficiently small $\eta$, we have almost surely that $\min_{i=0,...,t} \|\hat{G}_\eta(\omega_i)\| \leq \bar{c}/\sqrt{t+1}$ for a deterministic constant $\bar{c} > 0$, and that $\hat{G}_\eta(\omega_i) \to 0$ as $i \to \infty$. Using the fact that the orthogonal projection is 1-Lipschitz, (see, e.g., [10, Theorem 6.42]), we also have for all $\omega \in \Omega$:

$$ \|\hat{G}_{\eta}(\omega) - G_{\eta}(\omega)\| \leq \|\nabla \hat{\mathcal{F}}(\omega) - \nabla \mathcal{F}(\omega)\|. $$

The result follows by combining Proposition 3.1 and the uniform bound in Theorem 4.1. $\square$

## F Maximal Inequalities for Bounded and Lipschitz Family of Functions

Let $\mathcal{Z}$ be a subset of a Euclidean space and $\Omega$ be a compact subset of $\mathbb{R}^d$. Denote by $\otimes^k \mathcal{Z}$ the $k$-th tensor power of $\mathcal{Z}$, for any $k \geq 1$. Consider a parametric family $\mathcal{T}$ of real-valued functions defined over $\mathcal{Z}$ and indexed by a parameter $\omega \in \Omega$, i.e.,

$$ \mathcal{T} := \{\mathcal{Z} \ni z \mapsto t_{\omega}(z) \in \mathbb{R} \mid \omega \in \Omega\}. \quad (18) $$
For a given probability measure $\mu$ on $\mathcal{Z}$, denote by $L_2(\mu)$ the space of square $\mu$-integrable real-valued functions. We denote by $\|f\|_{\mathcal{D},2} := \mathbb{E}_{\mathcal{D}}[f(z)^2]^{\frac{1}{2}}$ the $L_2(\mu)$-norm of any function $f \in L_2(\mu)$. For any $\epsilon > 0$, we denote by $D(\epsilon, \mathcal{T}, L_2(\mu))$ the $\epsilon$-packing number of $\mathcal{T}$ w.r.t. $L_2(\mu)$. The next proposition provides a control on such a number under regularity conditions on the family $\mathcal{T}$.

**Proposition F.1 (Control on the packing number).** Assume that $\Omega$ is a compact subset of $\mathbb{R}^d$, that the parametric family $\mathcal{T}$ defined in Equation (18) is uniformly bounded by a positive constant $M$, and that there exists a positive constant $L$ so that, for any probability measure $\mu$ on $\mathcal{Z}$, $\omega \mapsto t_\omega(z)$ is $L$-Lipschitz for any $z \in \mathcal{Z}$. Then, there exists a positive constant $c(\Omega)$ greater than 1 that depends only on $\Omega$ and $d$ so that, for any probability measure $\mu$ on $\mathcal{Z}$, the following bound holds for any $0 < \epsilon \le M$:

$$D(\epsilon, \mathcal{T}, L_2(\mu)) \le c(\Omega) \left( \frac{\max(L \operatorname{diam}(\Omega), M)}{\epsilon} \right)^d.$$

*Proof.* First using [40, Lemma 9.18] and [40, Paragraph 8.1.2], we know that the $\epsilon$-packing number $D(\epsilon, \mathcal{T}, L_2(\mu))$ is smaller than the $\frac{\epsilon}{2}$-bracketing number $N_{\Box}(\frac{\epsilon}{2}, \mathcal{T}, L_2(\mu))$. Hence, we only need to control the bracketing number. To this end, we recall that the function $\omega \mapsto t_\omega(z)$ is $L$-Lipschitz for any $z \in \mathcal{Z}$, so that [71, Example 19.7] ensures the existence of a positive constant $c(\Omega)$ that depends only on $\Omega$ for which the following inequality holds for any $0 < \epsilon < L\operatorname{diam}(\Omega)$:

$$1 \le N_{\Box}(\epsilon, \mathcal{T}, L_2(\mu)) \le c(\Omega) \left( \frac{L \operatorname{diam}(\Omega)}{\epsilon} \right)^d.$$

Moreover, since the $\epsilon$-bracketing number is decreasing in $\epsilon$, it holds that:

$$N_{\Box}(\epsilon, \mathcal{T}, L_2(\mu)) \le N_{\Box}(\epsilon_-, \mathcal{T}, L_2(\mu)) \le c(\Omega) \left( \frac{L \operatorname{diam}(\Omega)}{\epsilon_-} \right)^d,$$

for any $\epsilon \ge L\operatorname{diam}(\Omega)$ and $\epsilon_- \le L\operatorname{diam}(\Omega)$. Taking the limit when $\epsilon_-$ approaches $L\operatorname{diam}(\Omega)$ yields $N_{\Box}(\epsilon, \mathcal{T}, L_2(\mu)) \le c(\Omega)$ for any $\epsilon \ge L\operatorname{diam}(\Omega)$. Hence, we have shown so far that for any $\epsilon > 0$:

$$N_{\Box}(\epsilon, \mathcal{T}, L_2(\mu)) \le c(\Omega) \max \left( 1, \left( \frac{L \operatorname{diam}(\Omega)}{\epsilon} \right)^d \right).$$

Moreover, by noticing that $\max(1, \frac{L \operatorname{diam}(\Omega)}{\epsilon}) \le \frac{\max(M, L \operatorname{diam}(\Omega))}{\epsilon}$ for any $\epsilon \le M$, we further have that:

$$N_{\Box}(\epsilon, \mathcal{T}, L_2(\mu)) \le c(\Omega) \left( \frac{\max(M, L \operatorname{diam}(\Omega))}{\epsilon} \right)^d.$$

Finally, recalling that $D(\epsilon, \mathcal{T}, L_2(\mu)) \le N_{\Box}(\frac{\epsilon}{2}, \mathcal{T}, L_2(\mu))$, we get that $D(\epsilon, \mathcal{T}, L_2(\mu)) \le 2^d c(\Omega) \left( \frac{\max(M, L \operatorname{diam}(\Omega))}{\epsilon} \right)^d$. The desired bound follows after redefining $c(\Omega)$ to include the factor $2^d$ (i.e., $c(\Omega) \to 2^d c(\Omega)$). $\square$

**Theorem F.2 (Maximal inequality for degenerate, bounded, and Lipschitz U-processes).** Let $k$ be either 1 or 2. Consider a parametric family $\mathcal{T} := \{\otimes^k \mathcal{Z} \ni (z_1, \dots, z_k) \mapsto t_\omega(z_1, \dots, z_k) \in \mathbb{R} \mid \omega \in \Omega\}$ of real-valued functions over $\otimes^k \mathcal{Z}$ indexed by a parameter $\omega \in \Omega$, where $\Omega$ is a compact subset of $\mathbb{R}^d$. For a given probability distribution $\mathbb{D}$ over $\mathcal{Z}$, assume that all elements $t_\omega$ are degenerate w.r.t. $\mathbb{D}$, meaning that:

$$\begin{cases} \mathbb{E}_{\bar{\mathcal{Z}} \sim \mathbb{D}}[t_\omega(\bar{z})] = 0, & \text{if } k=1 \\ \mathbb{E}_{\bar{\mathcal{Z}} \sim \mathbb{D}}[t_\omega(z, \bar{z})] = \mathbb{E}_{\bar{\mathcal{Z}} \sim \mathbb{D}}[t_\omega(\bar{z}, z)] = 0, & \text{if } k=2. \end{cases}$$

Furthermore, assume that all functions in $\mathcal{T}$ are uniformly bounded by a positive constant $M$ and that there exists a positive constant $L$ so that $\omega \mapsto t_\omega(z_1, \dots, z_k)$ is $L$-Lipschitz for any $(z_1, \dots, z_k) \in \otimes^k \mathcal{Z}$. Given i.i.d. samples $(z_i)_{1 \le i \le n}$ from $\mathbb{D}$, consider the following U-statistic $U_n^k$:

$$U_n^k t_\omega := \begin{cases} \frac{1}{n} \sum_{i=1}^n t_\omega(z_i), & \text{if } k=1 \\ \frac{1}{n(n-1)} \sum_{\substack{i,j=1 \\ i \ne j}}^n t_\omega(z_i, z_j), & \text{if } k=2. \end{cases}$$
Then, there exists a universal positive constant $c(\Omega)$ greater than 1 that depends only on $\Omega$ and $d$ such that for any $p \in \{1, 2\}$:

$$ \mathbb{E}_{\mathcal{D}} \left[ \sup_{\omega \in \Omega} |U_n^k t_\omega|^p \right]^{\frac{1}{p}} \le n^{-\frac{k}{2}} c(\Omega) \max (ML \operatorname{diam}(\Omega), M^2)^{\frac{1}{2}}. $$

**Proof.** **Maximal inequality for degenerate U-processes.** We will first apply the general result in [67, Maximal inequality] which controls $\mathbb{E}_{\mathcal{D}}[\sup_{\omega \in \Omega} |U_n^k t_\omega|]$ in terms of the packing number of $\mathcal{T}$. First note, by assumption, that the functions $t_\omega(z_1, \dots, z_k)$ are uniformly bounded by a positive constant $M$. Therefore, the constant function $T(z_1, \dots, z_k) := M$ is an envelope for $\mathcal{T}$, i.e., $T$ satisfies $T(z_1, \dots, z_k) \ge \sup_{\omega \in \Omega} |t_\omega(z_1, \dots, z_k)|$ for any $(z_1, \dots, z_k) \in \otimes^k \mathcal{Z}$. The envelope $T$ is, a fortiori, square $\mu$-integrable for any probability measure $\mu$ on $\otimes^k \mathcal{Z}$. Hence, we can apply [67, Maximal inequality] with the choice $T$ for the envelope function and set the integer $m$ appearing in the result to $m=d$ to get the following bound:

$$ \mathbb{E}_{\mathcal{D}} \left[ \sup_{\omega \in \Omega} |U_n^k t_\omega|^p \right]^{\frac{1}{p}} \le n^{-\frac{k}{2}} \Gamma \mathbb{E} \left[ \|T\|_{\mu_{n,2}} \int_0^{\delta_n} \left( D \left( \epsilon \|T\|_{\mu_{n,2}}, \mathcal{T}, L_2(\mu_n) \right) \right)^{\frac{1}{2\delta_p}} d\epsilon \right], \quad (19) $$

where $\Gamma$ is a positive universal constant<fnref n="3" /> that depends only on $d$ and that we choose to be greater than 1, while $\mu_n$ are suitably chosen probability measures on $\otimes^k \mathcal{Z}$ that possibly depend on the samples $z_1, \dots, z_n$ and other random variables, and $\delta_n \|T\|_{\mu_{n,2}} := \sup_{\omega \in \Omega} \|t_\omega\|_{\mu_{n,2}}$. Here, the expectation symbol in the right-hand side is over all randomness on which $\mu_n$ might depend. Note that the original result in [67, Maximal inequality] is stated using a slightly different definition of the packing number but which is still equivalent to the statement above in our setting<fnref n="4" />.

In our setting, the envelope function is constant and equal to $M$, and by definition $\delta_n \le 1$. Hence, the inequality in Equation (19) further becomes:

$$ \mathbb{E}_{\mathcal{D}} \left[ \sup_{\omega \in \Omega} |U_n^k t_\omega|^p \right]^{\frac{1}{p}} \le n^{-\frac{k}{2}} M \Gamma \mathbb{E} \left[ \int_0^1 \left( D \left( \epsilon \|T\|_{\mu_{n,2}}, \mathcal{T}, L_2(\mu_n) \right) \right)^{\frac{1}{2\delta_p}} d\epsilon \right]. \quad (20) $$

We simply need to control the packing number $D(\epsilon \|T\|_{\mu_{n,2}}, \mathcal{T}, L_2(\mu))$ independently of the probability measure $\mu$.

**Control on the packing number.** We have shown that the constant function $T(z_1, \dots, z_k) := M$ is an envelope for $\mathcal{T}$ which is, a fortiori, square $\mu$-integrable for any probability measure $\mu$ with $\|T\|_{\mu_{n,2}} = M < +\infty$. Moreover, the functions $\omega \mapsto t_\omega(z_1, \dots, z_k)$ are $L$-Lipschitz for any $(z_1, \dots, z_k) \in \otimes^k \mathcal{Z}$. We can therefore apply Proposition F.1 which ensures the existence of a positive constant $c(\Omega)$ greater than 1 and that depends only on $\Omega$ and $d$ so that the following estimate on the $\epsilon$-packing number of the class $\mathcal{T}$ w.r.t. $L_2(\mu)$ holds:

$$ D \left( \epsilon \|T\|_{\mu_{n,2}}, \mathcal{T}, L_2(\mu) \right) \le \underbrace{c(\Omega) \left( \max \left( \frac{L \operatorname{diam}(\Omega)}{M}, 1 \right) \right)^d \left( \frac{1}{\epsilon} \right)^d}_{A}, \quad \forall \epsilon \in (0, 1]. \quad (21) $$

Combining Equation (21) with Equation (20) yields:

$$ \begin{aligned} \mathbb{E}_{\mathcal{D}} \left[ \sup_{\omega \in \Omega} |U_n^k t_\omega|^p \right]^{\frac{1}{p}} &\le n^{-\frac{k}{2}} M \Gamma \mathbb{E} \left[ \int_0^1 (A \epsilon^{-d})^{\frac{1}{2\delta_p}} d\epsilon \right] = n^{-\frac{k}{2}} M \Gamma A^{\frac{1}{2\delta_p}} \underbrace{\int_0^1 \epsilon^{-\frac{1}{2p}} d\epsilon}_{\le 2} \\ &\le 2n^{-\frac{k}{2}} \Gamma c(\Omega)^{\frac{1}{2d}} \max (L \operatorname{diam}(\Omega), M^2)^{\frac{1}{2}}, \end{aligned} $$

<footnote n="3">³The constant $\Gamma$ appearing in [67, Maximal inequality] depends only on $k, p$ and $m$, i.e., $\Gamma := g(k, p, m)$. Since, we are only interested in $k \le 2$ and $p \le 2$ and $m$ is fixed to $d$, we choose $\Gamma$ to be $\max_{1 \le k, p \le 2} g(k, p, d)^{\frac{1}{p}}$, so that it is the same in all our cases.</footnote>

<footnote n="4">⁴In [67, Maximal inequality], the author considers a modified version of the $\epsilon$-packing number (call it $\tilde{D}(\epsilon, \mathcal{T}, L_2(\mu))$) associated to $L_2(\mu)$ but endowed with a normalized version of the standard norm on $L_2(\mu): \|f\|_{\mu} := \frac{\|f\|_{\mu,2}}{\|T\|_{\mu,2}}$. Both numbers are related by the following identity: $\tilde{D}(\epsilon, \mathcal{T}, L_2(\mu)) = D(\epsilon \|T\|_{\mu,2}, \mathcal{T}, L_2(\mu))$, thus making the statement (19) equivalent to the original statement in [67, Maximal inequality].</footnote>
Using the triangle inequality in the above equality and recalling that, by assumption, $t_\omega(z, z')$ is uniformly bounded by a positive constant $M$, it follows that:

$$|\tau_\omega| \le |U_n^2 s_\omega| + \frac{1}{n^2(n-1)} \sum_{\substack{i,j=1 \\ i \ne j}}^n |t_\omega(z_i, z_j)| + \frac{1}{n^2} \sum_{i=1}^n |t_\omega(z_i, z_i)| \le |U_n^2 s_\omega| + \frac{2M}{n}.$$

Furthermore, taking the supremum over $\omega$ followed by the expectation over samples yields:

$$\mathbb{E}_{\mathbb{D}} \left[ \sup_{\omega \in \Omega} |\tau_\omega| \right] \le \mathbb{E}_{\mathbb{D}} \left[ \sup_{\omega \in \Omega} |U_n^2 s_\omega| \right] + \frac{2M}{n}. \quad (23)$$

Hence, it only remains to control the first term in the above inequality. To this end, we will use a maximal inequality for degenerate U-processes due to [67].

**Maximal inequality for degenerate U-processes.** We will first check that $U_n^2 s_\omega$ is a degenerate statistic for a given $\omega \in \Omega$. Simple calculations show that for any $z$ in $\mathcal{Z}$:

$$\mathbb{E}_{\bar{z} \sim \mathbb{D}} [s_\omega(z, \bar{z})] = \mathbb{E}_{\bar{z} \sim \mathbb{D}} [s_\omega(\bar{z}, z)] = 0.$$

The above equalities precisely ensure that $U_n^2 s_\omega$ is a degenerate U-statistic for $\mathbb{D}$. Consider now the family $S := \{\mathcal{Z} \times \mathcal{Z} \ni (z, z') \mapsto s_\omega(z, z') \in \mathbb{R} \mid \omega \in \Omega\}$. We show that $S$ is uniformly bounded and Lipschitz which allows to directly apply the result stated in Theorem F.2, which is a special case of the more general result in [67, Maximal inequality]. First note, by assumption, that the functions $t_\omega(z, z')$ are uniformly bounded by a positive constant $M$. Hence, using Equation (22), it follows that $s_\omega(z, z')$ is uniformly bounded by $4M$. Moreover, the functions $\omega \mapsto t_\omega(z, z')$ are L-Lipschitz for any $z, z'$ in $\mathcal{Z}$. Hence, from Equation (22), we directly have that $\omega \mapsto s_\omega(z, z')$ is 4L-Lipschitz for any $z, z' \in \mathcal{Z}$. We can directly apply Theorem F.2 with $k=2$ and $p=1$ to $S$ and get the following maximal inequality:

$$\mathbb{E}_{\mathbb{D}} \left[ \sup_{\omega \in \Omega} |U_n^2 s_\omega| \right] \le 4n^{-1} c(\Omega) \max (ML \operatorname{diam}(\Omega), M^2)^{\frac{1}{2}}.$$

We obtain an upper bound on $\mathbb{E}_{\mathbb{D}}[\sup_{\omega \in \Omega} |\tau_\omega|]$ by combining the above inequality with Equation (23), then noticing that $2M \le 2c(\Omega) \max(L \operatorname{diam}(\Omega), M)$ so that:

$$\mathbb{E}_{\mathbb{D}} \left[ \sup_{\omega \in \Omega} |\tau_\omega| \right] \le 6n^{-1} c(\Omega) \max (ML \operatorname{diam}(\Omega), M^2)^{\frac{1}{2}}.$$

Finally, the desired result follows by redefining $c(\Omega)$ to include the factor 6 in the above inequality. $\square$

## G Differentiability Results

The proofs of Propositions B.1 and B.2 are direct applications of the following more general result.

**Proposition G.1.** Let $\mathcal{U}$ be an open non-trivial subset of $\mathbb{R}^d$. Consider a real-valued function $\ell: (\omega, v, y) \mapsto \ell(\omega, v, y)$ defined on $\mathcal{U} \times \mathbb{R} \times \mathcal{Y}$ that is of class $C^3$ jointly in $(\omega, v)$ and whose derivatives are jointly continuous in $(\omega, v, y)$. For a given probability distribution $\mathbb{D}$ over $\mathcal{X} \times \mathcal{Y}$, consider the following functional defined over $\mathcal{U} \times \mathcal{H}$:

$$L(\omega, h) := \mathbb{E}_{\mathbb{D}} [\ell(\omega, h(x), y)].$$

Under Assumptions (A) to (C), the following properties hold for $L$:

* $L$ admits finite values for any $(\omega, h) \in \mathcal{U} \times \mathcal{H}$.

* $(\omega, h) \mapsto L(\omega, h)$ is Fréchet differentiable with partial derivatives $\partial_\omega L(\omega, h)$ and $\partial_h L(\omega, h)$ at any point $(\omega, h) \in \mathcal{U} \times \mathcal{H}$ given by:

$$\partial_\omega L(\omega, h) = \mathbb{E}_{\mathbb{D}} [\partial_\omega \ell(\omega, h(x), y)] \in \mathbb{R}^d,$$

$$\partial_h L(\omega, h) = \mathbb{E}_{\mathbb{D}} [\partial_v \ell(\omega, h(x), y) K(x, \cdot)] \in \mathcal{H}.$$
* The map $(\omega, h) \mapsto \partial_h L(\omega, h)$ is differentiable. Moreover, for any $(\omega, h) \in \mathcal{U} \times \mathcal{H}$, its partial derivatives $\partial_{\omega,h}^2 L(\omega, h)$ and $\partial_h^2 L(\omega, h)$ at $(\omega, h)$ are Hilbert-Schmidt operators given by:

$$ 
\begin{aligned}
\partial_{\omega,h}^2 L(\omega, h) &= \mathbb{E}_{\mathcal{D}} [\partial_{\omega,v}^2 \ell(\omega, h(x), y) K(x, \cdot)] \in \mathcal{L}(\mathcal{H}, \mathbb{R}^d), \\
\partial_h^2 L(\omega, h) &= \mathbb{E}_{\mathcal{D}} [\partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot)] \in \mathcal{L}(\mathcal{H}, \mathcal{H}).
\end{aligned}
$$

**Proof. Finite values.** Fix $\omega \in \mathcal{U}$ and $h \in \mathcal{H}$. We will first show that $h$ is bounded on $\mathcal{X}$. By the reproducing property, we know that $|h(x)| \le \|h\|_{\mathcal{H}} \sqrt{K(x, x)}$ for any $x \in \mathcal{X}$. Moreover, the kernel $K$ is bounded by a constant $\kappa$ thanks to Assumption (B). Consequently, $|h(x)|$ is upper-bounded by $\|h\|_{\mathcal{H}} \sqrt{\kappa}$ for any $x \in \mathcal{X}$.

Denote by $\mathcal{I}$ the compact interval defined as $\mathcal{I} = [-\|h\|_{\mathcal{H}} \sqrt{\kappa}, \|h\|_{\mathcal{H}} \sqrt{\kappa}]$. By Assumption (C), the set $\mathcal{Y}$ is compact so that $\mathcal{I} \times \mathcal{Y}$ is also compact. Moreover, we know, by assumption on $\ell$, that $(\omega, v, y) \mapsto \ell(\omega, v, y)$ is continuous on $\mathcal{U} \times \mathbb{R} \times \mathcal{Y}$. Therefore, $(v, y) \mapsto \ell(\omega, v, y)$ must be bounded by some finite constant $C$ on the compact set $\mathcal{I} \times \mathcal{Y}$. This allows to deduce that $(x, y) \mapsto \ell(\omega, h(x), y)$ is bounded by $C$ for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$ and a fortiori $\mathbb{D}$-integrable, which shows that $L(\omega, h)$ is finite.

**Fréchet differentiability of L.** Let $(\omega, h) \in \mathcal{U} \times \mathcal{H}$. Consider $(\omega_j, h_j)_{j \ge 1}$ a sequence of elements in $\mathcal{U} \times \mathcal{H}$ converging to it, i.e., $(\omega_j, h_j) \to (\omega, h)$ with $(\omega_j, h_j) \ne (\omega, h)$ for any $j \ge 0$. Define the sequence of functions $r_j : \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$ as follows:

$$ r_j(x, y) = \frac{\ell(\omega_j, h_j(x), y) - \ell(\omega, h(x), y) - \langle \partial_v \ell(\omega, h(x), y) K(x, \cdot), h_j - h \rangle_{\mathcal{H}} - \langle \partial_{\omega} \ell(\omega, h(x), y), \omega_j - \omega \rangle}{\| (\omega_j, h_j) - (\omega, h) \|} \quad (24) $$

We will first show that $\mathbb{E}_{\mathcal{D}} [|r_j(x, y)|]$ converges to 0 by the dominated convergence theorem [62, Theorem 1.34]. By the reproducing property, note that $\ell(\omega, h(x), y) = \ell(\omega, \langle h, K(x, \cdot) \rangle_{\mathcal{H}}, y)$. Hence, since $\ell$ is jointly differentiable in $(\omega, v)$ for any $y$, it follows that $(\omega, h) \mapsto \ell(\omega, h(x), y)$ is also differentiable for any $(x, y)$ by composition with the evaluation map $(\omega, h) \mapsto (\omega, \langle h, K(x, \cdot) \rangle_{\mathcal{H}})$ which is differentiable. Hence, the sequence $r_j(x, y)$ converges to 0 for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$. Moreover, by the mean-value theorem, there exists $0 \le c_j \le 1$ such that:

$$ r_j(x, y) = \frac{\langle (\partial_v \ell(\bar{\omega}_j, \bar{h}_j(x), y) - \partial_v \ell(\omega, h(x), y)) K(x, \cdot), h_j - h \rangle_{\mathcal{H}} - \langle \partial_{\omega} \ell(\bar{\omega}_j, \bar{h}_j(x), y) - \partial_{\omega} \ell(\omega, h(x), y), \omega_j - \omega \rangle}{\| (\omega_j, h_j) - (\omega, h) \|} $$

where $(\bar{\omega}_j, \bar{h}_j) := (1 - c_j)(\omega, h) + c_j(\omega_j, h_j)$. We will show that $r_j(x, y)$ is bounded for $j$ large enough. We first construct a compact set that will contain all elements of the form $(\omega_j, h_j(x), y)$ and $(\bar{\omega}_j, \bar{h}_j(x), y)$ for all $j$ large enough. Since $\omega$ is an element in the open set $\mathcal{U}$, there exists a closed ball $\mathcal{B}(\omega, R)$ centered at $\omega$ and with some radius $R$ small enough so that $\mathcal{B}(\omega, R)$ is included in $\mathcal{U}$. For all $j$ large enough, $\omega_j$ and $\bar{\omega}_j$ belong to $\mathcal{B}(\omega, R)$ as these sequences converge to $\omega$. Moreover, $h_j$ and $\bar{h}_j$ are convergent sequences. Consequently, they must be bounded by some constant $B$. By the reproducing property, and recalling that the kernel $K$ is bounded by $\kappa$ by Assumption (B), it follows that $\max(|h_j(x)|, |\bar{h}_j(x)|) \le B\kappa$. Consider now the set $\mathcal{W} := \mathcal{B}(\omega, R) \times \mathcal{B}_1(0, B\kappa) \times \mathcal{Y}$ which is a product of compact sets (recalling that $\mathcal{Y}$ is compact by Assumption (C)), where $\mathcal{B}_1(0, B\kappa)$ is the closed ball in $\mathbb{R}$ centered at 0 and of radius $B\kappa$. For $j$ large enough, we have established that $(\omega_j, h_j(x), y)$ and $(\bar{\omega}_j, \bar{h}_j(x), y)$ belong to $\mathcal{W}$ for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$. Since, by assumption on $\ell$, $\partial_v \ell(\omega, v, y)$ and $\partial_{\omega} \ell(\omega, v, y)$ are continuous, they must be bounded on the compact set $\mathcal{W}$ by some constant $C$. This allows to deduce from the expression of $r_j(x, y)$ above that $r_j(x, y)$ is bounded, and a fortiori dominated by an integrable function (a constant function). We then deduce that $\mathbb{E}_{\mathcal{D}} [|r_j(x, y)|]$ converges to 0 by application of the dominated convergence theorem [62, Theorem 1.34].

Recalling Equation (24), $\mathbb{E}_{\mathcal{D}} [r_j(x, y)]$ admits the following expression:

$$ \mathbb{E}_{\mathcal{D}} [r_j(x, y)] = \frac{L(\omega_j, h_j) - L(\omega, h) - \mathbb{E}_{\mathcal{D}} [\langle \partial_v \ell(\omega, h(x), y) K(x, \cdot), h_j - h \rangle_{\mathcal{H}}] - \langle \mathbb{E}_{\mathcal{D}} [\partial_{\omega} \ell(\omega, h(x), y)], \omega_j - \omega \rangle}{\| (\omega_j, h_j) - (\omega, h) \|} \quad (25) $$
The convergence to 0 of the above expression precisely means that $L$ is differentiable at $(\omega, h)$ provided that the linear form $g \mapsto \mathbb{E}_{\mathbb{D}}[(\partial_v \ell(\omega, h(x), y) K(x, \cdot), g)_{\mathcal{H}}]$ is bounded. To establish this fact, consider the RKHS-valued function $(x, y) \mapsto \partial_v \ell(\omega, h(x), y) K(x, \cdot)$. This function is Bochner-integrable in the sense that $\mathbb{E}_{\mathbb{D}}[\|\partial_v \ell(\omega, h(x), y) K(x, \cdot)\|_{\mathcal{H}}]$ is finite [25, Definition 1, Chapter 2]. Indeed, we have the following:

$$ \mathbb{E}_{\mathbb{D}}[\|\partial_v \ell(\omega, h(x), y) K(x, \cdot)\|_{\mathcal{H}}] := \mathbb{E}_{\mathbb{D}}\left[|\partial_v \ell(\omega, h(x), y)| \sqrt{K(x, x)}\right] \\ \leq \sqrt{\kappa} \mathbb{E}_{\mathbb{D}}[|\partial_v \ell(\omega, h(x), y)|] < +\infty, $$

where, for the inequality, we used that $(x, y) \mapsto \partial_v \ell(\omega, h(x), y)$ is bounded as shown previously. Consequently, $\mathbb{E}_{\mathbb{D}}[\partial_v \ell(\omega, h(x), y) K(x, \cdot)]$ is an element in $\mathcal{H}$ satisfying:

$$ \langle \mathbb{E}_{\mathbb{D}}[\partial_v \ell(\omega, h(x), y) K(x, \cdot)], g \rangle_{\mathcal{H}} = \mathbb{E}_{\mathbb{D}}[\langle \partial_v \ell(\omega, h(x), y) K(x, \cdot), g \rangle_{\mathcal{H}}], \forall g \in \mathcal{H}. $$

The above property follows from [25, Theorem 6, Chapter 2] for Bochner-integrable functions that allows exchanging the integral and the application of a continuous linear map (here the scalar product with an element $g$). The above identity establishes that $g \mapsto \mathbb{E}_{\mathbb{D}}[\langle \partial_v \ell(\omega, h(x), y) K(x, \cdot), g \rangle_{\mathcal{H}}]$ is bounded and provides the desired expression for $\partial_h L(\omega, h)$. The expression for $\partial_\omega L(\omega, h)$ directly follows from the last term in Equation (25).

**Fréchet differentiability of $\partial_h L$.** We use the same proof strategy as for the differentiability of $L$.

Let $(\omega, h) \in \mathcal{U} \times \mathcal{H}$. Consider $(\omega_j, h_j)_{j \ge 1}$ a sequence of elements in $\mathcal{U} \times \mathcal{H}$ converging to it, i.e., $(\omega_j, h_j) \to (\omega, h)$ with $(\omega_j, h_j) \neq (\omega, h)$ for any $j \ge 0$. Define the sequence of functions $s_j : \mathcal{X} \times \mathcal{Y} \to \mathcal{H}$ as follows:

$$ \| (\omega_j, h_j) - (\omega, h) \|_{s_j(x, y)} = (\partial_v \ell(\omega_j, h_j(x), y) - \partial_v \ell(\omega, h(x), y)) K(x, \cdot) \\ - \partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot) (h_j - h) \\ - (\omega_j - \omega)^\top \partial_{\omega, v}^2 \ell(\omega, h(x), y) K(x, \cdot). \quad (26) $$

We will first show that $\mathbb{E}_{\mathbb{D}}[\|s_j(x, y)\|_{\mathcal{H}}]$ converges to 0 by the dominated convergence theorem for Bochner-integrable functions [25, Theorem 3, Chapter 2]. By the reproducing property, note that $\partial_v \ell(\omega, h(x), y) K(x, \cdot) = \partial_v \ell(\omega, \langle h, K(x, \cdot) \rangle_{\mathcal{H}}, y) K(x, \cdot)$. Hence, since $(\omega, v) \mapsto \partial_v \ell(\omega, v, y)$ is jointly differentiable in $(\omega, v)$ for any $y$, it follows that $(\omega, h) \mapsto \partial_v \ell(\omega, h(x), y) K(x, \cdot)$ is also differentiable for any $(x, y)$ by composition with the evaluation map $(\omega, h) \mapsto (\omega, \langle h, K(x, \cdot) \rangle_{\mathcal{H}})$ which is differentiable. Hence, the sequence $s_j(x, y)$ converges to 0 for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$. Moreover, by the mean-value theorem, there exists $0 \le c_j \le 1$ such that:

$$ \| (\omega_j, h_j) - (\omega, h) \|_{s_j(x, y)} = \partial_v^2 \ell (\bar{\omega}_j, \bar{h}_j(x), y) K(x, \cdot) \otimes K(x, \cdot) (h_j - h) \\ + (\omega_j - \omega)^\top \partial_{\omega, v}^2 \ell (\bar{\omega}_j, \bar{h}_j(x), y) K(x, \cdot) \\ - \partial_v^2 \ell (\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot) (h_j - h) \\ - (\omega_j - \omega)^\top \partial_{\omega, v}^2 \ell (\omega, h(x), y) K(x, \cdot), $$

where $(\bar{\omega}_j, \bar{h}_j) := (1 - c_j)(\omega, h) + c_j(\omega_j, h_j)$. Using the same construction as for the Fréchet differentiability, we find a compact set $\mathcal{W}$ containing all elements $(\omega_j, h_j(x), y)$ and $(\bar{\omega}_j, \bar{h}_j(x), y)$ for any $(x, y) \in \mathcal{X} \times \mathcal{Y}$ and all $j$ large enough. On such set, $\partial_v^2 \ell(\omega, v, y)$ and $\partial_{\omega, v}^2 \ell(\omega, v, y)$ are bounded by some constant $C$. Consequently, we can write:

$$ \| (\omega_j, h_j) - (\omega, h) \|_{s_j(x, y)\|_{\mathcal{H}}} \le 2C \| K(x, \cdot) \otimes K(x, \cdot) (h_j - h) \|_{\mathcal{H}} + 2C \| \omega_j - \omega \|_{K(x, \cdot)\|_{\mathcal{H}}} \\ \le 2C\kappa \| h_j - h \|_{\mathcal{H}} + 2C\sqrt{\kappa} \| \omega_j - \omega \|. $$

This already establishes that $s_j(x, y)$ is bounded so that $\mathbb{E}_{\mathbb{D}}[\|s_j(x, y)\|_{\mathcal{H}}]$ converges to 0 by application of the dominated convergence theorem. Recalling Equation (26), $\mathbb{E}_{\mathbb{D}}[s_j(x, y)]$ admits the following expression:

$$ \| (\omega_j, h_j) - (\omega, h) \|_{\mathbb{E}_{\mathbb{D}}[s_j(x, y)]} = \partial_h L(\omega_j, h_j) - \partial_h L(\omega, h) \\ - \mathbb{E}_{\mathbb{D}}[\partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot) (h_j - h)] \\ - (\omega_j - \omega)^\top \mathbb{E}_{\mathbb{D}}[\partial_{\omega, v}^2 \ell(\omega, h(x), y) K(x, \cdot)]. $$
The convergence to 0 of the above expression precisely means that $L$ is differentiable at $(\omega, h)$ provided that: (1) $\mathbb{E}_{\mathbb{D}}[\partial_{\omega,v}^2 \ell(\omega, h(x), y) K(x, \cdot)]$ is an element in $\mathcal{H}^d$, and (2) the linear map $g \mapsto \mathbb{E}_{\mathbb{D}}[\partial_v^2 \ell(\omega, h(x), y) (K(x, \cdot) \otimes K(x, \cdot))g]$ is bounded. Using the same strategy to establish Bochner's integrability of $(x, y) \mapsto \partial_v \ell(\omega, h(x), y) K(x, \cdot)$, we can show that $(x, y) \mapsto \partial_{\omega,v}^2 \ell(\omega, h(x), y) K(x, \cdot)$ is also Bochner-integrable so that $\mathbb{E}_{\mathbb{D}}[\partial_{\omega,v}^2 \ell(\omega, h(x), y) K(x, \cdot)]$ is indeed an element in $\mathcal{H}^d$. This also establishes the expression of $\partial_{\omega,h} L(\omega, h)$. Similarly, we consider the operator-valued function $\xi: (x, y) \mapsto \partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot)$ with values in the space of Hilbert-Schmidt operators on $\mathcal{H}$. The Hilbert-Schmidt (HS) norm of such function satisfies the following inequality:

$$ \mathbb{E}_{\mathbb{D}} [ \| \partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot) \|_{HS} ] := \mathbb{E}_{\mathbb{D}} [ | \partial_v^2 \ell(\omega, h(x), y) | K(x, x) ] \le \kappa C < +\infty. $$

Therefore, the function $\xi$ is Bochner-integrable, so that $\mathbb{E}_{\mathbb{D}}[\partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot)]$ is a Hilbert-Schmidt operator satisfying:

$$ \mathbb{E}_{\mathbb{D}}[\partial_v^2 \ell(\omega, h(x), y) K(x, \cdot) \otimes K(x, \cdot)] g = \mathbb{E}_{\mathbb{D}}[\partial_v^2 \ell(\omega, h(x), y) (K(x, \cdot) \otimes K(x, \cdot))g], \forall g \in \mathcal{H}. $$

The above property follows from [25, Theorem 6, Chapter 2] for Bochner-integrable functions that allows exchanging the integral and the application of a continuous linear map (here the scalar product with an element $g$). Hence, from the above identity, we deduce the desired expression for $\partial_h^2 L(\omega, h)$. $\square$

## H Auxiliary Technical Lemmas

**Lemma H.1.** Let $A$ and $A'$ be two bounded operators from $\mathcal{H}$ to $\mathbb{R}^d$, and $B$ and $B'$ be two bounded and invertible operators from $\mathcal{H}$ to itself. Assume that $B \ge \lambda \mathrm{Id}_{\mathcal{H}}$ and $B' \ge \lambda \mathrm{Id}_{\mathcal{H}}$. Then, the following inequalities hold:

$$ \begin{aligned} \|AB^{-1} - A'(B')^{-1}\|_{\mathrm{op}} &\le \frac{\|A\|_{\mathrm{op}}}{\lambda^2} \|B - B'\|_{\mathrm{op}} + \frac{1}{\lambda} \|A - A'\|_{\mathrm{op}}, \\ \|AB^{-1}\|_{\mathrm{op}} &\le \lambda^{-1} \|A\|_{\mathrm{op}}, \quad \|A'(B')^{-1}\|_{\mathrm{op}} \le \lambda^{-1} \|A'\|_{\mathrm{op}}. \end{aligned} $$

*Proof.* By the triangle inequality and the sub-multiplicative property of the operator norm $\| \cdot \|_{\mathrm{op}}$, we have:

$$ \begin{aligned} \|AB^{-1} - A'(B')^{-1}\|_{\mathrm{op}} &\le \|AB^{-1} - A(B')^{-1}\|_{\mathrm{op}} + \|A(B')^{-1} - A'(B')^{-1}\|_{\mathrm{op}} \\ &= \|A(B^{-1} - (B')^{-1})\|_{\mathrm{op}} + \|(A - A')(B')^{-1}\|_{\mathrm{op}} \\ &\le \|A\|_{\mathrm{op}} \|B^{-1} - (B')^{-1}\|_{\mathrm{op}} + \|A - A'\|_{\mathrm{op}} \|(B')^{-1}\|_{\mathrm{op}} \\ &= \|A\|_{\mathrm{op}} \|B^{-1}(B' - B)(B')^{-1}\|_{\mathrm{op}} + \|A - A'\|_{\mathrm{op}} \|(B')^{-1}\|_{\mathrm{op}} \\ &\le \|A\|_{\mathrm{op}} \|B^{-1}\|_{\mathrm{op}} \|B' - B\|_{\mathrm{op}} \|(B')^{-1}\|_{\mathrm{op}} + \|A - A'\|_{\mathrm{op}} \|(B')^{-1}\|_{\mathrm{op}}. \end{aligned} \quad (27) $$

Since $B \ge \lambda \mathrm{Id}_{\mathcal{H}}$ and $B' \ge \lambda \mathrm{Id}_{\mathcal{H}}$, we obtain:

$$ \|B^{-1}\|_{\mathrm{op}} \le \frac{1}{\lambda} \quad \text{and} \quad \|(B')^{-1}\|_{\mathrm{op}} \le \frac{1}{\lambda}. $$

Substituting these into Equation (27), we get:

$$ \|AB^{-1} - A'(B')^{-1}\|_{\mathrm{op}} \le \frac{\|A\|_{\mathrm{op}}}{\lambda^2} \|B - B'\|_{\mathrm{op}} + \frac{1}{\lambda} \|A - A'\|_{\mathrm{op}}. $$

This proves the first inequality. The remaining two inequalities follow directly from the sub-multiplicative property of the operator norm $\| \cdot \|_{\mathrm{op}}$ and the assumptions $B \ge \lambda \mathrm{Id}_{\mathcal{H}}$ and $B' \ge \lambda \mathrm{Id}_{\mathcal{H}}$. $\square$

**Lemma H.2.** Let $f: \mathcal{H} \to \mathbb{R}$ be a $\lambda$-strongly convex and Fréchet differentiable function. Denote by $h^* \in \mathcal{H}$ its minimizer. Then, for any $h \in \mathcal{H}$, the following holds:

$$ \|h - h^*\|_{\mathcal{H}} \le \frac{1}{\lambda} \|\partial_h f(h)\|_{\mathcal{H}}. $$
Proof. Let $h \in \mathcal{H}$.

**Case 1:** $h = h^*$. The proof is straightforward.

**Case 2:** $h \neq h^*$. Given that $f$ is $\lambda$-strongly convex, we have:

$$f(h) - f(h^*) \geq \langle \partial_h f(h^*), h - h^* \rangle_{\mathcal{H}} + \frac{\lambda}{2} \|h - h^*\|_{\mathcal{H}}^2,$$ 
and $f(h^*) - f(h) \geq \langle \partial_h f(h), h^* - h \rangle_{\mathcal{H}} + \frac{\lambda}{2} \|h - h^*\|_{\mathcal{H}}^2.$

After summing these two inequalities, noticing that $\partial_h f(h^*) = 0$, and rearranging the terms, we obtain:

$$\langle \partial_h f(h), h - h^* \rangle_{\mathcal{H}} \geq \lambda \|h - h^*\|_{\mathcal{H}}^2.$$

After using the Cauchy-Schwarz inequality, we get:

$$\|\partial_h f(h)\|_{\mathcal{H}} \|h - h^*\|_{\mathcal{H}} \geq \lambda \|h - h^*\|_{\mathcal{H}}^2.$$

Dividing by $\lambda \|h - h^*\|_{\mathcal{H}} \neq 0$ concludes the proof. $\square$

**Lemma H.3.** Let $\mathcal{X}$ be a subset of $\mathbb{R}^p$, $\mathcal{Y}$ be a subset of $\mathbb{R}^q$, and $\mathbb{D}$ be a probability distribution over $\mathcal{X} \times \mathcal{Y}$. Given i.i.d. samples $(x_i, y_i)_{1 \le i \le n}$ drawn from $\mathbb{D}$, consider a function $g: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ of class $C^1$ such that the operator $A: \mathcal{H} \to \mathcal{H}$ defined as:

$$A := \mathbb{E}_{(x,y) \sim \mathbb{D}} [g(x, y) K(x, \cdot) \otimes K(x, \cdot)] - \frac{1}{n} \sum_{i=1}^{n} g(x_i, y_i) K(x_i, \cdot) \otimes K(x_i, \cdot)$$

is Hilbert-Schmidt. Then, the following holds:

$$\|A\|_{\text{HS}}^2 = \mathbb{E}_{(x,y),(x',y') \sim \mathbb{D} \otimes \mathbb{D}} \left[ g(x,y)g(x',y')K^2(x,x') \right] + \frac{1}{n^2} \sum_{i,j=1}^{n} g(x_i,y_i)g(x_j,y_j)K^2(x_i,x_j) \\ - \frac{2}{n} \sum_{i=1}^{n} \mathbb{E}_{(x,y) \sim \mathbb{D}} \left[ g(x_i,y_i)g(x,y)K^2(x,x_i) \right].$$

Proof. Define $s := \mathbb{E}_{(x,y) \sim \mathbb{D}} [g(x, y) K(x, \cdot) \otimes K(x, \cdot)]$ and $\hat{s} := \frac{1}{n} \sum_{i=1}^{n} g(x_i, y_i) K(x_i, \cdot) \otimes K(x_i, \cdot)$. We have:

$$\|A\|_{\text{HS}}^2 = \|s - \hat{s}\|_{\text{HS}}^2 = \|s\|_{\text{HS}}^2 + \|\hat{s}\|_{\text{HS}}^2 - 2 \langle s, \hat{s} \rangle_{\text{HS}}. \qquad (28)$$
Next, we compute each of the following quantities: $\|s\|_{HS}^2$, $\|\hat{s}\|_{HS}^2$, and $\langle s, \hat{s} \rangle_{HS}$, separately. Simple calculations yield:

$$
\begin{align*}
\|s\|_{HS}^2 &= \left\| \mathbb{E}_{(x,y)\sim\mathbb{D}} [g(x,y)K(x,\cdot) \otimes K(x,\cdot)] \right\|_{HS}^2 \\
&= \left\langle \mathbb{E}_{(x,y)\sim\mathbb{D}} [g(x,y)K(x,\cdot) \otimes K(x,\cdot)], \mathbb{E}_{(x',y')\sim\mathbb{D}} [g(x',y')K(x',\cdot) \otimes K(x',\cdot)] \right\rangle_{HS} \\
&= \mathbb{E}_{(x,y),(x',y')\sim\mathbb{D}\otimes\mathbb{D}} \left[ g(x,y)g(x',y') \left\langle K(x,\cdot) \otimes K(x,\cdot), K(x',\cdot) \otimes K(x',\cdot) \right\rangle_{HS} \right] \\
&= \mathbb{E}_{(x,y),(x',y')\sim\mathbb{D}\otimes\mathbb{D}} \left[ g(x,y)g(x',y')K^2(x,x') \right],
\end{align*}
$$

$$
\begin{align*}
\|\hat{s}\|_{HS}^2 &= \frac{1}{n^2} \left\| \sum_{i=1}^{n} g(x_i, y_i) K(x_i, \cdot) \otimes K(x_i, \cdot) \right\|_{HS}^2 \\
&= \frac{1}{n^2} \left\langle \sum_{i=1}^{n} g(x_i, y_i) K(x_i, \cdot) \otimes K(x_i, \cdot), \sum_{j=1}^{n} g(x_j, y_j) K(x_j, \cdot) \otimes K(x_j, \cdot) \right\rangle_{HS} \\
&= \frac{1}{n^2} \sum_{i,j=1}^{n} g(x_i, y_i) g(x_j, y_j) \left\langle K(x_i, \cdot) \otimes K(x_i, \cdot), K(x_j, \cdot) \otimes K(x_j, \cdot) \right\rangle_{HS} \\
&= \frac{1}{n^2} \sum_{i,j=1}^{n} g(x_i, y_i) g(x_j, y_j) K^2(x_i, x_j),
\end{align*}
$$

$$
\begin{align*}
\langle s, \hat{s} \rangle_{\text{HS}} &= \left\langle \mathbb{E}_{(x,y)\sim\mathbb{D}} [g(x,y)K(x,\cdot) \otimes K(x,\cdot)], \frac{1}{n}\sum_{i=1}^{n} g(x_i, y_i)K(x_i,\cdot) \otimes K(x_i,\cdot) \right\rangle_{\text{HS}} \\
&= \frac{1}{n}\sum_{i=1}^{n} \mathbb{E}_{(x,y)\sim\mathbb{D}} [g(x_i, y_i)g(x,y) \langle K(x,\cdot) \otimes K(x,\cdot), K(x_i,\cdot) \otimes K(x_i,\cdot) \rangle_{\text{HS}}] \\
&= \frac{1}{n}\sum_{i=1}^{n} \mathbb{E}_{(x,y)\sim\mathbb{D}} [g(x_i, y_i)g(x,y) K^2(x, x_i)].
\end{align*}
$$

After substituting the obtained results into Equation (28) and rearranging, we obtain:

$$
\begin{align*}
\|A\|_{HS}^2 &= \mathbb{E}_{(x,y),(x',y')\sim\mathbb{D}\otimes\mathbb{D}} [g(x,y)g(x',y')K^2(x,x')] + \frac{1}{n^2} \sum_{i,j=1}^{n} g(x_i, y_i)g(x_j, y_j)K^2(x_i, x_j) \\
&\quad - \frac{2}{n} \sum_{i=1}^{n} \mathbb{E}_{(x,y)\sim\mathbb{D}} [g(x_i, y_i)g(x,y)K^2(x, x_i)].
\end{align*}
$$

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

# I Details on Experiments and Additional Numerical Results

In this section, we provide details on the experimental setting used to obtain Figure 2 and include additional numerical results in Appendix 1.5. We recall the formulation of the instrumental variable regression problem introduced in Section 2.2:

$$
\begin{align*}
\min_{\omega \in \mathbb{R}^d} \mathcal{F}(\omega) &:= L_{\text{out}}(\omega, h_{\omega}^{*}) = \frac{1}{2} \mathbb{E}_{(x,y)\sim\mathbb{Q}} \left[ |h_{\omega}^{*}(x) - y|^2 \right] \\
\text{s.t.} \quad & h_{\omega}^{*} = \arg \min_{h \in \mathcal{H}} L_{\text{in}}(\omega, h) = \frac{1}{2} \mathbb{E}_{(x,t)\sim\mathbb{P}} \left[ |h(x) - \omega^\top \phi(t)|^2 \right] + \frac{\lambda}{2} \|h\|_{\mathcal{H}}^2,
\end{align*}
$$

where $\phi(t) = (\phi_1(t), \dots, \phi_d(t))^{\top} \in \mathbb{R}^d$ is the feature map. We begin by deriving a closed-form expression for $\hat{h}_{\omega}$ (the empirical counterpart of $h_{\omega}^{*}$), which is key to obtaining closed-form expressions for $\hat{\mathcal{F}}(\omega)$ and $\overline{\nabla} \mathcal{F}(\omega)$, and thus accurate approximations of $\mathcal{F}(\omega)$ and $\nabla \mathcal{F}(\omega)$.