# On the Generalization Capacities of Neural Controlled Differential Equations

Linus Bleistein, Agathe Guilloux

► To cite this version:

Linus Bleistein, Agathe Guilloux. On the Generalization Capacities of Neural Controlled Differential Equations. Workshop on New Frontiers in Learning, Control, and Dynamical Systems at the International Conference on Machine Learning (ICML 2023), Jul 2023, Honolulu (Hawai), United States. hal-04876009

HAL Id: hal-04876009

https://hal.science/hal-04876009v1

Submitted on 9 Jan 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>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Generalization Capacities of Neural Controlled Differential Equations

Linus Bleistein ¹² Agathe Guilloux <fnref n="1" />

## Abstract

We consider a supervised learning setup in which the goal is to predict an outcome from a sample of irregularly sampled time series using Neural Controlled Differential Equations (Kidger et al., 2020). In our framework, the time series is a discretization of an unobserved continuous path, and the outcome depends on this path through a controlled differential equation with unknown vector field. Learning with discrete data thus induces a discretization bias, which we precisely quantify. Using theoretical results on the continuity of the flow of controlled differential equations, we show that the approximation bias is directly related to the approximation error of a Lipschitz function defining the generative model by a shallow neural network. By combining these result with recent work linking the Lipschitz constant of neural networks to their generalization capacities, we upper bound the generalization gap between the expected loss attained by the empirical risk minimizer and the expected loss of the true predictor.

## 1. Introduction

Time series are ubiquitous in many domains such as finance, agriculture, economics and healthcare. A common set of tasks consists in predicting an outcome $y \in \mathcal{Y}$, such as a scalar or a label, from a time-evolving set of features. This problem has been addressed with a great variety of methods, ranging from auto-regressive models, such as VAR, to deep learning models, such as Recurrent Neural Networks (RNN), Long-Short-Term-Memory Networks (LSTM) and many others. It has also been thoroughly studied through the lenses of stochastic processes, Gaussian processes and many more.

Both in practice and theory, most methods and theoretical setups often only tackle the case of regularly sampled time series. In this setup, a time series is seen as a collection $x = (x_{t_0}, \dots, x_{t_K})$ of $K+1$ datapoints in $\mathbb{R}^d$ for which $\Delta t := t_j - t_{j-1}$ is constant. In real-world scenarios, however, time series are often irregularly sampled. This is often the case when data collection is difficult or expensive. From a modelling perspective, it is natural to consider that the time series $x$ is a degraded version - through subsampling or missing values - of an unobserved underlying path $x = (x_t)_{t \in [0, 1]}$ which, in turn, determines the outcome $y \in \mathcal{Y} \subset \mathbb{R}$. Informally, we have that $y = F(x) + \varepsilon$, where $F$ is an operator that maps the space of paths $\{x : [0, 1] \to \mathbb{R}^d\}$ to $\mathcal{Y}$, and $\varepsilon$ is a bounded noise term. For instance, in healthcare, the value of a biomarker of interest of a patient is determined by the continuous and unobserved trajectory of her vitals, rather than by the discrete measurements made by a physician. Similarly, in agriculture, crop growth and yield is not determined by the punctual measurements of soil fertility, but by its continuous value through time.

As a consequence, learning from the discrete time series $x$ rather than from $x$ introduces a discretization bias. It is of great importance to better understand how this bias degrades the performance of learning algorithms. We tackle this question by studying how irregular sampling affects the generalization capacities of Neural Controlled Differential Equations (NCDE), a popular and state-of-the-art method for learning from irregular time dependent data introduced in the seminal work of Kidger et al. (2020).

**Out setup.** We restrict our attention to a supervised learning setup in which sampling times are irregularly spaced. Consider a sample $\{(y^i, x^{D,i})\}$, where $y^i$ is the label of data discretized on $D = \{t_0, \dots, t_K\}$. We will always require that this sampling grid is constant across individuals, and assume that $t_0 = 0$ and $t_K = 1$. Consider a predictor $\hat{f}^D$ obtained by empirical risk minimization on this sample such that

$$ \hat{f}^D(\mathbf{x}^D) \approx y. $$

Our goal is to upper bound the generalization gap

$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(\mathbf{f}^*) $$

where $\mathcal{R}^D(\hat{f}^D)$ refers to the generalization error of the learn predictor on discretized data, and $\mathcal{R}(\mathbf{f}^*)$ is the gener-

<footnote n="1">¹Inria Paris, F-75015 Paris, France ²LaMME, UEVE and UMR 8071, Paris Saclay University, F-91042, Evry, France. Correspondence to: Linus Bleistein <linus.bleistein@inria.fr>.</footnote>

<footnote>Workshop on New Frontiers in Learning, Control, and Dynamical Systems at the International Conference on Machine Learning (ICML), Honolulu, Hawaii, USA, 2023. Copyright 2023 by the author(s).</footnote>
alization error of the true predictor defining our generative model on continuous data (our framework is introduced in greater detail later on). Upper bounding this quantity crucially allows to quantify how sampling affects prediction performance. This is a central question both from a theoretical perspective and for practitioners, since it could, for instance, justify increasing sampling frequencies of ICU patients (Johnson et al., 2020) or intensifying monitoring of people suffering from sickle-cell anaemia (Bussy et al., 2019).

**Contribution.** Our contribution is threefold. First, building on a generative model proposed by Bleistein et al. (2023), we obtain an upper bound of the generalization gap which decomposes into two worst-case bounds, an approximation bias and a discretization bias. In a second time, we leverage the continuity of the flow of a CDE and approximation theorems for shallow neural networks to bound both the approximation and the discretization bias. Finally, by using the Lipschitz continuity of NCDE, we obtain an upper bound on the covering number of continuous and discretized NCDE, which explicitly depends on the discretization for the second class. This allows us to control both worst-case bounds and yields an upper bound on the generalization gap. As a byproduct, we also obtain a generalization bound for NCDEs.

**Overview.** Section 2 covers related works. Section 3 details our setup and assumptions. In Section 4, we state our main result, namely a generalization bound for NCDE and an upper bound on the generalization gap for a very general class of losses. Sections 5 and 6 give the two building blocks of this proof, namely upper bounds on the discretization and approximation biases, and upper bounds on the Rademacher complexity of our model.

## 2. Related works

**Hybridization of deep learning and differential equations.** The idea of combining differential systems and deep learning algorithms has been around for many years (Rico-Martinez et al., 1992; 1994). The recent work of Chen et al. (2018) has revived this idea by proposing a model which learns a representation of $x \in \mathbb{R}^d$ by using it to set the initial condition $z_0 = \varphi_{\xi}(x)$ of the ODE

$$dz_t = \mathbf{G}_{\psi}(z_t)dt$$

or equivalently in integral form

$$z_t = \varphi_{\xi}(x) + \int_{0}^{t} \mathbf{G}_{\psi}(z_s)ds$$

where $\mathbf{G}_{\psi}$ and $\varphi_{\xi}$ are neural networks. The terminal value $z_1$ of this ODE is then used as input of any classical machine

learning algorithm. Numeral contributions have built upon this idea and studied the connection between deep learning and differential equations in the recent years (Dupont et al., 2019; Chen et al., 2019b; 2020; Finlay et al., 2020). Mas-saroli et al. (2020) and Kidger (2022) offer comprehensive introductions to this topic. Links between deep learning and stochastic differential equations (SDE) have also been recently studied (Cohen et al., 2021; Marion et al., 2022; Hayou, 2022).

**Learning with irregular time series.** In the recent years, numerous competitive models have been proposed to handle time-dependent irregular data. Che et al. (2018) introduce a modified GRU with learn exponential decay of the hidden state between sampling times. Other models hybridizing classical deep learning architectures and ODE include GRU-ODE (De Brouwer et al., 2019) and RNN-ODE (Rubanova et al., 2019). Controlled Differential Equations have been introduced through two distinct frameworks, namely through Neural Rough Differential Equations (Morrill et al., 2021) and NCDEs (Kidger et al., 2020). The latter extends neural ODEs to sequential data by first interpolating a time series $\mathbf{x}^D$ with cubic splines. The first value of the time series and the interpolated path $\tilde{\mathbf{x}}$ are then used as initial condition and driving signal of the CDE

$$dz_t = \mathbf{G}_{\psi}(z_t)d\tilde{\mathbf{x}}_t$$

or equivalently

$$z_t = \varphi_{\xi}(x_{t_0}) + \int \mathbf{G}_{\psi}(z_s)d\tilde{\mathbf{x}}_s,$$

where the integral is to be understood as the Riemann-Stieltjes integral. As for neural ODEs, the terminal value of this NCDE is then used for classification or regression. Other methods for learning from irregular time series include Gaussian Processes (Li & Marlin, 2016) and the signature transform (Chevyrev & Kormilitzin, 2016; Kidger et al., 2019; Fermanian, 2021; Bleistein et al., 2023).

**Generalization bounds.** The generalization capacities of recurrent models have been studied since the end of the 1990s, mainly through their VC dimension (Dasgupta & Sontag, 1995; Koiran & Sontag, 1998). Bartlett et al. (2017) sparked a line of research connecting the Lipschitz constant of feedforward neural networks to their generalization capacities. Our central reference is the work of Che et al. (2018), which leverage the Lipschitz continuity of RNN, LSTM and GRU to derive generalization bounds of these models. This work assumes that the time series $\mathbf{x}^D$ is bounded and regularly sampled, which are unnecessary assumptions in our work. Fermanian et al. (2021) also obtain generalization bounds for RNN and a particular class of NCDE. However, they crucially impose restrictive regularity conditions which
then allow for linearization of the NCDE in the signature space. Generalization bounds are then derived using kernel learning theory. Hanson & Raginsky (2022) considers an encoder-decoder setup involving excess risk bounds on neural ODEs. Closest to our work is the recent work by Marion (2023), which leverages similar tools to study the generalization capacities of a particular class of neural ODE and deep ResNets. However, our work considers NCDE, which learn with time dependent data, and focuses on the effects of data irregularity on the generalization gap.

### 3. Setup

We first detail our generative model, before introducing NCDE and our learning problem.

#### 3.1. A CDE-based generative model

In order to fully track the approximation bias, we introduce a model that links the underlying path $x$ to the outcome $y$.

**Assumption 3.1.** There exists a vector field $\mathbf{G}^* : \mathbb{R}^p \to \mathbb{R}^{p \times d}$, a function $\varphi : \mathbb{R}^d \to \mathbb{R}^p$ and a vector $\alpha_* \in \mathbb{R}^p$ that governs a latent state $(z_t^*)_{t \in [0,1]}$ through the CDE

$$dz_t^* = \mathbf{G}^*(z_t^*)dx_t \quad (1)$$

with initial value $z_0^* = \varphi(x_0)$ such that the outcome $y$ is given by

$$y = \alpha_*^\top z_1^* + \varepsilon,$$

where $\varepsilon$ is a noise term bounded by $M_\varepsilon > 0$.

Put otherwise, the outcome $y$ associated to a path $x$ is the value of a linear transformation applied to the final value of a CDE with unknown vector field $\mathbf{G}^*$. This is a very general model, which encompasses ODE-based models used for instance in biology or physics. Indeed, by setting $x_t = t$ for all $t$, one recovers a generic ODE.

We call the map $\mathbf{f}^* : x \mapsto \alpha_*^\top z_1^*$ the true predictor. Uniqueness and existence to the CDE (1) is given by the Picard-Lindelöf theorem given in Appendix A.2. In line with this theorem, we make the following assumptions on the generative model and on the path $(x_t)_{t \in [0,1]}$.

**Assumption 3.2.** The parameters of the generative model satisfy $\|\alpha_*\| \le B_\alpha$, $\varphi \in C(\mathbb{R}^d, \mathbb{R}^p)$ and $\mathbf{G}^* \in \text{Lip}(\mathbb{R}^p, \mathbb{R}^{p \times d})$ with Lipschitz constant $L_{\mathbf{G}^*}$.

**Assumption 3.3.** The path $(x_t)_{t \in [0,1]}$ is $L_x$-Lipschitz continuous, that is $\|x_t - x_s\| \le L_x |t - s|$ for all $t, s \in [0, 1]$.

Assumption 3.3 crucially implies that for all $t \in [0, 1]$,

$$\|x\|_{1-\text{var},[0,t]} := \sup_D \sum \|x_{t_{i+1}} - x_{t_i}\| \le L_x t \le L_x,$$

where the supremum is taken on all finite discretizations $D = \{0 = t_0 < t_1 < \cdots < t_N = 1\}$ of $[0, t]$, for $N \in \mathbb{N}^*$.

We make a last assumption on the starting point of the considered time series, which we will always assume to be true when speaking of paths in the following.

**Assumption 3.4.** There exists a constant $B_x > 0$ such that $\|x_0\| \le B_x$.

#### 3.2. Neural Controlled Differential Equations

Neural Controlled Differential Equations, introduced by Kidger et al. (2020), are a particular form of CDE in which the vector field is chosen as a neural network. They can be seen as a natural extension of neural ODEs (Chen et al., 2018) to sequential data. Informally, NCDEs learn representations of time series by using them as drivers of controlled differential equations, in combination with a neural vector field. This vector field is optimized with respect to the considered learning task as in Neural ODEs. They bear close resemblance with ResNets and Recurrent Neural Networks (Fermanian et al., 2021).

In our setup, the outcome $y$ is determined by the endpoint of an unknown CDE. Our goal is to approximate the dynamics of this CDE using a NCDE, such that the endpoint of this neural model matches the outcome.

**Definition 3.5.** Let $(x_t)_{t \in [0,1]}$ be a continuous path of bounded variation and let $\mathbf{G}_\psi : \mathbb{R}^p \to \mathbb{R}^{p \times d}$ be a neural vector field parametrized by $\psi$. Consider the solution $z := (z_t)_{t \in [0,1]}$ of the controlled differential equation

$$dz_t = \mathbf{G}_\psi(z_t)dx_t,$$

with initial condition $z_0 = \varphi_\xi(x_0)$, where $\varphi_\xi : \mathbb{R}^d \to \mathbb{R}^d$ is a neural network parametrized by $\xi$. We call the latent space trajectory and

$$\Phi^\top z_1 =: f_\theta(x)$$

the prediction of the NCDE with parameters $\theta = (\Phi, \psi, \xi)$.

The vector field $\mathbf{G}_\psi$ can be any common neural network, since these architectures are Lipschitz continuous (Virmaux & Scaman, 2018). This assures that the solution of the NCDE is well defined. We restrict our attention to neural vector fields of the form

$$\mathbf{G}_{\psi}(z) = \sigma(\mathbf{A}z + \mathbf{b}), \quad (2)$$

with $L_\sigma$-Lipschitz activation function $\sigma : \mathbb{R} \to \mathbb{R}$ that verifies $\sigma(0) = 0$. $\mathbf{A} : \mathbb{R}^p \to \mathbb{R}^{p \times d}$ is a linear operator and $\mathbf{b} \in \mathbb{R}^{p \times d}$ a matrix. The activation function is evaluated entry-wise. We also restrict ourselves to initializations of the form

$$z_0 = \varphi_\xi(x_0) = \sigma(\mathbf{U}x_0 + v), \quad (3)$$

where $\mathbf{U} \in \mathbb{R}^{p \times d}$ and $v \in \mathbb{R}^p$. We use the notation

$$\text{NN}_{\mathbf{U},v} : u \mapsto \sigma(\mathbf{U}u + v)$$
to refer to this initialization network. The activation is assumed to be identical for the initialization layer and the neural vector field, but such an assumption can be relaxed. The learnable parameters of the model thus are $\theta = (\Phi, \mathbf{A}, \mathbf{b}, \mathbf{U}, v)$.

The restriction on the depth of the neural vector field is mainly made for the sake of simplicity. Indeed, our proof only relies on the Lipschitz constant of $\mathbf{G}_{\psi}$. Such a constant can be upper bounded for deeper neural vector fields at the price of heavier notations. We refer to Bartlett et al. (2017) and Virmaux & Scaman (2018) for an in-depth discussion of the Lipschitz continuity of deeper neural networks. Also note that in practice, the vector fields used in real-world applications of NCDE have typically only a few hidden layers (Kidger et al., 2020).

**Learning with time series.** If one has access to the continuous path $x$ the time series $\mathbf{x}^D$ is sampled from, one can directly use this path as an input of a NCDE. However, to apply NCDEs to time series $\mathbf{x}^D$, one needs to embed $\mathbf{x}^D$ into the space of paths of bounded variation. This can be done through any reasonable embedding mapping

$$ \rho : \mathbf{x}^D \in (\mathbb{R}^d)^{K+1} \mapsto (\tilde{\mathbf{x}}_t)_{t \in [0,1]} $$

such as splines, polynomials or linear interpolation.

In this work, we focus on the fill-forward embedding, which simply defines the value of $\tilde{x}$ between two consecutive points of $D$ as the last observed value of $\mathbf{x}^D$.

**Definition 3.6.** The fill-forward embedding of a time series $\mathbf{x}^D = (x_{t_0}, \dots, x_{t_K})$ sampled on $D = \{0 = t_0 < \dots < t_K = 1\}$ is the piecewise constant path $(\tilde{\mathbf{x}}_t)_{t \in [0,1]}$ defined for every $t \in [t_k, t_{k+1}$ [as $\tilde{\mathbf{x}}_t = \mathbf{x}^D_{t_k}$ for $k \in \{0, \dots, K-1\}$, and $\tilde{\mathbf{x}}_1 = x_{t_K} = x_1$.

Using this embedding in a NCDE with parameters $\theta = (\Phi, \mathbf{A}, \mathbf{b}, \mathbf{U}, v)$ recovers a piecewise constant latent space trajectory $z^D := (z_t^D)_{t \in [0,1]}$ recursively defined for every $t \in [t_k, t_{k+1}$ [by

$$ z_t^D = z_{t_{k-1}}^D + \sigma(\mathbf{A}z_{t_{k-1}} + \mathbf{b})(x_{t_k} - x_{t_{k-1}}) \quad (4) $$

for $k = 1, \dots, K-1$ and initialized as $z_0^D := \sigma(\mathbf{U}x_0 + v)$. The terminal value is equal to

$$ z_1 = z_{t_K} = z_{t_{K-1}} + \sigma(\mathbf{A}z_{t_{K-1}} + \mathbf{b})(x_{t_K} - x_{t_{K-1}}). $$

Formally, the prediction $\Phi^\top z_1^D$ is equal to

$$ \Phi^\top z_1^D = f_\theta \circ \rho_{\text{FF}}(\mathbf{x}^D), $$

where $\rho_{\text{FF}}$ is the fill-forward operator. To lighten notations, we simply write $\Phi^\top z_1^D = f_\theta(\mathbf{x}^D)$.

We highlight that the latent state $z$ is more informative about the outcome $y$ than $z^D$. Indeed, it embeds the full trajectory

of $x$, which determines $y$ through the CDE (1), while $z^D$ embeds a version of $x$ degraded through sampling.

This recursive architecture has been studied with random $\mathbf{A}, \mathbf{b}$ by Cirone et al. (2023) under the name of homogenous controlled ResNet because of its resemblance with the popular ResNet (He et al., 2015).

**Restrictions on the parameter space.** In order to obtain generalization bounds, we must furthermore restrict the size of parameter space by requiring that the parameters $\theta$ lie in a bounded set $\Theta$. This means that there exist $B_{\mathbf{A}}, B_{\mathbf{b}}, B_{\mathbf{U}}, B_{\mathbf{v}}$ such that

$$ \|\mathbf{A}\| \le B_{\mathbf{A}}, \|\mathbf{b}\| \le B_{\mathbf{b}}, \|\mathbf{U}\| \le B_{\mathbf{U}}, \|\mathbf{v}\| \le B_{\mathbf{v}} \quad (5) $$

and

$$ \|\Phi\| \le B_{\alpha} \quad (6) $$

for all NCDEs considered, where $\|W\|$ is the Frobenius norm. We recall that $B_\alpha$ is defined in Assumption 3.2. Such an restriction is classical for deriving generalization bounds (Bartlett et al., 2017; Bach, 2021; Fermanian et al., 2021). We call

$$ \mathcal{F}_{\Theta} = \{ f_{\theta} : (x_t)_{t \in [0,1]} \mapsto f_{\theta}(x) \in \mathbb{R} \text{ s.t. } \theta \in \Theta \} $$

this class of predictors.

We now state an important lemma. It is a direct consequence of Gronwall's Lemma, stated in Lemma A.4, and of our generative model and the restrictions on the parameter space. First define

$$ \|\mathbf{G}^*(0)\|_{\text{op}} := \max_{\|u\|=1} \|\mathbf{G}^*(0)u\|, $$

which is finite since $\mathbf{G}^*$ is continuous.

**Lemma 3.7.** Let $y$ be generated from the CDE (1) from a $L_x$-Lipschitz path $x$. One has

$$ |y| \le B_{\alpha} \left( B_{\varphi} + \|\mathbf{G}^*(0)\|_{\text{op}} L_x \right) \exp \left( L_{\mathbf{G}^*} L_x \right) + M_{\varepsilon}, $$

where $B_\varphi := \max_{\|u\| \le B_x} \|\varphi(u)\|$. Also, both $f_\theta(x)$ and $f_\theta(\mathbf{x}^D)$ are upper bounded by

$$ M_{\Theta} := B_{\alpha} L_{\sigma} \exp(L_{\mathbf{A}} L_{\sigma} L_x) (B_{\mathbf{U}} B_x + B_{\mathbf{v}} + B_{\mathbf{b}} L_x) $$

for all $f_\theta \in \mathcal{F}_{\Theta}$.

This lemma ensures that the predictions as well as the outcome of our generative model remain bounded. This, in turn, means that the loss function is bounded, which is crucial to leverage results from Bartlett et al. (2017).
### 3.3. The learning problem

We now detail our learning setup. We consider an i.i.d. sample $\{(y^i, \mathbf{x}^{D,i})\}_{i=1}^n \sim y, x$ with given discretization $D$.

For a given predictor $f_\theta \in \mathcal{F}_\Theta$, define

$$ \mathcal{R}_n^D(f_\theta) = \frac{1}{n} \sum_{i=1}^n \ell(y^i, f_\theta(\mathbf{x}^{D,i})) $$

and

$$ \mathcal{R}^D(f_\theta) = \mathbb{E}_{x,y} [\ell(y, f_\theta(\mathbf{x}^D))] $$

as the empirical risk and the expected risk on the discretized data. Similarly, define

$$ \mathcal{R}_n(f_\theta) = \frac{1}{n} \sum_{i=1}^n \ell(y^i, f_\theta(x^i)) $$

and

$$ \mathcal{R}(f_\theta) = \mathbb{E}_{x,y} [\ell(y, f_\theta(x))] $$

the empirical risk and expected risk on the continuous data. We stress that $\mathcal{R}_n(f_\theta)$ cannot be optimized, since we do not have access to the continuous data. Let

$$ \hat{f}^D \in \arg \min_{\theta \in \Theta} \mathcal{R}_n^D(f_\theta) $$

be an optimal predictor obtained by empirical risk minimization on the discretized data. In order to obtain generalization bounds, the following technical assumption on the loss are necessary (Mohri et al., 2018).

**Assumption 3.8.** The loss $\ell: \mathbb{R} \times \mathbb{R} \to \mathbb{R}_+$ is Lipschitz with respect to its second variable, that is there exists $L_\ell$ such that for all $u, u' \in \mathcal{Y}$ and $y \in \mathcal{Y}$,

$$ |\ell(y, u) - \ell(y, u')| \le L_\ell |u - u'|. $$

This hypothesis is satisfied for most classical losses, such as the mean squared error, as long as the outcome and the predictions are bounded. This is the case in our setup, as stated in Lemma 3.7.

### 4. A bound on the generalization gap

We first decompose the difference between the expected risk of the learnt predictor $\hat{f}^D$ learn from the sample $S$ and the expected risk of the true predictor $\mathbf{f}^*$.

**Lemma 4.1.** For all $f_\theta \in \mathcal{F}_\theta$, the generalization gap

$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(\mathbf{f}^*) $$

is almost surely bounded from above by

$$ \sup_{g \in \mathcal{F}_\Theta} \left[ \mathcal{R}^D(g) - \mathcal{R}_n^D(g) \right] + \sup_{g \in \mathcal{F}_\Theta} \left[ \mathcal{R}(g) - \mathcal{R}_n(g) \right] + \mathcal{R}_n^D(f_\theta) - \mathcal{R}_n(f_\theta) + \mathcal{R}(f_\theta) - \mathcal{R}(\mathbf{f}^*). $$

We briefly comment on this inequality. The two first suprema are common worst-case type bounds and can be controlled using the Rademacher complexity of the class $\mathcal{F}_\Theta$. The third term $\mathcal{R}_n^D(f_\theta) - \mathcal{R}_n(f_\theta)$ corresponds to a discretization bias, which arises because one can only minimize the empirical risk on the discretized data, and not on the continuous data. It will be bounded using the continuity of the flow of a CDE with respect to its driving signal. The last term is an approximation bias.

We let $M_\ell$ be a bound on the loss function, which is finite since the loss is continuous and the values of $y$ and $f_\theta(x)$ are bounded as stated in Lemma 3.7. We also let

$$ |D| := \max_{j=1,...,K-1} |t_{j+1} - t_j| $$

be the greatest gap between two sampling times. The following theorems are our central results.

**Theorem 4.2.** One has, with probability at least $1 - \delta$, that

$$ \mathcal{R}^D(\hat{f}^D) \le \mathcal{R}_n(\hat{f}^D) + \frac{24M_\Theta L_\ell}{\sqrt{n}} \sqrt{2pU_1 + dp(p+2)U_2} + M_\ell \sqrt{\frac{\log 1/\delta}{2n}} $$

with $U_1 := \log 20 \sqrt{n} p K_1$ and $U_2 := \log 20 \sqrt{n} d p K_2$ and $K_1, K_2$ two constants depending on $\Theta, L_x, B_x$ and $L_\sigma$.

Turning to the generalization gap, we get the following result.

**Theorem 4.3.** The generalization gap

$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(\mathbf{f}^*) $$

is bounded from above for any $f_\theta \in \mathcal{F}_\Theta$ with parameters $\theta = (\Phi, \psi, \mathbf{U}, v)$ such that $\Phi = \alpha$ by

$$ \frac{16}{n} + \frac{96M_\Theta L_\ell}{\sqrt{n}} \sqrt{2pU_1 + dp(p+2)U_2} + L_\ell C_2(\Theta, L_x, L_\sigma) |D| + L_\ell B_\alpha \exp(L_G \cdot L_x) \left[ L_x \|G_\psi - G^*\|_{\infty, \Omega(\Theta, L_x, B_x)} + \|\varphi - \text{NN}_{\mathbf{U},v}\|_{\infty, B_x} \right]. $$

The constant $M_\Theta$ is given in Lemma 3.7. The constant $C_2(\Theta, L_x, L_\sigma)$ is given in Proposition 5.1. The quantity $\|G_\psi - G^*\|_{\infty, \Omega(\Theta, L_x, B_x)}$ is defined as

$$ \max_{u \in \Omega(\Theta, L_x, B_x)} \|G_\psi(u) - G^*(u)\|, $$

and

$$ \|\varphi - \text{NN}_{\mathbf{U},v}\|_{\infty, B_x} := \max_{\|u\| \le B_x} \|\varphi(u) - \text{NN}_{\mathbf{U},v}(u)\|. $$
The diameter of $\Omega(\Theta, L_x, B_x)$ is upper bounded in Proposition 5.2.

Not assuming that $\Phi = \alpha$ introduces a supplementary source of bias directly proportional to the difference $\|\Phi - \alpha\|$. However, since we assume that $\alpha$ and $\Phi$ lie in a common ball of radius $B_\alpha$, we make this assumption which simplifies the presentation of our results.

Let us comment the three terms of the bound of Theorem 4.3.

1. The first term is a common worst case bound when deriving generalization bounds. It depends on the Rademacher complexity of $\mathcal{F}_\Theta$, and not on the discretization. We show in our proofs that an identical upper bounds hold for discrete and continuous input.

2. The second term corresponds to the discretization bias, and is directly proportional to $|D|$. It therefor vanishes at the same speed than $|D|$. For instance, if we consider the sequence $D_K$ of equidistant discretizations of $[0, 1]$ with $K$ points, the discretization bias vanishes at linear speed.

3. Finally, observe that the approximation bias writes as the sum of approximation errors on the vector fields and the initial condition, rather than an general approximation error on the true predictor $\mathbf{f}^*$.

**Outline of the proof.** The proof schematically works in two times. We first bound the two sources of bias by leveraging the continuity of the flow of a CDE. Informally, this theorem states that the difference $\Delta(\text{terminal value})$ between the terminal value of two CDEs decomposes as

$$ \Delta(\text{vector fields}) + \Delta(\text{initial conditions}) + \Delta(\text{driving paths}). $$

This neat decomposition allows us to bound the discretization bias since it depends on the terminal value of two CDEs whose driving paths differ. This is done in Proposition 5.1. It also allows to control the approximation bias, which depends on the terminal value of two CDEs with identical driving path but whose vector fields and initial conditions differ, which we do in Proposition 5.2. Continuity of the flow is stated in full generality in Theorem A.7. This property is illustrated in Figure 1.

In a second time, we proceed to show that NCDEs are Lipschitz with respect to their parameters $\theta$. Building on Bartlett et al. (2017), we first control the covering number of $\mathcal{F}_\Theta$ in Proposition 6.1. This gives us a bound on the Rademacher complexity of class of predictors, stated in 6.2. Theorem 4.2 then follows using standard arguments from Mohri et al. (2018). Combining these results with bounds on both sources of bias gives Theorem 4.3.

# 5. Bounding the discretization and the approximation bias

We bound the two sources of bias, namely the approximation bias and the discretization bias.

**Proposition 5.1.** For any $f_\theta \in \mathcal{F}_\Theta$, one has

$$ \mathcal{R}_n^D(f_\theta) - \mathcal{R}_n(f_\theta) \le L_\ell C_2(\Theta, L_x, L_\sigma)|D|, \\ \text{where } C_2(\Theta, L_x, L_\sigma) \text{ is equal to} $$

$$ B_\alpha \left[ L_\sigma B_{\mathbf{A}} \left( L_\sigma (B_{\mathbf{U}} B_x + B_v) + L_\sigma B_{\mathbf{b}} L_x \right) \times \exp(L_\sigma B_{\mathbf{A}} L_x) + L_\sigma B_{\mathbf{b}} \right] \exp(L_\sigma B_{\mathbf{A}} L_x). $$

The proof is given in Appendix B.1. Turning to the approximation bias, we get the following bound, using the same technique as for bounding the discretization bias.

**Proposition 5.2.** For any $f_\theta \in \mathcal{F}_\Theta$ with parameters $\theta = (\Phi, \psi, \mathbf{U}, v)$ such that $\Phi = \alpha$, the approximation bias $\mathcal{R}(f_\theta) - \mathcal{R}(\mathbf{f}^*)$ is bounded from above by

$$ L_\ell B_\alpha \exp(L_{\mathbf{G}^*} L_x) L_x \|\mathbf{G}_\psi - \mathbf{G}^*\|_{\infty, \Omega(\Theta, L_x, B_x)} \\ + L_\ell B_\alpha \exp(L_{\mathbf{G}^*} L_x) \|\varphi - NN_{\mathbf{U},v}\|_{\infty, B_x}, $$

where

$$ \Omega(\Theta, L_x, B_x) \subset \{u \in \mathbb{R}^p \mid \|u\| \le (L_\sigma(B_{\mathbf{U}}B_x + B_v) \\ + L_\sigma B_{\mathbf{b}}L_x) \exp(L_\sigma L_{\mathbf{A}}L_x)\}. $$

Notice that this bound has a double dependence on the unknown vector field $\mathbf{G}^*$. The first one is through its Lipschitz constant $L_{\mathbf{G}^*}$. The second one comes from the error made when approximating $\mathbf{G}^*$ by $\mathbf{G}_\psi$ and $\varphi$ by the shallow initialization neural network. Controlling this second term using approximation results for neural networks is left for future work.

# 6. Bounding the Rademacher complexity of NCDE

We first upper bound the covering number of $\mathcal{F}_\Theta$. This result is obtained by showing that NCDE are Lipschitz with respect to their parameters, such that covering each parameter class yields a covering of the whole class of predictors.

**Proposition 6.1.** The covering number $\mathcal{N}(\mathcal{F}_\Theta, \beta)$ of $\mathcal{F}_\Theta$ is bounded from above by

$$ \left(1 + \frac{10\sqrt{p}K_1}{\beta}\right)^{2p} \left(1 + \frac{10\sqrt{dp}K_3}{\beta}\right)^{dp(p+2)} $$

where $K_1, K_2$ are two constants depending on $\Theta, L_x, B_x$ and $L_\sigma$.
<figure><img src="image_2.png" /><figcaption>Figure 1: On the **left**, we plot in bold the solutions of two NCDEs who only differ in their vector field $A_1$ and $A_2$. We then plot the solutions of the NCDEs with interpolated vector field $\delta A_1 + (1 - \delta)A_2$ for $\delta \in [0, 1]$ in red-blue gradient. In the **middle**, we consider a given NCDE and interpolate linearly between two initial conditions $x_0^1$ and $x_0^2$. On the **right**, we consider a given NCDE and drive it with the linear interpolation of two paths $(x_t^1)_t$ and $(x_t^2)_t$. In all three cases, the solutions evolve continuously as we interpolate between the models.</figcaption></figure>

In combination with arguments borrowed from Chen et al. (2019a) and Bartlett et al. (2017), Proposition 6.1 allows to upper bound the empirical Rademacher complexity, which we denote by $\text{Rad}(\mathcal{H})$ for a function class $\mathcal{H}$.

**Proposition 6.2.** The empirical Rademacher complexity $\text{Rad}(\mathcal{F}_{\Theta})$ associated to an i.i.d. sample of size $n$ is upper bounded by

$$ \frac{4}{n} + \frac{24M_{\Theta}}{\sqrt{n}} \sqrt{2pU_1 + dp(p+2)U_2} $$

where $U_1 := \log 20\sqrt{npK_1}$ and $U_2 := \log 20\sqrt{ndpK_2}$, and $K_1, K_2$ are two constants depending on $\Theta, L_x, B_x$ and $L_\sigma$. The constant $M_\Theta$ is given in Lemma 3.7.

## 7. Experiments

Our bounds rely crucially on restrictions on our parameter space: the vaster the parameter space we consider - i.e. the bigger the norm of the model's parameters, the looser our bounds.

In a first experiment, we analyse how the norm of the NCDE's parameters evolve during training in a teacher-student setup. Paths are not downsampled in this experiment to neutralize the discretization bias. The target is generated from a teacher NCDE: the approximation bias is thus null. Figure 2 considers a single training run. First, despite only seeing the endpoint of the true red trajectory, the trained model achieves almost perfect interpolation of the whole trajectory on the training set. Secondly, the model's parameters stay within close range of their initialized value. Since our model is initialized with Pytorch's standard scaled Gaussian initialization, this means that the parameter's norms remain

<figure><img src="image_3.png" /><figcaption>Figure 2: On the **left**, we display the evolution of the latent state trajectory $(\alpha^\top z_t)_t$ and the true latent state trajectory $(\alpha^\top z_t^*)_t$ through training. On the **right**, we display the evolution of the parameter's norms normalized by their value at initialization. Training is performed with Adam (Kingma & Ba, 2014), standard $\alpha, \beta$ parameters and a learning rate of $9 \times 10^{-4}$ for 3000 iterations on a training set of size $n = 100$.</figcaption></figure>

well behaved during training and that our bounds do not explode.

We then iterate this experiment and display the distribution of the parameter's norms - this time in absolute value - after training in Figure 3. Since the norms of the parameters scale with their dimensions, greater values are observed for the parameters $A$ and $b$ who are of dimension $dp^2$ and $dp$ respectively. The parameters norm remain relatively small and do not disperse excessively over training runs.

Finally, we conduct a last experiment to analyse the interaction between training error, generalization error and discretization. We first fix a test and train dataset. For
<figure><img src="image_4.png" /><figcaption>Figure 3: Absolute value of the parameters norm after training for 1000 iterations on a training set of size $n = 100$. Training is performed with Adam (Kingma & Ba, 2014), standard $\alpha, \beta$ parameters and a learning rate of $9 \times 10^{-4}$. Each dot corresponds to one random initialization and training run of an NCDE.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Train and test error with paths downsampled to 5 points. Every point corresponds to a random downsampling of the paths $(x^i)$, initialization and training of an NCDE model. The colorbar indicates the coarsness of the sampling, that is the biggest time between two sampling points. Training is performed with Adam (Kingma & Ba, 2014), standard $\alpha, \beta$ parameters and a learning rate of $1 \times 10^{-3}$ for 1200 iterations.</figcaption></figure>

every random initialization and training of an NCDE, we randomly downsample the paths $(x^i)$ on a 5 point grid. We then compute the loss on the downsampled training and the downsampled test dataset. The results are displayed in Figure 4. The darker the points, the greater the sampling gap $|D|$. First, one can see that the training error is positively correlated with the generalization error, as predicted by Theorem 4.2. Interestingly, the sampling gap $|D|$ is also positively correlated with the training and generalization error. Indeed, as more information is missed between sampling times when $|D|$ increases, the prediction of the label becomes less precise.

## 8. Conclusion

Several perspectives are of high interest for extending our work. First, a deeper empirical assessment of our claims is needed. Concerning the approximation bias, we believe that approximation theorems for deeper neural networks can be used for quantifying its dependence in $\Theta$ and in the size of the latent space $p$, leading to a trade-off between complexity of the model and approximation capacities.

**Acknowledgements.** We thank Pierre Marion for insightful comments and clarifications, and Cristopher Salvi for discussions on Controlled ResNets. LB conducted this research while visiting MBZUAI's Machine Learning Department (UAE).
References

Bach, F. Learning theory from first principles. *Online version*, 2021.

Bartlett, P. L., Foster, D. J., and Telgarsky, M. J. Spectrally-normalized margin bounds for neural networks. *Advances in neural information processing systems*, 30, 2017.

Bleistein, L., Fermanian, A., Jannot, A.-S., and Guilloux, A. Learning the dynamics of sparsely observed interacting systems. *arXiv preprint arXiv:2301.11647*, 2023.

Bussy, S., Veil, R., Looten, V., Burgun, A., Gaïffas, S., Guilloux, A., Ranque, B., and Jannot, A.-S. Comparison of methods for early-readmission prediction in a high-dimensional heterogeneous covariates and time-to-event outcome framework. *BMC medical research methodology*, 19:1–9, 2019.

Che, Z., Purushotham, S., Cho, K., Sontag, D., and Liu, Y. Recurrent neural networks for multivariate time series with missing values. *Scientific reports*, 8(1):6085, 2018.

Chen, M., Li, X., and Zhao, T. On generalization bounds of a family of recurrent neural networks. *arXiv preprint arXiv:1910.12947*, 2019a.

Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. Neural ordinary differential equations. *Advances in neural information processing systems*, 31, 2018.

Chen, R. T., Behrmann, J., Duvenaud, D. K., and Jacobsen, J.-H. Residual flows for invertible generative modeling. *Advances in Neural Information Processing Systems*, 32, 2019b.

Chen, R. T., Amos, B., and Nickel, M. Learning neural event functions for ordinary differential equations. *arXiv preprint arXiv:2011.03902*, 2020.

Chevyrev, I. and Kormilitzin, A. A primer on the signature method in machine learning. *arXiv preprint arXiv:1603.03788*, 2016.

Cirone, N. M., Lemercier, M., and Salvi, C. Neural signature kernels as infinite-width-depth-limits of controlled resnets. *arXiv preprint arXiv:2303.17671*, 2023.

Clark, D. S. Short proof of a discrete gronwall inequality. *Discrete applied mathematics*, 16(3):279–281, 1987.

Cohen, A.-S., Cont, R., Rossier, A., and Xu, R. Scaling properties of deep residual networks. In *International Conference on Machine Learning*, pp. 2039–2048. PMLR, 2021.

Dasgupta, B. and Sontag, E. Sample complexity for learning recurrent perceptron mappings. *Advances in Neural Information Processing Systems*, 8, 1995.

De Brouwer, E., Simm, J., Arany, A., and Moreau, Y. Gru-ode-bayes: Continuous modeling of sporadically-observed time series. *Advances in neural information processing systems*, 32, 2019.

Dupont, E., Doucet, A., and Teh, Y. W. Augmented neural odes. *Advances in neural information processing systems*, 32, 2019.

Fermanian, A. Embedding and learning with signatures. *Computational Statistics & Data Analysis*, 157:107148, 2021.

Fermanian, A., Marion, P., Vert, J.-P., and Biau, G. Framing rnn as a kernel method: A neural ode approach. *Advances in Neural Information Processing Systems*, 34: 3121–3134, 2021.

Finlay, C., Jacobsen, J.-H., Nurbekyan, L., and Oberman, A. How to train your neural ode: the world of jacobian and kinetic regularization. In *International conference on machine learning*, pp. 3154–3164. PMLR, 2020.

Friz, P. K. and Victoir, N. B. *Multidimensional stochastic processes as rough paths: theory and applications*, volume 120. Cambridge University Press, 2010.

Hanson, J. and Raginsky, M. Fitting an immersed submanifold to data via sussmann’s orbit theorem. In *2022 IEEE 61st Conference on Decision and Control (CDC)*, pp. 5323–5328. IEEE, 2022.

Hayou, S. On the infinite-depth limit of finite-width neural networks. *arXiv preprint arXiv:2210.00688*, 2022.

He, K., Zhang, X., Ren, S., and Sun, J. Deep residual learning for image recognition. *arXiv preprint arXiv:1512.03385*, 14, 2015.

Holte, J. M. Discrete gronwall lemma and applications. In *MAA-NCS meeting at the University of North Dakota*, volume 24, pp. 1–7, 2009.

Johnson, A., Bulgarelli, L., Pollard, T., Horng, S., Celi, L. A., and Mark, R. Mimic-iv. *PhysioNet. Available online at: https://physionet.org/content/mimiciv/1.0/(accessed August 23, 2021)*, 2020.

Kidger, P. On neural differential equations. *arXiv preprint arXiv:2202.02435*, 2022.

Kidger, P., Bonnier, P., Perez Arribas, I., Salvi, C., and Lyons, T. Deep signature transforms. *Advances in Neural Information Processing Systems*, 32, 2019.

Kidger, P., Morrill, J., Foster, J., and Lyons, T. Neural controlled differential equations for irregular time series. *Advances in Neural Information Processing Systems*, 33: 6696–6707, 2020.
Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. *arXiv preprint arXiv:1412.6980*, 2014.

Koiran, P. and Sontag, E. D. Vapnik-chervonenkis dimension of recurrent neural networks. *Discrete Applied Mathematics*, 86(1):63–79, 1998.

Li, S. C.-X. and Marlin, B. M. A scalable end-to-end gaussian process adapter for irregularly sampled time series classification. *Advances in neural information processing systems*, 29, 2016.

Marion, P. Generalization bounds for neural ordinary differential equations and deep residual networks, 2023.

Marion, P., Fermanian, A., Biau, G., and Vert, J.-P. Scaling resnets in the large-depth regime. *arXiv preprint arXiv:2206.06929*, 2022.

Massaroli, S., Poli, M., Park, J., Yamashita, A., and Asama, H. Dissecting neural odes. *Advances in Neural Information Processing Systems*, 33:3952–3963, 2020.

Mohri, M., Rostamizadeh, A., and Talwalkar, A. *Foundations of machine learning*. MIT press, 2018.

Morrill, J., Salvi, C., Kidger, P., and Foster, J. Neural rough differential equations for long time series. In *International Conference on Machine Learning*, pp. 7829–7838. PMLR, 2021.

Rico-Martinez, R., Krischer, K., Kevrekidis, I., Kube, M., and Hudson, J. Discrete-vs. continuous-time nonlinear signal processing of cu electrodissolution data. *Chemical Engineering Communications*, 118(1):25–48, 1992.

Rico-Martinez, R., Anderson, J., and Kevrekidis, I. Continuous-time nonlinear signal processing: a neural network based approach for gray box identification. In *Proceedings of IEEE Workshop on Neural Networks for Signal Processing*, pp. 596–605. IEEE, 1994.

Rubanova, Y., Chen, R. T., and Duvenaud, D. K. Latent ordinary differential equations for irregularly-sampled time series. *Advances in neural information processing systems*, 32, 2019.

Virmaux, A. and Scaman, K. Lipschitz regularity of deep neural networks: analysis and efficient estimation. *Advances in Neural Information Processing Systems*, 31, 2018.
The appendix is structured as follows. Appendix A gives preliminary results on CDEs and covering numbers. Appendix B details the proof of the bounds of the approximation and discretization biases. Appendix C collects the proofs of the generalization bound, and the bound on the generalization gap. Proofs of Appendix A are given in Appendix D for the sake of clarity.

# A. Mathematical background

## A.1. Riemann-Stieltjes integral

**Lemma A.1.** Let $f : [0, 1] \to \mathbb{R}^{p \times d}$ be a continuous function, and $g : [0, 1] \to \mathbb{R}^d$ be of finite total variation. Then Riemann-Stieltjes integral

$$ \int_{0}^{t} f_{s} d g_{s} \in \mathbb{R}^{p}, $$

where the product is understood in the sense of matrix-vector multiplication, is well-defined and finite for every $t \in [0, 1]$.

We refer to Friz & Victoir (2010) for a thorough introduction to the Riemann-Stieltjes integral.

## A.2. Picard-Lindelöf Theorem

**Theorem A.2.** Let $x$ be a continuous path of bounded variation, and assume that $G : \mathbb{R}^p \to \mathbb{R}^{p \times d}$ is Lipschitz continuous. Then the CDE

$$ dy_{t} = G(y_{t}) dx_{t} $$

with initial condition $y_0 = \varphi(x_0)$ has a unique solution.

A full proof can be found in Fermanian et al. (2021). Remark that since in our setting, NCDE are Lipschitz since the neural vector fields are Lipschitz (Virmaux & Scaman, 2018). This ensures that the solutions to NCDEs are well defined.

We also need the following variation, which ensures that the CDE driven by $\tilde{x}$ is well-defined.

**Lemma A.3.** Let $x$ be a piecewise constant and right-continuous path taking finite values in $\mathbb{R}^d$, with a finite number of discontinuities at $0 < t_1, \dots, t_K = 1$, i.e.

$$ \lim_{t \to t_i^+} x_t = x_{t_i} $$

and

$$ x_t = x_{t_i} $$

for all $t \in [t_i, t_{i+1}]$, for all $i = 1, \dots, K-1$. Assume that $G : \mathbb{R}^p \to \mathbb{R}^{p \times d}$ is Lipschitz continuous. Then the CDE

$$ dy_{t} = G(y_{t}) dx_{t} $$

with initial condition $y_0 = \varphi(x_0)$ has a unique solution.

This result can be obtained by first remarking that since $x$ is piecewise constant, the solution to this CDE will also be piecewise constant, with discontinuities at $0 < t_1 < \dots < t_K$: indeed, the variations of $x$ between two points of discontinuity being null, the variations of the solution will also be null between these two points. The solution can then be recursively obtained by seeing that for all $t \in [t_i, t_{i+1}]$

$$ y_{t} = y_{t_{i}} = y_{t_{i-1}} + G(y_{t_{i-1}})(x_{t_{i}} - x_{t_{i-1}}) $$

with $y_t = \varphi(x_0)$ for $t \in [t_0, t_1]$, where $t_0 = 0$, and $y_{t_K} = y_1 = y_{t_{K-1}} + G(y_{t_{K-1}})(x_{t_K} - x_{t_{K-1}})$.

## A.3. Gronwall's Lemmas

**Lemma A.4 (Gronwall's Lemma for CDEs).** Let $x : [0, 1] \to \mathbb{R}^d$ be a continuous path of bounded variations, and $\phi : [0, 1] \to \mathbb{R}^d$ be a bounded measurable function. If

$$ \phi(t) \leq K_t + L_x \int_0^t \phi(s) \|dx_s\| $$
for all $t \in [0, 1]$, where $K_t$ is a time dependent constant, then

$$ \phi(t) \leq K_t \exp (L_x \|x\|_{1-\text{var},[0,t]}) $$

for all $t \geq 0$.

See Friz & Victoir (2010) for a proof. Remark that this Lemma does not require $\phi(\cdot)$ to be continuous. We also state a variant of the discrete Gronwall Lemma, which will allow us to obtain the bound for discrete inputs $x^D$.

**Lemma A.5 (Gronwall's Lemma for sequences).** Let $(y_k)_{k \geq 0}$, $(b_k)_{k \geq 0}$ and $(f_k)_{k \geq 0}$ be positive sequences of real numbers such that

$$ y_n \leq f_n + \sum_{l=0}^{n-1} b_l y_l $$

for all $n \geq 0$. Then

$$ y_n \leq f_n + \sum_{l=0}^{n-1} f_l b_l \prod_{j=l+1}^{n-1} (1 + b_j) $$

for all $n \geq 0$.

A proof can be found in Holte (2009) and Clark (1987).

### A.4. An upper bound on the total variation of the solution of a CDE

We have the following bound on the total variation of the solution to a CDE.

**Proposition A.6.** Let $\mathbf{F}: \mathbb{R}^{p \times d} \to \mathbb{R}^d$ be a Lipschitz vector field with Lipschitz constant $L_{\mathbf{F}}$ and $(x_t)_{t \in [0,1]}$ be a continuous path of total variation bounded by $L_x$. Let $(z_t)_{t \in [0,1]}$ be the solution of the CDE

$$ dz_t = \mathbf{F}(z_t)dx_t $$

with initial condition $z_0 \in \mathbb{R}^p$. Then for all $t \in [0, 1]$, one has

$$ \|z\|_{1-\text{var},[0,t]} \leq C_1(L_{\mathbf{F}}, \mathbf{F}, L_x) \|x\|_{1-\text{var},[0,t]} $$

with

$$ C_1(L_{\mathbf{F}}, \mathbf{F}, L_x) := \left[ L_{\mathbf{F}}(\|z_0\| + \|\mathbf{F}(0)\|_{\text{op}} L_x) \exp(L_{\mathbf{F}}L_x) + \|\mathbf{F}(0)\|_{\text{op}} \right] \exp(L_{\mathbf{F}}L_x). $$

For $f \in \mathcal{F}_\Theta$, let $z$ be the trajectory of the latent state associated to an input $x$, and $z^D$ be the trajectory of the latent state associated to an input $x^D$. Both $\|z_1\|_{1-\text{var}}$ and $\|z_1^D\|_{1-\text{var}}$ are upper bounded by

$$ C_1(L_{G_\psi}, G_\psi, L_x)L_x $$

where $C_1(L_{G_\psi}, G_\psi, L_x)$ is defined as

$$ \left[ L_{\sigma} B_{\mathbf{A}} L_{\sigma} (B_{\mathbf{U}} B_{x} + B_{v} + B_{\mathbf{b}} L_{x}) \exp(L_{\sigma} B_{\mathbf{A}} L_{x}) + L_{\sigma} B_{\mathbf{b}} \right] \exp(L_{\sigma} B_{\mathbf{A}} L_{x}). $$

The proof is given in Appendix D.2.

### A.5. Continuity of the flow of a CDE

**Theorem A.7.** Let $\mathbf{F}, \mathbf{G}: \mathbb{R}^p \to \mathbb{R}^{p \times d}$ be two Lipschitz vector fields with Lipschitz constants $L_{\mathbf{F}}, L_{\mathbf{G}}$. Let $x, r$ be either continuous or piecewise constant paths of total variation bounded by $L_x$ and $L_r$. Consider the controlled differential equations

$$ dw_t = \mathbf{F}(w_t)dx_t \quad \text{and} \quad dv_t = \mathbf{G}(v_t)dr_t $$
with initial conditions $w_0 \in \mathbb{R}^p$ and $v_0 \in \mathbb{R}^p$ respectively.

One has that

$$ \|w_t - v_t\| \le \left( \begin{aligned} & \|w_0 - v_0\| + \|x_0 - r_0\| \\ & + \|x - r\|_{\infty, [0, t]} (1 + L_F L_r C_1(L_F, F, L_x)) + \max_{v \in \Omega(G)} \|F(v) - G(v)\|_{\text{op}} L_r \end{aligned} \right) \times \exp(L_F L_x). $$

and symmetrically that

$$ \|w_t - v_t\| \le \left( \begin{aligned} & \|w_0 - v_0\| + \|x_0 - r_0\| \\ & + \|x - r\|_{\infty, [0, t]} (1 + L_G L_x C_1(L_G, G, L_r)) + \max_{v \in \Omega(F)} \|F(v) - G(v)\|_{\text{op}} L_x \end{aligned} \right) \times \exp(L_G L_r), $$

where the constant $C_1(\cdot, \cdot, \cdot)$ is given in Proposition A.6 and the sets $\Omega(G)$ and $\Omega(F)$ are defined as

$$ \Omega(G) = \{u \in \mathbb{R}^p \mid \|u\| \le (\|v_0\| + \|G(0)\| L_r) \exp(L_G L_r)\}. $$

and

$$ \Omega(F) = \{u \in \mathbb{R}^p \mid \|u\| \le (\|w_0\| + \|F(0)\| L_x) \exp(L_F L_x)\}. $$

The proof is given in Appendix D.3. We stress that any combination of continuous and piecewise constant paths can be used.

Using this result, we obtain the following theorem.

**Theorem A.8.** Let $f_{\theta_1}, f_{\theta_2} \in \mathcal{F} \oplus$ two predictors with respective parameters $\theta_1 = (\mathbf{A}_1, \mathbf{b}_1, \mathbf{U}_1, v_1, \Phi_1)$ and $\theta_2 = (\mathbf{A}_2, \mathbf{b}_2, \mathbf{U}_2, v_2, \Phi_2)$.

Both

$$ |f_{\theta_1}(x) - f_{\theta_2}(x)| $$

and

$$ |f_{\theta_1}(\mathbf{x}^D) - f_{\theta_2}(\mathbf{x}^D)| $$

are upper bounded by

$$ L_1 \|\Phi_1 - \Phi_2\| + L_2 \|\mathbf{A}_1 - \mathbf{A}_2\| + L_3 \|\mathbf{b}_1 - \mathbf{b}_2\| + L_4 \|\mathbf{U}_1 - \mathbf{U}_2\| + L_5 \|v_1 - v_2\| $$

where $L_1, L_2, L_3, L_4, L_5$ are explicit Lipschitz constants.

The proof is in Appendix D.4.

## A.6. Covering numbers

First, we recall the definition of the covering number of a class of functions.

**Definition A.9.** Let $\mathcal{H} = \{h : \mathcal{U} \subset \mathbb{R}^e \to \mathbb{R}^f\}$ be a class of functions. The covering number $\mathcal{N}(\mathcal{H}, \beta)$ of $\mathcal{H}$ is the minimal cardinality of a subset $C \subset \mathcal{H}$ such that for all $h \in \mathcal{H}$, there exists $\hat{h} \in C$ such that

$$ \sup_{x \in \mathcal{U}} \|h(x) - \hat{h}(x)\| \le \beta. $$
We state a Lemma from Chen et al. (2019a).

**Lemma A.10.** Let $\mathcal{G} = \{A \in \mathbb{R}^{d_1 \times d_2} \text{ s.t. } \|A\| \le \lambda\}$ for a given $\lambda > 0$. The covering number $\mathcal{N}(\mathcal{G}, \beta)$ is upper bounded by

$$ \left( 1 + \frac{2 \min(d_1^{\frac{1}{2}}, d_2^{\frac{1}{2}}) \lambda}{\varepsilon} \right)^{d_1 d_2} . $$
We now proceed to prove the bounds on the approximation and discretization biases, and the generalization inequalities. Let us precisely detail the different steps of our proofs.

1. The bound of the discretization bias, stated in Proposition 5.1), follows from the continuity of the flow given in Theorem A.7. Indeed, the discretization error depends on the difference between $z_1$ and $z_1^D$, which are solutions to identical CDE, but with different driving paths $x$ and embedded path $\tilde{x}$.

2. The bound of the approximation bias is obtained in a similar fashion. Indeed, it directly depends on the difference between $z_1^*$ and $z_1$. However, this time, only the initial condition and the vector fields are different, while the driving paths are identical.

3. Turning to the generalization inequalities, the first step consists in showing that NCDEs are Lipschitz with respect to their parameters, which we do in Theorem A.8. This leverages once again the continuity of the flow, since we aim at bounding the difference between two NCDEs with different parameters but identical driving paths.

4. We then use the central idea of Chen et al. (2019a) which consists in obtaining an upper bound on the covering number of a parameterized class by covering each of its parameters in Proposition 6.1. While the authors of this article use such a bound for RNN, the same technique applies for NCDE.

5. We then proceed to connect the covering number of $\mathcal{F}_\Theta$ to its Rademacher complexity building on arguments made by Bartlett et al. (2017) in Proposition 6.2.

6. Once this last result is obtained, we can obtain the generalization bound by resorting to classical techniques (Mohri et al., 2018).

7. The bound on the generalization gap is obtained by combining all the precedent elements and using a symmetrization argument to upper bound the two worst case bounds.

# B. Proof of bias bounds

## B.1. Proof of Proposition 5.1

Take $f_\theta \in \mathcal{F}_\Theta$. One has

$$ \mathcal{R}_n^D(f_\theta) - \mathcal{R}_n(f_\theta) = \frac{1}{n} \sum_{i=1}^n \left[ \ell(y^i, f_\theta(\mathbf{x}^{D,i})) - \ell(y^i, f_\theta(x^i)) \right]. \quad (7) $$

Considering a single individual $i \in \{1, \dots, n\}$, one has

$$ \ell(y^i, f_\theta(\mathbf{x}^{D,i})) - \ell(y^i, f_\theta(x^i)) \le L_\ell |f_\theta(\mathbf{x}^{D,i}) - f_\theta(x^i)| \quad (8) $$

using the Lipschitz continuity of the loss function with respect to its second argument. From the Cauchy-Schwarz inequality, it follows that

$$ |f_\theta(\mathbf{x}^{D,i}) - f_\theta(x^i)| = |\Phi^\top(z_1^D - z_1)| \le B_\alpha \|z_1^{D,i} - z_1^i\| \quad (9) $$

where $z_1^{D,i}$ and $z_1^i$ refer to the endpoint of the latent space trajectory of the NCDE, resp. with discrete and continuous input, associated to the predictor $f_\theta$.

$z_1^D$ and $z_1$ correspond to the endpoint of two CDEs with identical vector field and identical initial condition, but whose driving path differ. Using the continuity of the flow stated in Theorem A.7, and using the fact that the total variation of the piecewise constant path corresponding to the fill-forward embedding of the time series $\mathbf{x}^{D,i}$ is bounded by $L_x$, one has the inequality

$$ \|z_1^D - z_1\| \le \|x - \tilde{x}^i\|_{\infty, [0,t]} (1 + L_\sigma L_A L_x C_1(L_{G_\psi}, G_\psi, L_x)) \exp(L_\sigma B_A L_x) \quad (10) $$
where $\tilde{x}^i$ is the piecewise constant path corresponding to the fill-forward embedding of the time series $x^{D,i}$ and we recall that $C_1(L_{G_\psi}, G_\psi, L_x)$ is equal to

$$ \left[ L_{\sigma} B_{\mathbf{A}} L_{\sigma} (B_{\mathbf{U}} B_{x} + B_{v} + B_{\mathbf{b}} L_{x}) \exp(L_{\sigma} B_{\mathbf{A}} L_{x}) + L_{\sigma} B_{\mathbf{b}} \right] \exp(L_{\sigma} B_{\mathbf{A}} L_{x}). $$

One also has that

$$ \|\mathbf{x}^{D,i} - \mathbf{x}^i\|_{\infty, [0,1]} = \max_{j=1,...,K} \max_{s \in [t_j, t_{j+1}]} \|x_{t_j}^i - x_s^i\| \le L_x |D| \quad (11) $$

since the discretization is identical between individuals. Putting everything together, this yields for all $i = 1, ..., n$ that

$$ |f_{\theta}(\mathbf{x}^{D,i}) - f_{\theta}(\mathbf{x}^i)| \le \underbrace{B_{\alpha}(1 + L_{\sigma} L_{\mathbf{A}} L_x C_1(L_{G_{\psi}}, G_{\psi}, L_x)) \exp(L_{\sigma} B_{\mathbf{A}} L_x) L_x}_{C_2(\Theta, L_x, L_{\sigma})} |D| \quad (12) $$

and finally

$$ \mathcal{R}_n^D(f_\theta) - \mathcal{R}_n(f_\theta) \le L_\ell C_2(\Theta, L_x, L_\sigma) |D|. \quad (13) $$

This concludes the proof.

## B.2. Proof of Proposition 5.2

For any $f_\theta \in \mathcal{F}_\Theta$ with continuous input, the approximation bias writes

$$ \mathcal{R}(f_\theta) - \mathcal{R}(\mathbf{f}^*) = \mathbb{E}[\ell(y, f_\theta(x)) - \ell(y, \mathbf{f}^*(x))]. \quad (14) $$

Using the Lipschitz continuity of $\ell$ and the Cauchy-Schwarz inequality, one gets

$$ \ell(y, f_\theta(x)) - \ell(y, \mathbf{f}^*(x)) \le L_\ell \|\alpha_* - \Phi\| \|z_1\| + L_\ell B_\alpha \|z_1^* - z_1\|. \quad (15) $$

Since $z_1^*$ and $z_1$ are the solution to two CDEs with different vector fields and different initial conditions, the continuity of the flow stated in Theorem A.7 yields

$$ \|z_1^* - z_1\| \le \exp(L_{G^*} L_x) \left[ L \max_{u \in \Omega(G_\psi)} \|G_\psi(u) - G^*(u)\| + \max_{\|u\| \le B_x} \|\varphi(u) - \text{NN}_{U,v}(u)\| \right] \quad (16) $$

where

$$ \Omega(\mathbf{G}_\psi) = \{u \in \mathbb{R}^p \mid \|u\| \le (\|z_0\| + L_\sigma B_\mathbf{b} L_x) \exp(L_\sigma L_\mathbf{A} L_x)\} \\ \subset \{u \in \mathbb{R}^p \mid \|u\| \le (L_\sigma (B_\mathbf{U} B_x + B_v) + L_\sigma B_\mathbf{b} L_x) \exp(L_\sigma L_\mathbf{A} L_x)\}. $$

This means that $\text{diam}(\Omega(\mathbf{G}_\psi))$ is a function of $\Theta, L_x$ and $B_x$. To clarify this, we now write $\Omega(\Theta, L_x, B_x)$.

Putting everything together and taking expectations yields

$$ \mathcal{R}(f_\theta) - \mathcal{R}(\mathbf{f}^*) \le L_\ell B_\alpha \exp(L_{G^*} L_x) \left[ L_x \max_{u \in \Omega(\Theta, L_x, B_x)} \|\mathbf{G}_\psi(u) - \mathbf{G}^*(u)\| \right] \quad (17) $$

$$ + \max_{\|u\| \le B_x} \|\varphi(u) - \text{NN}_{U,v}(u)\| \quad (18) $$

$$ + L_\ell \|\alpha_* - \Phi\| \|z_1\|. \quad (19) $$

We now turn to a predictor for which one has $\Phi = \alpha_*$. This yields

$$ \mathcal{R}(f_\theta) - \mathcal{R}(\mathbf{f}^*) \quad (21) $$

$$ \le L_\ell B_\alpha \exp(L_{G^*} L_x) \left[ L_x \|\mathbf{G}_\psi - \mathbf{G}^*\|_{\infty, \Omega(\Theta, L_x, B_x)} + \|\varphi - \text{NN}_{U,v}\|_{\infty, B_x} \right]. \quad (22) $$

This concludes the proof.
## C. Proof of generalization bound

### C.1. Proof of Lemma 4.1

For all $f \in \mathcal{F}_{\Theta}$, one has the decomposition

$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(\mathbf{f}^*) = \mathcal{R}^D(\hat{f}^D) - \mathcal{R}_n^D(\hat{f}^D) \quad (23) $$

$$ + \mathcal{R}_n^D(\hat{f}^D) - \mathcal{R}_n^D(f) \quad (24) $$

$$ + \mathcal{R}_n^D(f) - \mathcal{R}_n(f) \quad (25) $$

$$ + \mathcal{R}_n(f) - \mathcal{R}(f) \quad (26) $$

$$ + \mathcal{R}(f) - \mathcal{R}(\mathbf{f}^*). \quad (27) $$

By optimality of $\hat{f}^D$, one has almost surely

$$ \mathcal{R}_n^D(\hat{f}^D) - \mathcal{R}_n^D(f) \le 0. \quad (28) $$

One is then left with the inequality

$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(\mathbf{f}^*) \le \mathcal{R}^D(\hat{f}^D) - \mathcal{R}_n^D(\hat{f}^D) \quad (29) $$

$$ + \mathcal{R}_n^D(f) - \mathcal{R}_n(f) \quad (30) $$

$$ + \mathcal{R}_n(f) - \mathcal{R}(f) \quad (31) $$

$$ + \mathcal{R}(f) - \mathcal{R}(\mathbf{f}^*). \quad (32) $$

Taking the supremum on the two differences between the empirical and expected risk yields that

$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(\mathbf{f}^*) \le \sup_{g \in \mathcal{F}_{\Theta}} [\mathcal{R}^D(g) - \mathcal{R}_n^D(g)] + \sup_{g \in \mathcal{F}_{\Theta}} [\mathcal{R}(g) - \mathcal{R}_n(g)] \quad (33) $$

$$ + \mathcal{R}_n^D(f) - \mathcal{R}_n(f) + \mathcal{R}(f) - \mathcal{R}(\mathbf{f}^*). \quad (34) $$

This concludes the proof.

### C.2. Proof of Proposition 6.1

We follow the proof strategy of Chen et al. (2019a). Starting from Theorem A.8, one can see that since

$$ \sup_x \|f_{\theta_1}(x) - f_{\theta_2}(x)\| \le L_1 \|\Phi_1 - \Phi_2\| + L_2 \|\mathbf{A}_1 - \mathbf{A}_2\| + L_3 \|\mathbf{b}_1 - \mathbf{b}_2\| \quad (35) $$

$$ + L_4 \|\mathbf{U}_1 - \mathbf{U}_2\| + L_5 \|v_1 - v_2\|, \quad (36) $$

where the supremum is taken on all $x$ such that $\|x\|_{1-\text{var},[0,1]} \le L_x$ and $\|x_0\| \le B_x$, it is sufficient to have $\hat{\theta} = (\hat{\Phi}, \hat{\mathbf{A}}, \hat{\mathbf{b}}, \hat{\mathbf{U}}, \hat{v})$ such that

$$ \|\Phi - \hat{\Phi}\| \le \frac{\beta}{5L_1}, \quad \|\mathbf{A} - \hat{\mathbf{A}}\| \le \frac{\beta}{5L_2}, \quad \|\mathbf{b} - \hat{\mathbf{b}}\| \le \frac{\beta}{5L_3}, \quad \|\mathbf{U} - \hat{\mathbf{U}}\| \le \frac{\beta}{5L_4} \quad (37) $$

and

$$ \|v - \hat{v}\| \le \frac{\beta}{5L_5} \quad (38) $$

to obtain an $\beta$ covering of $\mathcal{F}_{\Theta}$, since in this case we get that

$$ \sup_x \|f_{\theta_1}(x) - f_{\hat{\theta}}(x)\| \le \beta. \quad (39) $$
Using Lemma A.10 and denoting by $\mathcal{G}_{\Phi} = \{h: z \mapsto \Phi^{\top}z, \|\Phi\| \le B_{\alpha}\}$, and using the corresponding definitions for $\mathcal{G}_{\mathbf{A}}$, $\mathcal{G}_{\mathbf{b}}$, $\mathcal{G}_{\mathbf{U}}$ and $\mathcal{G}_{v}$, we get that

$$ \mathcal{N}(\mathcal{G}_{\Phi}, \frac{\beta}{5L_1}) \le \left(1 + \frac{10\sqrt{p}B_{\alpha}L_1}{\beta}\right)^p, \mathcal{N}(\mathcal{G}_{\mathbf{A}}, \frac{\beta}{5L_2}) \le \left(1 + \frac{10\sqrt{dp}B_{\mathbf{A}}L_2}{\beta}\right)^{dp} \quad (40) $$

$$ \mathcal{N}(\mathcal{G}_{\mathbf{b}}, \frac{\beta}{5L_3}) \le \left(1 + \frac{10 \min\{\sqrt{p}, \sqrt{d}\} B_{\mathbf{b}} L_3}{\beta}\right)^{dp}, \quad (41) $$

$$ \mathcal{N}(\mathcal{G}_{\mathbf{U}}, \frac{\beta}{5L_4}) \le \left(1 + \frac{10 \min\{\sqrt{p}, \sqrt{d}\} B_{\mathbf{U}} L_4}{\beta}\right)^{dp} \quad (42) $$

and

$$ \mathcal{N}(\mathcal{G}_{v}, \frac{\beta}{5L_5}) \le \left(1 + \frac{10\sqrt{p}B_vL_5}{\beta}\right)^p. \quad (43) $$

The covering number of $\mathcal{F}_{\Theta}$ is obtained by multiplying the covering number of each functional class (Chen et al., 2019a). Defining

$$ K_1 := \max\{B_{\alpha}L_1, B_vL_5\}, \quad (44) $$

$$ K_2 := \max\{B_{\mathbf{b}}L_3, B_{\mathbf{U}}L_4\}, \quad (45) $$

$$ K_3 := \max\{B_{\mathbf{b}}L_3, B_{\mathbf{U}}L_4, B_{\mathbf{A}}L_2\} \quad (46) $$

and using the fact that $\min\{\sqrt{d}, \sqrt{p}\} \le \sqrt{dp}$, we finally get that

$$ \mathcal{N}(\mathcal{F}_{\Theta}, \beta) \le \left(1 + \frac{10\sqrt{p}K_1}{\beta}\right)^{2p} \left(1 + \frac{10 \min\{\sqrt{p}, \sqrt{d}\} K_2}{\beta}\right)^{2dp} \left(1 + \frac{10\sqrt{dp}B_{\mathbf{A}}L_2}{\beta}\right)^{dp^2} \quad (47) $$

$$ \le \left(1 + \frac{10\sqrt{p}K_1}{\beta}\right)^{2p} \left(1 + \frac{10\sqrt{dp}K_3}{\beta}\right)^{dp(2+p)}. \quad (48) $$

The proof is identical for inputs $\mathbf{x}^D$. This concludes the proof.

### C.3. Proof of Proposition 6.2

We use the following Lemma from Bartlett et al. (2017).

**Lemma C.1.** Let $\mathcal{H}$ be a class of real-valued functions taking values in $[-M, M]$, for $M > 0$, and assume that $0 \in \mathcal{H}$. Then the empirical Rademacher complexity associated to a sample of $n$ datapoints verifies

$$ \mathrm{Rad}(\mathcal{H}) \le \inf_{\beta > 0} \left[ \frac{4\beta}{\sqrt{n}} + \frac{12}{n} \int_{\beta}^{2M\sqrt{n}} \sqrt{\log \mathcal{N}(\mathcal{H}, \beta)} d\beta \right] \quad (49) $$

We apply this Lemma to the class $\mathcal{F}_{\Theta}$. We trivially have that $0 \in \mathcal{F}$, since this function is recovered by taking $\theta = 0$. By Lemma 3.7, the value of $f_{\theta}$ is bounded by

$$ M_{\Theta} = B_{\alpha}L_{\sigma} \exp(L_{\mathbf{A}}L_{\sigma}L_{x}) (B_{\mathbf{U}}B_{x} + B_{v} + B_{\mathbf{b}}L_{x}). \quad (50) $$
In our setup, we get, from Proposition 6.1, that

$$ \int_{\beta}^{2M_{\Theta}\sqrt{n}} \sqrt{\log \mathcal{N}(\mathcal{F}_{\Theta}, \beta)} d\beta \qquad (51) $$

$$ \leq \int_{\beta}^{2M_{\Theta}\sqrt{n}} \sqrt{2p \log \left(1 + \frac{10\sqrt{p}K_1}{\beta}\right) + dp(2+p) \log \left(1 + \frac{10\sqrt{dp}K_2}{\beta}\right)} d\beta \qquad (52) $$

$$ \leq 2M_{\Theta}\sqrt{n}\sqrt{2p \log \left(1 + \frac{10\sqrt{p}K_1}{\beta}\right) + dp(2+p) \log \left(1 + \frac{10\sqrt{dp}K_2}{\beta}\right)}. \qquad (53) $$

Since for $x > 1$, $\log(1 + x) \leq \log(2x)$, we have, for $\beta$ small enough to ensure that both

$$ \frac{10\sqrt{dp}K_2}{\beta} \geq \frac{10\sqrt{p}K_1}{\beta} > 1 $$

so that

$$ \int_{\beta}^{2M_{\Theta}\sqrt{n}} \sqrt{\log \mathcal{N}(\mathcal{F}_{\Theta}, \beta)} d\beta \qquad (54) $$

$$ \leq 2M_{\Theta}\sqrt{n}\sqrt{2p \log \left(\frac{20\sqrt{p}K_1}{\beta}\right) + dp(p+2) \log \left(\frac{20\sqrt{dp}K_2}{\beta}\right)} \qquad (55) $$

Taking $\beta = \frac{1}{\sqrt{n}}$, one gets

$$ \mathrm{Rad}(\mathcal{F}_{\Theta}) \leq \frac{4}{n} + \frac{24M_{\Theta}}{\sqrt{n}} \sqrt{2p \log 20\sqrt{np}K_1 + dp(p+2) \log 20\sqrt{ndp}K_2}. \qquad (56) $$

yielding the same upper bound on the complexity.

### C.4. Proof of Theorems 4.2 and 4.3

The first part of the theorem is a straightforward application of Mohri et al. (2018), Theorem 11.3. Since our loss is $L_\ell$-Lipschitz and bounded by $M_\ell$, this gets us immediately that with probability at least $1 - \delta$, one has

$$ \mathcal{R}^D(\hat{f}^D) \leq \mathcal{R}_n(\hat{f}^D) + \frac{24M_{\Theta}L_{\ell}}{\sqrt{n}} \sqrt{2pU_1 + dp(p+2)U_2} + M_{\ell} \sqrt{\frac{\log 1/\delta}{2n}} \qquad (57) $$

with $U_1 := \log 20\sqrt{np}K_1$ and $U_2 := \log 20\sqrt{ndp}K_2$.

Note that this also gets us

$$ \mathcal{R}(\hat{f}) \leq \mathcal{R}_n(\hat{f}) + \frac{24M_{\Theta}L_{\ell}}{\sqrt{n}} \sqrt{2pU_1 + dp(p+2)U_2} + M_{\ell} \sqrt{\frac{\log 1/\delta}{2n}}, \qquad (58) $$

with

$$ \hat{f} \in \arg \min_{\theta \in \Theta} \mathcal{R}_n(f_{\theta}) $$

We now turn to the generalization gap. Using the generalization gap decomposition given in Lemma 4.1, it is clear that for any predictor $f_\theta$ parametrized by $\theta = (\Phi, \psi, \mathbf{U}, v)$ such that $\Phi = \alpha$,
$$ \mathcal{R}^D(\hat{f}^D) - \mathcal{R}(f^*) \qquad (59) $$

$$ \leq \sup_{g \in \mathcal{F}_{\Theta}} [\mathcal{R}^D(g) - \mathcal{R}_n^D(g)] + \sup_{g \in \mathcal{F}_{\Theta}} [\mathcal{R}(g) - \mathcal{R}_n(g)] \qquad (60) $$

$$ + L_{\ell} C_2(\Theta, L_x, L_{\sigma}) |D| \qquad (61) $$

$$ + L_{\ell} B_{\alpha} \exp(L_{\mathbf{G}} \cdot L_x) [L_x \|\mathbf{G}_{\psi} - \mathbf{G}^*\|_{\infty, \Omega(\Theta, L_x, B_x)} + \|\varphi - \mathbf{NN}_{\mathbf{U},v}\|_{\infty, B_x}] \qquad (62) $$

Now, by a classical symmetrization argument - see for instance Bach (2021), Proposition 4.2 - the obtained bounds on the Rademacher complexity imply that

$$ \sup_{g \in \mathcal{F}_{\Theta}} [\mathcal{R}(g) - \mathcal{R}_n(g)] \leq 2\mathrm{Rad}(\mathcal{F}_{\Theta}) \qquad (63) $$

and

$$ \sup_{g \in \mathcal{F}_{\Theta}} [\mathcal{R}^D(g) - \mathcal{R}_n^D(g)] \leq 2\mathrm{Rad}(\mathcal{F}_{\Theta}). \qquad (64) $$

Combining all these elements gets us the final result.
# D. Supplementary proofs

## D.1. Proof of Lemma 3.7

We recall that $||W||_{op} := \max_{||x||=1} ||Wx||$, and that $||W||_{op} \le ||W||$, which means that $||Wx|| \le ||W|| ||x||$ for all $x$. We first prove this result for a general CDE

$$dz_t = \mathbf{F}(z_t)dx_t$$

with initial condition $z_0 \in \mathbb{R}^p$, where the driving path $x$ is supposed to be continuous of bounded variation (or piecewise constant with a finite number of discontinuities.)

By definition,

$$z_t = z_0 + \int_0^t \mathbf{F}(z_t)dx_t. \qquad (65)$$

Taking norms, this yields

$$\|z_t\| \le \|z_0\| + \int_0^t \|\mathbf{F}(z_t)\|_{op} \|dx_t\|. \qquad (66)$$

Notice that since we have assumed $\mathbf{F}$ to be Lipschitz, one has for all $z \in \mathbb{R}^p$

$$\|\mathbf{F}(z)\|_{op} \le \|\mathbf{F}(z) - \mathbf{F}(0)\|_{op} + \|\mathbf{F}(0)\|_{op} \qquad (67)$$

$$\le \|\mathbf{F}(z) - \mathbf{F}(0)\| + \|\mathbf{F}(0)\|_{op} \qquad (68)$$

$$\le L_F \|z\| + \|\mathbf{F}(0)\|_{op}, \qquad (69)$$

where the last inequality follows from the fact that $\mathbf{F}$ is Lipschitz. It follows that

$$\|z_t\| \le \|z_0\| + \int_0^t (L_F \|z_t\| + \|\mathbf{F}(0)\|_{op}) \|dx_t\|. \qquad (70)$$

Using the fact that $\int_0^t \|dx_s\| = \|x\|_{1-\text{var},[0,t]}$, one gets

$$\|z_t\| \le \|z_0\| + \|\mathbf{F}(0)\|_{op} \|x\|_{1-\text{var},[0,t]} + L_F \int_0^t \|z_s\| \|dx_s\|. \qquad (71)$$

Applying Gronwall's Lemma for CDEs yields

$$\|z_t\| \le \left( \|z_0\| + \|\mathbf{F}(0)\|_{op} \|x\|_{1-\text{var},[0,t]} \right) \exp \left( L_F \|x\|_{1-\text{var},[0,t]} \right). \qquad (72)$$

Now, turning to the generative CDE (1), we obtain as a consequence that

$$|y| \le B_{\alpha} \left( B_{\varphi} + \|\mathbf{G}^{*}(0)\|_{op} L_x \right) \exp \left( L_{\mathbf{G}^{*}} L_x \right) + M_{\varepsilon}, \qquad (73)$$

since $y$ is the sum of a linear transformation of the endpoint of a CDE and a noise term $\varepsilon$ bounded by $M_{\varepsilon}$.

Turning now to $f_\theta \in \mathcal{F}_\Theta$, one has that

$$|f_\theta(x)| = |\Phi^\top z_1| \le \|\Phi\| \|z_1\| \le B_\alpha \|z_1\| \qquad (74)$$

and since $z_1$ is the solution of a NCDE, we can directly leverage the previous result to bound $\|z_1\|$. From the definition of $\mathcal{F}_\Theta$, it follows that

$$\|z_0\| \le L_{\sigma} B_{\mathbf{U}} B_x + L_{\sigma} B_v, \quad \|\mathbf{F}(0)\|_{op} \le L_{\sigma} B_{\mathbf{b}} \text{ and } L_{\mathbf{F}} \le L_{\sigma} B_{\mathbf{A}}. \qquad (75)$$
As a direct consequence, one has

$$|f_{\theta}(x)| \leq B_{\alpha} L_{\sigma} \exp(L_{\mathbf{A}} L_{\sigma} L_{x}) (B_{\mathbf{U}} B_{x} + B_{v} + B_{\mathbf{b}} L_{x}). \quad (76)$$

We now turn to $f_{\theta}(\mathbf{x}^D)$. Bounding the value of $z_1^D$ can be done by applying this lemma to the values of $z^D$ at points in $D$, since it is constant between those points. The proof leverages the discrete version of Gronwall's Lemma, stated in Lemma A.5. One has

$$\|z_{t_k}\| \le \|z_0\| + \sum_{l=0}^{k-1} \|\sigma(\mathbf{A}z_{t_l} + \mathbf{b})\| \|\Delta x_{t_{l+1}}\| \le \sum_{l=0}^{k-1} \|\sigma(\mathbf{A}z_{t_l} + \mathbf{b})\| \|\Delta x_{t_{l+1}}\|, \quad (77)$$

for all $k = 1, \dots, K$. Now, since

$$\|\sigma(\mathbf{A}z_{t_l} + \mathbf{b})\| \le L_{\sigma} B_{\mathbf{A}} \|z_l\| + L_{\sigma} B_{\mathbf{b}} \quad (78)$$

one has

$$\|z_{t_k}\| \le \|z_0\| + L_{\sigma} B_{\mathbf{b}} \sum_{l=0}^{k-1} \|\Delta x_{t_{l+1}}\| + L_{\sigma} B_{\mathbf{A}} \sum_{l=0}^{k-1} \|z_l\| \|\Delta x_{l+1}\| \quad (79)$$

$$\le \|z_0\| + L_{\sigma} B_{\mathbf{b}} L_x + L_{\sigma} B_{\mathbf{A}} \sum_{l=0}^{k-1} \|z_l\| \|\Delta x_{l+1}\| \quad (80)$$

and one gets by the discrete Gronwall Lemma that

$$\|z_{t_k}\| \le (\|z_0\| + L_{\sigma} B_{\mathbf{b}} L_x) \prod_{l=0}^{k-1} (1 + L_{\sigma} B_{\mathbf{A}} \|\Delta x_{t_{l+1}}\|) \quad (81)$$

$$\le (\|z_0\| + L_{\sigma} B_{\mathbf{b}} L_x) \exp(L_{\sigma} B_{\mathbf{A}} \|x\|_{1-\text{var},[0,1]}) \quad (82)$$

$$\le L_{\sigma} \exp(L_{\mathbf{A}} L_{\sigma} L_x) (B_{\mathbf{U}} B_x + B_v + B_{\mathbf{b}} L_x). \quad (83)$$

From this we get that

$$|f_{\theta}(\mathbf{x}^D)| \le B_{\alpha} L_{\sigma} \exp(L_{\mathbf{A}} L_{\sigma} L_x) (B_{\mathbf{U}} B_x + B_v + B_{\mathbf{b}} L_x) \quad (84)$$

This concludes the proof.

## D.2. Proof of Proposition A.6

We first consider a general CDE

$$dz_t = \mathbf{F}(z_t)dx_t \quad (85)$$

with initial condition $z_0 \in \mathbb{R}^p$. Take $s, r \in [0, t]$. By definition,

$$\|z_r - z_s\| = \left\| \int_s^r \mathbf{F}(z_u) dx_u \right\| \quad (86)$$

$$= \left\| \int_s^r \mathbf{F}(z_u - z_s + z_s) dx_u \right\| \quad (87)$$

$$\le \int_s^r (L_{\mathbf{F}} \|z_u - z_s + z_s\| + \|\mathbf{F}(0)\|_{\text{op}}) \|dx_u\| \quad (88)$$

$$\le \int_s^r (L_{\mathbf{F}} \|z_u - z_s\| + L_{\mathbf{F}} \|z_s\| + \|\mathbf{F}(0)\|_{\text{op}}) \|dx_u\| \quad (89)$$

$$= (L_{\mathbf{F}} \|z_s\| + \|\mathbf{F}(0)\|_{\text{op}}) \|x\|_{1-\text{var},[s,r]} + L_{\mathbf{F}} \int_s^r \|z_u - z_s\| \|dx_u\|. \quad (90)$$
Now since $z_s$ is the solution of a CDE evaluated at $s$, it can be bounded by Lemma 3.7 by

$$ \|z_s\| \le ( \|z_0\| + \|\mathbf{F}(0)\|_{\text{op}} \|x\|_{\text{1-var},[0,s]} ) \exp(L_{\mathbf{F}} \|x\|_{\text{1-var},[0,s]}). \quad (91) $$

This means that

$$ \|z_r - z_s\| \le \left[ L_{\mathbf{F}}(\|z_0\| + \|\mathbf{F}(0)\|_{\text{op}} \|x\|_{\text{1-var},[0,s]}) \exp(L_{\mathbf{F}} \|x\|_{\text{1-var},[0,s]}) + \|\mathbf{F}(0)\|_{\text{op}} \|x\|_{\text{1-var},[s,r]} \right] \|x\|_{\text{1-var},[s,r]} \quad (92) $$

$$ + L_{\mathbf{F}} \int_{s}^{r} \|z_u - z_s\| \|dx_u\|. \quad (94) $$

We may now use Gronwall's Lemma and the fact that $\|x\|_{\text{1-var},[0,s]} \le L_x$ for all $s \in [0, 1]$ to obtain that

$$ \|z_r - z_s\| \le \left[ L_{\mathbf{F}}(\|z_0\| + \|\mathbf{F}(0)\|_{\text{op}} L_x) \exp(L_{\mathbf{F}} L_x) + \|\mathbf{F}(0)\|_{\text{op}} \right] \exp(L_{\mathbf{F}} L_x) \|x\|_{\text{1-var},[s,r]}. \quad (95) $$

and thus

$$ \|z_r - z_s\| \le C_1(L_{\mathbf{F}}, \mathbf{F}, L_x) \|x\|_{\text{1-var},[s,r]}. \quad (96) $$

This means that the variations of $z$ on arbitrary intervals $[r, s]$ are bounded by the variations of $x$ on the same interval, times an interval independent constant.

From this, we may immediately conclude that

$$ \|z\|_{\text{1-var},[0,t]} \le C_1(L_{\mathbf{F}}, \mathbf{F}, L_x) \|x\|_{\text{1-var},[0,t]}. \quad (97) $$

which concludes the proof for a general CDE.

We now turn to $f_\theta \in \mathcal{F}_\Theta$. The previous result allows to bound the total variation of a CDE and thus of the trajectory of the latent state. Using this proposition with

$$ C_1(L_{\mathbf{G}_\psi}, \mathbf{G}_\psi, L_x) := \left[ L_\sigma B_{\mathbf{A}} L_\sigma (B_{\mathbf{U}} B_x + B_v + B_{\mathbf{b}} L_x) \exp(L_\sigma B_{\mathbf{A}} L_x) + L_\sigma B_{\mathbf{b}} \right] \exp(L_\sigma B_{\mathbf{A}} L_x) $$

yields

$$ \|z\|_{\text{1-var}} \le C_1(L_{\mathbf{G}_\psi}, \mathbf{G}_\psi, L_x)L_x. \quad (98) $$

This concludes the proof for $f_\theta$ with continuous input. For a discrete input, the proof directly transfers by remarking that

$$ \|z^D\|_{\text{1-var},[0,1]} = \sum_{k=1}^{K-1} \|z_{t_k}^D - z_{t_{k-1}}^D\| \quad (99) $$

and by using the discrete Gronwall Lemma.
### D.3. Proof of Theorem A.7

For all $t \in [0, 1]$, one has the decomposition

$$w_t - v_t = w_0 - v_0 + \int_0^t \mathbf{F}(w_s)dx_s - \int_0^t \mathbf{G}(v_s)dr_s \quad (100)$$

$$= w_0 - v_0 + \int_0^t (\mathbf{F}(w_s) - \mathbf{F}(v_s))dx_s + \int_0^t \mathbf{F}(v_s)dx_s - \int_0^t \mathbf{G}(v_s)dr_s \quad (101)$$

$$= w_0 - v_0 \quad (102)$$

$$+ \int_0^t (\mathbf{F}(w_s) - \mathbf{F}(v_s))dx_s \quad (103)$$

$$+ \int_0^t \mathbf{F}(v_s)d(x_s - r_s) \quad (104)$$

$$+ \int_0^t (\mathbf{F}(v_s) - \mathbf{G}(v_s))dr_s. \quad (105)$$

We control every one of these terms separately and conclude by applying Gronwall's Lemma. Writing

$$A := \left\| \int_0^t (\mathbf{F}(w_s) - \mathbf{F}(v_s))dx_s \right\|,$$

$$B := \left\| \int_0^t \mathbf{F}(v_s)d(x_s - r_s) \right\|,$$

$$C := \left\| \int_0^t (\mathbf{F}(v_s) - \mathbf{G}(v_s))dr_s \right\|$$

it is clear that

$$\|w_t - z_t\| \le \|w_0 - z_0\| + A + B + C. \quad (106)$$

**Control of the term A.** One has

$$A \le \int_0^t \|\mathbf{F}(w_s) - \mathbf{F}(v_s)\|_{\text{op}} \|dx_s\|. \quad (107)$$

Since $F$ is $L_F$-Lipschitz, this gives

$$A \le L_F \int_0^t \|w_s - v_s\| \|dx_s\|. \quad (108)$$

**Control of the term B.** One gets using integration by parts

$$\int_0^t \mathbf{F}(v_s)d(x_s - r_s) = (x_t - r_t) - (x_0 - r_0) - \int_0^t (x_s - r_s)^\top d\mathbf{F}(v_s). \quad (109)$$

This gives

$$B \le \|x_0 - r_0\| + \|x - r\|_{\infty, [0, t]} (1 + \|\mathbf{F}(v)\|_{1-\text{var}, [0, t]}). \quad (110)$$

Since $F$ is $L_F$-Lipschitz, one gets

$$\|\mathbf{F}(v)\|_{1-\text{var}, [0, t]} \le L_F \|v\|_{1-\text{var}, [0, t]}. \quad (111)$$

Since $v$ is the solution to a CDE, we can resort to Proposition A.6 to bound its total variation. This yields

$$\|v\|_{1-\text{var}, [0, t]} \le C_1(L_F, \mathbf{F}, L_r) \|r\|_{1-\text{var}, [0, t]} \quad (112)$$

and finally

$$B \le \|x_0 - r_0\| + \|x - r\|_{\infty, [0, t]} (1 + L_F C_1 (L_F, \mathbf{F}, L_x) \|r\|_{1-\text{var}, [0, t]}). \quad (113)$$
**Control of the term C.** Finally, we get that

$$C \le \int_0^t \|\mathbf{F}(v_s) - \mathbf{G}(v_s)\|_{\text{op}} \|dr_s\| \quad (114)$$

$$\le \max_{v \in \Omega(\mathbf{G})} \|\mathbf{F}(v) - \mathbf{G}(v)\|_{\text{op}} L_r \quad (115)$$

where we recall that $\Omega(\mathbf{G})$ is the closed ball

$$\Omega(\mathbf{G}) = \{u \in \mathbb{R}^p \mid \|u\| \le (\|v_0\| + \|\mathbf{G}(0)\| L_r) \exp(L_{\mathbf{G}} L_r)\}.$$

Indeed, since $v$ is the solution of a CDE, its norm at any time $t \in [0, 1]$ is bounded as stated in Lemma 3.7. One can therefore bound the difference $\|\mathbf{F}(v_s) - \mathbf{G}(v_s)\|_{\text{op}}$ by considering all possible values of $v$.

**Putting everything together.** Combining the obtained bounds on all terms, one is left with

$$\|w_t - v_t\| \le \|w_0 - v_0\| \quad (116)$$

$$+ \|x_0 - r_0\| + \|x - r\|_{\infty, [0, t]} (1 + L_{\mathbf{F}} L_r C_1(L_{\mathbf{F}}, \mathbf{F}, L_x)) \quad (117)$$

$$+ L_r \max_{v \in \Omega(\mathbf{G})} \|\mathbf{F}(v) - \mathbf{G}(v)\|_{\text{op}} \quad (118)$$

$$+ L_{\mathbf{F}} \int_0^t \|w_s - v_s\| \|dx_s\|. \quad (119)$$

The proof is concluded by using Gronwall's Lemma A.4. If the path $(x_t)$ is continuous, we resort to its continuous version. If it is piecewise constant, we use its discrete version. The path $(\|w_t - v_t\|)_t$ only needs to be measurable and bounded. There are no assumptions on its continuity. This finally yields

$$\|w_t - v_t\| \le \left( \|w_0 - v_0\| + \|x_0 - r_0\| \right) \quad (120)$$

$$+ \|x - r\|_{\infty, [0, t]} (1 + L_{\mathbf{F}} L_r C_1(L_{\mathbf{F}}, \mathbf{F}, L_x)) + \max_{v \in \Omega(\mathbf{G})} \|\mathbf{F}(v) - \mathbf{G}(v)\|_{\text{op}} L_r \quad (121)$$

$$\times \exp(L_{\mathbf{F}} L_x). \quad (122)$$

Now, notice that in our proof, the two CDEs are exchangeable. This means that we immediately get the alternative bound

$$\|w_t - v_t\| \le \left( \|w_0 - v_0\| + \|x_0 - r_0\| \right) \quad (123)$$

$$+ \|x - r\|_{\infty, [0, t]} (1 + L_{\mathbf{G}} L_x C_1(L_{\mathbf{G}}, \mathbf{G}, L_r)) + \max_{v \in \Omega(\mathbf{F})} \|\mathbf{F}(v) - \mathbf{G}(v)\|_{\text{op}} L_x \quad (124)$$

$$\times \exp(L_{\mathbf{G}} L_r) \quad (125)$$

where we recall that

$$\Omega(\mathbf{F}) = \{u \in \mathbb{R}^p \mid \|u\| \le (\|w_0\| + \|\mathbf{F}(0)\| L_x) \exp(L_{\mathbf{F}} L_x)\}.$$

### D.4. Proof of Theorem A.8

One has

$$\|f_{\theta_1}(x) - f_{\theta_2}(x)\| = \|\Phi_1^\top z_1^1 - \Phi_2^\top z_1^2\| \le \|\Phi_1 - \Phi_2\| \|z_1^1\| + B_\alpha \|z_1^1 - z_1^2\|. \quad (126)$$

Using Theorem A.7, one gets that

$$\|z_1^1 - z_1^2\| \le \exp(L_\sigma B_{\mathbf{A}} L_x) (\|z_0^1 - z_0^2\| + L_x \|\mathbf{G}_1 - \mathbf{G}_2\|_{\infty, \Omega}), \quad (127)$$
where

$$ \Omega \subset \{u \in \mathbb{R}^p \mid \|u\| \le (L_{\sigma}(B_{\mathbf{U}}B_x + B_v) + L_{\sigma}B_{\mathbf{b}}L_x) \exp(L_{\sigma}L_{\mathbf{A}}L_x)\}. $$

We first bound

$$ \|z_0^1 - z_0^2\| \le L_{\sigma} \|\mathbf{U}_1 - \mathbf{U}_2\| B_x + L_{\sigma} \|v_1 - v_2\|. \quad (128) $$

One also has that

$$ \|\mathbf{G}_1 - \mathbf{G}_2\|_{\infty, \Omega} \le L_{\sigma} \|\mathbf{A}_1 - \mathbf{A}_2\| \times \max_{u \in \Omega} \|u\| + L_{\sigma} \|\mathbf{b}_1 - \mathbf{b}_2\| \quad (129) $$

It follows from the inclusion of $\Omega$ that

$$ \max_{u \in \Omega} \|u\| \le (L_{\sigma}(B_{\mathbf{U}}B_x + B_v) + L_{\sigma}B_{\mathbf{b}}L_x) \exp(L_{\sigma}L_{\mathbf{A}}L_x). \quad (130) $$

Since $z_1^1$ is bounded as the endpoint of a NCDE, one gets using Theorem A.7 that

$$ \|\Phi_1 - \Phi_2\| \|z_1^1\| \le L_{\sigma} \left[ B_{\mathbf{U}}B_x + B_v + (B_{\mathbf{A}} + B_{\mathbf{b}})L \right] \exp(L_{\sigma}L_x B_{\mathbf{A}}) \|\Phi_1 - \Phi_2\| \quad (131) $$

Putting everything together, one gets that

$$ \|f_1(x) - f_2(x)\| \le L_1 \|\Phi_1 - \Phi_2\| \quad (132) $$

$$ + L_2 \|\mathbf{A}_1 - \mathbf{A}_2\| + L_3 \|\mathbf{b}_1 - \mathbf{b}_2\| \quad (133) $$

$$ + L_4 \|\mathbf{U}_1 - \mathbf{U}_2\| + L_5 \|v_1 - v_2\| \quad (134) $$

with

$$ L_1 := L_{\sigma} \left[ B_{\mathbf{U}} B_x + B_v + (B_{\mathbf{A}} + B_{\mathbf{b}}) L \right] \exp(L_{\sigma} L_x B_{\mathbf{A}}), \quad (135) $$

$$ L_2 := B_{\alpha} \exp(2L_{\sigma}B_{\mathbf{A}}L_x)L_xL_{\sigma}(L_{\sigma}(B_{\mathbf{U}}B_x + B_v) + L_{\sigma}B_{\mathbf{b}}L_x), \quad (136) $$

$$ L_3 := B_{\alpha} \exp(L_{\sigma}B_{\mathbf{A}}L_x)L_xL_{\sigma} \quad (137) $$

$$ L_4 := B_{\alpha}B_x \exp(L_{\sigma}B_{\mathbf{A}}L_x)L_{\sigma} \quad (138) $$

$$ L_5 := B_{\alpha} \exp(L_{\sigma}B_{\mathbf{A}}L_x)L_{\sigma}. \quad (139) $$

which concludes our proof. The proof for $|f_{\theta_1}(\mathbf{x}^D) - f_{\theta_2}(\mathbf{x}^D)|$ is identical.