From NeurODEs to AutoencODEs: a mean-field control framework for
  width-varying Neural Networks
The connection between Residual Neural Networks (ResNets) and continuous-time
control systems (known as NeurODEs) has led to a mathematical analysis of
neural networks which has provided interesting results of both theoretical and
practical significance. However, by construction, NeurODEs have been limited to
describing constant-width layers, making them unsuitable for modeling deep
learning architectures with layers of variable width. In this paper, we propose
a continuous-time Autoencoder, which we call AutoencODE, based on a
modification of the controlled field that drives the dynamics. This adaptation
enables the extension of the mean-field control framework originally devised
for conventional NeurODEs. In this setting, we tackle the case of low Tikhonov
regularization, resulting in potentially non-convex cost landscapes. While the
global results obtained for high Tikhonov regularization may not hold globally,
we show that many of them can be recovered in regions where the loss function
is locally convex. Inspired by our theoretical findings, we develop a training
method tailored to this specific type of Autoencoders with residual
connections, and we validate our approach through numerical experiments
conducted on various examples.
Cristina Cipriani, Massimo Fornasier, Alessandro Scagliotti
# From NeurODEs to AutoencODEs:

###### Abstract

The connection between Residual Neural Networks (ResNets) and continuous-time control systems (known as NeurODEs) has led to a mathematical analysis of neural networks which has provided interesting results of both theoretical and practical significance. However, by construction, NeurODEs have been limited to describing constant-width layers, making them unsuitable for modeling deep learning architectures with layers of variable width. In this paper, we propose a continuous-time Autoencoder, which we call AutoencODE, based on a modification of the controlled field that drives the dynamics. This adaptation enables the extension of the mean-field control framework originally devised for conventional NeurODEs. In this setting, we tackle the case of low Tikhonov regularization, resulting in potentially non-convex cost landscapes. While the global results obtained for high Tikhonov regularization may not hold globally, we show that many of them can be recovered in regions where the loss function is locally convex. Inspired by our theoretical findings, we develop a training method tailored to this specific type of Autoencoders with residual connections, and we validate our approach through numerical experiments conducted on various examples.

_Keywords--_ Machine Learning, Optimal Control, Gradient Flow, Minimising Movement Scheme, Autoencoders

## 1 Introduction

In recent years, the field of artificial intelligence has witnessed remarkable progress across diverse domains, including computer vision and natural language processing. In particular, neural networks have emerged as a prominent tool, revolutionizing numerous machine learning tasks. Consequently, there is an urgent demand for a robust mathematical framework to analyze their intricate characteristics. A deep neural network can be seen as map \(\phi:\mathbb{R}^{d_{\mathrm{in}}}\to\mathbb{R}^{d_{\mathrm{out}}}\), obtained as the composition of \(L\gg 1\) applications \(\phi=\phi_{L}\circ\ldots\circ\phi_{1}\), where, for every \(n=1,\ldots,L\), the function \(\phi_{n}:\mathbb{R}^{d_{n}}\to\mathbb{R}^{d_{n+1}}\) (also referred as _the \(n\)-th layer_ of the network) depends on a _trainable_ parameter \(\theta_{a}\in\mathbb{R}^{m_{n}}\). The crucial process of choosing the values of the parameters \(\theta_{1},\ldots,\theta_{L}\) is known as the _training of the network_. For a complete survey on the topic, we recommend the textbook [24].

Recent advancements have explored the link between dynamical systems, optimal control, and deep learning, proposing a compelling perspective. In the groundbreaking work [29], it was highlighted how the problem of training very deep networks can be alleviated by the introduction of a new type layer called "Residual Block". This consists in using the identity map as skip connection, and after-addition activations. In other words, every layer has the following form:

\[X_{n+1}=\phi_{n}(X_{n})=X_{n}+\mathcal{F}(X_{n},\theta_{n}), \tag{1.1}\]

where \(X_{n+1}\) and \(X_{n}\) are, respectively, the output and the input of the \(n\)-th layer. This kind of architecture is called _Residual Neural Network_ (or ResNet). It is important to observe that, in order to give sense to the sum in (1.1), in each layer the dimension of the input should coincide with the dimension of the output. In the practice of Deep
Learning, this novel kind of layer has turned out to be highly beneficial, since it is effective in avoiding the "vanishing of the gradients" during the training [5], or the saturation of the network's accuracy [28]. Indeed, before [29], these two phenomena had limited for long time the large-scale application of deep architectures.

Despite the original arguments in support of residual blocks being based on empirical considerations, their introduction revealed nevertheless a more mathematical and rigorous bridge between residual deep networks and controlled dynamical systems. Indeed, what makes Residual Neural Networks particularly intriguing is that they can be viewed as discretized versions of continuous-time dynamical systems. This dynamical approach was proposed independently in [18] and [26], and it was greatly popularized in the machine learning community under the name of NeurODEs by [13]. This connection with dynamical systems relies on reinterpreting the iteration (1.1) as a step of the forward-Euler approximation of the following dynamical system:

\[\dot{X}(t)=\mathcal{F}(X(t),\theta(t)), \tag{1.2}\]

where \(t\mapsto\theta(t)\) is the map that, instant by instant, specifies the value of the parameter \(\theta\). Moreover, the training of these neural networks, typically formulated as empirical risk minimization, can be reinterpreted as an optimal control problem. Given a labelled dataset \(\{(X_{0}^{i},Y_{0}^{i})\}_{i=1}^{N}\) of size \(N\geq 1\), the depth of the time-continuous neural network (1.2) is denoted by \(T>0\). Then, training this network amounts to learning the control signals \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) in such a way that the terminal output \(X_{T}^{i}\) of (1.2) is close to it corresponding label \(Y_{0}^{i}\) for all \(i=1,\ldots,N\), with respect to some distortion measure \(\ell(\cdot,\cdot)\in C^{1}\). A typical choice is \(\ell(x,y):=|x-y|^{2}\), which is often referred as the _squared loss function_ in the machine learning literature. Therefore, it is possible to formulate the following optimal control problem

\[\inf_{\theta\in L^{2}([0,T];\mathbb{R}^{m})}J^{N}(\theta):=\begin{cases} \frac{1}{N}\sum_{i=1}^{N}\ell\big{(}X^{i}(T),Y^{i}(T)\big{)}+\lambda \int_{0}^{T}|\theta(t)|^{2}\,dt,\\ \text{s.t.}\ \begin{cases}\dot{X}^{i}(t)=\mathcal{F}(t,X^{i}(t),\theta(t)), \hskip 28.452756pt\dot{Y}^{i}(t)=0,\\ \big{(}X^{i}(t),Y^{i}(t)\big{)}\big{|}_{t=0}=(X_{0}^{i},Y_{0}^{i}),\end{cases}&i \in\{1,\ldots,N\},\end{cases}\]

where, differently from (1.2), we admit here the explicit dependence of the dynamics on the time variable. Notice that the objective function also comprises of Tikhonov regularization, tuned by the parameter \(\lambda\), which plays a crucial role in the analysis of this control problem. The benefit of interpreting the training process in this manner results in the possibility of exploiting established results from the branch of mathematical control theory, to better understand this process. A key component of optimal control theory is a set of necessary conditions, known as Pontryagin Maximum Principle (PMP), that must be satisfied by any (local) minimizer \(\theta\). These conditions were introduced in [37] and have served as inspiration for the development of innovative algorithms [33] and network structures [12] within the machine learning community.

This work specifically addresses a variant of the optimal control problem presented above, in which the focus is on the case of an infinitely large dataset. This formulation gives rise to what is commonly known as a _mean-field optimal control problem_, where the term "mean-field" emphasizes the description of a multiparticle system through its averaged effect. In this context, the focus is on capturing the collective behavior of the system rather than individual particle-level dynamics, by considering the population as a whole. As a consequence, the parameter \(\theta\) is shared by the entire population of input-target pairs, and the optimal control is required to depend on the initial distribution \(\mu_{0}(x,y)\in\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d})\) of the input-target pairs. Therefore, the optimal control problem needs to be defined over spaces of probability measures, and it is formulated as follows:

\[\inf_{\theta\in L^{2}([0,T];\mathbb{R}^{m})}J(\theta):=\begin{cases}\int_{ \mathbb{R}^{2d}}\ell(x,y)\,d\mu_{T}(x,y)+\lambda\int_{0}^{T}|\theta(t)|^{2} \,dt,\\ \text{s.t.}\ \begin{cases}\partial_{t}\mu_{t}(x,y)+\nabla_{x}\cdot( \mathcal{F}(t,x,\theta)\mu_{t}(x,y))=0&t\in[0,T],\\ \mu_{t}|_{t=0}(x,y)=\mu_{0}(x,y),\end{cases}\end{cases}\]

This area of study has gained attention in recent years, and researchers have derived the corresponding Pontryagin Maximum Principle in various works, such as [19] and [7]. It is worth mentioning that there are other types of mean-field analyses of neural networks, such as the well-known work [36], which focus on mean-field at the parameter level, where the number of parameters is assumed to be infinitely large. However, our approach in this work takes a different viewpoint, specifically focusing on the control perspective in the case of an infinitely large dataset.

One of the contributions of this paper is providing a more accessible derivation of the necessary conditions for optimality, such as the well-known Pontryagin Maximum Principle. Namely, we characterize the stationary points of the cost functional, and we are able to recover the PMP that was deduced in [7] under the assumption of large values of the regularization parameter \(\lambda\), and whose proof relied on an infinite-dimensional version of the Lagrange multiplier rule. This alternative perspective offers a clearer and more intuitive understanding of the PMP, making it easier to grasp and apply it in practical scenarios.
In addition, we aim at generalizing the applicability of the results presented in [7] by considering a possibly non-convex regime, corresponding to small values of the parameter \(\lambda>0\). As mentioned earlier, the regularization coefficient \(\lambda\) plays a crucial role in determining the nature of the cost function. Indeed, when \(\lambda\) is sufficiently large, the cost function is convex on the sub-level sets, and it is possible to prove the existence and uniqueness of the solution of the optimal control problem that arises from training NeurODEs. Additionally, in this highly-regularized scenario, desirable properties of the solution, such as its continuous dependence on the initial data and a bound on the generalization capabilities of the networks, have been derived in [7].

However, in practical applications, a large regularization parameter may cause a poor performance of the trained NeurODE on the task. In other words, in the highly-regularized case, the cost functional is unbalanced towards the \(L^{2}\)-penalization, at the expenses of the term that promotes that each datum \(X_{0}^{i}\) is driven as close as possible to the corresponding target \(Y_{0}^{i}\). This motivated us to investigate the case of low Tikhonov regularization. While we cannot globally recover the same results as in the highly-regularized regime, we find interesting results concerning local minimizers. Moreover, we also show that the (mean field) optimal control problem related to the training of the NeurODE induces a gradient flow in the space of admissible controls. The perspective of the gradient flow leads us to consider the well-known minimizing movement scheme, and to introduce a proximal stabilization term to the cost function in numerical experiments. This approach effectively addresses the well-known instability issues (see [14]) that arise when solving numerically optimal control problems (or when training NeurODEs) with iterative methods based on the PMP. It is important to note that our stabilization technique differs from previous methods, such as the one introduced in [33].

From NeurODEs to AutoencODEs.Despite their huge success, it should be noted that NeurODEs (as well as ResNets, their discrete-time counterparts) in their original form face a limitation in capturing one of the key aspects of modern machine learning architectures, namely the discrepancy in dimensionality between consecutive layers. As observed above, the use of skip connections with identity mappings requires a "rectangular" shape of the network, where the width of the layers are all identical and constant with respect to the input's dimension. This restriction poses a challenge when dealing with architectures that involve layers with varying dimensions, which are common in many state-of-the-art models. Indeed, the inclusion of layers with different widths can enhance the network's capacity to represent complex functions and to capture intricate patterns within the data. In this framework, Autoencoders have emerged as a fundamental class of models specifically designed to learn efficient representations of input data by capturing meaningful features through an encoder-decoder framework. More precisely, the encoder compresses the input data into a lower-dimensional latent space, while the decoder reconstructs the original input from the compressed representation. The concept of Autoencoders was first introduced in the 1980s in [39], and since then, it has been studied extensively in various works, such as [30], among many others. Nowadays, Autoencoders have found numerous applications, including data compression, dimensionality reduction, anomaly detection, and generative modeling. Their ability to extract salient features and capture underlying patterns in an unsupervised manner makes them valuable tools in scenarios where labeled training data is limited or unavailable. Despite their huge success in practice, there is currently a lack of established theory regarding the performance guarantees of these models.

Prior works, such as [20], have extended the control-theoretic analysis of NeurODEs to more general width-varying neural networks. Their model is based on an integro-differential equation that was first suggested in [34] in order to study the continuum limit of neural networks with respect to width and depth. In such an equation the state variable has a dependency on both time and space since the changing dimension over time is viewed as an additional spatial variable. In [20, Section 6] the continuous space-time analog of residual neural networks proposed in [34] has been considered and discretized in order to model variable width ResNets of various types, including convolutional neural networks. The authors assume a simple time-dependent grid, and use forward difference discretization for the time derivative and Newton-Cotes for discretizing the integral term, but refer to more sophisticated moving grids in order to possibly propose new types of architectures. In this setting, they are also able to derive some stability estimates and generalization properties in the overparametrized regime, making use of turnpike theory in optimal control [22]. In principle, there could be several different ways to model width-varying neural networks with dynamical systems, e.g., forcing some structure on the control variables, or formulating a viability problem. In this last case, a possibility could be to require admissible trajectories to visit some lower-dimensional subsets during the evolution. For an introduction to viability theory, we recommend the monograph [4], while we refer to [8, 9] for recent results on viability theory for differential inclusions in Wasserstein spaces.

In contrast, our work proposes a simpler extension of the control-theoretical analysis. It is based on a novel design of the vector field that drives the dynamics, allowing us to develop a continuous-time model capable of accommodating various types of width-varying neural networks. This approach has the advantage of leveraging insights and results obtained from our previous work [7]. Moreover, the simplicity of our model facilitates the implementation of residual networks with variable width and allows us to test their performance in machine learning tasks. In order to capture
width-varying neural networks, we need to extend the previous control-theoretical framework to a more general scenario, in particular we need to relax some of the assumptions of [7]. This is done in Subsection 2.2, where we introduce a discontinuous-in-time dynamics that can describe a wider range of neural network architectures. By doing so, we enable the study of Autoencoders (and, potentially, of other width-varying architectures) from a control-theoretic point of view, with the perspective of getting valuable insights into their behavior.

Furthermore, we also generalize the types of activation functions that can be employed in the network. The previous work [7] primarily focused on sigmoid functions, which do not cover the full range of activations commonly employed in practice. Our objective is to allow for unbounded activation functions, which are often necessary for effectively solving certain tasks. By considering a broader set of activation functions, we aim at enhancing the versatility and applicability of our model.

Furthermore, in contrast to [7], we introduce a stabilization method to allow the numerical resolution of the optimal control problem in the low-regularized regime, as previously discussed. This stabilization technique provides the means to test the architecture with our training approach on various tasks: from low-dimensional experiments, which serve to demonstrate the effectiveness of our method, to more sophisticated and high-dimensional tasks such as image reconstruction. In Section 5, we present all the experiments and highlight noteworthy behaviors that we observe. An in-depth exploration of the underlying reasons for these behaviors is postponed to future works.

The structure of the paper is the following: Section 2 discusses the dynamical model of NeurODEs and extends it to the case of width-varying neural networks, including Autoencoders, which we refer to as AutoencODEs. In Section 3, we present our mean-field analysis, focusing on the scenario of an infinitely large dataset. We formulate the mean-field optimal control problem, we derive a set of necessary optimality conditions, and we provide a convergence result for the finite-particles approximation. At the end of this section, we compare our findings with the ones previously obtained in [7]. Section 4 covers the implementation and the description of the training procedure, and we compare it with other methods for NeurODEs existing in the literature. Finally, in Section 5, we present the results of our numerical experiments, highlighting interesting properties of the AutoencODEs that we observe.

### Measure-theoretic preliminaries

Given a metric space \((X,d_{X})\), we denote by \(\mathcal{M}(X)\) the space of signed Borel measures in \(X\) with finite total variation, and by \(\mathcal{P}(X)\) the space of probability measures, while \(\mathcal{P}_{c}(X)\subset\mathcal{P}(X)\) represents the set of probability measures with compact support. Furthermore, \(\mathcal{P}_{c}^{N}(X)\subset\mathcal{P}_{c}(X)\) denotes the subset of empirical or atomic probability measures. Given \(\mu\in\mathcal{P}(X)\) and \(f:X\to Y\), with \(f\)\(\mu-\)measurable, we denote with \(f_{\#}\mu\in\mathcal{P}(Y)\) the push-forward measure defined by \(f_{\#}\mu(B)=\mu(f^{-1}(B))\) for any Borel set \(B\subset Y\). Moreover, we recall the change-of-variables formula

\[\int_{Y}g(y)\,d\big{(}f_{\#}\mu\big{)}(y)=\int_{X}g\circ f(x)\,d\mu(x) \tag{1.3}\]

whenever either one of the integrals makes sense.

We now focus on the case \(X=\mathbb{R}^{d}\) and briefly recall the definition of the Wasserstein metrics of optimal transport in the following definition, and refer to [2, Chapter 7] for more details.

**Definition 1.1**.: _Let \(1\leq p<\infty\) and \(\mathcal{P}_{p}(\mathbb{R}^{d})\) be the space of Borel probability measures on \(\mathbb{R}^{d}\) with finite \(p\)-moment. In the sequel, we endow the latter with the \(p\)-Wasserstein metric_

\[W_{p}^{p}(\mu,\nu):=\inf\left\{\int_{\mathbb{R}^{2d}}|z-\hat{z}|^{p}\ d\pi(z, \hat{z})\ \big{|}\ \pi\in\Pi(\mu,\nu)\right\},\]

_where \(\Pi(\mu,\nu)\) denotes the set of transport plan between \(\mu\) and \(\nu\), that is the collection of all Borel probability measures on \(\mathbb{R}^{d}\times\mathbb{R}^{d}\) with marginals \(\mu\) and \(\nu\) in the first and second component respectively._

It is a well-known result in optimal transport theory that when \(p=1\) and \(\mu,\nu\in\mathcal{P}_{c}(\mathbb{R}^{d})\), then the following alternative representation holds for the Wasserstein distance

\[W_{1}(\mu,\nu)=\sup\left\{\int_{\mathbb{R}^{d}}\varphi(x)\,d\big{(}\mu-\nu \big{)}(x)\,\big{|}\,\varphi\in\mathrm{Lip}(\mathbb{R}^{d}),\ \mathrm{Lip}(\varphi)\leq 1 \right\}\,, \tag{1.4}\]

by Kantorovich's duality [2, Chapter 6]. Here, \(\mathrm{Lip}(\mathbb{R}^{d})\) stands for the space of real-valued Lipschitz continuous functions on \(\mathbb{R}^{d}\), and \(\mathrm{Lip}(\varphi)\) is the Lipschitz constant of a mapping \(\varphi\) defined ad

\[Lip(\varphi):=\sup_{x,y\in\mathbb{R}^{d},x\neq y}\frac{\|\varphi(x)-\varphi(y )\|}{\|x-y\|}\]
## 2 Dynamical Model of NeurODEs

### Notation and basic facts

In this paper, we consider controlled dynamical systems in \(\mathbb{R}^{d}\), where the velocity field is prescribed by a function \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) that satisfies these basic assumptions.

**Assumption 1**.: _The vector field \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) satisfies the following:_

1. _For every_ \(x\in\mathbb{R}^{d}\) _and every_ \(\theta\in\mathbb{R}^{m}\)_, the map_ \(t\mapsto\mathcal{F}(t,x,\theta)\) _is measurable in_ \(t\)_._
2. _For every_ \(R>0\) _there exists a constant_ \(L_{R}>0\) _such that, for every_ \(\theta\in\mathbb{R}^{m}\)_, it holds_ \[|\mathcal{F}(t,x_{1},\theta)-\mathcal{F}(t,x_{2},\theta)|\leq L_{R}(1+|\theta |)|x_{1}-x_{2}|,\quad\text{ for a.e. }t\in[0,T]\text{ and every }x_{1},x_{2}\in B_{R}(0),\] _from which it follows that_ \(|\mathcal{F}(t,x,\theta)|\leq L_{R}(1+|x|)(1+|\theta|)\) _for a.e._ \(t\in[0,T]\)_._
3. _For every_ \(R>0\) _there exists a constant_ \(L_{R}>0\) _such that, for every_ \(\theta_{1},\theta_{2}\in\mathbb{R}^{m}\)_, it holds_ \[|\mathcal{F}(t,x,\theta_{1})-\mathcal{F}(t,x,\theta_{2})|\leq L_{R}(1+|\theta _{1}|+|\theta_{2}|)|\theta_{1}-\theta_{2}|,\quad\text{ for a.e. }t\in[0,T]\text{ and every }x\in B_{R}(0).\]

The control system that we are going to study is

\[\begin{cases}\dot{x}(t)=\mathcal{F}(t,x(t),\theta(t)),\ \ \text{ a.e. in }[0,T],\\ x(0)=x_{0},\end{cases} \tag{2.1}\]

where \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) is the control that drives the dynamics. Owing to Assumption 1, the classical Caratheodory Theorem (see [27, Theorem 5.3]) guarantees that, for every \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) and for every \(x_{0}\in\mathbb{R}^{d}\), the Cauchy problem (2.1) has a unique solution \(x:[0,T]\to\mathbb{R}^{d}\). Hence, for every \((t,\theta)\in[0,T]\times L^{2}([0,T],\mathbb{R}^{m})\), we introduce the flow map \(\Phi^{\theta}_{(0,t)}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) defined as

\[\Phi^{\theta}_{(0,t)}(x_{0}):=x(t), \tag{2.2}\]

where \(t\mapsto x(t)\) is the absolutely continuous curve that solves (2.1), with Cauchy datum \(x(0)=x_{0}\) and corresponding to the admissible control \(t\mapsto\theta(t)\). Similarly, given \(0\leq s<t\leq T\), we write \(\Phi^{\theta}_{(s,t)}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) to denote the flow map obtained by prescribing the Cauchy datum at the more general instant \(s\geq 0\). We now present the properties of the flow map defined in (2.2) that describes the evolution of the system: we show that is well-posed, and we report some classical properties.

**Proposition 2.1**.: _For every \(t\in[0,T]\) and for every \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\), let \(\mathcal{F}\) satisfy Assumption 1. Then, the flow \(\Phi^{\theta}_{(0,t)}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) is well-defined for any \(x_{0}\in\mathbb{R}^{d}\) and it satisfies the following properties._

* _For every_ \(R>0\) _and_ \(\rho>0\)_, there exists a constant_ \(\bar{R}>0\) _such that_ \[|\Phi^{\theta}_{(0,t)}(x)|\leq\bar{R}\] _for every_ \(x\in B_{R}(0)\) _and every_ \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) _such that_ \(||\theta||_{L^{2}}\leq\rho\)_._
* _For every_ \(R>0\) _and_ \(\rho>0\)_, there exists a constant_ \(\bar{L}>0\) _such that, for every_ \(t\in[0,T]\)_, it holds_ \[|\Phi^{\theta}_{(0,t)}(x_{1})-\Phi^{\theta}_{(0,t)}(x_{2})|\leq\bar{L}|x_{1}-x _{2}|\] _for every_ \(x_{1},x_{2}\in B_{R}(0)\) _and every_ \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) _such that_ \(||\theta||_{L^{2}}\leq\rho\)_._
Even though the framework introduced in Assumption 1 is rather general, in this paper we specifically have in mind the case where the mapping \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) represents the feed-forward dynamics associated to residual neural networks. In this scenario, the parameter \(\theta\in\mathbb{R}^{m}\) encodes the _weights_ and _shifts_ of the network, i.e., \(\theta=(W,b)\), where \(W\in\mathbb{R}^{d\times d}\) and \(b\in\mathbb{R}^{d}\). Moreover, the mapping \(\mathcal{F}\) has the form:

\[\mathcal{F}(t,x,\theta)=\sigma(Wx+b),\]

where \(\sigma:\mathbb{R}^{d}\to\mathbb{R}^{d}\) is a nonlinear function acting component-wise, often called in literature _activation function_. In this work, we consider sigmoidal-type activation functions, such as the hyperbolic tangent function:

\[\sigma(x)=\tanh(x),\]

as well as smooth approximations of the Rectified Linear Unit (ReLU) function, which is defined as:

\[\sigma(x)=\max\{0,x\}. \tag{2.3}\]

We emphasize the need to consider smoothed versions of the ReLU function due to additional differentiability requirements on \(\mathcal{F}\), which will be further clarified in Assumption 2. Another useful activation function covered by Assumption 2 is the Leaky Rectified Linear Unit (Leaky ReLU) function:

\[\sigma(x)=\max\{0,x\}-\max\{-\alpha x,0\} \tag{2.4}\]

where \(\alpha\in[0,1]\) is a predetermined parameter that allows the output of the function to have negative values. The smooth approximations of (2.3) and (2.4) that we consider will be presented in Section 4.

### From NeurODEs to AutoencODEs

As explained in the Introduction, NeurODEs and ResNets -their discrete-time counterparts- face the limitation of a "rectangular" shape of the network because of formulas (1.2) and (1.1), respectively. To overcome this fact, we aim at designing a continuous-time model capable of describing width-varying neural networks, with a particular focus on Autoencoders, as they represent the prototype of neural networks whose layers operate between spaces of different dimensions. Indeed, Autoencoders consist of an _encoding phase_, where the layers' dimensions progressively decrease until reaching the "latent dimension" of the network. Subsequently, in the _decoding phase_, the layers' widths are increased until the same dimensionality as the input data is restored. For this reason, Autoencoders are prominent examples of width-varying neural networks, since the changes in layers' dimensions lie at the core of their functioning. Sketches of encoders and Autoencoders are presented in Figure 1. Finally, we insist on the fact that our model can encompass as well other types of architectures. In this regard, in Remark 2.2 we discuss how our approach can be extended to U-nets.

EncoderOur goal is to first model the case of a network which sequentially reduces the dimensionality of the layers' outputs. For this purpose, we artificially force some of the components not to evolve anymore, while we let the others be active part of the dynamics. More precisely, given an input variable \(x_{0}\in\mathbb{R}^{d}\), we denote with \((\mathcal{I}_{j})_{j=0,\dots,r}\) an increasing filtration, where each element \(\mathcal{I}_{j}\) contains the sets of indices whose corresponding components are _inactive_, i.e., they are constant and do not contribute to the dynamics. Clearly, since the layers' width will decrease sequentially, the filtration of inactive components \(\mathcal{I}_{j}\) will increase, i.e.

\[\emptyset=:\mathcal{I}_{0}\subsetneq\mathcal{I}_{1}\subsetneq...\subsetneq \mathcal{I}_{r}\subsetneq\{1,\dots,d\},\quad r<d,\quad j=0,\dots,r.\]

Figure 1: Left: network with an encoder structure. Right: Autoencoder.
On the other hand, the sets of indices of _active_ components define an decreasing filtration \(\mathcal{A}_{j}:=\{1,\ldots,d\}\setminus\mathcal{I}_{j}\) for \(j=0,\ldots,r\). As opposed to before, the sets of active components \((\mathcal{A}_{j})_{j=0,\ldots,r}\) satisfy

\[\{1,\ldots,d\}=:\mathcal{A}_{0}\supseteq\mathcal{A}_{1}\supseteq... \supseteq\mathcal{A}_{r}\supseteq\emptyset,\quad r<d,\quad j=0,\ldots,r.\]

We observe that, for every \(j=0,\ldots,r\), the sets \(\mathcal{A}_{j}\) and \(\mathcal{I}_{j}\) provide a partition of \(\{1,\ldots,d\}\). A visual representation of this model for encoders is presented on the left side of Figure 2.

Now, in the time interval \([0,T]\), let us consider \(r+1\) nodes \(0=t_{0}<t_{1}<...<t_{r}<t_{r+1}=T\). For \(j=0,\ldots,r\), we denote with \([t_{j},t_{j+1}]\) the sub-interval and, for every \(x\in\mathbb{R}^{d}\), we use the notation \(x_{\mathcal{I}_{j}}:=(x_{i})_{i\in\mathcal{I}_{j}}\) and \(x_{\mathcal{A}_{j}}:=(x_{i})_{i\in\mathcal{A}_{j}}\) to access the components of \(x\) belonging to \(\mathcal{I}_{j}\) and \(\mathcal{A}_{j}\), respectively. Hence, the controlled dynamics for any \(t\in[t_{j},t_{j+1}]\) can be described by

\[\begin{cases}\dot{x}_{\mathcal{I}_{j}}(t)=0,\\ \dot{x}_{\mathcal{A}_{j}}(t)=\mathcal{G}_{j}(t,x_{\mathcal{A}_{j}}(t),\theta( t)),\end{cases} \tag{2.5}\]

where \(\mathcal{G}_{j}:[t_{j},t_{j+1}]\times\mathbb{R}^{|\mathcal{A}_{j}|}\times \mathbb{R}^{m}\rightarrow\mathbb{R}^{|\mathcal{A}_{j}|}\), for \(j=0,\ldots,r\), and \(x(0)=x_{\mathcal{A}_{0}}(0)=x_{0}\). Furthermore, the dynamical system describing the encoding part is

\[\begin{cases}\dot{x}(t)=\mathcal{F}(t,x(t),\theta(t)),\quad\text{a.e }t\in[0,T],\\ x(0)=x_{0}\end{cases}\]

where, for \(t\in[t_{j},t_{j+1}]\), we define the discontinuous vector field as follows

\[\left(\mathcal{F}(t,x,\theta)\right)_{k}=\begin{cases}\left(\mathcal{G}(t,x_{ \mathcal{A}_{j}},\theta)\right)_{k},&\text{if }k\in\mathcal{A}_{j},\\ 0,&\text{if }k\in\mathcal{I}_{j}.\end{cases}\]

**Remark 2.1**.: _Notice that \(\theta(t)\in\mathbb{R}^{m}\) for every \(t\in[0,T]\), according to the model that we have just described. However, it is natural to expect that, since \(x\) has varying active components, in a similar way the controlled dynamics \(\mathcal{F}(t,x,\theta)\) shall not explicitly depend at every \(t\in[0,T]\) on every component of \(\theta\)._

AutoencoderWe now extend the previous model to the case of networks which not only decrease the dimensionality of the layers, but they are also able to increase the layers' width in order to restore the original dimension of the input data. Here we denote by \(z_{0}\in\mathbb{R}^{\bar{d}}\) the input variable, and we fictitiously augment the input's dimension, so that we consider the initial datum \(x_{0}=(z_{0},\underline{0})\in\mathbb{R}^{d}=\ \mathbb{R}^{\bar{d}}\times \mathbb{R}^{\bar{d}}\), where \(\underline{0}\in\mathbb{R}^{\bar{d}}\). We make use of the following notation for every \(x\in\mathbb{R}^{d}\):

\[x=\left((z_{i})_{i=1,\ldots,\bar{d}},(z^{H}_{i})_{i=1,\ldots,\bar{d}})\right)\]

where \(z^{H}\) is the augmented (or _shadow_) part of the vector \(x\). In this model, the time horizon \([0,T]\) is splitted using the following time-nodes:

\[0=t_{0}\leq t_{1}\leq...\leq t_{r}\leq...\leq t_{2r}\leq t_{2r+1}:=T\]

where \(t_{r}\), which was the end of the encoder in the previous model, is now the instant corresponding to the bottleneck of the autoencoder. Similarly as before, we introduce two families of partitions of \(\{1,\ldots,\bar{d}\}\) modeling the active and

Figure 2: Left: Embedding of an encoder into a dynamical system. Right: model for an Autoencoder.
non-active components of, respectively, \(z\) and \(z^{H}\). The first filtrations are relative to the encoding phase and they involve the component of \(z\):

\[\begin{cases}\mathcal{I}_{j-1}\subsetneq I_{j}&\text{if }1\leq j\leq r,\\ \mathcal{I}_{j}=\mathcal{I}_{j-1}&\text{if }j>r,\end{cases}\qquad\qquad \qquad\qquad\begin{cases}\mathcal{A}_{j-1}\supsetneq\mathcal{A}_{j}&\text{if }1\leq j \leq r,\\ \mathcal{A}_{j}=\mathcal{A}_{j-1}&\text{if }j>r.\end{cases}\]

where \(\mathcal{I}_{0}:=\emptyset\), \(\mathcal{I}_{r}\subsetneq\{1,\ldots,\tilde{d}\}\) and \(\mathcal{A}_{0}=\{1,\ldots,\tilde{d}\}\), \(\mathcal{A}_{r}\supsetneq\emptyset\). The second filtrations, that aim at modeling the decoder, act on the shadow part of \(x\), i.e., they involve the components of \(z^{H}\):

\[\begin{cases}\mathcal{I}_{j-1}^{H}=\{1,\ldots,d\}&\text{if }1\leq j\leq r,\\ \mathcal{I}_{j}^{H}\subsetneq\mathcal{I}_{j-1}^{H}&\text{if }r<j\leq 2r,\end{cases} \qquad\qquad\qquad\begin{cases}\mathcal{A}_{j-1}^{H}=\emptyset& \text{if }1\leq j\leq r,\\ \mathcal{A}_{j}^{H}\supsetneq\mathcal{A}_{j-1}^{H}&\text{if }r<j\leq 2r.\end{cases}\]

While the encoder structure acting on the input data \(z_{0}\) is the same as before, in the decoding phase we aim at activating the components that have been previously turned off during the encoding. However, since the information contained in the original input \(z_{0}\) should be first compressed and then decompressed, we should not make use of the values of the components that we have turned off in the encoding and hence, we cannot re-activate them. Therefore, in our model the dimension is restored by activating components of \(z^{H}\), the shadow part of \(x\), which we recall was initialized equal to \(\underline{0}\in\mathbb{R}^{\tilde{d}}\). This is the reason why we introduce sets of active and inactive components also for the shadow part of the state variable. A sketch of this type of model is presented on the right of Figure 2. Moreover, in order to be consistent with the classical structure of an autoencoder, the following identities must be satisfied:

1. \(\mathcal{A}_{j}\cap\mathcal{A}_{j}^{H}=\emptyset\) for every \(j=1,\ldots,2r\),
2. \(\mathcal{A}_{2r}\cup\mathcal{A}_{2r}^{H}=\{1,\ldots,\tilde{d}\}\).

The first identity formalizes the constraint that the active component of \(z\) and those of \(z^{H}\) cannot overlap and must be distinct, while the second identity imposes that, at the end of the evolution, the active components in \(z\) and \(z^{H}\) should sum up exactly to \(1,\ldots,\tilde{d}\). Furthermore, from the first identity we derive that \(\mathcal{A}_{j}\subseteq(\mathcal{A}_{j}^{H})^{C}=\mathcal{I}_{j}^{H}\) and, similarly, \(\mathcal{A}_{j}^{H}\subseteq\mathcal{I}_{j}\) for every \(j=1,\ldots,2r\). Moreover, \(\mathcal{A}_{r}\) satisfies the inclusion \(\mathcal{A}_{r}\subseteq\mathcal{A}_{j}\) for every \(j=1,\ldots,2r\), which is consistent with the fact that layer with the smallest width is located in the bottleneck, i.e., in the interval \([t_{r},t_{r+1}]\). Finally, from the first and the second assumption, we obtain that \(\mathcal{A}_{2r}^{H}=\mathcal{I}_{2r}\), i.e., the final active components of \(z^{H}\) coincide with the inactive components of \(z\), and, similarly, \(\mathcal{I}_{2r}^{H}=\mathcal{A}_{2r}\). Finally, to access the active components of \(x=(z,z^{H})\), we make use of the following notation:

\[x_{\mathcal{A}_{j}}=(z_{k})_{k\in\mathcal{A}_{j}},\quad x_{\mathcal{A}_{j}^{H }}=(z_{k}^{H})_{k\in\mathcal{A}_{j}^{H}}\quad\text{and}\quad x_{\mathcal{A}_{j },\mathcal{A}_{j}^{H}}=(z_{\mathcal{A}_{j}},z_{\mathcal{A}_{j}^{H}}^{H}),\]

and we do the same for the inactive components:

\[x_{\mathcal{I}_{j}}=(z_{k})_{k\in\mathcal{I}_{j}},\quad x_{\mathcal{I}_{j}^{H }}=(z_{k}^{H})_{k\in\mathcal{I}_{j}^{H}}\quad\text{and}\quad x_{\mathcal{I}_{ j},\mathcal{I}_{j}^{H}}=(z_{\mathcal{I}_{j}},z_{\mathcal{I}_{j}^{H}}^{H}).\]

We are now in position to write the controlled dynamics in the interval \(t_{j}\leq t\leq t_{j+1}\):

\[\begin{cases}\dot{x}_{\mathcal{I}_{j},\mathcal{I}_{j}^{H}}(t)=0,\\ \dot{x}_{\mathcal{A}_{j},\mathcal{A}_{j}^{H}}(t)=\mathcal{G}_{j}(t,x_{\mathcal{ A}_{j},\mathcal{A}_{j}^{H}}(t),\theta(t)),\end{cases} \tag{2.6}\]

where \(\mathcal{G}_{j}:[t_{j},t_{j+1}]\times\mathbb{R}^{|\mathcal{A}_{j}|+|\mathcal{A }_{j}^{H}|}\times\mathbb{R}^{m}\to\mathbb{R}^{|\mathcal{A}_{j}|+|\mathcal{A} _{j}^{H}|}\,,\) for \(j=0,\ldots,2r\), and \(x_{\mathcal{I}_{0}}^{H}(0)=\underline{0}\), \(x_{\mathcal{A}_{0}}(0)=x_{0}\). As before, we define the discontinuous vector field \(\mathcal{F}\) for \(t\in[t_{j},t_{j+1}]\), as follows

\[\big{(}\mathcal{F}(t,x,\theta)\big{)}_{k}=\begin{cases}\big{(}\mathcal{G}(t,x _{\mathcal{A}_{j}},\theta)\big{)}_{k}\,,&\text{if }k\in\mathcal{A}_{j}\cup\mathcal{A}_{j}^{H}\\ 0,&\text{if }k\in\mathcal{I}_{j}\cup\mathcal{I}_{j}^{N}.\end{cases}\]

Hence, we are now able to describe any type of width-varying neural network through a continuous-time model depicted by the following dynamical system

\[\begin{cases}\dot{x}(t)=\mathcal{F}(t,x(t),\theta(t))\quad\text{a.e. in }[0,T],\\ x(0)=x_{0}.\end{cases}\]

It is essential to highlight the key difference between the previous NeurODE model in (2.6) and the current model: the vector field \(\mathcal{F}\) now explicitly depends on the time variable \(t\) to account for sudden dimensionality drops, where certain components are forced to remain constant. As a matter of fact, the resulting dynamics exhibit high discontinuity in
the variable \(t\). To the best of our knowledge, this is the first attempt to consider such discontinuous dynamics in NeurODEs. Previous works, such as [26, 18], typically do not include an explicit dependence on the time variable in the right-hand side of NeurODEs, or they assume a continuous dependency on time, as in [7]. Furthermore, it is worth noting that the vector field \(\mathcal{F}\) introduced to model autoencoders satisfies the general assumptions outlined in Assumption 1 at the beginning of this section.

**Remark 2.2**.: _The presented model, initially designed for Autoencoders, can be easily extended to accommodate various types of width-varying neural networks, including architectures with long skip-connections such as U-nets [38]. While the specific details of U-nets are not discussed in detail, their general structure is outlined in Figure 3. U-nets consist of two main components: the contracting path (encoder) and the expansive path (decoder). These paths are symmetric, with skip connections between corresponding layers in each part. Within each path, the input passes through a series of convolutional layers, followed by a non-linear activation function (often ReLU), and other operations (e.g., max pooling) which are not encompassed by our model. The long skip-connections that characterize U-nets require some modifications to the model of autoencoder described above. If we denote with \(\bar{d}_{i}\) for \(i=0,\ldots,r\) the dimensionality of each layer in the contracting path, we have that \(\bar{d}_{2r-i}=\bar{d}_{i}\) for every \(i=0,\ldots,r\). Then, given an initial condition \(z_{0}\in\mathbb{R}^{\bar{d}_{0}}\), we embed it into the augmented state variable_

\[x_{0}=(z_{0},\underline{0}),\text{ where }\underline{0}\in\mathbb{R}^{\bar{d}_{ 1}+\ldots+\bar{d}_{r}}.\]

_As done in the previous model for autoencoder, we consider time-nodes \(0=t_{0}<\ldots<t_{2r}=T\), and in each sub-interval we introduce a controlled dynamics with the scheme of active/inactive components depicted in Figure 3._

## 3 Mean-field Analysis

In this section, we extend the dynamical model introduced in Section 2 to its mean-field limit, which corresponds to the scenario of an infinitely large dataset. Within this framework, we formulate the training of NeurODEs and AutoencODEs as a mean-field optimal control problem and provide the associated necessary optimality conditions. It is worth noting that our analysis covers both the high-regularized regime, as studied in previous work [7], as well as the low-regularized regime, which has not been extensively addressed before. In this regard, we dedicate a subsection to a detailed comparison with the results obtained in [7]. Additionally, we investigate the case of finite-particles approximation and we establish a quantitative bound on the generalization capabilities of these networks.

### Mean-field dynamical model

In this section, we employ the same view-point as in [7], and we consider the case of a dataset with an infinite number of observations. In our framework, each datum is modeled as a point \(x_{0}\in\mathbb{R}^{d}\), and it comes associated to its corresponding label \(y_{0}\in\mathbb{R}^{d}\). Notice that, in principle, in Machine Learning applications the label (or _target_) datum

Figure 3: Embedding of the U-net into a higher-dimensional dynamical system.
\(y_{0}\) may have dimension different from \(d\). However, the labels' dimension is just a matter of notation and does not represent a limit of our model. Following [7], we consider the curve \(t\mapsto(x(t),y(t))\) which satisfies

\[\dot{x}(t)=\mathcal{F}(t,x(t),\theta(t))\qquad\text{and}\qquad\dot{y}(t)=0 \tag{3.1}\]

for a.e. \(t\in[0,T]\), and \((x(0),y(0))=(x_{0},y_{0})\). We observe that the variable \(y\) corresponding to the labels is not changing, nor it is affecting the evolution of the variable \(x\). We recall that the flow associated to the dynamics of the variable \(x\) is denoted by \(\Phi^{\theta}_{(0,t)}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) for every \(t\in[0,T]\), and it has been defined in (2.2). Moreover, in regards to the full dynamics prescribed by (3.1), for every admissible control \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) we introduce the extended flow \(\boldsymbol{\Phi}^{\theta}_{(0,t)}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to \mathbb{R}^{d}\times\mathbb{R}^{d}\), which reads

\[\boldsymbol{\Phi}^{\theta}_{(0,t)}(x_{0},y_{0})=(\Phi^{\theta}_{(0,t)}(x_{0} ),y_{0}) \tag{3.2}\]

for every \(t\in[0,T]\) and for every \((x_{0},y_{0})\in\mathbb{R}^{d}\times\mathbb{R}^{d}\). We now consider the case of an infinite number of labeled data \((X^{i}_{0},Y^{i}_{0})_{i\in I}\), where \(I\) is an infinite set of indexes. In our mathematical model, we understand this data distribution as a compactly-supported probability measure \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d})\). Moreover, for every \(t\in[0,T]\), we denote by \(t\mapsto\mu_{t}\) the curve of probability measures in \(\mathcal{P}_{c}(\mathbb{R}^{d}\times\mathbb{R}^{d})\) that models the evolution of the solutions of (3.1) corresponding to the Cauchy initial conditions \((X^{i}_{0},Y^{i}_{0})_{i\in I}\). In other words, the curve \(t\mapsto\mu_{t}\) satisfies the following continuity equation:

\[\partial_{t}\mu_{t}(x,y)+\nabla_{x}\cdot\big{(}\mathcal{F}(t,x,\theta_{t})\mu _{t}(x,y)\big{)}=0,\qquad\mu_{t}|_{t=0}(x,y)=\mu_{0}(x,y), \tag{3.3}\]

understood in the sense of distributions, i.e.

**Definition 3.1**.: _For any given \(T>0\) and \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\), we say that \(\mu\in\mathcal{C}([0,T],\mathcal{P}_{c}(\mathbb{R}^{2d}))\) is a weak solution of (3.3) on the time interval \([0,T]\) if_

\[\int_{0}^{T}\int_{\mathbb{R}^{2d}}\big{(}\partial_{t}\psi(t,x,y)+\nabla_{x} \psi(t,x,y)\cdot\mathcal{F}(t,x,\theta_{t})\big{)}\,d\mu_{t}(x,y)\,dt=0, \tag{3.4}\]

_for every test function \(\psi\in\mathcal{C}^{1}_{c}((0,T)\times\mathbb{R}^{2d})\)._

Let us now discuss the existence and the characterisation of the solution.

**Proposition 3.2**.: _Under Assumptions 1, for every \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\) we have that (3.3) admits a unique solution \(t\mapsto\mu_{t}\) in the sense of Definition 3.1. Moreover, we have that for every \(t\in[0,T]\)_

\[\mu_{t}=\boldsymbol{\Phi}^{\theta}_{(0,t)\#}\mu_{0}. \tag{3.5}\]

Proof.: Existence and uniqueness of the measure solution of (3.3) follow from [1, Proposition 2.1, Theorem 3.1 and Remark 2.1]. 

From the characterisation of the solution of (3.3) provided in (3.5), it follows that the curve \(t\mapsto\mu_{t}\) inherits the properties of the flow map \(\Phi^{\theta}\) described in Proposition 2.1. These facts are collected in the next result.

**Proposition 3.3**.: _Let us fix \(T>0\) and \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\), and let us consider \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) satisfying Assumption 1. Let \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) be an admissible control, and let \(t\mapsto\mu_{t}\) be the corresponding solution of (3.3). Then, the curve \(t\mapsto\mu_{t}\) satisfies the properties listed below._

* _For every_ \(R>0\) _and_ \(\rho>0\)_, there exists_ \(\bar{R}>0\) _such that, for every_ \(t\in[0,T]\)_, it holds that_ \[\operatorname{supp}(\mu_{t})\subset B_{R}(0)\] _for every_ \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) _such that_ \(\|\theta\|_{L^{2}}\leq\rho\)_, and for every_ \(\mu_{0}\) _such that_ \(\operatorname{supp}(\mu_{0})\subset B_{R}(0)\)_._
* _For every_ \(R>0\) _and_ \(\rho>0\)_, there exists_ \(\bar{L}>0\) _such that, for every_ \(t\in[0,T]\)_, it holds that_ \[W_{1}(\mu_{t},\nu_{t})\leq\bar{L}W_{1}(\mu_{0},\nu_{0})\] _for every_ \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) _such that_ \(\|\theta\|_{L^{2}}\leq\rho\)_, and for every initial conditions_ \(\mu_{0},\nu_{0}\) _such that the supports satisfy_ \(\operatorname{supp}(\mu_{0}),\operatorname{supp}(\nu_{0})\subset B_{R}(0)\)_, where_ \(\mu_{t}=\boldsymbol{\Phi}^{\theta}_{(0,t)\#}\mu_{0}\) _and_ \(\nu_{t}=\boldsymbol{\Phi}^{\theta}_{(0,t)\#}\nu_{0}\)_._
* _For every_ \(R>0\) _and_ \(\rho>0\)_, there exists_ \(\bar{L}>0\) _such that, for every_ \(t_{1},t_{2}\in[0,T]\)_, it holds that_ \[W_{1}(\mu_{t_{1}},\mu_{t_{2}})\leq\bar{L}\cdot|t_{1}-t_{2}|^{\frac{1}{2}}\] _for every_ \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) _such that_ \(\|\theta\|_{L^{2}}\leq\rho\)_, and for every_ \(\mu_{0}\) _such that_ \(\operatorname{supp}(\mu_{0})\subset B_{R}(0)\)_._
* _For every_ \(R>0\) _and_ \(\rho>0\)_, there exists_ \(\bar{L}>0\) _such that, for every_ \(t\in[0,T]\)_, it holds that_ \[W_{1}(\mu_{t},\nu_{t})\leq\bar{L}\|\theta_{1}-\theta_{2}\|_{L^{2}}\] _for every_ \(\theta_{1},\theta_{2}\in L^{2}([0,T],\mathbb{R}^{m})\) _such that_ \(\|\theta\|_{L^{2}},\|\theta_{2}\|_{L^{2}}\leq\rho\)_, and for every initial condition_ \(\mu_{0}\) _such that_ \(\operatorname{supp}(\mu_{0})\subset B_{R}(0)\)_, where_ \(\mu_{t}=\boldsymbol{\Phi}_{(0,t)\#}^{\theta_{1}}\mu_{0}\) _and_ \(\nu_{t}=\boldsymbol{\Phi}_{(0,t)\#}^{\theta_{2}}\mu_{0}\)_._

Proof.: All the results follow from Proposition 3.2 and from the properties of the flow map presented in Proposition 2.1, combined with the Kantorovich duality (1.4) for the distance \(W_{1}\), and the change-of-variables formula (1.3). Since the argument is essentially the same for all the properties, we detail the computations only for the second point, i.e., the Lipschitz-continuous dependence on the initial distribution. Owing to (1.4), for any \(t\in[0,T]\), for any \(\varphi\in\operatorname{Lip}(\mathbb{R}^{2d})\) such that its Lipschitz constant \(\text{Lip}(\varphi)\leq 1\), it holds that

\[W_{1}(\mu_{t},\nu_{t})\leq\int_{\mathbb{R}^{2d}}\varphi(x,y)\,d(\mu_{t}-\nu_{ t})(x,y)=\int_{\mathbb{R}^{2d}}\varphi(\Phi_{(0,t)}^{\theta}(x),y)\,d(\mu_{0}- \nu_{0})(x,y)\leq\bar{L}W_{1}(\mu_{0},\nu_{0}),\]

where the equality follows from the definition of push-forward and from (3.2), while the constant \(\bar{L}\) in the second inequality descends from the local Lipschitz estimate of \(\Phi_{(0,t)}^{\theta}\) established in Proposition 2.1. 

### Mean-field optimal control

Using the transport equation (3.3), we can now formulate the mean-field optimal control problem that we aim to address. To this end, we introduce the functional \(J:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\), defined as follows:

\[J(\theta)=\left\{\begin{aligned} &\int_{\mathbb{R}^{2d}}\ell(x,y)\,d\mu_{T}(x, y)+\lambda\int_{0}^{T}|\theta(t)|^{2}\,dt,\\ &\text{s.t.}\ \left\{\begin{aligned} &\partial_{t}\mu_{t}(x,y)+ \nabla_{x}\cdot(\mathcal{F}(t,x,\theta_{t})\mu_{t}(x,y))=0\quad t\in[0,T],\\ &\mu_{t}|_{t=0}(x,y)=\mu_{0}(x,y),\end{aligned}\right. \end{aligned}\right. \tag{3.6}\]

for every admissible control \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\). The objective is to find the optimal control \(\theta^{*}\) that minimizes \(J(\theta^{*})\), subject to the PDE constraint (3.3) being satisfied by the curve \(t\mapsto\mu_{t}\). The term "mean-field" emphasizes that \(\theta\) is shared by an entire population of input-target pairs, and the optimal control must depend on the distribution of the initial data. We observe that when the initial measure \(\mu_{0}\) is empirical, i.e.

\[\mu_{0}:=\mu_{0}^{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{(X_{\delta}^{i},Y_{\delta }^{i})},\]

then minimization of (3.6) reduces to a classical finite particle optimal control problem with ODE constraints.

We now state the further regularity hypotheses that we require, in addition to the one contained in Assumption 1.

**Assumption 2**.: _For any given \(T>0\), the vector field \(\mathcal{F}\) satisfies the following._

* _For every_ \(R>0\) _there exists a constant_ \(L_{R}>0\) _such that, for every_ \(x_{1},x_{2}\in B_{R}(0)\)_, it holds_ \[|\nabla_{x}\mathcal{F}(t,x_{1},\theta)-\nabla_{x}\mathcal{F}(t,x_{2},\theta)| \leq L_{R}(1+|\theta|^{2})|x_{1}-x_{2}|,\quad\text{ for a.e. }t\in[0,T]\text{ and every }\theta\in\mathbb{R}^{m}.\]
* _There exists another constant_ \(L_{R}>0\) _such that, for every_ \(\theta_{1},\theta_{2}\in\mathbb{R}^{m}\)_, it holds_ \[|\nabla_{\theta}\mathcal{F}(t,x,\theta_{1})-\nabla_{\theta}\mathcal{F}(t,x, \theta_{2})|\leq L_{R}|\theta_{1}-\theta_{2}|,\quad\text{ for a.e. }t\in[0,T]\text{ and every }x\in B_{R}(0).\] _From this, it follows that_ \(|\nabla_{\theta}\mathcal{F}(t,x,\theta)|\leq L_{R}(1+|\theta|)\) _for every_ \(x\in B_{R}(0)\) _and for every_ \(\theta\in\mathbb{R}^{m}\)_._
* _There exists another constant_ \(L_{R}>0\) _such that, for every_ \(\theta_{1},\theta_{2}\in\mathbb{R}^{m}\)_, it holds_ \[|\nabla_{x}\mathcal{F}(t,x,\theta_{1})-\nabla_{x}\mathcal{F}(t,x,\theta_{2})| \leq L_{R}(1+|\theta_{1}|+|\theta_{2}|)|\theta_{1}-\theta_{2}|,\quad\text{ for a.e. }t\in[0,T]\text{ and every }x\in B_{R}(0).\] _From this, it follows that_ \(|\nabla_{x}\mathcal{F}(t,x,\theta)|\leq L_{R}(1+|\theta|^{2})\) _for every_ \(x\in B_{R}(0)\) _and for every_ \(\theta\in\mathbb{R}^{m}\)_._
* _There exists another constant_ \(L_{R}>0\) _such that_ \[|\nabla_{\theta}\mathcal{F}(t,x_{1},\theta)-\nabla_{\theta}\mathcal{F}(t,x_{2}, \theta)|\leq L_{R}(1+|\theta|)|x_{1}-x_{2}|,\quad\text{ for a.e. }t\in[0,T]\text{ and every }x_{1},x_{2}\in B_{R}(0).\]

Additionally, it is necessary to specify the assumptions on the function \(\ell\) that quantifies the discrepancy between the output of the network and its corresponding label.
**Assumption 3**.: _The function \(\ell:\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}_{+}\) is \(C^{1}\)-regular and non-negative. Moreover, for every \(R>0\), there exists a constant \(L_{R}>0\) such that, for every \(x_{1},x_{2}\in B_{R}(0)\), it holds_

\[|\nabla_{x}\ell(x_{1},y_{1})-\nabla_{x}\ell(x_{2},y_{2})|\leq L_{R}\left(|x_{1} -x_{2}|+|y_{1}-y_{2}|\right). \tag{3.7}\]

Let us begin by establishing a regularity result for the reduced final cost, which refers to the cost function without the regularization term.

**Lemma 3.4** (Differentiability of the cost).: _Let \(T,R>0\) and \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\) be such that \(\operatorname{supp}(\mu_{0})\subset B_{R}(0)\), and let us consider \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) and \(\ell:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\) that satisfy, respectively, Assumptions 1-2 and Assumption 3. Then, the reduced final cost_

\[J_{t}:\theta\in L^{2}([0,T];\mathbb{R}^{m})\mapsto\left\{\begin{array}{l} \int_{\mathbb{R}^{2d}}\ell(x,y)\,d\mu_{T}^{\theta}(x,y),\\ \mathrm{s.t.}\begin{cases}\partial_{t}\mu_{t}^{\theta}(x,y)+\nabla_{x}\big{(} \mathcal{F}(t,x,\theta_{t})\mu_{t}^{\theta}(x,y)\big{)}=0,\\ \mu_{t}^{\theta}|_{t=0}(x,y)=\mu_{0}(x,y),\end{cases}\end{array}\right. \tag{3.8}\]

_is Frechet-differentiable. Moreover, using the standard Hilbert space structure of \(L^{2}([0,T],\mathbb{R}^{m})\), we can represent the differential of \(J_{t}\) at the point \(\theta_{0}\) as the function:_

\[\nabla_{\theta}J_{\ell}(\theta):t\mapsto\int_{\mathbb{R}^{2d}}\nabla_{\theta} \mathcal{F}^{\top}\big{(}t,\Phi_{(0,t)}^{\theta_{0}}(x),\theta(t)\big{)}\cdot \mathcal{R}_{(t,T)}^{\theta}(x)^{\top}\cdot\nabla_{x}\ell^{\top}\big{(}\Phi_{ (0,T)}^{\theta}(x),y\big{)}\,d\mu_{0}(x,y) \tag{3.9}\]

_for a.e. \(t\in[0,T]\)._

Before proving the statement, we need to introduce the linear operator \(\mathcal{R}_{\tau,s}^{\theta}(x):\mathbb{R}^{d}\to\mathbb{R}^{d}\) with \(\tau,s\in[0,T]\), that is related to the linearization along a trajectory of the dynamics of the control system (2.1), and that appears in (3.9). Given an admissible control \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\), let us consider the corresponding trajectory curve \(t\mapsto\Phi_{(0,t)}^{\theta}(x)\) for \(t\in[0,T]\), i.e., the solution of (2.1) starting at the point \(x\in\mathbb{R}^{d}\) at the initial instant \(t=0\). Given any \(\tau\in[0,T]\), we consider the following linear ODE in the phase space \(\mathbb{R}^{d\times d}\):

\[\begin{cases}\frac{d}{ds}\mathcal{R}_{(\tau,s)}^{\theta}(x)=\nabla_{x} \mathcal{F}(s,\Phi_{(0,s)}^{\theta}(x),\theta(s))\cdot\mathcal{R}_{(\tau,s)}^ {\theta}(x)\quad\text{for a.e. }s\in[0,T],\\ \mathcal{R}_{(\tau,\tau)}^{\theta}(x)=\mathrm{Id}.\end{cases} \tag{3.10}\]

We insist on the fact that, when we write \(\mathcal{R}_{(\tau,s)}^{\theta}(x)\), \(x\) denotes the starting point of the trajectory along which the dynamics has been linearized. We observe that, using Assumption 2-\((iv)-(vi)\) and Caratheodory Theorem, it follows that (3.10) admits a unique solution, for every \(x\in\mathbb{R}^{d}\) and for every \(\tau\in[0,T]\). Since it is an elementary object in Control Theory, the properties of \(\mathcal{R}^{\theta}\) are discussed in the Appendix (see Proposition A.7). We just recall here that the following relation is satisfied:

\[\mathcal{R}_{\tau,s}^{\theta}(x)=\nabla_{x}\Phi_{(\tau,s)}^{\theta}|_{\Phi_{(0, \tau)}^{\theta}(x)}^{\theta} \tag{3.11}\]

for every \(\tau,s\in[0,T]\) and for every \(x\in\mathbb{R}^{d}\) (see, e.g., [10, Theorem 2.3.2]). Moreover, for every \(\tau,s\in[0,T]\) the following identity holds:

\[\mathcal{R}_{\tau,s}^{\theta}(x)\cdot\mathcal{R}_{s,\tau}^{\theta}(x)=\mathrm{ Id},\]

i.e., the matrices \(\mathcal{R}_{\tau,s}^{\theta}(x),\mathcal{R}_{s,\tau}^{\theta}(x)\) are one the inverse of the other. From this fact, it is possible to deduce that

\[\frac{\partial}{\partial\tau}\mathcal{R}_{\tau,s}^{\theta}(x)=-\mathcal{R}_{ \tau,s}^{\theta}(x)\cdot\nabla_{x}\mathcal{F}(\tau,\Phi_{(0,\tau)}^{\theta}(x ),\theta(\tau)) \tag{3.12}\]

for almost every \(\tau,s\in[0,T]\) (see, e.g., [10, Theorem 2.2.3] for the details).

Proof of Lemma 3.4.: Let us fix an admissible control \(\theta\in L^{2}([0,T];\mathbb{R}^{m})\) and let \(\mu_{\cdot}^{\theta}\in\mathcal{C}^{0}([0,T];\mathcal{P}_{c}(\mathbb{R}^{2d}))\) be the unique solution of the continuity equation (3.3), corresponding to the control \(\theta\) and satisfying \(\mu^{\theta}|_{t=0}=\mu_{0}\). According to Proposition 3.2, this curve can be expressed as \(\mu_{t}^{\theta}=\mathbf{\Phi}_{(0,t)\#}^{\theta}\mu_{0}\) for every \(t\in[0,T]\), where the map \(\mathbf{\Phi}_{(0,t)}^{\theta}=(\Phi_{(0,t)}^{\theta},\mathrm{Id}):\mathbb{R}^{2 d}\to\mathbb{R}^{2d}\) has been introduced in (3.2) as the flow of the extended control system (3.1). In particular, we can rewrite the terminal cost \(J_{\ell}\) defined in (3.8) as

\[J_{\ell}(\theta)=\int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta}(x),y)d\mu_{0}(x,y).\]
In order to compute the gradient \(\nabla_{\theta}J_{\ell}\), we preliminarily need to focus on the differentiability with respect to \(\theta\) of the mapping \(\theta\mapsto\ell(\Phi^{\theta}_{(0,T)}(x),y)\), when \((x,y)\) is fixed. Indeed, given another control \(\vartheta\in L^{2}([0,T];\mathbb{R}^{m})\) and \(\varepsilon>0\), from Proposition A.6 it descends that

\[\begin{split}\Phi^{\theta+\varepsilon\vartheta}_{(0,T)}(x)& =\Phi^{\theta}_{(0,T)}(x)+\varepsilon\xi^{\theta}(T)+o_{\theta}( \varepsilon)\\ &=\Phi^{\theta}_{(0,T)}(x)+\varepsilon\int_{0}^{T}\mathcal{R}^{ \theta}_{(s,T)}(x)\nabla_{\theta}\mathcal{F}(s,\Phi^{\theta}_{(0,s)}(x), \theta(s))\vartheta(s)ds+o_{\theta}(\varepsilon)\end{split}\qquad \text{ as }\varepsilon\to 0, \tag{3.13}\]

where \(o_{\theta}(\varepsilon)\) is uniform for every \(x\in B_{R}(0)\subset\mathbb{R}^{d}\), and as \(\vartheta\) varies in the unit ball of \(L^{2}\). Owing to Assumption 3, for every \(x,y,v\in B_{R}(0)\) we observe that

\[|\ell(x+\varepsilon v+o(\varepsilon),y)-\ell(x,y)-\varepsilon\nabla_{x}\ell (x,y)\cdot v|\leq|\nabla_{x}\ell(x,y)|o(\varepsilon)+\frac{1}{2}L_{R}| \varepsilon v+o(\varepsilon)|^{2}\qquad\text{as }\varepsilon\to 0. \tag{3.14}\]

Therefore, combining (3.13) and (3.14), we obtain that

\[\ell(\Phi^{\theta+\varepsilon\vartheta}_{(0,T)}(x),y)-\ell(\Phi^{\theta}_{(0,T)}(x),y)=\varepsilon\int_{0}^{T}\big{(}\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}( x),y)\cdot\mathcal{R}^{\theta}_{(s,T)}(x)\cdot\nabla_{\theta}\mathcal{F}(s, \Phi^{\theta}_{(0,s)}(x),\theta(s))\big{)}\cdot\vartheta(s)ds+o_{\theta}( \varepsilon).\]

Since the previous expression is uniform for \(x,y\in B_{R}(0)\), then if we integrate both sides of the last identity with respect to \(\mu_{0}\), we have that

\[J_{\ell}(\theta+\varepsilon\vartheta)-J_{\ell}(\theta)=\varepsilon\int_{ \mathbb{R}^{2d}}\int_{0}^{T}\big{(}\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}(x),y) \cdot\mathcal{R}^{\theta}_{(s,T)}(x)\cdot\nabla_{\theta}\mathcal{F}(s,\Phi^{ \theta}_{(0,s)}(x),\theta(s))\big{)}\cdot\vartheta(s)\,ds\,d\mu_{0}(x,y)+o_{ \theta}(\varepsilon). \tag{3.15}\]

This proves the Frechet differentiability of the functional \(J_{\ell}\) at the point \(\theta\). We observe that, from Proposition 2.1, Proposition A.7 and Assumption 2, it follows that the function \(s\mapsto\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}(x),y)\cdot\mathcal{R}^{\theta}_{ (s,T)}(x)\cdot\nabla_{\theta}\mathcal{F}(s,\Phi^{\theta}_{(0,s)}(x),\theta(s))\) is uniformly bounded in \(L^{2}\), as \(x,y\) vary in \(B_{R}(0)\subset\mathbb{R}^{d}\). Then, using Fubini Theorem, the first term of the expansion (3.15) can be rewritten as

\[\int_{0}^{T}\left(\int_{\mathbb{R}^{2d}}\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}( x),y)\cdot\mathcal{R}^{\theta}_{(s,T)}(x)\cdot\nabla_{\theta}\mathcal{F}(s, \Phi^{\theta}_{(0,s)}(x),\theta(s))\,d\mu_{0}(x,y)\right)\cdot\vartheta(s)\,ds.\]

Hence, from the previous asymptotic expansion and from Riesz Representation Theorem, we deduce (3.9). 

We now prove the most important result of this subsection, concerning the Lipschitz regularity of the gradient \(\nabla_{\theta}J_{\ell}\).

**Proposition 3.5**.: _Under the same assumptions and notations as in Lemma 3.4, we have that the gradient \(\nabla_{\theta}J_{\ell}:L^{2}([0,T],\mathbb{R}^{m})\to L^{2}([0,T],\mathbb{R}^ {m})\) is Lipschitz-continuous on every bounded set of \(L^{2}\). More precisely, given \(\theta_{1},\theta_{2}\in L^{2}([0,T];\mathbb{R}^{m})\), there exists a constant \(\mathcal{L}(T,R,\|\theta_{1}\|_{L^{2}},\|\theta_{2}\|_{L^{2}})>0\) such that_

\[\big{\|}\nabla_{\theta}J_{\ell}(\theta_{1})-\nabla_{\theta}J_{\ell}(\theta_{2}) \big{\|}_{L^{2}}\leq\mathcal{L}(T,R,\|\theta_{1}\|_{L^{2}},\|\theta_{2}\|_{L^{ 2}})\,\big{\|}\theta_{1}-\theta_{2}\big{\|}_{L^{2}}.\]

Proof.: Let us consider two admissible controls \(\theta_{1},\theta_{2}\in L^{2}([0,T],\mathbb{R}^{m})\) such that \(\|\theta_{1}\|_{L^{2}},\|\theta_{1}\|_{L^{2}}\leq C\). In order to simplify the notations, given \(x\in B_{R}(0)\subset\mathbb{R}^{d}\), we define the curves \(x_{1}:[0,T]\to\mathbb{R}^{d}\) and \(x_{2}:[0,T]\to\mathbb{R}^{d}\) as

\[x_{1}(t):=\Phi^{\theta_{1}}_{(0,t)}(x),\quad x_{2}(t):=\Phi^{\theta_{2}}_{(0,t)}(x)\]

for every \(t\in[0,T]\), where the flows \(\Phi^{\theta_{1}},\Phi^{\theta_{2}}\) where introduced in (2.2). We recall that, in virtue of Proposition 2.1, \(x_{1}(t),x_{2}(t)\in B_{R}(0)\) for every \(t\in[0,1]\). Then, for every \(y\in B_{R}(0)\), we observe that

\[\begin{split}\Big{|}\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1} (t),&\theta_{1}(t)\big{)}\mathcal{R}^{\theta_{1}}_{(t,T)}(x^{\top} )^{\top}\nabla_{x}\ell^{\top}\big{(}x_{1}(T),y\big{)}-\nabla_{\theta}\mathcal{F} ^{\top}\big{(}t,x_{1}(t),\theta_{2}(t)\big{)}\mathcal{R}^{\theta_{2}}_{(t,T)}(x ^{\top})^{\top}\nabla_{x}\ell^{\top}\big{(}x_{2}(T),y\big{)}\Big{|}\\ &\leq\Big{|}\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t), \theta_{1}(t)\big{)}\Big{|}\Big{|}\mathcal{R}^{\theta_{1}}_{(t,T)}(x^{\top})^{ \top}\Big{|}\,\Big{|}\nabla_{x}\ell^{\top}\big{(}x_{1}(T),y\big{)}-\nabla_{x} \ell^{\top}\big{(}x_{2}(T),y\big{)}\Big{|}\\ &+\Big{|}\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t),\theta_{ 1}(t)\big{)}\Big{|}\mathcal{R}^{\theta_{1}}_{(t,T)}(x)^{\top}-\mathcal{R}^{ \theta_{2}}_{(t,T)}(x)^{\top}\Big{|}\,\Big{|}\nabla_{x}\ell^{\top}\big{(}x_{2}(T),y \big{)}\Big{|}\\ &+\Big{|}\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t),\theta_{ 1}(t)\big{)}-\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{2}(t),\theta_{2}(t) \big{)}\Big{|}\,\Big{|}\mathcal{R}^{\theta_{2}}_{(t,T)}(x)^{\top}\Big{|}\, \Big{|}\nabla_{x}\ell^{\top}\big{(}x_{2}(T),y\big{)}\Big{|}\end{split} \tag{3.16}\]

for a.e. \(t\in[0,T]\). We bound separately the three terms at the right-hand side of (3.16). As regards the first addend, from Assumption 2-(\(v\)), Assumption 3, Proposition A.7 and Lemma A.4, we deduce that there exists a positive constant \(C_{1}>0\) such that

\[\begin{split}\Big{|}\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1} (t),\theta_{1}(t)\big{)}\Big{|}\,\Big{|}\mathcal{R}^{\theta_{1}}_{(t,T)}(x)^{ \top}\Big{|}&\nabla_{x}\ell^{\top}\big{(}x_{1}(T),y\big{)}- \nabla_{x}\ell^{\top}\big{(}x_{2}(T),y\big{)}\Big{|}\\ &\leq C_{1}\left(1+|\theta_{1}(t)|\right)\|\theta_{1}-\theta_{2}\|_{L^ {2}}\end{split} \tag{3.17}\]
for a.e. \(t\in[0,T]\). Similarly, using again Assumption 2-\((v)\), Assumption 3, and Proposition A.7 on the second addend at the right-hand side of (3.16), we obtain that there exists \(C_{2}>0\) such that

\[\left|\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t),\theta_{1}(t)\big{)} \right|\Big{|}\mathcal{R}_{(t,T)}^{\theta_{1}}(x)^{\top}-\mathcal{R}_{(t,T)}^{ \theta_{2}}(x)^{\top}\Big{|}\left|\nabla_{x}\ell^{\top}\big{(}x_{2}(T),y\big{)} \right|\leq C_{2}\left(1+|\theta_{1}(t)|\right)\|\theta_{1}-\theta_{2}\|_{L^{ 2}} \tag{3.18}\]

for a.e. \(t\in[0,T]\). Moreover, the third term can be bounded as follows:

\[\left|\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t),\theta_{1}(t)\big{)} \right|-\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{2}(t),\theta_{2}(t) \big{)}\Big{|}\left|\mathcal{R}_{(t,T)}^{\theta_{2}}(x)^{\top}\right|\Big{|} \nabla_{x}\ell^{\top}\big{(}x_{2}(T),y\big{)}\Big{|} \tag{3.19}\]

for a.e. \(t\in[0,T]\), where we used Assumption 2-\((v)-(vii)\), Proposition A.7 and Lemma A.4. Therefore, combining (3.16)-(3.19), we deduce that

\[\left|\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t),\theta_{1}(t)\big{)} \mathcal{R}_{(t,T)}^{\theta_{1}}(x)^{\top}\nabla_{x}\ell^{\top}\big{(}x_{1}( T),y\big{)}-\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,x_{1}(t),\theta_{2}(t) \big{)}\mathcal{R}_{(t,T)}^{\theta_{2}}(x)^{\top}\nabla_{x}\ell^{\top}\big{(} x_{2}(T),y\big{)}\right| \tag{3.20}\]

for a.e. \(t\in[0,T]\). We observe that the last inequality holds for every \(x,y\in B_{R}(0)\). Therefore, if we integrate both sides of (3.20) with respect to the probability measure \(\mu_{0}\), recalling the expression of the gradient of \(J_{\ell}\) reported in (3.9), we have that

\[\left|\nabla_{\theta}J_{\ell}(\theta_{1})[t]-\nabla_{\theta}J_{\ell}(\theta_ {1})[t]\right|\leq\bar{C}\Big{[}(1+|\theta_{1}(t)|)\|\theta_{1}-\theta_{2}\|_ {L^{2}}+|\theta_{1}(t)-\theta_{2}(t)|\Big{]} \tag{3.21}\]

for a.e. \(t\in[0,T]\), and this concludes the proof. 

From the previous result we can deduce that the terminal cost \(J_{\ell}:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) is locally semi-convex.

**Corollary 3.6** (Local semiconvexity of the cost functional).: _Under the same assumptions and notations as in Lemma 3.4, let us consider a bounded subset \(\Gamma\subset L^{2}([0,T];\mathbb{R}^{m})\). Then, \(\nabla_{\theta}J:L^{2}([0,T])\to L^{2}([0,T])\) is Lipschitz continuous on \(\Gamma\). Moreover, there exists a constant \(\mathcal{L}(T,R,\Gamma)>0\) such that the cost functional \(J:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) defined in (3.6) satisfies the following semiconvexity estimate:_

\[J\big{(}(1-\zeta)\theta_{1}+\zeta\theta_{2}\big{)}\leq(1-\zeta)J(\theta_{1})+ \zeta J(\theta_{2})-(2\lambda-\mathcal{L}(T,R,\Gamma))\tfrac{\zeta(1-\zeta)} {2}\|\theta_{1}-\theta_{2}\|_{L^{2}}^{2} \tag{3.22}\]

_for every \(\theta_{1},\theta_{2}\in\Gamma\) and for every \(\zeta\in[0,1]\). In particular, if \(\lambda>\frac{1}{2}\mathcal{L}(T,R,\Gamma)\), the cost functional \(J\) is strictly convex over \(\Gamma\)._

Proof.: We recall that \(J(\theta)=J_{\ell}(\theta)+\lambda\|\theta\|_{L^{2}}^{2}\), where \(J_{\ell}\) has been introduced in (3.8). Owing to Proposition 3.5, it follows that \(\nabla_{\theta}J_{\ell}\) is Lipschitz continuous on \(\Gamma\) with constant \(\mathcal{L}(T,R,\Gamma)\). This implies that \(J\) is Lipschitz continuous as well on \(\Gamma\). Moreover, it descends that

\[J_{\ell}\big{(}(1-\zeta)\theta_{1}+\zeta\theta_{2}\big{)}\leq(1-\zeta)J_{\ell} (\theta_{1})+\zeta J_{\ell}(\theta_{2})+\mathcal{L}(T,R,\Gamma)\tfrac{\zeta(1- \zeta)}{2}\|\theta_{1}-\theta_{2}\|_{L^{2}}^{2}\]

for every \(\theta_{1},\theta_{2}\in\Gamma\) and for every \(\zeta\in[0,1]\). On the other hand, recalling that

\[\|(1-\zeta)\theta_{1}+\zeta\theta_{2}\|_{L^{2}}^{2}=(1-\zeta)\|\theta_{1}\|_{L ^{2}}^{2}+\zeta\|\theta_{2}\|_{L^{2}}^{2}-\zeta(1-\zeta)\|\theta_{1}-\theta_{2} \|_{L^{2}}^{2}\]

for every \(\theta_{1},\theta_{2}\in L^{2}\), we immediately deduce (3.22). 

**Remark 3.1**.: _When the parameter \(\lambda>0\) that tunes the \(L^{2}\)-regularization is large enough, we can show that the functional \(J\) defined by (3.6) admits a unique global minimizer. Indeed, since the control identically \(0\) is an admissible competitor, we have that_

\[\inf_{\theta\in L^{2}}J(\theta)\leq J(0)=J_{\ell}(0),\]

_where we observe that the right-hand side is not affected by the value of \(\lambda\). Hence, recalling that \(J(\theta)=J_{\ell}(\theta)+\lambda\|\theta\|_{L^{2}}^{2}\), we have that the sublevel set \(\{\theta:J(\theta)\leq J_{\ell}(0)\}\) is included in the ball \(B_{\lambda}:=\{\theta:\|\theta\|_{L^{2}}^{2}\leq\frac{1}{\lambda}J_{\ell}(0)\}\). Since these balls are decreasing as \(\lambda\) increases, owing to Corollary 3.6, we deduce that there exists a parameter \(\bar{\lambda}>0\) such that the cost functional \(J\) is strongly convex when restricted to \(B_{\bar{\lambda}}\). Then, Lemma 3.4 guarantees that the functional \(J:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) introduced in (3.6) is continuous with respect to the strong topology of \(L^{2}\), while the convexity implies that it is weakly lower semi-continuous as well. Being the ball \(B_{\bar{\lambda}}\) weakly compact, we deduce that the restriction to \(B_{\bar{\lambda}}\) of the functional \(J\) admits a unique minimizer \(\theta^{*}\). However, since \(B_{\bar{\lambda}}\) includes the sublevel set \(\{\theta:J(\theta)\leq J_{\ell}(0)\}\), it follows that \(\theta^{*}\) is actually the unique global minimizer. It is interesting to observe that, even though \(\lambda\) is chosen large enough to ensure existence (and uniqueness) of the global minimizer, it is not possible to conclude that the functional \(J\) is globally convex. This is essentially due to the fact that Corollary 3.6 holds only on bounded subsets of \(L^{2}\)._
Taking advantage of the representation of the gradient of the terminal cost \(J_{\ell}\) provided by (3.9), we can formulate the necessary optimality conditions for the cost \(J\) introduced in (3.6). In order to do that, we introduce the function \(p:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d}\) as follows:

\[p_{t}(x,y):=\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}(x),y)\cdot\mathcal{R}^{\theta }_{(t,T)}(x), \tag{3.23}\]

where \(\mathcal{R}^{\theta}_{(t,T)}(x)\) is defined according to (3.10). We observe that \(p\) (as well as \(\nabla_{x}\ell\)) should be understood as a row vector. Moreover, using (3.12), we deduce that, for every \(x,y\in\mathbb{R}^{d}\), the \(t\mapsto p_{t}(x,y)\) is solving the following backward Cauchy problem:

\[\frac{\partial}{\partial t}p_{t}(x,y)=-p_{t}(x,y)\cdot\nabla_{x}\mathcal{F}(t, \Phi^{\theta}_{(0,t)}(x),\theta(t)),\qquad p_{T}(x,y)=\nabla_{x}\ell(\Phi^{ \theta}_{(0,T)}(x),y). \tag{3.24}\]

Hence, we can equivalently rewrite \(\nabla_{\theta}J_{\ell}\) using \(p\):

\[\nabla_{\theta}J_{\ell}(\theta)[t]=\int_{\mathbb{R}^{2d}}\nabla_{\theta} \mathcal{F}^{\top}\big{(}t,\Phi^{\theta}_{(0,t)}(x),\theta(t)\big{)}\cdot p_ {t}^{\top}(x,y)\,d\mu_{0}(x,y) \tag{3.25}\]

for almost every \(t\in[0,T]\). Therefore, recalling that \(J(\theta)=J_{\ell}(\theta)+\lambda\|\theta\|_{L^{2}}^{2}\), we deduce that the stationary condition \(\nabla_{\theta}J(\theta^{*})=0\) can be rephrased as

\[\begin{cases}\partial_{t}\mu_{t}^{*}(x,y)+\nabla_{x}\cdot\big{(}\mathcal{F}(t,x,\theta^{*}(t))\mu_{t}^{*}(x,y)\big{)}=0,&\mu_{t}^{*}|_{t=0}(x,y)=\mu_{0}(x,y ),\\ \partial_{t}p_{t}^{*}(x,y)=-p_{t}^{*}(x,y)\cdot\nabla_{x}\mathcal{F}(t,\Phi^{ \theta^{*}}_{(0,t)}(x),\theta^{*}(t)),&p_{t}^{*}|_{t=T}(x,y)=\nabla_{x}\ell( \Phi^{\theta^{*}}_{(0,T)}(x),y),\\ \theta^{*}(t)=-\frac{1}{2\lambda}\int_{\mathbb{R}^{2d}}\nabla_{\theta} \mathcal{F}^{\top}\big{(}t,\Phi^{\theta^{*}}_{(0,t)}(x),\theta^{*}(t)\big{)} \cdot p_{t}^{*\top}(x,y)\,d\mu_{0}(x,y).\end{cases} \tag{3.26}\]

**Remark 3.2**.: _The computation of \(p\) through the backward integration of (3.24) can be interpreted as the control theoretic equivalent of the "back-propagation of the gradients". We observe that, in order to check whether (3.26) is satisfied, it is sufficient to evaluate \(p^{*}\) only on \(\operatorname{supp}(\mu_{0})\). Moreover, the evaluation of \(p^{*}\) on different points \((x_{1},y_{1}),(x_{2},y_{2})\in\operatorname{supp}(\mu_{0})\) involves the resolution of two uncoupled backward ODEs. This means that, when dealing with a measure \(\mu_{0}\) that charges only finitely many points, we can solve the equation (3.24) in parallel for every point in \(\operatorname{supp}(\mu_{0})\)._

In virtue of Proposition 3.5, we can study the gradient flow induced by the cost functional \(J:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) on its domain. More precisely, given an admissible control \(\theta_{0}\in L^{2}([0,T],\mathbb{R}^{m})\), we consider the gradient flow equation:

\[\begin{cases}\dot{\theta}(\omega)=-\nabla_{\theta}J(\theta(\omega))\quad \text{for $\omega\geq 0$},\\ \theta(0)=\theta_{0}.\end{cases} \tag{3.27}\]

In the next result we show that the gradient flow equation (3.27) is well-posed and that the solution is defined for every \(\omega\geq 0\). In the particular case of linear-control systems, the properties of the gradient flow trajectories has been investigated in [42].

**Lemma 3.7**.: _Let \(T,R>0\) and \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\) be a probability measure such that \(\operatorname{supp}(\mu_{0})\subset B_{R}(0)\), and let us consider \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) and \(\ell:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\) that satisfy, respectively, Assumptions 1-2 and Assumption 3. Then, for every \(\theta_{0}\in L^{2}([0,T],\mathbb{R}^{m})\), the gradient flow equation (3.27) admits a unique solution \(\omega\mapsto\theta(\omega)\) of class \(C^{1}\) that is defined for every \(\omega\in[0,+\infty)\)._

Proof.: Let us consider \(\theta_{0}\in L^{2}([0,T],\mathbb{R}^{m})\), and let us introduce the sub-level set

\[\Gamma:=\{\theta\in L^{2}([0,T],\mathbb{R}^{m}):J(\theta)\leq J(\theta_{0})\},\]

where \(J\) is the functional introduced in (3.6) defining the mean-field optimal control problem. Using the fact that the end-point cost \(\ell:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}\) is non-negative, we deduce that \(\Gamma\subset\{\theta\in L^{2}([0,T],\mathbb{R}^{m}):\|\theta\|_{L^{2}}^{2}\leq \frac{1}{\lambda}J(\theta_{0})\}\). Hence, from Proposition 3.5 it follows that the gradient field \(\nabla_{\theta}J\) is Lipschitz (and bounded) on \(\Gamma\). Hence, using a classical result on ODE in Banach spaces (see, e.g., [31, Theorem 5.1.1]), it follows that the initial value problem (3.27) admits a unique small-time solution \(\omega\mapsto\theta(\omega)\) of class \(C^{1}\) defined for \(\omega\in[-\delta,\delta]\), with \(\delta>0\). Moreover, we observe that

\[\frac{d}{d\omega}J(\theta(\omega))=\langle\nabla_{\theta}J(\theta(\omega)), \dot{\theta}(\omega))\rangle=-\|\nabla_{\theta}J(\theta(\omega))\|_{L^{2}}\leq 0,\]

and this implies that \(\theta(\omega)\in\Gamma\) for every \(\omega\in[0,\delta]\). Hence, it is possible to recursively extend the solution to every interval of the form \([0,M]\), with \(M>0\)
We observe that, under the current working assumptions, we cannot provide any convergence result for the gradient flow trajectories. This is not surprising since, when the regularization parameter \(\lambda>0\) is small, it is not even possible to prove that the functional \(J\) admits minimizers. Indeed, the argument presented in Remark 3.1 requires the regularization parameter \(\lambda\) to be sufficiently large.

We conclude the discussion with an observation on a possible discretization of (3.27). If we fix a sufficiently small parameter \(\tau>0\), given an initial guess \(\theta_{0}\), we can consider the sequence of controls \((\theta_{k}^{\tau})_{k\geq 0}\subset L^{2}([0,T],\mathbb{R}^{m})\) defined through the Minimizing Movement Scheme:

\[\theta_{0}^{\tau}=\theta_{0},\qquad\theta_{k+1}^{\tau}\in\arg\min_{\theta} \left[J(\theta)+\frac{1}{2\tau}\|\theta-\theta_{k}^{\tau}\|_{L^{2}}^{2}\right] \quad\text{for every }k\geq 0. \tag{3.28}\]

**Remark 3.3**.: _We observe that the minimization problems in (3.28) are well-posed as soon as the functionals \(\theta\mapsto J_{\theta_{k}}^{\tau}(\theta):=J(\theta)+\frac{1}{2\tau}\| \theta-\theta_{k}^{\tau}\|_{L^{2}}^{2}\) are strictly convex on the bounded sublevel set \(K_{\theta_{0}}:=\{\theta:J(\theta)\leq J(\theta_{0}^{\tau})\}\), for every \(k\geq 0\). Hence, the parameter \(\tau>0\) can be calibrated by means of the estimates provided by Corollary 3.6, considering the bounded set \(K_{\theta_{0}}\). Then, using and inductive argument, it follows that, for every \(k\geq 0\), the functional \(J_{\theta_{k}}^{\tau}:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) admits a unique global minimizer \(\theta_{k+1}^{\tau}\). Also for \(J_{\theta_{k}}^{\tau}\) we can derive the necessary conditions for optimality satisfied by \(\theta_{k+1}\), that are analogous to the ones formulated in (3.26), and that descend as well from the identity \(\nabla_{\theta}J_{\theta_{k}}^{\tau}(\theta_{k+1}^{\tau})=0\):_

\[\begin{cases}\partial_{t}\mu_{t}(x,y)+\nabla_{x}\cdot\big{(}\mathcal{F}(t,x, \theta_{k+1}^{\tau}(t))\mu_{t}(x,y)\big{)}=0,&\mu_{t}|_{t=0}(x,y)=\mu_{0}(x,y), \\ \partial_{t}p_{t}(x,y)=-p_{t}(x,y)\cdot\nabla_{x}\mathcal{F}(t,\Phi_{(0,t)}^{ \sigma_{t+1}^{\tau}}(x),\theta_{k+1}^{\tau}(t)),&p_{t}|_{t=T}(x,y)=\nabla_{x} \ell(\Phi_{(0,T)}^{\theta_{k+1}^{\tau}}(x),y),\\ \theta_{k+1}^{\tau}(t)=-\frac{1}{1+2\lambda\tau}\left(\theta_{k}^{\tau}(t)- \tau\int_{\mathbb{R}^{2d}}\nabla_{\theta}\mathcal{F}^{\top}\big{(}t,\Phi_{(0, t)}^{\theta_{k+1}^{\tau}}(x),\theta_{k+1}^{\tau}(t)\big{)}\cdot p_{t}^{ \tau}(x,y)\,d\mu_{0}(x,y)\right).\end{cases} \tag{3.29}\]

_Finally, we observe that the mapping \(\Lambda^{\tau}:L^{2}([0,T],\mathbb{R}^{m})\to L^{2}([0,T],\mathbb{R}^{m})\) defined for a.e. \(t\in[0,T]\) as_

\[\Lambda_{\theta_{k}^{\tau}}^{\tau}(\theta)[t]:=-\frac{1}{1+2\lambda\tau}\left( \theta_{k}^{\tau}(t)-\tau\int_{\mathbb{R}^{2d}}\nabla_{\theta}\mathcal{F}^{ \top}\big{(}t,\Phi_{(0,t)}^{\theta}(x),\theta(t)\big{)}\cdot p_{t}^{\top}(x,y )\,d\mu_{0}(x,y)\right) \tag{3.30}\]

_is a contraction on \(K_{\theta_{0}}\) as soon as_

\[\frac{\tau}{1+2\lambda\tau}\mathrm{Lip}\left(\nabla_{\theta}J_{\ell}|_{K_{ \theta_{0}}}\right)<1.\]

For every \(\tau>0\) such that the sequence \((\theta_{k}^{\tau})_{k\geq 0}\) is defined, we denote with \(\tilde{\theta}^{\tau}:[0,+\infty)\to L^{2}([0,T],\mathbb{R}^{m})\) the piecewise affine interpolation obtained as

\[\tilde{\theta}^{\tau}(\omega)=\theta_{k}^{\tau}+\frac{\theta_{k+1}^{\tau}- \theta_{k}^{\tau}}{\tau}(\omega-k\tau)\quad\text{for }\omega\in[k\tau,(k+1)\tau]. \tag{3.31}\]

We finally report a classical result concerning the convergence of the piecewise affine interpolation \(\tilde{\theta}^{\tau}\) to the gradient flow trajectory solving (3.27).

**Proposition 3.8**.: _Under the same assumptions and notations as in Lemma 3.7, let us consider an initial point \(\theta_{0}\in L^{2}([0,T],\mathbb{R}^{m})\) and a sequence \((\tau_{j})_{j\in\mathbb{N}}\) such that \(\tau_{j}\to 0\) as \(j\to\infty\), and let \((\tilde{\theta}^{\tau_{j}})_{j\in\mathbb{N}}\) be the sequence of piecewise affine curves defined by (3.31). Then, for every \(\Omega>0\), there exists a subsequence \((\tilde{\theta}^{\tau_{j_{k}}})_{k\in\mathbb{N}}\) converging uniformly on the interval \([0,\Omega]\) to the solution of (3.27) starting from \(\theta_{0}\)._

Proof.: The proof follows directly from [41, Proposition 2.3]. 

### Finite particles approximation

In this section, we study the stability of the mean-field optimal control problem (3.6) with respect to finite-samples distributions. More precisely, assume that we are given samples \(\{(X_{0}^{i},Y_{0}^{i})\}_{i=1}^{N}\) of size \(N\geq 1\) independently and identically distributed according to \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\), and consider the empirical loss minimization problem

\[\inf_{\theta\in L^{2}([0,T];\mathbb{R}^{m})}J^{N}(\theta):=\begin{cases}\frac{1 }{N}\sum_{i=1}^{N}\ell\big{(}X^{i}(T),Y^{i}(T)\big{)}+\lambda\int_{0}^{T}| \theta(t)|^{2}\,dt\\ \text{s.t.}\begin{cases}\dot{X}^{i}(t)=\mathcal{F}(t,X^{i}(t),\theta(t)),&\dot {Y}^{i}(t)=0,\\ \big{(}X^{i}(t),Y^{i}(t)\big{)}\big{|}_{t=0}=(X_{0}^{i},Y_{0}^{i}),&i\in\{1, \ldots,N\}.\end{cases}\end{cases} \tag{3.32}\]
By introducing the empirical measure \(\mu_{0}^{N}\in\mathcal{P}_{c}^{N}(\mathbb{R}^{2d})\), defined as

\[\mu_{0}^{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{(X_{0}^{i},Y_{0}^{i})},\]

the cost function in (3.32) can be rewritten as

\[J^{N}(\theta)=\int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta}(x),y)\,d\mu_{0} ^{N}(x,y)+\lambda\|\theta\|_{L^{2}}^{2} \tag{3.33}\]

for every \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\), and the empirical loss minimization problem in (3.32) can be recast as a mean-field optimal control problem with initial datum \(\mu_{0}^{N}\). In this section we are interested in studying the asymptotic behavior of the functional \(J^{N}\) as \(N\) tends to infinity. More precisely, we consider a sequence of probability measures \((\mu_{0}^{N})_{N\geq 1}\) such that \(\mu_{0}^{N}\) charges uniformly \(N\) points, and such that

\[W_{1}(\mu_{0}^{N},\mu_{0})\ \underset{N\to+\infty}{\longrightarrow}\ 0.\]

Then, in Proposition 3.9 we study the uniform convergence of \(J^{N}\) and of \(\nabla_{\theta}J^{N}\) to \(J\) and \(\nabla_{\theta}J^{N}\), respectively, where \(J:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) is the functional defined in (3.6) and corresponding to the limiting measure \(\mu_{0}\). Moreover, in Theorem 3.10, assuming the existence of a region where the functionals \(J^{N}\) are uniformly strongly convex, we provide an estimate of the so called _generalization error_ in terms of the distance \(W_{1}(\mu_{0}^{N},\mu_{0})\).

**Proposition 3.9** (Uniform convergence of \(J^{N}\) and \(\nabla_{\theta}J^{N}\)).: _Let us consider a probability measure \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\) and a sequence \((\mu_{0}^{N})_{N\geq 1}\) such that \(\mu_{0}^{N}\in\mathcal{P}_{c}^{N}(\mathbb{R}^{2d})\) for every \(N\geq 1\). Let us further assume that \(W_{1}(\mu_{0}^{N},\mu_{0})\to 0\) as \(N\to\infty\), and that there exists \(R>0\) such that \(\operatorname{supp}(\mu_{0}),\operatorname{supp}(\mu_{0}^{N})\subset B_{R}(0)\) for every \(N\geq 1\). Given \(T>0\), let \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) and \(\ell:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\) satisfy, respectively, Assumptions 1-2 and Assumption 3, and let \(J,J^{N}:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) be the cost functionals defined in (3.6) and (3.32), respectively. Then, for every bounded subset \(\Gamma\subset L^{2}([0,T],\mathbb{R}^{m})\), we have that_

\[\lim_{N\to\infty}\ \sup_{\theta\in\Gamma}|J^{N}(\theta)-J(\theta)|=0 \tag{3.34}\]

_and_

\[\lim_{N\to\infty}\ \sup_{\theta\in\Gamma}\|\nabla_{\theta}J^{N}(\theta)-\nabla_{ \theta}J(\theta)\|_{L^{2}}=0, \tag{3.35}\]

_where \(J\) was introduced in (3.6), and \(J^{N}\) is defined as in (3.33)._

Proof.: Since we have that \(J(\theta)=J_{\ell}(\theta)+\lambda\|\theta\|_{L^{2}}^{2}\) and \(J^{N}(\theta)=J_{\ell}^{N}(\theta)+\lambda\|\theta\|_{L^{2}}^{2}\), it is sufficient to prove (3.34)-(3.35) for \(J_{\ell}\) and \(J_{\ell}^{N}\), where we set

\[J_{\ell}^{N}(\theta):=\int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta}(x),y)\, d\mu_{0}^{N}(x,y)\]

for every \(\theta\in L^{2}\) and for every \(N\geq 1\). We first observe that, for every \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) such that \(\|\theta\|_{L^{2}}\leq\rho\), from Proposition 3.3 it follows that \(\operatorname{supp}(\mu_{0}),\operatorname{supp}(\mu_{0}^{N})\subset B_{R}(0)\), for some \(\bar{R}>0\). Then, denoting with \(t\mapsto\mu_{t}^{N}\) and \(t\mapsto\mu_{t}\) the solutions of the continuity equation (3.3) driven by the control \(\theta\) and with initial datum, respectively, \(\mu_{0}^{N}\) and \(\mu_{0}\), we compute

\[\begin{split}|J_{\ell}^{N}(\theta)-J_{\ell}(\theta)|& =\left|\int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta}(x),y)\, \big{(}d\mu_{0}^{N}-d\mu_{0}\big{)}(x,y)\right|=\left|\int_{\mathbb{R}^{2d}} \ell(x,y)\,\big{(}d\mu_{T}^{N}-d\mu_{T}\big{)}(x,y)\right|\\ &\leq\bar{L}_{1}\bar{L}_{2}W_{1}(\mu_{0}^{N},\mu_{0}),\end{split} \tag{3.36}\]

where we have used (3.2) and Proposition 3.2 in the second identity, and we have indicated with \(\bar{L}_{1}\) the Lipschitz constant of \(\ell\) on \(B_{\bar{R}}(0)\), while \(\bar{L}_{2}\) descends from the continuous dependence of solutions of (3.3) on the initial datum (see Proposition 3.3). We insist on the fact that both \(\bar{L}_{1},\bar{L}_{2}\) depend on \(\rho\), i.e., the upper bound on the \(L^{2}\)-norm of the controls.

We now address the uniform converge of \(\nabla_{\theta}J_{\ell}^{N}\) to \(\nabla_{\theta}J_{\ell}\) on bounded sets of \(L^{2}\). As before, let us consider an admissible control \(\theta\) such that \(\|\theta\|_{L^{2}}\leq\rho\). Hence, using the representation provided in (3.9), for a.e. \(t\in[0,T]\) we have:

\[\begin{split}\big{|}\nabla_{\theta}J_{\ell}^{N}(\theta)[t]& -\nabla_{\theta}J_{\ell}(\theta)[t]\big{|}=\\ &\left|\int_{\mathbb{R}^{2d}}\nabla_{\theta}\mathcal{F}^{\top} \big{(}t,\Phi_{(0,t)}^{\theta_{0}}(x),\theta_{0}(t)\big{)}\cdot\mathcal{R}_{(t,T)}^{\theta_{0}}(x)^{\top}\cdot\nabla_{x}\ell^{\top}\big{(}\Phi_{(0,T)}^{\theta_ {0}}(x),y\big{)}\,\big{(}d\mu_{0}^{N}-d\mu_{0}\big{)}(x,y)\right|,\end{split} \tag{3.37}\]
In order to prove uniform convergence in \(L^{2}\) norm, we have to show that the integrand is Lipschitz-continuous in \((x,y)\) for a.e. \(t\in[0,T]\), where the Lipschitz constant has to be \(L^{2}\)-integrable as a function of the \(t\) variable. First of all, combining Assumption 2\(-(v)\) and Lemma A.2, we can prove that there exists constants \(C_{1},\tilde{L}_{3}>0\) (depending on \(\rho\)) such that

\[\begin{split}|\nabla_{\theta}\mathcal{F}(t,\Phi^{\theta}_{(0,t)} (x),\theta(t))|\leq& C_{1}(1+|\theta(t)|),\\ |\nabla_{\theta}\mathcal{F}(t,\Phi^{\theta}_{(0,t)}(x_{1}), \theta(t))-\nabla_{\theta}\mathcal{F}(t,\Phi^{\theta}_{(0,t)}(x_{2}),\theta( t))|\leq&\bar{L}_{3}\bar{L}_{2}(1+|\theta(t)|)|x_{1}-x_{2}|\end{split} \tag{3.38}\]

for a.e. \(t\in[0,T]\). We recall that the quantity \(\bar{L}_{2}>0\) (that already appeared in (3.36)) represents the Lipschitz constant of the flow \(\Phi_{(0,t)}\) with respect to the initial datum. Moreover, from Proposition A.7, it descends that

\[\begin{split}|\mathcal{R}^{\theta}_{(t,T)}(x)|\leq& C_{2},\\ |\mathcal{R}^{\theta}_{(t,T)}(x_{1})-\mathcal{R}^{\theta}_{(t,T)}( x_{2})|\leq&\bar{L}_{4}|x_{1}-x_{2}|\end{split} \tag{3.39}\]

for every \(t\in[0,T]\), where the constants \(C_{2},\bar{L}_{4}\) both depend on \(\rho\). Finally, owing to Assumption 3 and Proposition 2.1, we deduce

\[\begin{split}|\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}(x),y)|\leq& C_{3},\\ |\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}(x_{1}),y_{1})-\nabla_{x} \ell(\Phi^{\theta}_{(0,T)}(x_{2}),y_{2})|\leq&\bar{L}_{5}(\bar{L} _{2}|x_{1}-x_{2}|+|y_{1}-y_{2}|)\end{split} \tag{3.40}\]

for every \(x,y\in B_{R}(0)\), where the constants \(C_{3},\bar{L}_{2}\) and the Lipschitz constant \(\bar{L}_{5}\) of \(\nabla_{x}\ell\) depend, once again, on \(\rho\). Combining (3.38), (3.39) and (3.40), we obtain that there exists a constant \(\tilde{L}_{\rho}>0\) such that

\[|\nabla_{\theta}J^{N}[t]-\nabla_{\theta}J[t]|\leq\tilde{L}_{\rho}(1+|\theta( t)|)W_{1}(\mu_{0}^{N},\mu_{0}),\]

for a.e. \(t\in[0,T]\). Observing that the right-hand side is \(L^{2}\)-integrable in \(t\), the previous inequality yields

\[\|\nabla_{\theta}J^{N}-\nabla_{\theta}J\|_{L^{2}}\leq\tilde{L}_{\rho}(1+\rho) W_{1}(\mu_{0}^{N},\mu_{0}),\]

and this concludes the proof. 

In the next result we provide an estimate of the _generalization error_ in terms of the distance \(W_{1}(\mu_{0}^{N},\mu_{0})\). In this case, the important assumption is that there exists a sequence \((\theta^{\star,N})_{N\geq 1}\) of local minimizers of the functionals \((J^{N})_{N\geq 1}\), and that it is contained in a region where \((J^{N})_{N\geq 1}\) are uniformly strongly convex.

**Theorem 3.10**.: _Under the same notations and hypotheses as in Proposition 3.9, let us further assume that the functional \(J\) admits a local minimizer \(\theta^{\star}\) and, similarly, that, for every \(N\geq 1\), \(\theta^{\star,N}\) is a local minimizer for \(J^{N}\). Moreover, we require that there exists a radius \(\rho>0\) such that, for every \(N\geq\bar{N}\), \(\theta^{\star,N}\in B_{\rho}(\theta^{\star})\) and the functional \(J^{N}\) is \(\eta-\)strongly convex in \(B_{\rho}(\theta^{\star})\), with \(\eta>0\). Then, there exists a constant \(C>0\) such that, for every \(N\geq\bar{N}\), we have_

\[\bigg{|}\int_{\mathbb{R}^{2d}}\ell(x,y)\,d\mu_{T}^{\theta^{\star,N}}(x,y)- \int_{\mathbb{R}^{2d}}\ell(x,y)\,d\mu_{T}^{\theta^{\star}}(x,y)\bigg{|}\leq C \left(W_{1}(\mu_{0}^{N},\mu_{0})+\frac{1}{\sqrt{\eta}}\sqrt{W_{1}(\mu_{0}^{N}, \mu_{0})}\right). \tag{3.41}\]

Proof.: According to our assumptions, the control \(\theta^{\star,N}\in B_{\rho}(\theta^{\star})\) is a local minimizer for \(J^{N}\), and, being \(J^{N}\) strongly convex on \(B_{\rho}(\theta^{\star})\) for \(N\geq\bar{N}\), we deduce that \(\{\theta^{\star,N}\}=\arg\min_{B_{\rho}(\theta^{\star})}J^{N}\). Furthermore, from the \(\eta\)-strong convexity of \(J^{N}\), it follows that for every \(\theta_{1},\theta_{2}\in B_{\rho}(\theta^{\star})\), it holds

\[\langle\nabla_{\theta}J^{N}(\theta_{1})-\nabla_{\theta}J^{N}(\theta_{2}),\, \theta_{1}-\theta_{2}\rangle\geq\eta||\theta_{1}-\theta_{2}||_{L^{2}}^{2}.\]

According to Proposition 3.9, we can pass to the limit in the latter and deduce that

\[\langle\nabla_{\theta}J(\theta_{1})-\nabla_{\theta}J(\theta_{2}),\,\theta_{1}- \theta_{2}\rangle\geq\eta||\theta_{1}-\theta_{2}||_{L^{2}}^{2}\]

for every \(\theta_{1},\theta_{2}\in B_{\rho}(\theta^{\star})\). Hence, \(J\) is \(\eta\)-strongly convex in \(B_{\rho}(\theta^{\star})\) as well, and that \(\{\theta^{\star}\}=\arg\min_{B_{\rho}(\theta^{\star})}J\). Therefore, from the \(\eta\)-strong convexity of \(J^{N}\) and \(J\), we obtain

\[J^{N}(\theta^{\star})-J^{N}(\theta^{\star,N})\geq\frac{\eta}{2}|| \theta^{\star,N}-\theta^{\star}||_{L^{2}}^{2}\] \[J(\theta^{\star,N})-J^{N}(\theta^{\star})\geq\frac{\eta}{2}|| \theta^{\star,N}-\theta^{\star}||_{L^{2}}^{2}.\]

Summing the last two inequalities, we deduce that

\[\eta||\theta^{\star,N}-\theta^{\star}||_{L^{2}}^{2}\leq\left(J^{N}(\theta^{\star})- J(\theta^{\star})\right)+\left(J^{N}(\theta^{\star,N})-J(\theta^{\star,N})\right) \leq 2C_{1}W_{1}(\mu_{0}^{N},\mu_{0}), \tag{3.42}\]
where the second inequality follows from the local uniform convergence of Proposition 3.9. We are now in position to derive a bound on the generalization error:

\[\begin{split}\left|\int_{\mathbb{R}^{2d}}\ell(x,y)\left(d\mu_{T}^{ \theta^{*,N}}-d\mu_{T}^{\theta^{*}}\right)(x,y)\right|&=\left| \int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta^{*,N}}(x),y)\,d\mu_{0}^{N}(x,y) -\int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta^{*}}(x),y)\,d\mu_{0}(x,y) \right|\\ &\leq\int_{\mathbb{R}^{2d}}\left|\ell(\Phi_{(0,T)}^{\theta^{*,N}}( x),y)-\ell(\Phi_{(0,T)}^{\theta^{*}}(x),y)\right|\,d\mu_{0}^{N}(x,y)\\ &\quad+\left|\int_{\mathbb{R}^{2d}}\ell(\Phi_{(0,T)}^{\theta^{*}} (x),y)\big{(}d\mu_{0}^{N}(x,y)-d\mu_{0}(x,y)\big{)}\right|\\ &\leq\bar{L}\sup_{x\in\mathrm{supp}(\mu_{0}^{N})}\left|\Phi_{(0, T)}^{\theta^{*,N}}(x)-\Phi_{(0,T)}^{\theta^{*}}(x)\right|+\bar{L}_{R}W_{1}( \mu_{0}^{N},\mu_{0}),\end{split} \tag{3.43}\]

where \(\bar{L}\) and \(\bar{L}_{R}\) are constants coming from Assumption 3 and Proposition 2.1. Then, we combine Proposition 2.1 with the estimate in (3.42), in order to obtain

\[\sup_{x\in\mathrm{supp}(\mu_{0}^{N})}\left|\Phi_{(0,T)}^{\theta^{*,N}}(x)- \Phi_{(0,T)}^{\theta^{*}}(x)\right|\leq C_{2}\|\theta^{*,N}-\theta^{*}\|_{L^{2 }}\leq C_{2}\sqrt{\frac{2C_{1}}{\eta}W_{1}(\mu_{0}^{N},\mu_{0})}.\]

Finally, from the last inequality and (3.43), we deduce (3.41). 

**Remark 3.4**.: _Since the functional \(J:L^{2}([0,T],\mathbb{R}^{m})\to\mathbb{R}\) defined in (3.6) is continuous (and, in particular, lower semi-continuous) with respect to the strong topology of \(L^{2}\), the locally uniform convergence of the functionals \(J^{N}\) to \(J\) (see Proposition 3.9) implies that \(J^{N}\) is \(\Gamma\)-converging to \(J\) with respect to the strong topology of \(L^{2}\). However, this fact is of little use, since the functionals \(J,J^{N}\) are not strongly coercive. On the other hand, if we equip \(L^{2}\) with the weak topology, in general the functional \(J\) is not lower semi-continuous. In our framework, the only circumstance where one can hope for \(\Gamma\)-convergence with respect to the weak topology corresponds to the highly-regularized scenario, i.e., when the parameter \(\lambda>0\) is sufficiently large. Therefore, in the situations of practical interest when \(\lambda\) is small, we cannot rely on this tool, and the crucial aspect is that the dynamics (2.1) is nonlinear with respect to the control variable. Indeed, in the case of affine-control systems considered in [44], it is possible to establish \(\Gamma\)-convergence results in the \(L^{2}\)-weak topology (see [43] for an application to diffeomorphisms approximation). Finally, we report that in [46], in order to obtain the \(L^{2}\)-strong equi-coercivity of the functionals, the authors introduced in the cost the \(H^{1}\)-seminorm of the controls._

### Convex regime and previous result

In order to conclude our mean-field analysis, we now compare our results with the ones obtained in the similar framework of [7], where the regularization parameter \(\lambda\) was assumed to be _sufficiently large_, leading to a convex regime in the sublevel sets (see Remark 3.1). We recall below the main results presented in [7].

**Theorem 3.11**.: _Given \(T,R,R_{T}>0\), and an initial datum \(\mu_{0}\in\mathcal{P}_{c}(\mathbb{R}^{2d})\) with \(\mathrm{supp}(\mu_{0})\subset B_{R}(0)\), let us consider a terminal condition \(\psi_{T}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\) such that \(\mathrm{supp}(\psi_{T})\subset B_{R_{T}}(0)\) and \(\psi_{T}(x,y)=\ell(x,y)\ \forall x,y\in B_{R}(0)\). Let \(\mathcal{F}\) satisfy [7, Assumptions 1-2] and \(\ell\in C^{2}(\mathbb{R}^{d}\times\mathbb{R}^{d},\mathbb{R})\). Assume further that \(\lambda>0\) is large enough. Then, there exists a triple \((\mu^{*},\theta^{*},\psi^{*})\in\mathcal{C}([0,T],\mathcal{P}_{c}(\mathbb{R}^ {2d}))\times Lip([0,T],\mathbb{R}^{m})\times\mathcal{C}^{1}([0,T],\mathcal{C} _{c}^{2}(\mathbb{R}^{2d}))\) solution of_

\[\begin{cases}\partial_{t}\mu_{t}^{*}(x,y)+\nabla_{x}\cdot(\mathcal{F}(t,x, \theta^{*}(t))\mu_{t}^{*}(x,y))=0,&\mu_{t}^{*}|_{t=0}(x,y)=\mu_{0}(x,y),\\ \partial_{t}\psi_{t}^{*}(x,y)+\nabla_{x}\psi_{t}^{*}(x,y)\cdot\mathcal{F}(t,x,\theta^{*}(t))=0,&\psi_{t}^{*}|_{t=T}(x,y)=\ell(x,y),\\ \theta^{*^{\top}}(t)=-\frac{1}{2\lambda}\int_{\mathbb{R}^{2d}}\nabla_{x}\psi_{t }^{*}(x,y)\cdot\nabla_{\theta}\mathcal{F}(t,x,\theta^{*}(t))\,d\mu_{t}^{*}(x,y ),\end{cases} \tag{3.44}\]

_where \(\psi^{*}\in\mathcal{C}^{1}([0,T],\mathcal{C}_{c}^{2}(\mathbb{R}^{2d}))\) is in characteristic form. Moreover, the control solution \(\theta^{*}\) is unique in a ball \(\Gamma_{C}\subset L^{2}([0,T],\mathbb{R}^{m})\) and continuously dependent on the initial datum \(\mu_{0}\)._

We observe that the condition on \(\lambda>0\) to be large enough is crucial to obtain local convexity of the cost functional and, consequently, existence and uniqueness of the solution. However, in the present paper we have not done assumptions on the magnitude of \(\lambda\), hence, as it was already noticed in Remark 3.1, we might end up in a non-convex regime. Nevertheless, in Proposition 3.12 we show that, in the case of \(\lambda\) sufficiently large, the previous approach and the current one are "equivalent".

**Proposition 3.12**.: _Under the same hypotheses as in Theorem 3.11, let \(J:L^{2}([0,T])\to\mathbb{R}\) be the functional defined in (3.6). Then, \(\theta^{*}\) satisfies (3.44) if and only if it is a critical point for \(J\)._
Proof.: According to Lemma (3.4), the gradient of the functional \(J\) at \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) is defined for a.e. \(t\in[0,T]\) as

\[\nabla_{\theta}J(\theta)[t]=\int_{\mathbb{R}^{2d}}\nabla_{\theta}\mathcal{F}^{ \top}\big{(}t,\Phi^{\theta}_{(0,t)}(x),\theta(t)\big{)}\cdot\mathcal{R}^{ \theta}_{(t,T)}(x)^{\top}\cdot\nabla_{x}\ell^{\top}\big{(}\Phi^{\theta}_{(0,T )}(x),y\big{)}\,d\mu_{0}(x,y)+2\lambda\theta(t).\]

Hence, if we set the previous expression equal to zero, we obtain the characterization of the critical point

\[\theta(t)=-\frac{1}{2\lambda}\int_{\mathbb{R}^{2d}}\nabla_{\theta}\mathcal{F}^ {\top}\big{(}t,\Phi^{\theta}_{(0,t)}(x),\theta(t)\big{)}\cdot\mathcal{R}^{ \theta}_{(t,T)}(x)^{\top}\cdot\nabla_{x}\ell^{\top}\big{(}\Phi^{\theta}_{(0,T )}(x),y\big{)}\,d\mu_{0}(x,y) \tag{3.45}\]

for a.e. \(t\in[0,T]\). On the other hand, according to Theorem 3.11, the optimal \(\theta\) satisfies for a.e. \(t\in[0,T]\) the following

\[\theta(t) =-\frac{1}{2\lambda}\int_{\mathbb{R}^{2d}}\big{(}\nabla_{x}\psi_ {t}(x,y)\cdot\nabla_{\theta}\mathcal{F}(t,x,\theta(t))\big{)}^{\top}\,d\mu_{t} (x,y) \tag{3.46}\] \[=-\frac{1}{2\lambda}\int_{\mathbb{R}^{2d}}\nabla_{\theta} \mathcal{F}^{\top}(t,\Phi^{\theta}_{(0,t)}(x),\theta(t))\cdot\nabla_{x}\psi_{ t}^{\top}(\Phi^{\theta}_{(0,t)}(x),y)\,d\mu_{0}(x,y).\]

Hence, to conclude that \(\nabla_{\theta}J=0\) is equivalent to condition stated in Theorem 3.11, we are left to show that

\[\mathcal{R}^{\theta}_{(t,T)}(x)^{\top}\cdot\nabla_{x}\ell^{\top}\big{(}\Phi^ {\theta}_{(0,T)}(x),y\big{)}=\nabla_{x}\psi_{t}^{\top}(\Phi^{\theta}_{(0,t)}( x),y), \tag{3.47}\]

where the operator \(\mathcal{R}^{\theta}_{(t,T)}(x)\) is defined as the solution of (3.10). First of all, we recall that \((t,x,y)\mapsto\psi(t,\Phi^{\theta}_{(0,t)}(x),y)\) is defined as the characteristic solution of the second equation in (3.44) and, as such, it satisfies

\[\psi_{t}(x,y)=\ell(\Phi^{\theta}_{(t,T)}(x),y),\]

for every \(t\in[0,T]\) and for every \(x,y\in B_{R_{T}}(0)\). By taking the gradient with respect to \(x\), we obtain that

\[\nabla_{x}\psi_{t}(x,y)=\nabla_{x}\ell(\Phi^{\theta}_{(t,T)}(x),y))\cdot \nabla_{x}\Phi^{\theta}_{(t,T)}\big{|}_{x},\]

for all \(x,y\in B_{R_{T}}(0)\). Hence, using (3.11), we deduce that

\[\nabla_{x}\psi_{t}(\Phi^{\theta}_{(0,t)}(x),y)=\nabla_{x}\ell(\Phi^{\theta}_{ (0,t)}\circ\Phi^{\theta}_{(0,t)}(x),y))\cdot\nabla_{x}\Phi^{\theta}_{(t,T)} \big{|}_{\Phi^{\theta}_{(0,t)}(x)}=\nabla_{x}\ell(\Phi^{\theta}_{(0,T)}(x),y) \cdot\mathcal{R}^{\theta}_{(t,T)}(x)\]

which proves (3.47). 

## 4 Algorithm

In this section, we present our training procedure, which is derived from the necessary optimality conditions related to the minimizing movement scheme (see (3.29)). Since the mean-field optimal control problem as presented in (3.6) is numerically intractable (especially in high-dimension), in the practice we always consider the functional corresponding to the finite particles approximation (see (3.32)). For its resolution, we employ an algorithm belonging to the family of shooting methods, which consists in the forward evolution of the trajectories, in the backward evolution of the adjoint variables, and in the update of the controls. Variants of this method have already been employed in different works, e.g. [33, 26, 6, 7, 14], with the name of _method of successive approximations_, and they have been proven to be an alternative way of performing the training of NeurODEs for a range of tasks, including high-dimensional problems.

In our case, we start with a random guess for the control parameter \(\theta_{0}\in L^{2}([0,T],\mathbb{R}^{m})\). Subsequently, we solve the necessary optimality conditions specified in equation (3.29) for a suitable \(\tau>0\) to obtain an updated control parameter \(\theta_{1}\). More precisely, since the last identity in (3.29) has the form \(\theta_{1}=\Lambda^{\tau}_{\theta_{0}}(\theta_{1})\), the computation of \(\theta_{1}\) is performed via fixed-point iterations of the mapping \(\Lambda^{\tau}_{\theta_{0}}\), which is defined as in (3.30). In this regard, we recall that \(\Lambda^{\tau}_{\theta_{0}}\) is a contraction if \(\tau\) is small enough. The scheme that we implemented is presented is Algorithm 1.

**Remark 4.1**.: _It is interesting to observe that, in the highly-regularized regime considered in [7], the authors managed to obtain a contractive map directly from the necessary conditions for optimality, and they did not need to consider the minimizing movements scheme. This is rather natural since, when the parameter \(\lambda>0\) that tunes the \(L^{2}\)-penalization is large enough, the functional associated to the optimal control problem is strongly convex in the sublevel set corresponding to the control \(\theta\equiv 0\), as discussed in Remark 3.1. However, as reported in [7], determining the appropriate value for \(\lambda\) in each application can be challenging. On the other hand, from the practitioners' perspective, dealing with high regularization is not always desirable, since the machine learning task that the system should learn is encoded in the final-time cost. The authors highlighted the complexity involved in selecting a regularization parameter that is large enough to achieve contractivity, while ensuring that the resulting controls are not excessively small (due to high regularization) and of little use._
These considerations motivated us to consider a scenario where the regularization parameter does not need to be set sufficiently large. From a numerical perspective, the parameter \(\tau\) in Equation (3.29) (coming from the minimizing movement scheme) plays the role of the _learning rate_, and it provides the lacking amount of convexity, addressing the stability issues related to the resolution of optimal control problems in non-convex regime. These kinds of instabilities were already known in the Soviet literature of numerical optimal control (see the review paper [14]), and various solutions have been proposed to address them. For example, in [40] the authors proposed an iterative method based on the Maximum Principle and on an augmented Hamiltonian, with an approach that is somehow reminiscent of minimizing movements. More recently, in the framework of NeurODEs, in [33] it was proposed another stabilization strategy, which is different from ours since it enforces similarity between the evolution of state and co-state variables after the control update. Implicitly, the approach of [33] leads to a penalization of significant changes in the controls. On the other hand, in our approach this penalization is more explicit, and it is enforced via the memory term of the minimizing movement scheme. To the best of our knowledge, this is the first instance where a regularization based on the minimizing movement scheme is employed for training NeurODEs.

**Remark 4.2**.: _Although we formulate and analyze theoretically our problem within the mean-field framework, it is not advantageous to numerically solve the forward equation as a partial differential equation. In [7], various numerical methods for solving PDEs were employed and compared. However, these methods encounter limitations when applied to high-dimensional data, which is often the case in Machine Learning scenarios. Therefore, in this study, we employ a particle method to solve both the forward partial differential equation and the backward dynamics. This particle-based approach involves reformulating the PDE as a system of ordinary differential equations in which particles represent mathematical collocation points that discretize the continuous fields. By employing this particle method, we address the challenges associated with high-dimensional data, enabling efficient numerical solutions for the forward and backward dynamics._

To conclude this section, we briefly present the forward and the backward systems that are solved during the execution of the method. For the sake of simplicity, we will focus on the case of an encoder. The objective is to minimize the following function:

\[J(\theta)=\frac{1}{N}\sum_{i=1}^{N}\ell(X_{\mathcal{A}_{r}}^{i}(T),Y^{i}(0))+ \frac{\lambda}{2}\left\|\theta\right\|_{2}^{2}, \tag{4.1}\]

where \(\mathcal{A}_{r}\) denotes the active indices in the bottleneck, i.e. at \(t_{r}=T\), of the state-vector \(X^{i}(T)\). The latter denotes the encoded output at time \(T\) for the \(i\)-th particle, while \(Y^{i}(0)\) represents the corresponding target at time \(0\) (which we recall is the same at time \(T\), being \(\dot{Y}^{i}\equiv 0\) for every \(i=1,\ldots,N\)). For each \(i\)-th particle and every \(t\) such that
\(t_{j}\leq t\leq t_{j+1}\), the forward dynamics can be described as follows:

\[\begin{cases}\dot{X}_{Z_{j}}^{i}(t)=0,\\ \dot{X}_{\mathcal{A}_{j}}^{i}(t)=\mathcal{G}_{j}\big{(}t,X_{\mathcal{A}_{j}}^{i} (t),\theta(t)\big{)},\end{cases} \tag{4.2}\]

subject to the initial condition \(X^{i}(0)=X_{\mathcal{A}_{0}}^{i}(0)=X_{0}^{i}\in\mathbb{R}^{d}\). In the same interval \(t_{j}\leq t\leq t_{j+1}\), the backward dynamics reads

\[\begin{cases}\dot{P}_{\mathcal{I}_{j}}^{i}(t)=0,\\ \dot{P}_{\mathcal{A}_{j}}^{i}(t)=-P_{\mathcal{A}_{j}}^{i}(t)\cdot\nabla_{x_{ \mathcal{A}_{j}}}\mathcal{G}_{j}\big{(}t,X_{\mathcal{A}_{j}}^{i}(t),\theta(t) \big{)},\end{cases} \tag{4.3}\]

where the final co-state is

\[P^{i}(T)=\begin{cases}-\partial_{t}\ell(X_{\mathcal{A}_{r}}^{i}(T),Y^{i}(0)), &\text{if }k\in\mathcal{A}_{r},\\ 0,&\text{if }k\notin\mathcal{A}_{r}.\end{cases}\]

We notice that, for \(t_{j}\leq t\leq t_{j+1}\) and every \(i\in\{0,\dots,N\}\), we have

\[\mathcal{F}(t,X^{i}(t),\theta(t))=\mathcal{F}(t,(X_{\mathcal{A}_{r}}^{i},X_{ \mathcal{I}_{j}}^{i})(t),\theta(t))=\begin{pmatrix}\mathcal{G}_{j}(t,X_{ \mathcal{A}_{j}}^{i}(t),\theta(t))\\ 0\end{pmatrix}, \tag{4.4}\]

and, consequently, we deduce that

\[\nabla_{x}\mathcal{F}(t,X^{i}(t),\theta(t))=\begin{pmatrix}\nabla_{x_{ \mathcal{A}_{j}}}\mathcal{G}_{j}(t,X_{\mathcal{A}_{j}}^{i}(t),\theta(t))&0\\ 0&0\end{pmatrix}, \tag{4.5}\]

where the null blocks are due to the fact that, for \(t_{j}\leq t\leq t_{j+1}\), \(\nabla_{x}\mathcal{F}_{k}(t,x,\theta)=0\) if \(k\in\mathcal{I}_{j}\), and \(\nabla_{x_{\mathcal{I}_{j}}}\mathcal{G}_{j}(t,x,\theta)=0\). In the case of an Autoencoder, the structure of the forward and backward dynamics is analogous.

**Remark 4.3**.: _From the calculations reported above it is evident that the matrices and the vectors involved in our forward and backward dynamics are quite sparse (see (4.5) and (4.4)), and that the state and co-state variables contain components that are constant in many sub-intervals (see (4.2) and (4.3)). Hence, in the practical implementation, especially when dealing with an Autoencoder, we do not actually need to double the original state variables and to introduce the shadow ones, but we can simply overwrite those values and, in this way, we obtain a more memory-efficient code. A similar argument holds as well for the co-state variables. Moreover, we expect the control variable \(\theta\) to have several null components during the evolution. This relates to Remark 2.1 and descends from the fact that, even though in our model \(\theta\in\mathbb{R}^{m}\) for every \(t\in[0,T]\), in the internal sub-intervals \([t_{j},t_{j+1}]\) only few of its components are influencing the dynamics. Hence, owing to the \(L^{2}\)-squared regularization on \(\theta\), it results that, if in an interval \([t_{j},t_{j+1}]\) a certain component of \(\theta\) is not affecting the velocity, then it is convenient to keep it null._

## 5 Numerical Experiments

In this section, we present a series of numerical examples to illustrate the practical application of our approach. We consider datasets of varying dimensions, ranging from low-dimensional data to a more typical Machine Learning dataset such as MNIST. Additionally, we provide justifications and insights into some of the choices made in our theoretical analysis. For instance, we examine the process of choosing the components to be deactivated during the modeling phase, and we investigate whether this hand-picked selection can lead to any issues or incorrect results. In this regard, in our first experiment concerning a classification task, we demonstrate that this a priori choice does not pose any problem, as the network effectively learns to separate the dataset into two classes before accurately classifying them. Furthermore, as we already pointed out, we have extended some of the assumptions from [7] to accommodate the use of a smooth approximation of the ReLU function. This extension is not merely a theoretical exercise, since in our second numerical example we show how valuable it is to leverage unbounded activation functions. While both of these examples involve low-dimensional data and may not be representative of typical tasks for an Autoencoder architecture, we address this limitation in our third experiment by performing a reconstruction task on the MNIST dataset. Lastly, we present noteworthy results obtained from analyzing the performance on MNIST, highlighting specific behaviors that warrant further investigation in future research.

The layers of the networks that we employ in all our experiments have the form:

\[\mathbb{R}^{d}\ni X=\big{(}X_{\mathcal{A}_{j}},X_{\mathcal{I}_{j}}\big{)}^{ \top}\mapsto\phi_{n}^{W,b}(X)=\big{(}X_{\mathcal{A}_{j}},X_{\mathcal{I}_{j}} \big{)}^{\top}+h\Big{(}\sigma\left(W_{\mathcal{A}_{j}}\cdot X_{\mathcal{A}_{j }}+b_{\mathcal{A}_{j}}\right),0\Big{)}^{\top},\]

where \(\mathcal{A}_{j},\mathcal{I}_{j}\) are, respectively, the sets of active and inactive components at the layer \(n\), \(b_{\mathcal{A}_{j}}\) are the components of \(b\in\mathbb{R}^{d}\) belonging to \(\mathcal{A}_{j}\), while \(W_{\mathcal{A}_{j}}\) is the square sub-matrix of \(W\in\mathbb{R}^{d\times d}\) corresponding to the active components. Finally, the activation function \(\sigma\) will be specified case by case.
### Bidimensional Classification

In our initial experiment, we concentrate on a bidimensional classification task that has been extensively described in [7]. Although this task deviates from the typical application of Autoencoders, where the objective is data reconstruction instead of classification, we believe it gives valuable insights on how our model works. The objective is to classify particles sampled from a standard Gaussian distribution in \(\mathbb{R}^{2}\) based on the sign of their first component. Given an initial data point \(x_{0}\in\mathbb{R}^{2}\), denoted by \(x_{0}[i]\) with \(i=1,2\) representing its \(i-\)th component, we assign a positive label \(+1\) to it if \(x_{0}[1]>0\), and a negative label \(-1\) otherwise. To incorporate the labels into the Autoencoder framework, we augment the labels to obtain a positive label \([1,0]\) and a negative one \([-1,0]\). In such a way, we obtain target vectors in \(\mathbb{R}^{2}\), i.e., with the same dimension as the input data-points in the first layer.

The considered architecture is an Autoencoder comprising twenty-one layers, corresponding to \(T=2\) and \(dt=0.05\). The first seven layers maintain a constant active dimension equal to \(2\), followed by seven layers of active dimension \(1\). Finally, the last seven layers, representing the prototype of a decoder, have again constant active dimension \(2\), restoring the initial one.

A sketch of the architecture is presented on the right side of Figure 4. We underline that we make use of the observation presented in Remark 4.3 to construct the implemented network, and we report that we employ the hyperbolic tangent as activation function.

The next step is to determine which components to deactivate, i.e., we have to choose the sets \(\mathcal{I}_{j}\) for \(j=1,\ldots,2r\): the natural choice is to deactivate the second component, since the information which the classification is based on is contained in the first component (the sign) of the input data-points. Since we use the memory-saving regime of Remark 4.3, we observe that, in the decoder, the particles are "projected" onto the \(x\)-axis, as their second component is deactivated and set equal to \(0\). Then, in the decoding phase, both the components have again the possibility of evolving. This particular case is illustrated on the left side of Figure 4.

Now, let us consider a scenario where the network architecture remains the same, but instead of deactivating the second component, we turn-off the first component. This has the effect of "projecting" the particles onto the \(y\)-axis in the encoding phase.

The results are presented in Figure 5, where an interesting effect emerges. In the initial phase (left), where the particles can evolve in the whole space \(\mathbb{R}^{2}\), the network is capable of rearranging the particles in order to separate them. More precisely, in this part, the relevant information for the classification (i.e., the sign of the first

Figure 4: Left: Classification task performed when the turned off component is the natural one.

Right: sketch of the AutoencODE architecture considered.

Figure 5: Left: Initial phase, i.e., separation of the data along the \(y\)-axis. Center: Encoding phase, i.e., only the second component is active. Right: Decoding phase and classification result after the “unnatural turn off”. Notice that, for a nice clustering of the classified data, we have increased the number of layers from \(20\) to \(40\). However, we report that the network accomplishes the task even if we use the same structure as in Figure 4.
component), is transferred to the second component, that will not be deactivated. Therefore, once the data-points are projected onto the \(y\)-axis in the bottleneck (middle), two distinct clusters are already formed, corresponding to the two classes of particles. Finally, when the full dimension is restored, the remaining task consists in moving these clusters towards the respective labels, as demonstrated in the plot on the right of Figure 5.

This numerical evidence confirms that our a priori choice (even when it is very unnatural) of the components to be deactivated does not affect the network's ability to learn and classify the data. Finally, while studying this low-dimensional numerical example, we test one of the assumptions that we made in the theoretical setting. In particular, we want to check if it is reasonable to assume that the cost landscape is convex around local minima, as assumed in Theorem 3.10. In Table 1, we report the smallest and highest eigenvalues of the Hessian matrix of the loss function recorded during the training process, i.e., starting from a random initial guess, until the convergence to an optimal solution.

### Parabola Reconstruction

In our second numerical experiment, we focus on the task of reconstructing a two-dimensional parabola. To achieve this, we sample points from the parabolic curve and we use them as the initial data for our network. The network architecture consists of a first block of seven layers with active dimension 2, followed by seven additional layers with active dimension 1. Together, these two blocks represent the encoding phase in which the set of active components are \(\mathcal{A}_{j}=\{0\}\) for \(j=7,\ldots,14\). Similarly as in the previous example, the points at the 7-th layer are "projected" onto the \(x-\)axis, and for the six subsequent layers they are constrained to stay in this subspace. After the 14-th layer, the original active dimension is restored, and the particles can move in the whole space \(\mathbb{R}^{2}\), aiming at reaching their original positions. Despite the low dimensionality of this task, it provides an interesting application that allows us to observe the distinct phases of our mode, which are presented in Figure 6.

Notably, in the initial seven layers, the particles show quite tiny movements (top left of Figure 6). This is since the relevant information to reconstruct the position is encoded in the first component, which is kept active in the bottleneck. On the other hand, if in the encoder we chose to deactivate the first component instead of the second one, we would expect that the points need to move considerably before the projection takes place, as it was the case in the previous classification task. During the second phase (top right of Figure 6), the particles separate along the \(x\)-axis, preparing for the final decoding phase, which proves to be the most challenging to learn (depicted in the

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline Epochs & 0 & 80 & 160 & 240 & 320 & 400 & 480 & 560 & 640 & 720 \\ \hline Min Eingenvalue & -1.72e-2 & -1.19e-2 & -1.09e-2 & -8.10e-3 & -3.44e-3 & -6.13e-3 & 6.80e-4 & 7.11e-4 & 7.25e-4 & 7.33e-4 \\ Max Eigenvalue & 3.78e-2 & 2.84e-1 & 7.30e-1 & 9.34e-1 & 1.11 & 1.18 & 1.22 & 1.25 & 1.26 & 1.27 \\ \hline \end{tabular}
\end{table}
Table 1: Minimum and maximum eigenvalues of the Hessian matrix across epochs.

Figure 6: Top Left: Initial phase. Top Rights: Encoding phase. Bottom Left: Decoding phase. Bottom Right: network’s reconstruction with alternative architecture.
bottom left of Figure 6). Based on our theoretical knowledge and the results from initial experiments, we attempt to improve the performance of the AutoencODE network by modifying its structure. One possible approach is to design the network in a way that allows more time for the particles to evolve during the decoding phase, while reducing the time spent in the initial and bottleneck phases. Indeed, we try to use 40 layers instead of 20, and most of the new ones are allocated in the decoding phase. The result is illustrated in the bottom right of Figure 6, where we observe that changing the network's structure has a significant positive impact on the reconstruction quality, leading to better results. This result is inspired by the heuristic observation that the particles "do not need to move" in the first two phases. On this point, a more theoretical analysis of the network's structure will be further discussed in the next paragraph, where we perform a sanity check, and we relate the need for extra layers to the Lipschitz constant of the trained network.

This experiment highlights an important observation regarding the choice of activation functions. Specifically, it becomes evident that certain bounded activation functions, such as the hyperbolic tangent, are inadequate for moving the particles back to their original positions during the decoding phase. The bounded nature of these activation functions limits their ability to move a sufficiently large range of values, which can lead to the points getting stuck at suboptimal positions and failing to reconstruct the parabolic curve accurately. To overcome this limitation and achieve successful reconstruction, it is necessary to employ unbounded activation functions that allow for a wider range of values, in particular the well-known Leaky Relu function. An advantage of our approach is that our theory permits the use of smooth approximations for well-known activation functions, such as the Leaky ReLU (2.4). Specifically, we employ the following smooth approximation of the Leaky ReLU function:

\[\sigma_{smooth}(x)=\alpha x+\frac{1}{s}\log\big{(}1+e^{sx}\big{)}, \tag{5.1}\]

where \(s\) approaching infinity ensures convergence to the original Leaky ReLU function. While alternative approximations are available, we employed (5.1) in our study. This observation emphasizes the importance of considering the characteristics and properties of activation functions when designing and training neural networks, and it motivates our goal in this work to encompass unbounded activation functions in our working assumptions.

### MNIST Reconstruction

In this experiment, we apply the AutoencODE architecture and our training method to the task of reconstructing images from the MNIST dataset. The MNIST dataset contains 70000 grayscale images of handwritten digits ranging from zero to nine. Each image has a size of \(28\times 28\) pixels and has been normalized. This dataset is commonly used as a benchmark for image classification tasks or for evaluating image recognition and reconstruction algorithms. However, our objective in this experiment is not to compare our reconstruction error with state-of-the-art results, but rather to demonstrate the applicability of our method to high-dimensional data, and to highlight interesting phenomena that we encounter. In general, when performing an autoencoder reconstruction task, the goal is to learn a lower-dimensional representation of the data that captures its essential features. On the other hand, determining the dimension of the lower-dimensional representation, often referred to as the _latent dimension_, requires setting a hyperparameter, i.e., the width of the bottleneck's layers, which might depend on the specific application.

We now discuss the architecture we employed and the choice we made for the latent dimension. Our network consists of twenty-three layers, with the first ten layers serving as encoder, where the dimension of the layers is gradually reduced from the initial value \(d_{0}=784\) to a latent dimension of \(d_{r}=32\). Then, this latent dimension is kept in the bottleneck for three layers, and the last ten layers act as decoder, and, symmetrically to the encoder, it increases the width of the layers from 32 back to \(d_{2r}=784\). Finally, for each layer we employ a smooth version of the Leaky Relu, see (5.1), as activation function. The architecture is visualized in Figure 7, while the achieved reconstruction results are presented in Figure 8. We observe that, once again, we made use of Remark 4.3 for the implementation of the AutoencODE-based model.

Latent dimensionality in the bottleneck:One of the first findings that we observe in our experiments pertains to the latent dimension of the network and to the intrinsic dimension of the dataset. The problem of determining the intrinsic dimension has been object of previous studies such as [47, 16, 17], where it was estimated to be approximately equal to 13 in the case of MNIST dataset. On this interesting topic, we also report the paper [32], where a maximum likelihood estimator was proposed and datasets of images were considered, and the recent contribution [35]. Finally, the model of the _hidden manifold_ has been formulated and studied in [23].

Notably, our network exhibits an interesting characteristic in which, starting from the initial guess of weights and biases initialized at 0, the training process automatically identifies an intrinsic dimensionality of 13. Namely, we observe that the latent vectors of dimension 32 corresponding to each image in the dataset are sparse vectors with 13 non-zero components, forming a consistent support across all latent vectors derived from the original images. To further analyze this phenomenon, we compute the means of all the latent vectors for each digit and we compare them,
as depicted in the left and middle of Figure 9. These mean vectors always have exactly the same support of dimension 13, and, interestingly, we observe that digits that share similar handwritten shapes, such as the numbers 4 and 9 or digits 3 and 5, actually have latent means that are close to each other. Additionally, we explore the generative capabilities of our network by allowing the latent means to evolve through the decoding phase, aiming to generate new images consistent with the mean vector. On the right of Figure 9, we present the output of the network when using a latent vector corresponding to the mean of all latent vectors representing digit 3.

This intriguing behavior of our network warrants further investigation into its ability to detect the intrinsic dimension of the input data, and into the exploration of its generative potential. Previous studies have demonstrated that the ability of neural networks to converge to simpler solutions is significantly influenced by the initial parameter values (see e.g. [15]). Indeed, in our case we have observed that this phenomenon only occurs when initializing the parameters with zeros. Moreover, it is worth mentioning that this behavior does not seem to appear in standard Autoencoders without residual connections.

Sanity check of the network's architecture.An advantage of interpreting neural networks as discrete approximations of dynamical systems is that we can make use of typical results of numerical resolutions of ODEs in

Figure 8: Reconstruction of some numbers achieved by the Autoencoder.

Figure 7: Architecture used for the MNIST reconstruction task. The inactive nodes are marked in green.
order to better analyze our results. Indeed, we notice that, according to well-known results, in order to solve a generic ODEs we need to take as discretization step-size \(dt\) a value smaller than the inverse of the lipschitz constant of the vector field driving the dynamics. We recall that the quantity \(dt\) is related to the number of layers of the network through the relation \(n_{\text{layers}}=\frac{T}{dt}\), where \(T\) is the right-extreme of the evolution interval \([0,T]\).

In our case, we choose _a priori_ the amplitude of \(dt\), we train the network and, once we have computed \(\theta^{*}\), we can compare _a posteriori_ the discretization step-size chosen at the beginning with the quantity \(\Delta=\frac{1}{Lip(\mathcal{F}(t,X,\theta^{*})}\) for each time-node \(t\) and every datum \(x\). In Figure 5, we show the time discretization \(dt\) in orange and in blue the quantity \(\Delta\),

for the case of a wrongly constructed autoencoder (on the left) and the correct one (on the right). From this plots, we can perform a "sanity check" and we can make sure that the number of layers that we chose is sufficient to solve the task. Indeed, in the wrong autoencoder on the left, we see that in the last layer the quantity \(\Delta\) is smaller than \(dt\), and this violates the condition that guarantees the stability of the explicit Euler discretization.

Indeed, the introduction of two symmetric layers to the network (corresponding to the plot on the right of Figure 5) allows the network to satisfy everywhere the relation \(\Delta>dt\). Moreover, we also notice that during the encoding phase the inverse of the Lipschitz constant of \(\mathcal{F}\) is quite high, which means that the vector field does not need to move a lot the points. This suggests that we could get rid of some of the layers in the encoder and only keep the necessary ones, i.e., the ones in the decoder where \(\Delta\) is small and a finer discretization step-size is required. We report that this last observation is consistent with the results recently obtained in [11]. Finally, we also draw attention to the work [45], which shares a similar spirit with our experiments, since the Lipschitz constant of the layers is the main subject of investigation. In their study, the authors employ classical results on the numerical integration of ordinary differential equations in order to understand how to constrain the weights of the network with the aim of designing stable architectures. This approach leads to networks with non-expansive properties, which is highly advantageous for mitigating instabilities in various scenarios, such as testing adversarial examples [25], training generative adversarial networks [3], or solving inverse problems using deep learning.

Entropy across layers.We present our first experiments on the study of the information propagation within the network, where some intriguing results appear. This phenomenon is illustrated in Figure 5, where we examine the entropy across the layers after the network has been trained. We introduce two different measures of entropy, depicted in the two graphs of the figure. In first place, we consider the well-known Shannon entropy, denoted as \(H(E)\), which quantifies the information content of a discrete random variable \(E\), distributed according to a discrete

Figure 10: Left: wrong autoencoder detected with \(\Delta\). Right: correct version of the same autoencoder.

Figure 9: Left: Comparing two similar latent means. Center: again two similar latent means. Right: Output of the decoding of one latent mean.
probability measure \(p:\Omega\to[0,1]\) such that \(p(e)=p(E=e)\). The Shannon entropy is computed as follows:

\[H(E)=\mathbb{E}\big{[}-\log(p(E)\big{]}=\sum_{e\in E}-p(e)\log(p(e))\]

In our context, the random variable of interest is \(E=\sum_{j=1}^{N}\mathds{1}_{|X_{0}^{i}-X_{0}^{j}|\leq\epsilon}\), where \(X_{0}^{i}\) represents a generic image from the MNIST dataset. Additionally, we introduce another measure of entropy, denoted as \(\mathcal{E}\), which quantifies the probability that the dataset can be partitioned into ten clusters corresponding to the ten different digits. This quantity has been introduced in [21] and it is defined as

\[\mathcal{E}=\mathbb{P}\left(X\in\bigcup_{i=1}^{k}B_{\varepsilon}(X_{0}^{i}) \right),\]

where \(\varepsilon>0\) is a small radius, and \(X_{0}^{1},\ldots,X_{0}^{k}\) are samplings from the dataset. Figure 5 suggests the existence of a distinct pattern in the variation of information entropy across the layers, which offers a hint for further investigations.

Let us first focus on the Shannon entropy: as the layers' dimensionality decreases in the encoding phase, there is an expected decrease of entropy, reflecting the compression and reduction of information in the lower-dimensional representation. The bottleneck layer, where the dimension is kept constant, represents a critical point where the entropy reaches a minimum. This indicates that the information content is highly concentrated and compressed in this latent space. Then, during the decoding phase, the Shannon entropy do not revert to its initial value but instead exhibit a slower increase. This behavior suggests that the network retains some of the learned structure and information from the bottleneck layer. Something similar happens for the second measure of entropy: at the beginning, the data is unlikely to be highly clustered, since two distinct images of the same digit may be quite distant one from the other. In the inner layers, this probability increases until it reaches its maximum (rather close to 1) in the bottleneck, where the data can then be fully partitioned into clusters of radius \(\epsilon\). As for the Shannon entropy, the information from the bottleneck layer is retained during the decoding phase, which is why the entropy remains constant for a while and then decreases back in a slower manner.

It is worth noticing that in both cases the entropy does not fully return to its initial level. This might be attributed to the phenomenon of mode collapse, where the network fails to capture the full variability in the input data and instead produces similar outputs for different inputs, hence inducing some sort of _implicit bias_. Mode collapse is often considered undesirable in generative models, as it hinders the ability to generate diverse and realistic samples. However, in the context of understanding data structure and performing clustering, the network's capability to capture the main modes or clusters of the data can be seen as a positive aspect. The network learns to extract salient features and represent the data in a compact and informative manner, enabling tasks such as clustering and classification. Further investigation is needed to explore the relationship between the observed entropy patterns, mode collapse, and the overall performance of the network on different tasks.

## Acknowledgments

The authors would like to thank Prof. Giuseppe Savare for the fruitful discussions during his permanence in Munich. Moreover, the authors are grateful to Dr. Oleh Melnyk for the suggestion on the extension of the dynamical model to the U-net architecture.

This work has been funded by the German Federal Ministry of Education and Research and the Bavarian State Ministry for Science and the Arts. C.C. and M.F. acknowledge also the partial support of the project "Online Firestorms And Resentment Propagation On Social Media: Dynamics, Predictability and Mitigation" of the TUM

Figure 11: Left: Shannon entropy across layers. Right: Our measure of entropy across layers.
Institute for Ethics in Artificial Intelligence and of the DFG Project "Implicit Bias and Low Complexity Networks" within the DFG SPP 2298 "Theoretical Foundations of Deep Learning". A.S. acknowledges the partial support from INdAM-GNAMPA.

## Appendix A Appendix

**Lemma A.1** (Boundedness of trajectories).: _Let us consider the controlled system_

\[\dot{x}=\mathcal{F}(t,x,\theta),\quad x(0)=x_{0},\]

_where \(\mathcal{F}:[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}\) satisfies Assumptions 1, and \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\). Then, for every \(R>0\) and any \(x_{0}\in B_{R}(0)\), we have that \(x(t)\in B_{R}(0)\) for every \(t\in[0,T]\), where \(\bar{R}=(R+L_{R}(1+\|\theta\|_{L^{1}}))e^{L_{R}(1+\|\theta\|_{L^{1}})}\)._

Proof.: According to Assumption 1\(-(ii)\) on \(\mathcal{F}\), the trajectories can be bounded as follows:

\[|x(t)|\leq|x_{0}|+\int_{0}^{t}|\mathcal{F}(s,x(s),\theta(s))|\,ds\leq|x_{0}|+L _{R}\int_{0}^{t}(1+|x(s)|)(1+|\theta(s)|)\,ds\]

for every \(t\in[0,T]\). Using Gronwall's lemma, it follows that

\[|x(t)|\leq\Big{(}|x_{0}|+L_{R}(1+\|\theta\|_{L^{1}})\Big{)}e^{L_{R}(1+\| \theta\|_{L^{1}})}.\]

**Lemma A.2** (Flow's dependency on initial datum).: _For every \(t\in[0,T]\), let us consider the flow mapping \(\Phi^{\theta}_{(0,t)}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) defined in (2.2) and driven by the control \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\). Let us assume that the controlled dynamics \(\mathcal{F}:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{d}\) satisfies Assumption 1. Then, for every \(R>0\), and every \(x_{1},x_{2}\in B_{R}(0)\), it follows that_

\[|\Phi^{\theta}_{(0,t)}(x_{1})-\Phi^{\theta}_{(0,t)}(x_{2})|\leq e^{L_{R}(1+\| \theta\|_{L^{1}})}|x_{1}-x_{2}|,\]

_where \(\bar{R}\) is defined as in Lemma A.1, and \(L_{R}\) is prescribed by Assumption 1-(ii)._

Proof.: Let us denote with \(t\mapsto x_{1}(t),t\mapsto x_{2}(t)\) the solutions of (2.1) driven by \(\theta\) and starting, respectively, from \(x_{1}(0)=x_{1},x_{2}(0)=x_{2}\). Then, for every \(t\in[0,T]\), we have

\[|x_{1}(t)-x_{2}(t)| \leq|x_{1}-x_{2}|+\int_{0}^{t}|\mathcal{F}(s,x_{1}(s),\theta(s))- \mathcal{F}(s,x_{2}(s),\theta(s))|\,ds\] \[\leq|x_{1}-x_{2}|+L_{R}\int_{0}^{t}(1+|\theta(s)|)|x_{1}(s)-x_{2} (s))|\,ds,\]

by using Assumption 1\(-(ii)\). As before, the statement follows from Gronwall's Lemma. 

**Lemma A.3** (Flow's dependency on time).: _Under the same assumptions and notations as in Lemma A.2, for every \(R>0\), for every \(x\in B_{R}(0)\) and for every \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\), we have that_

\[|\Phi^{\theta}_{(0,t_{2})}(x)-\Phi^{\theta}_{(0,t_{1})}(x)|\leq L_{R}(1+\bar{ R})(1+\|\theta\|_{L^{2}})|t_{2}-t_{1}|^{\frac{1}{2}}\]

_for every \(0\leq t_{1}<t_{2}\leq T\), where \(\bar{R}\) is defined as in Lemma A.1, and \(L_{R}\) is prescribed by Assumption 1-(\(ii)\). Moreover, if \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\cap L^{\infty}([0,T],\mathbb{R}^{m})\), then, for every \(0\leq t_{1}<t_{2}\leq T\), it holds:_

\[|\Phi^{\theta}_{(0,t_{2})}(x)-\Phi^{\theta}_{(0,t_{1})}(x)|\leq L_{\bar{R}}(1+ \bar{R})(1+\|\theta\|_{L^{2}})|t_{2}-t_{1}|.\]

Proof.: If we denote by \(t\mapsto x(t)\) the solution of (2.1) driven by the control \(\theta\), then

\[|x(t_{2})-x(t_{1})|\leq\int_{t_{1}}^{t_{2}}|\mathcal{F}(s,x(s),\theta(s)|\,ds \leq\int_{t_{1}}^{t_{2}}L_{\bar{R}}(1+\bar{R})(1+|\theta(s)|)\,ds.\]

The thesis follows by using Cauchy-Schwarz for \(\theta\in L^{2}\), or from basic estimates if \(\theta\in L^{\infty}\). 

**Lemma A.4** (Flow's dependency on controls).: _For every \(t\in[0,T]\), let \(\Phi^{\theta_{1}}_{(0,t)},\Phi^{\theta_{2}}_{(0,t)}:\mathbb{R}^{d}\to\mathbb{R}^{d}\) be the flows defined in (2.2) and driven, respectively, by \(\theta_{1},\theta_{2}\in L^{2}([0,T],\mathbb{R}^{m})\). Let us assume that the controlled dynamics \(\mathcal{F}:[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}\) satisfies Assumption 1. Then, for every \(R>0\) and for every \(x\in B_{R}(0)\), it holds that_

\[|\Phi^{\theta_{1}}_{(0,t)}(x)-\Phi^{\theta_{2}}_{(0,t)}(x)|\leq L_{\bar{R}}(1+ \|\theta_{1}\|_{L^{2}}+\|\theta_{2}\|_{L^{2}})e^{L_{R}(1+\|\theta_{1}\|_{L^{1}}) }\left\|\theta_{1}-\theta_{2}\right\|_{L^{2}},\]

_where \(\bar{R}\) is defined as in Lemma A.1, and \(L_{\bar{R}}\) is prescribed by Assumption 1-(ii)._
Proof.: By using Assumption 1\(-(ii),(iii)\) and the triangle inequality, we obtain that

\[|\Phi^{\theta_{1}}_{(0,t)}(x)-\Phi^{\theta_{2}}_{(0,t)}(x)| \leq\int_{0}^{t}|\mathcal{F}(s,x_{1}(s),\theta_{1}(s))-\mathcal{F}( s,x_{2}(s),\theta_{2}(s))|\,ds\] \[\leq\int_{0}^{t}|\mathcal{F}(s,x_{1}(s),\theta_{1}(s))-\mathcal{F }(s,x_{2}(s),\theta_{1}(s))|\,ds\] \[\quad+\int_{0}^{t}|\mathcal{F}(s,x_{2}(s),\theta_{1}(s))-\mathcal{ F}(s,x_{2}(s),\theta_{2}(s))|\,ds\] \[\leq L_{R}\int_{0}^{t}(1+\theta_{1}(s))|x_{1}(s)-x_{2}(s)|\,ds+L_ {\bar{R}}(1+\left\|\theta_{1}\right\|_{L^{2}}+\left\|\theta_{2}\right\|_{L^{2} })\left\|\theta_{1}-\theta_{2}\right\|_{L^{2}}.\]

The statement follows again by applying Gronwall's Lemma. 

**Proposition A.5** (Differentiability with respect to trajectories perturbations).: _Let us assume that the controlled dynamics \(\mathcal{F}\) satisfies Assumptions 1-2. Given an admissible control \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) and a trajectory \(t\mapsto x(t)=\Phi^{\theta}_{(0,t)}(x_{0})\) with \(x_{0}\in B_{R}(0)\), let \(\xi:[0,T]\to\mathbb{R}^{d}\) be the solution of the linearized problem_

\[\begin{cases}\dot{\xi}(t)=\nabla_{x}\mathcal{F}(t,x(t),\theta(t))\xi(t),\\ \xi(t)=v,\end{cases}\] (A.1)

_where \(\bar{t}\in[0,T]\) is the instant of perturbation and \(v\) is the direction of perturbation of the trajectory. Then, for every \(t\in(\bar{t},T)\), it holds_

\[|\Phi^{\theta}_{(\bar{t},t)}(x(\bar{t})+\epsilon v)-\Phi^{\theta}_{(\bar{t},t )}(x(\bar{t}))-\epsilon\xi(t)|\leq C|v|^{2}\epsilon^{2}\]

_where \(C\) is a constant depending on \(T,R,\left\|\theta\right\|_{L^{2}}\)._

Proof.: For \(t\geq\bar{t}\), let us denote with \(t\mapsto y(t):=\Phi^{\theta}_{(\bar{t},t)}(x(\bar{t})+\epsilon v)\) the solution of the modified problem, obtained by perturbing the original trajectory with \(\epsilon v\) at instant \(\bar{t}\). Then, since \(\xi\) solves (A.1), we can write

\[|y(t)-x(t)-\epsilon\xi(t)| =|\Phi^{\theta}_{(\bar{t},t)}(x(\bar{t})+\epsilon v)-\Phi^{ \theta}_{(\bar{t},t)}(x(\bar{t}))-\epsilon\xi(t)|\] \[\leq\int_{\bar{t}}^{t}|\mathcal{F}(s,y(s),\theta(s))-\mathcal{F} (s,x(s),\theta(s))-\epsilon\nabla_{x}\mathcal{F}(s,x(s),\theta(s))\xi(s)|ds\] \[\leq\int_{\bar{t}}^{t}|\mathcal{F}(s,y(s),\theta(s))-\mathcal{F} (s,x(s),\theta(s))-\nabla_{x}\mathcal{F}(s,x(s),\theta(s))(y(s)-x(s))|ds\] \[\quad+\int_{\bar{t}}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s) )||y(s)-x(s)-\epsilon\xi(s)|ds\] \[\leq\int_{\bar{t}}^{t}\left[\int_{0}^{1}|\nabla_{x}\mathcal{F}(s,x(s)+\tau(y(s)-x(s)),\theta(s)-\nabla_{x}\mathcal{F}(s,x(s),\theta(s)||y(s)- x(s)|d\tau\right]ds\] \[\quad+\int_{\bar{t}}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s)) ||y(s)-x(s)-\epsilon\xi(s)|ds\]

for every \(t\geq\bar{t}\). We now address the two integrals separately. Using Assumption 2\(-(iv)\) and the result of Lemma A.4, we obtain the following bound

\[\int_{\bar{t}}^{t}\left[\int_{0}^{1}|\nabla_{x}\mathcal{F}(s,x(s )+\tau(y(s)-x(s))),\theta(s)-\nabla_{x}\mathcal{F}(s,x(s),\theta(s))||y(s)-x( s)|d\tau\right]ds\] \[\leq\int_{\bar{t}}^{t}L_{R}(1+|\theta(s)|^{2})\frac{1}{2}|y(s)-x (s)|^{2}ds\] \[\leq\frac{1}{2}L_{\bar{R}}\left(1+\|\theta\|_{L^{2}}^{2}\right)e^{ 2L_{R}(1+\|\theta\|_{L^{1}})}|\epsilon v|^{2}\]

Similarly, for the second integral, owing to Assumption 2\(-(iv)\), we can compute:

\[\int_{\bar{t}}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s))||y(s)-x(s)- \epsilon\xi(s)|ds\leq\int_{\bar{t}}^{t}L_{R}(1+|\theta(s)|^{2})(1+\bar{R})|y(s )-x(s)-\epsilon\xi(s)|ds\]

Finally, by combining the two results together and using Gronwall's Lemma, we prove the statement.
**Proposition A.6** (Differentiability with respect to control perturbations).: _Consider the solution \(\xi\) of the linearized problem_

\[\begin{cases}\dot{\xi}(t)=\nabla_{x}\mathcal{F}(t,x^{\theta}(t),\theta(t))\xi(t) +\nabla_{\theta}\mathcal{F}(t,x^{\theta}(t),\theta(t))\nu(t)\\ \xi(0)=0\end{cases}\] (A.2)

_where the control \(\theta\) is perturbed at the initial time with \(\theta+\epsilon\nu\), when starting with an initial datum \(x_{0}\in B_{R}(0)\). Then,_

\[|\Phi^{\theta+\epsilon\nu}_{(0,t)}(x_{0})-\Phi^{\theta}_{(0,t)}(x_{0})- \epsilon\xi(t)|\leq C||\nu||_{L^{2}}^{2}\epsilon^{2}\] (A.3)

_where \(C\) is a constant depending on \(T,\bar{R},L_{R},\left\|\theta\right\|_{L^{1}}\). Moreover, we have that for every \(t\in[0,T]\)_

\[\xi(t)=\int_{0}^{t}\mathcal{R}^{\theta}_{(s,t)}(x_{0})\cdot\nabla_{\theta} \mathcal{F}(s,x^{\theta}(s),\theta(s))\nu(s)\,ds,\] (A.4)

_where \(\mathcal{R}^{\theta}_{(s,t)}(x_{0})\) has been defined in (3.10)._

Proof.: We first observe that the dynamics in (A.2) is affine in the \(\xi\) variable. Moreover, Assumptions 1-2 guarantee that the coefficients are \(L^{1}\)-regular in time. Hence, from the classical Caratheodory Theorem we deduce the existence and the uniqueness of the solution of (A.2). Finally, the identity (A.4) follows as a classical application of the resolvent map (see,e.g., in [10, Theorem 2.2.3]).

Let us denote with \(t\mapsto x(t)\) and \(t\mapsto y(t)\) the solutions of Cauchy problem (2.1) corresponding, respectively, to the admissible controls \(\theta\) and \(\theta+\epsilon\nu\). We observe that, in virtue of Lemma A.1, we have that there exists \(\bar{R}>0\) such that \(x(t),y(t)\in B_{\bar{R}}(0)\) for every \(t\in[0,T]\). Then, recalling the definition of the flow map provided in (2.2), we compute

\[|y(t)-x(t)-\epsilon\xi(t)| =|\Phi^{\theta+\epsilon\nu}_{(0,t)}(x_{0})-\Phi^{\theta}_{(0,t)}( x_{0})-\epsilon\xi(t)|\] \[\leq\int_{0}^{t}|\mathcal{F}(s,y(s),\theta(s)+\epsilon\nu(s))- \mathcal{F}(s,x(s),\theta(s))-\epsilon\dot{\xi}(s)|\,ds\] \[\leq\int_{0}^{t}|\mathcal{F}(s,y(s),\theta(s)+\epsilon\nu(s))- \mathcal{F}(s,x(s),\theta(s)+\epsilon\nu(s))\] \[\qquad\qquad-\epsilon\nabla_{x}\mathcal{F}(s,x(s),\theta(s)+ \epsilon\nu(s))\cdot(y(s)-x(s))|\,ds\] \[\quad+\int_{0}^{t}|\mathcal{F}(s,x(s),\theta(s)+\epsilon\nu(s))- \mathcal{F}(s,x(s),\theta(s))-\epsilon\nabla_{\theta}\mathcal{F}(s,x(s), \theta(s))\cdot\nu(s)|\,ds\] \[\quad+\int_{0}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s)+ \epsilon\nu(s))-\nabla_{x}\mathcal{F}(s,x(s),\theta(s))||y(s)-x(s)|ds\] \[\quad+\int_{0}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s))||y(s)- x(s)-\epsilon\xi(s)|\,ds.\]

We now handle each term separately:

\[\int_{0}^{t}|\mathcal{F}(s,y(s),\theta(s)+\epsilon\nu(s))- \mathcal{F}(s,x(s),\theta(s)+\epsilon\nu(s))-\epsilon\nabla_{x}\mathcal{F}(s, x(s),\theta(s)+\epsilon\nu(s))(y(s)-x(s))|\,ds\] (A.5) \[\leq\int_{0}^{t}\left[\int_{0}^{1}L_{\bar{R}}(1+|\theta(s)+ \epsilon\nu(s)|^{2})\tau|y(s)-x(s)|^{2}d\tau\right]ds\] \[\leq L_{R}^{3}(1+\left\|\theta\right\|_{L^{2}}+\epsilon\left\| \nu\right\|_{L^{2}})^{4}e^{2L_{R}(1+\left\|\theta\right\|_{L^{1}})}\left\|\nu \right\|_{L^{2}}^{2}\epsilon^{2}\]

where we used Assumption 2\(-(iv)\) and Lemma A.4. By using Assumption 2\(-(v)\), we obtain the following bounds for the second integral:

\[\int_{0}^{t}|\mathcal{F}(s,x(s),\theta(s)+\epsilon\nu(s))- \mathcal{F}(s,x(s),\theta(s))-\nabla_{\theta}\mathcal{F}(s,x(s),\theta(s)) \cdot\epsilon\nu(s)|\,ds\] (A.6) \[\leq\int_{0}^{t}\left[\int_{0}^{1}L_{R}|\nu(s)|^{2}\epsilon^{2} \tau d\tau\right]ds=\frac{1}{2}L_{R}\left\|\nu\right\|_{L^{2}}^{2}\epsilon^{2}.\]

Similarly, the third integral can be bounded by using Assumption 2\(-(vi)\) and Lemma A.4, and it yields

\[\int_{0}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s)+\epsilon\nu(s ))-\nabla_{x}\mathcal{F}(s,x(s),\theta(s))||y(s)-x(s)|\,ds\] (A.7) \[\leq\int_{0}^{t}L_{R}(1+|\theta(s)|+\epsilon|\nu(s)|)\epsilon|y(s)- x(s)||\nu(s)|\,ds\] \[\leq L_{R}^{2}(1+\left\|\theta\right\|_{L^{2}}+\epsilon\left\| \nu\right\|_{L^{2}})^{2}e^{L_{R}(1+\left\|\theta\right\|_{L^{1}})}\left\|\nu \right\|_{L^{2}}^{2}\epsilon^{2}.\]
Finally, the fourth integral can be bounded using Assumption 2\(-(iv)\) as follows:

\[\int_{0}^{t}|\nabla_{x}\mathcal{F}(s,x(s),\theta(s))||y(s)-x(s)-\epsilon\xi(s)|ds \leq\int_{0}^{t}L_{R}(1+\tilde{R})(1+|\theta(s)|^{2})|y(s)-x(s)-\epsilon\xi(s)| \,ds.\] (A.8)

Hence, by combining (A.5), (A.6), (A.7) and (A.8), the thesis follows from Gronwall Lemma. 

**Proposition A.7** (Properties of the resolvent map).: _Let us assume that the controlled dynamics \(\mathcal{F}\) satisfies Assumptions 1-2. Given an admissible control \(\theta\in L^{2}([0,T],\mathbb{R}^{m})\) and a trajectory \(t\mapsto x(t)=\Phi^{\theta}_{(0,t)}(x)\) with \(x\in B_{R}(0)\), for every \(\tau\in[0,T]\) the resolvent map \(\mathcal{R}^{\theta}_{(\tau,\cdot)}(x):[0,T]\to\mathbb{R}^{d\times d}\) is the curve \(s\mapsto\mathcal{R}^{\theta}_{(\tau,s)}(x_{0})\) that solves_

\[\begin{cases}\frac{d}{ds}\mathcal{R}^{\theta}_{(\tau,s)}(x)=\nabla_{x} \mathcal{F}(s,\Phi^{\theta}_{(0,s)}(x),\theta(s))\cdot\mathcal{R}^{\theta}_{( \tau,s)}(x)\quad\text{for $a.e.$ $s\in[0,T]$},\\ \mathcal{R}^{\theta}_{(\tau,\tau)}(x)=\mathrm{Id}.\end{cases}\] (A.9)

_Then for every \(\tau,s\in[0,T]\), there exists a constant \(C_{1}\) depending on \(T,R,\left\|\theta\right\|_{L^{2}}\) such that_

\[|\mathcal{R}^{\theta}_{(\tau,s)}(x)|:=\sup_{v\neq 0}\frac{|\mathcal{R}^{ \theta}_{(\tau,s)}(x)\cdot v|}{|v|}\leq C_{1}.\] (A.10)

_Moreover, for every \(x,y\in B_{R}(0)\), there exists a constant \(C_{2}\) depending on \(T,R,\left\|\theta\right\|_{L^{2}}\) such that_

\[|\mathcal{R}^{\theta}_{(\tau,s)}(x)-\mathcal{R}^{\theta}_{(\tau,s)}(y)|:=\sup _{v\neq 0}\frac{|\mathcal{R}^{\theta}_{(\tau,s)}(x)\cdot v-\mathcal{R}^{\theta }_{(\tau,s)}(y)\cdot v|}{|v|}\leq C_{2}|x-y|.\] (A.11)

_Finally, if \(\theta_{1},\theta_{2}\) satisfy \(\left\|\theta_{1}\right\|,\left\|\theta_{2}\right\|\leq\rho\), then there exists a constant \(C_{3}\) depending on \(T,R,\rho\) such that_

\[|\mathcal{R}^{\theta_{1}}_{(\tau,s)}(x)-\mathcal{R}^{\theta_{2}}_{(\tau,s)}(x )|:=\sup_{v\neq 0}\frac{|\mathcal{R}^{\theta_{1}}_{(\tau,s)}(x)\cdot v- \mathcal{R}^{\theta_{2}}_{(\tau,s)}(x)\cdot v|}{|v|}\leq C_{3}\|\theta_{1}- \theta_{2}\|_{L^{2}}.\] (A.12)

Proof.: We first prove the boundedness of the resolvent map. Let us fix \(v\in\mathbb{R}^{d}\) with \(v\neq 0\), and let us define \(\xi(s):=\mathcal{R}^{\theta}_{(\tau,s)}(x)\cdot v\) for every \(s\in[0,T]\). Then, in virtue of Assumption 2\(-(vi)\), we have:

\[|\xi(s)|\leq|\xi(\tau)|+\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma, \Phi^{\theta}_{(0,\sigma)}(x),\theta(\sigma))\right|\left|\xi(\sigma)\right|d \sigma\leq|v|+L_{R}\int_{0}^{t}(1+\theta(\sigma)^{2})|\xi(\sigma)|\,d\sigma,\]

and, by Gronwall's Lemma, we deduce (A.10). Similarly as before, given \(x,y\in B_{R}(0)\) and \(v\neq 0\), let us define \(\xi^{x}(s):=\mathcal{R}^{\theta}_{(\tau,s)}(x)\cdot v\) and \(\xi^{y}(s):=\mathcal{R}^{\theta}_{(\tau,s)}(y)\cdot v\) for every \(s\in[0,T]\). Then, we have that

\[|\xi^{x}(s)-\xi^{y}(s)| \leq\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma,\Phi^{ \theta}_{(0,\sigma)}(x),\theta(\sigma))\xi^{x}(\sigma)-\nabla_{x}\mathcal{F}( \sigma,\Phi^{\theta}_{(0,\sigma)}(y),\theta(\sigma))\xi^{y}(\sigma)\right|\,d\sigma\] \[\leq\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma,\Phi^{ \theta}_{(0,\sigma)}(x),\theta(\sigma))-\nabla_{x}\mathcal{F}(\sigma,\Phi^{ \theta}_{(0,\sigma)}(y),\theta(\sigma))\right|\left|\xi^{y}(\sigma)\right|d\sigma\] \[\quad+\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma,\Phi^{ \theta}_{(0,\sigma)}(x),\theta(\sigma))\right|\left|\xi^{x}(\sigma)-\xi^{y}( \sigma)\right|d\sigma\] \[\leq C_{1}|v|\int_{\tau}^{s}L_{R}(1+\theta(\sigma)^{2})\left|\Phi^ {\theta}_{(0,\sigma)}(x)-\Phi^{\theta}_{(0,\sigma)}(y)\right|\,d\sigma\] \[\quad+\int_{\tau}^{s}L_{R}(1+\theta(\sigma)^{2})|\xi^{x}(\sigma)- \xi^{y}(\sigma)|\,d\sigma,\]

where we used (A.10) and Assumption 2\(-(iv)\). Hence, combining Lemma A.2 with Gronwall's Lemma, we deduce (A.11). Finally, we prove the dependence of the resolvent map on different controls \(\theta_{1},\theta_{2}\in L^{2}([0,T];\mathbb{R}^{m})\). Given \(x\in B_{R}(0)\) and \(v\neq 0\), let us define \(\xi^{\theta_{1}}(s):=\mathcal{R}^{\theta_{1}}_{(\tau,s)}(x)\cdot v\) and \(\xi^{\theta_{2}}(s):=\mathcal{R}^{\theta_{2}}_{(\tau,s)}(x)\cdot v\) for every \(s\in[0,T]\). Then, we
compute

\[|\xi^{\theta_{1}}(s)-\xi^{\theta_{2}}(s)| \leq\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma,\Phi_{(0, \sigma)}^{\theta_{1}}(x),\theta_{1}(\sigma))\xi^{\theta_{1}}(\sigma)-\nabla_{x} \mathcal{F}(\sigma,\Phi_{(0,\sigma)}^{\theta_{2}}(x),\theta_{2}(\sigma))\xi^{ \theta_{2}}(\sigma)\right|\,d\sigma\] \[\leq\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma,\Phi_{(0, \sigma)}^{\theta_{1}}(x),\theta_{1}(\sigma))-\nabla_{x}\mathcal{F}(\sigma,\Phi _{(0,\sigma)}^{\theta_{2}}(x),\theta_{2}(\sigma))\right|\left|\xi^{\theta_{1}} (\sigma)\right|d\sigma\] \[\quad+\int_{\tau}^{s}\left|\nabla_{x}\mathcal{F}(\sigma,\Phi_{(0, \sigma)}^{\theta_{2}}(y),\theta(\sigma))\right|\left|\xi^{\theta_{1}}(\sigma) -\xi^{\theta_{2}}(\sigma)\right|d\sigma\] \[\leq C_{1}|v|\int_{\tau}^{s}L_{R}(1+\theta_{1}(\sigma)^{2})\left| \Phi_{(0,\sigma)}^{\theta_{1}}(x)-\Phi_{(0,\sigma)}^{\theta_{2}}(x)\right|\,d\sigma\] \[\quad+C_{1}|v|\int_{\tau}^{s}L_{R}(1+|\theta_{1}(\sigma)|+| \theta_{2}(\sigma)|)|\theta_{1}(\sigma)-\theta_{2}(\sigma)|\,d\sigma\] \[\quad+\int_{\tau}^{s}L_{R}(1+\theta(\sigma)^{2})|\xi^{\theta_{1} }(\sigma)-\xi^{\theta_{2}}(\sigma)|\,d\sigma,\]

where we used Assumption 2\(-(iv)\)-\((vi)\). 
