{'fraction_non_alphanumeric': 0.04354445387548334, 'fraction_numerical': 0.011088207159505814, 'mean_word_length': 4.077228264587181, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural ODE control for classification, approximation and transport

We analyze Neural Ordinary Differential Equations (NODEs) from a control theoretical perspective to address some of the main properties and paradigms of Deep Learning (DL), in particular, data classification and universal approximation. These objectives are tackled and achieved from the perspective of the simultaneous control of systems of NODEs. For instance, in the context of classification, each item to be classified corresponds to a different initial datum for the control problem of the NODE, to be classified, all of them by the same common control, to the location (a subdomain of the euclidean space) associated to each label. Our proofs are genuinely nonlinear and constructive, allowing us to estimate the complexity of the control strategies we develop. The nonlinear nature of the activation functions governing the dynamics of NODEs under consideration plays a key role in our proofs, since it allows deforming half of the phase space while the other half remains invariant, a property that classical models in mechanics do not fulfill. This very property allows to build elementary controls inducing specific dynamics and transformations whose concatenation, along with properly chosen hyperplanes, allows achieving our goals in finitely many steps. The nonlinearity of the dynamics is assumed to be Lipschitz. Therefore, our results apply also in the particular case of the ReLU activation function. We also present the counterparts in the context of the control of neural transport equations, establishing a link between optimal transport and deep neural networks.


Introduction
In this paper, we prove that the dynamical control properties of Neural Ordinary Differential Equations (NODEs) allow for the understanding of some of the main Deep Neural Networks' properties. In particular, we will focus on data classification and universal approximation.
Given a dataset {(x i , y i )} N i=1 ⊂ R d × R d , x i being the input and y i the corresponding label, Supervised Learning (SL) aims to find a map φ : R d → R d that associates x to y for later use on the allocation of a label to new unknown data. In classification, the labels take values in a finite set of R d (typically associated with a subset of the natural numbers), while in regression, the labels correspond to continuous values of R d .
Neural networks (NN) are widely used in this context. A shallow neural network can be written as a linear combination of the form: where σ : R d → R d is the so called activation function (more details will be given later on), and W i , A i ∈ R d×d , b i ∈ R d , i ∈ {1, ..., n}.
In the seminal paper by Cybenko [10], it was proved that, under certain assumptions on σ, any continuous function in a compact set can be approximated by NN of the form (1.1) in the supremum norm as n → ∞. The product dn is the width of the network. This density result is the so called Universal Approximation Theorem and it has been the object of intensive study and several improvements and generalizations (see [19,22,29]).
A deep neural network is obtained from the composition of shallow neural networks of the type (1.1) and in recent years they have played an important role in diverse applications, [21]. Each of this composition steps is refereed to as a layer and the depth of the network is the total number of compositions operated. Deep neural networks can also be interpreted as discrete dynamical systems and can be employed, in particular, to prove universal approximation theorems, [11]. The parameters W, A and b entering in the design of the networks are often found via optimization, using, for instance, a least-squares error type functionals [8], in the so-called learning process.
Residual neural networks (ResNets) (see [18]) refer to the variant in which an "inertia" term is added in each layer, leading to a discrete dynamical system of the form (1.2) x l+1 = x l + hW l σ(A l x l + b l ), where h is a positive real number. ResNets can also be understood as an Euler discretization scheme of an ODE, the so called Neural ODEs (NODEs) ( [6,17,42]): This paper aims to develop a dynamical control theoretical analysis of NODEs (1.3). The time-dependent coefficients A(t), W (t) and b(t) are not chosen using least squares and optimization approaches (as in [14] and the references therein). Rather, we adopt a dynamic controllability approach according to which the time-dependent parameters A(t), W (t) and b(t) are viewed as controls and are built in a constructive way, being defined in a piecewiseconstant switching manner, adapted to the distribution of data to be classified. This allows for algorithmic constructions that permit to measure the complexity of the controls needed to achieve the main goals, in particular, classification and universal approximation. Our approach does not require either the regularity of the nonlinearity that classical ODE control methods need, often based on linearization and Lie brackets. In this way our results apply also in the case of the Rectified Linear Unit (ReLU) activation function which is Lipschitz but fails to be C 1 .
The properties and problems we study are addressed in an ordered and systematic manner, according to their complexity.
The first one we address is the classification problem. From a control perspective, it can be formulated as a simultaneous control problem. The data to be classified are viewed as the initial data of the Cauchy problem associated to (1.3). Classification requires to build controls A(t), W (t) and b(t) that simultaneously drive all initial data to their corresponding final destination, allocated depending on their label. The final targets that we prefix are disjoint sets corresponding to a partition of the space R d . The choice of the partition is not unique, and the specific control problem to be addressed depends on it. Although our methods allow to handle arbitrary partitions, for the sake of simplicity, we consider the particular one in which the space is partitioned in parallel strips, which suffices to achieve the goal of classification. In Theorem 1, Section 3, we develop our new constructive proof of classification. As mentioned above, our construction allows to master the complexity of the required classifying controls, in terms of the distribution and structure of the dataset to be classified.
Theorem 1 assures that each initial data can be driven simultaneously to the strip corresponding to its label. But, in fact, our constructions allow driving the system to much more precise targets so that each trajectory approaches arbitrarily a prescribed target point in R d . This is proved in Section 4, Theorem 2.
Both results hold in an arbitrary time horizon T > 0. This can be easily achieved by scaling arguments, the general case T > 0 being reduced to T = 1. Obviously, the L ∞ -norm of the controls employed depends then monotonically and linearly on 1/T , T being the length of time-horizon.
From a control perspective, these results can be viewed as particular instances of simultaneous or ensemble controllability properties. This topic has been considered in the control community for a variety of systems [24,25,39]. From that perspective, the results in the present paper are new since the ODE system under consideration is the same for all initial data and the number N of the controlled initial data is arbitrarily large, while the number of available controls is limited: 2d 2 + d. This kind of simultaneous control result would be impossible for linear finite or infinite-dimensional systems, where the simultaneous control never occurs when different trajectories solve the same system. The specific nature of the nonlinearity of the NODEs and, as explained above, the very structure of the activation function is the key of our construction.
Our proofs are genuinely nonlinear and strongly rely on the structure of the activation function σ that models the nonlinearity of the NODE (1.3) under consideration. The key ingredient of our proofs is that the activation function allows splitting the phase-space into two half-spaces, so that one remains invariant under the action of the control, while the other one evolves in the direction we wish, approaching the prescribed final target. Such construction is not possible for the classical nonlinear ODE systems in Mechanics in which all trajectories evolve simultaneously, driven by the nonlinear flow. The proofs are constructive, allowing us to estimate the complexity of the needed piecewise-constant controls for achieving such results. Moreover, our techniques are of dynamic nature, and they are not based on optimization or linearization arguments, that would not allow to obtain the global results we prove.
Once the classification problem is solved and its upgraded simultaneous controllability version is proved in Theorems 1 and 2, our second main objective is to show that these techniques, when refined, allow also to prove the Universal Approximation Theorem (Section 5). More precisely, in Theorem 4 we prove that, given a simple function in f : Ω → R d (a linear finite combination of characteristic functions of measurable sets), under the assumption that characteristic sets have finite perimeter, for any bounded set Ω ⊂ R d and > 0 arbitrary being given, we can find controls A(t), W (t) and b(t) (depending in particular on ) such that the associated flow φ T of (1.3) satisfies: The proof, which is of constructive nature, yields a connection between the complexity of the function to be approximated and the cost of control in terms of the box-counting dimension (see [15,Chapter 2]). The control complexity is measured in terms of its L ∞ -norm and the number of switches and it depends on, from one side, the level of approximation we desire, and from the other, on the complexity of the target function. In fact, any simple function can be approximated and its control cost approximated provided that the box-counting dimension of the boundaries of its characteristic sets is smaller than the dimension of the ambient space.
Our third main objective is to formulate these results in terms of control of transport equations, where the dynamical control perspective is made more evident. In Section 6 (see Theorems 5 and 6), we prove that the Neural Transport equation is approximately controllable. Furthermore, we will see that one can achieve a simultaneous control result for a Neural Transport system of equations of the form (1.4) under certain conditions on the initial data and targets (roughly, having disjoint supports), which mimic the configuration of classification problems.
1.1. Roadmap and related work. This paper is organized as follows. In Section 2, we introduce the notation and formulate the main problems more precisely, presenting some of the fundamental tools that will systematically be used in our proofs. Section 3 is devoted to the classification problem and in Section 4 we reinterpret and further develop it in the simultaneous control context. Section 5 is devoted to proving the Universal Approximation Theorem, while Section 6 to analyze neural transport equations. Section 7 is devoted to present some open problems and further directions of research.
For the sake of simplicity of the exposition most of our results are presented in the particular case of the ReLU activation function, although the methods we develop apply for larger classes of activation functions that will be made clear below.
We conclude this introduction with a discussion of some of the related literature. Control issues for NODEs have been addressed from different perspectives. In [35], for instance, the authors analyze the classical control problem of NODEs. Note however that this is done in the context in which one typically considers one single trajectory and not the simultaneous control problems that classification requires.
In the context of simultaneous control of Neural ODEs, in [9] a Lie bracket technique is developed for certain smooth activation functions that allow to achieve universal interpolation results, similar to those in Section 4. We also refer to [1] for further developments on the geometric Lie bracket interpretation of these concepts in the context of classification, also using smooth vector fields. Note however that our approach can be implement for Lipschitz continuous activation functions as the ReLU as in Figure 2.1.
In [14], the authors show a simultaneous control result for sufficiently regular activations under smallness conditions on the targets, and employing linearization arguments, to later derive Turnpike results for NODEs. Our techniques, of constructive nature, are of independent interest and pave the path to analyze Turnpike properties of NODEs in a much broader context.
Recently, in [23], universal approximation results were proved by means of NODEs, see also [37,38]. In Section 5 we provide a simpler and more comprehensive proof of those results which allows, in particular, to handle the complexity of the controls in terms of the data to be classified or the target functions.
The results in Section 6, which provide a transport formulation of Deep Learning, making use of the classical link between ODEs and transport equations, can be seen as bilinear-type simultaneous control results for Neural transport equations and systems (Theorems 5 and 6). The bilinear control of transport equations has also been considered (see, for instance, [12,13]) with controls localized in space motivated by collective behavior problems [28]. The nature of the controls we employ here and the dynamics that the systems under consideration, exploiting the very features of activation functions, are completely different. Furthermore, the use of Neural transport equations allows a probabilistic formulation of the classification problem.
Our results for Neural Transport equations can also be related to optimal transport, that can be formulated as a minimal kinetic energy controllability problem as in [3], [2] and the references therein.
1.2. Notation. Throughout this article we will use the following notation.
-We understand by E q F the existence of a constant C(q) depending on q such that E ≤ C(q)F . Analogously, we denote by E ∼ q F the existence of C(q) such that E = C(q)F . -Let x ∈ R d , we denote by x (k) the k-th component of the vector x.
-Given ω ⊂ R d , we denote by diam the quantity diam(ω) := sup x1,x2∈ω Consider now the Neural ODE The coefficients A, W and b play the role of controls, With the notations above (2.3) is a system of d coupled differential equations. The first question we address is the classification problem, which, in the context of the NODEs (2.3), can be reformulated as follows. We fix a time horizon 3). In fact {x i }, the data to be classified, will be taken to be the initial data at time t = 0 for the Cauchy problem (2.3), while {y i } are the corresponding labels that will allow us to define the targets at the final time t = T .
To reformulate this problem we consider a partition of R d into M (the number of distinct labels) disjoint open sets with nonempty interior. We allocate one of these sets to each label. From a control theoretical perspective the classification problem can be reformulated as follows: to find controls W, A, b bringing each initial data to the corresponding set along the trajectories of (2.3).
Note that this problem can be considered as a simultaneous or ensemble control one, in the sense that the same controls are aimed to control all trajectories.
The Partition of the Euclidean Space. Our techniques can be applied for an arbitrary partition of the euclidean space. But the complexity of the control maps will depend on the structure, geometry and topology of the partition. To simplify the construction of the controls and our presentation we shall consider a particular partition of the space into parallel strips that we introduce now.
Let M be the number of different labels y i , i.e.
x (1) −axis Partition of the Euclidean plane (d = 2) into three parallel strips (M = 3). Nine points to be classified are denoted by three different colors and shapes corresponding to their labels (black triangles, red circles and green diamonds) and they are distributed in an unclassified manner. The goal of the approximation process, that we aim to conduct by controlling the flow of the Neural ODE, is to drive each subclass of three points to the corresponding label, i.e. to the strip corresponding to its color.
We split the space R d into M disjoint sets S := {S m } M m=1 , each set associated with one label. For simplicity of the presentation we set the partition into strips The Control Process. Given controls W, A, b we denote the solution of (2.3) by φ T (x; A, W, b).
where m(i) ∈ {1, ..., M } denotes the index of the label corresponding to the input x i .
The Prediction. Once this map φ T has been built out of the corresponding controls A, W, b, it can be used for prediction purposes, assigning to each new datum x the corresponding label y m associated with the subset S m of the partition satisfying φ T (x; A, W, b) ∈ S m .
2.2. The fundamental operations: Place, freeze, compress/expand and translate. All our results will be proved by the same methodology. It consists of building piece-wise constant controls which, concatenated, will assure that the NODE fulfills the needed requirements. Each of these constant actions exploits the activation function's properties, the ReLU, in this presentation. In each step the available control parameters W, A, and b are chosen so that the ODE induces a transformation of the Euclidean space R d which leaves invariant half of the space while the other evolves towards the desired direction. Obviously, the choice of the correct controls is essential to assure that the invariant and the moving half-spaces are the appropriate ones, together with the direction of motion. For simplicity of the proofs we will also choose separating hyperplanes parallel to the Cartesian axes. The algorithm could be optimized by choosing oblique hyperplanes adapted to the datasets to be classified (in the context of the classical problem) or the function to be approximated (in the context of universal approximation). But, for simplicity of the presentation, we will only consider hyperplanes parallel to the Cartesian axes.
(1) Consider a hyperplane of the form: This hyperplane can be represented as Ax + b = 0 with the choice of the controls A, b as below The hyperplane divides the Euclidean space in two half-spaces. Here and in the sequel δ jk stands for the Kronecker delta. (2) When applying the vector-valued activation function σ to Ax + b, taking into account that the ReLU activation function σ vanishes for negative inputs, the choice of the sign of A and b is crucial to determine which half-space remains invariant; See Figure 2.3.A. (3) Choosing W to be a rotation matrix, when applied to σ(Ax + b), we can change the direction of the vector field in the direction we desire, see Figure 2.3.B. Note that, in particular, we can choose W so that the hyperplane is attractive or repulsive with respect to the moving half-space, but also we can generate a movement parallel to the hyperplane. Throughout the proofs, we will only consider controls W such that the resulting vector fields point in some of the canonical euclidean directions. As it occurs for the choice of the hyperplanes, proofs could be implemented with a lower number of concatenations if W were chosen in oblique directions, depending on the dataset. But, for simplicity of the exposition, we shall only employ these elementary W 's. Summarizing, we can: (1) Place and orient the dividing hyperplane arbitrarily.
(2) Choose the invariant and the active half-spaces. In particular we may choose the halfspace that we aim to freeze. (3) Choose the orientation of the vector-field σ(Ax + b) so that either -it induces an expansion or contraction of the active half-space. or -it induces a translation on the active half-space parallel to the dividing hyperplane. As we shall see, the concatenation of these elementary transformations will suffice to achieve our classification and approximation goals.

Classification
This section is devoted to present Theorem 1 and its proof, which allows to classify an arbitrary finite dataset in strips, according to their label. The proof is based on combining the essential movements explained in subsection 2.2. For simplicity we state and prove the Theorem for the ReLU (2.1) activation function. However our arguments apply to more general activation functions σ, as we will make precise later in Remark 3.2.
Theorem 1. Let d ≥ 2 and M ≥ 2 be natural numbers and let σ be as in (2.1).
Then, for every T > 0, there exist control functions A, W ∈ L ∞ (0, T ); R d×d and b ∈ L ∞ (0, T ), R d such that the flow associated to (2.3), when applied to all initial data {x i } N i=1 , classifies the trajectories by strips (2.4), according to their labels S mi being the subset corresponding to the label y i , with i ∈ {1, ..., M }. Furthermore, controls are piecewise constant with a finite number of switches of the order of O(N ). Therefore, they also lie in BV .
Remark 3.1 (d = 1). Theorem 1 requires that d ≥ 2. If d = 1 and x 1 < x 2 , then for all choices of the controls W, A, b, the associated flow will fulfill that φ T (x 1 ; A, W, b) < φ T (x 2 ; A, W, b) for all T > 0. This is due to the ODE formulation and the uniqueness of solutions for the corresponding Cauchy problem. Note that this restriction does not necessarily limit time-discrete NNs as, for instance, the one in (1.2), since it does not necessarily fulfill this monotonicity property when the step-size is large. Our proofs apply in an ODE setting and therefore are limited to d ≥ 2. As a corollary, they apply as well for the time-discrete ResNets (1.2) provided the step-size h is small enough.
x (k) −axis Proof. As we mentioned above, the actual value of the final time T > 0 is irrelevant since, by scaling, it can be set to be, in particular, T = 1.
Therefore we present a strategy that in a finite number of steps allows to achieve the classification in a finite time T * > 0. The scheme of the proof is the following: (1) Preparation of the dataset. The goal of this first step is to assure that the data to be classified fulfill the condition: j . Obviously, for that to hold, one needs to exploit the dynamics of the system by properly choosing the controls. This condition is relevant to assure that data can be efficiently separated by hyperplanes parallel to the Cartesian axes.
(2) Classification of the points so that the corresponding trajectories reach the allocated strips.
(3) Time rescaling allows to assure that the goal is achieved in the given final time T > 0.
We now present in more detail the process of each of these steps.
(1) Preparation of the dataset. We need the data to be classified to fulfill the condition do not fulfill this property, by a proper choice of the controls A, W, b in a time τ > 0, we may assure that The construction of the controls guaranteeing that this separation property is fulfilled can be made by induction.
• Suppose we have only two points, x i and x j , sharing the first component x (1) .
j . We then take the mean value r = (x j )/2 and consider the matrix A and the vector b defined by .., d}.
x (1) −axis The hyperplane {Ax + b = 0} separates the two points x i and x j . We choose W such that the resulting vector field σ(Ax + b) can change the x (1) -component of one of the elements. For doing so, we set acting as a translation in the half-space. Solving the NODE in a sufficiently small time τ > 0, we can assure that as shown in Figure 3.1.qquad Note that τ > 0 needs to be taken small enough so to guarantee that when assuring that (3.2) holds, the other points do not lose the distinction of all their coordinates.
• This argument can be applied recursively on all pairs of data sharing the component x (1) . Once the separation property is guaranteed we can reset t = 0 to proceed to classification. This first step shows that, without loss of generality, we can assume that the dataset to be classified fulfills the separation property (3.1).
(2) Classification. By the previous step we have a set of points that, in particular, do not share the x (1) -component.
Recall also that to each initial datum x i , i ∈ {1, ..., N }, it corresponds the strip S m(i) as a target, associated to its label m(i) ∈ {1, ..., M }. Whenever a datum x j lies in the set corresponding to its target, i.e. x j ∈ S m(j) , no action is needed, which corresponds to simply taking W = 0, i.e. the trivial dynamics. But for each j ∈ {1, ..., N } such that x j / ∈ S m(j) one needs to force the allocation of the points to the corresponding strip by a suitable choice of the controls.
We do it iteratively as follows. We start with any point x j / ∈ S m(j) . To simplify the notation, let us assume that this is the point x 1 corresponding to the index j = 1. Then: In this manner the green diamond that was originally located in the strip S 3 ends up lying in the one corresponding to its label (green color) S 1 .
(a) we choose the hyperplane {Ax + b = 0} with A and b of the form i . We choose W so that the field is parallel to the hyperplane, for simplicity: Figure 3.2. We solve (2.3) up to a time T 1 so that (b) Then, we set r to be such that x and we choose our controls A and b as In this way σ(Ax + b) yields a null field (= 0) in all points that satisfy x Now we choose W = W as before and we solve (2.3) up to a time T 2 so that: (c) Then, we can choose a hyperplane represented by with r so that the image of x 1 under the previous transformations, the one to be properly classified, lies alone to one side of the hyperplane. We can then choose W so that the corresponding vector field pushes this point towards S m (1) . Note that this procedure not only assures that the image of x 1 gets into S m(1) but also allows to determine exactly the location of the first coordinate of that point. We choose a target location z 1 . Moreover, note that his process keeps the x (1) coordinate of all points i = 1 unaltered while allocating x 1 to S m(1) (see Figure 3.3 for the final result). Observe that the process above will not change the proper allocation of the points that already lie in the corresponding strip. Applying this construction recursively, in finitely many steps, we can guarantee the correct collocation all points that initially did not fulfill the classification criterion, i. e. such that x j / ∈ S m(j) . (3) Time rescaling. The argument above assures that the complete classification can be achieved in some time T * > 0. Classification can be assured in any time T > 0, by time-scaling, s = (tT )/T * . Indeed, we then have and the multiplicative term T * /T can be absorbed within the control W , keeping A and b unchanged (other than time rescaling).

Remark 3.2 (On the activation function).
Although the proof was presented in the particular case of the ReLU, the properties of the activation function σ that we used in the proof are the following: (1) σ is Lipschitz continuous, σ is positive in (0, +∞).
The above properties of σ suffice to assure that the main steps of our proof apply as represented in Figure 3.2. In fact, these conditions could also be further relaxed to consider activation functions that do not necessarily vanish in (−∞, 0] but simply decay as s tends to −∞. But this would increase the cost of control. We shall not pursue this case here but the main steps of the proof would be essentially the same. Our proof shows that the role of the control b is very much related to the distribution and spreading of the data along the Euclidean space. On the other hand, the controls W and A are of norm one in L ∞ . Their complexity can be estimated in terms of the number of switchings, which is reflected on the BV -norm.
In our proof, the cost of the control is measured in terms of the needed time of control, which is the accumulation of the time spent in each iteration. The length of the n-th time interval T n for which the control is constant is larger when points whose ordering needs to be swapped are very close to each other, and this leads to a higher cost. Note also that in our proof, we did not exploit the possible structuring of data into clusters. Indeed, when a given point to be classified is embedded into a cluster corresponding to the same label, one may take advantage of that fact to classify the whole cluster simultaneously and, in this manner, reduce the number of switchings.
By the rescaling step 3 of Theorem 1, we have seen that classification can be achieved in any time horizon, in particular in T = 1. A natural way of measuring the cost of control is to fix T = 1 and then to estimate the L ∞ and BV norm of the controls.
A complete analysis of the cost of classification in terms of the complexity of the data-set is a challenging problem that requires substantial further work. However, we will give an insight of this issue relating it with the control cost in Section 5, Remark 5.1.

Simultaneous Control
The proof of the previous classification Theorem shows that, other than classifying the data to the corresponding strips, we can also, for instance, completely control the first coordinate x (1) of all points to be classified at the final time, with the same controls.
Based on that observation, in this section, we explore to which extent the classification can be enhanced to allocate to each datum to be classified a given point of destination. As we shall see, this can be done in an approximate manner, finding controls driving the ensemble of data to be classified as close as we wish from the given final locations chosen. This is an approximate simultaneous or ensemble controllability result. In fact, as we shall see, the exact simultaneous control can be achieved when initial data and targets are well separated. • Note that it is impossible to bring two different points to the same target by means of a Lipschitz vector field. Therefore, the control to distinct targets is the sharpest result one can achieve. • By continuity, this also implies the approximate simultaneous control when the initial data or targets are not completely distinct. Indeed, let us assume for instance that z i = z j . We can then slightly modify these targets to new points points z j = z j + , with | | ≤ arbitrarily small. The exact simultaneous control of the system to the targets {z i } N i=1 ⊂ R d is then feasible and assures the control to a distance of the original In other words, even when the targets Proof. We proceed in several steps.
(1) Preparing the target. The backward NODE is also a solution of the forward NODE (2.3) by the controls . Therefore, we can apply the same argument done in Step 1 of Theorem 1 for preparing the targets. At the price of N extra switches, we can find controls W, A and b such that the targets fulfill (2) Control the first component. Applying the arguments in the proof of Theorem 1, we may control the first component x (1) of each trajectory. In this way we can deduce there existence of A, W and b such that: (3) Controlling to the targets. Now, we consider the hyperplanes {x (1) = c}, for c negative enough. Let us take We use W to rotate the field in the direction to control the point φ T (x 1 ; A, W, b) to its target z 1 with time T > 0. Then we set c = z 1 + δ for δ small enough and we consider the hyperplane x (1) = c and apply the same argument, finding W that controls φ T (φ T (x 2 ; A, W, b); A , W , b ) to z 2 . Note that z 1 remains static along this deformation (see Figure 4.2).
Since z for all i, j holds, we can proceed inductively. After a time-rescaling, as in Theorem 1, we conclude the proof of the simultaneous controllability.   Remark 4.4 (Complexity of the controls II). Note that, with respect to the classification result and Remark 3.3, the number of switches in this simultaneous control result has increased by 2N , remaining of the order of N . However, its norm will depend also on the chosen targets. The reason is that the gap for the choice of a hyperplane in step 3 of the proof of Theorem 2 can be small, depending on the target configuration. This will imply the smallness of the vector field chosen, and therefore the needed time will be larger.
Note however, that this discussion is specific for the controls we have constructed. The analysis of optimal control strategies is a challenging open problem.
Remark 4.5 (Nonlinearities with drift). Note that the simultaneous control result we have proved would also hold in an approximate manner for systemṡ with f bounded. Indeed, the control W can enhance the effect of the nonlinearity making f negligible. Whether an exact simultaneous control result in the spirit of Theorem 2 can be obtained in this setting is an interesting open problem. To be more precise, let us consider a hyperplane with normal vector a and displaced by b from the origin. Then, the sign of scalar product a, x − b determines whether x is located at one side or the other of the hyperplane. Therefore, applying the scalar function σ to a, x − b , In the expression above only 3d controls are activated. But, taking into account that a, b is a scalar value, the total number of active controls is 2d + 1.
As a conclusion of this discussion we see that our proofs hold with controls W, A and b with the following particular structure: . . .
, involving a total of 2d + 1 controls. The form of dynamics (2.3) is richer since it allows for more general matrices W, A. But the arguments above show that the actual dimension of the needed controls can be drastically diminished. This is a manifestation of the fact that the results own this paper can be achieved assuring further sparsity or polarization conditions on the controls.

Universal Approximation
In this section we prove the Universal Approximation Theorem using the methods developed to treat simultaneous controllability.
We mainly consider simple functions of the form α m χ Ωm , as in Figure 5.1. Our goal is to show how, controlling a NODE, such a function can be approximated in L 2 loc . Obviously, by density, the result extends immediately to any bounded measurable function. But, as we shall see, the complexity of the controls we will build depends on the structure of the simple function and, more precisely, on the box-counting dimension of the boundaries of the sets Ω m . In the context of NODEs the goal of approximating a function can be interpreted as a simultaneous control problem for an uncountable ensemble of trajectories. We do it by viewing the values f at each point x as the target or label that each x have to be associated with along the trajectories of the NODE. Consequently, the proofs of the main results of this section will employ the ideas of the previous ones, but with substantial added work. For simplicity, we assume that all the boundaries of the sets Ω m constituting the simple target function f , have a finite perimeter. Later on in Remark 5.1 we shall show that the same proof works under weaker conditions, to the price of increasing the complexity of the needed controls. We will see how the box-counting dimension of the boundary of the characteristic sets arises in the control cost. The approximation can only be guaranteed over bounded sets Ω, as we shall see in Section 6 (Remark 6.8). Otherwise the use of unbounded controls would be required.
More precisely, the main result of this section is as follows: Theorem 3 (Universal Approximation). Let d ≥ 2 and T > 0.
Let Ω ⊂ R d be a bounded set, and σ be the ReLU activation function in (2.1). Then, for any f ∈ L 2 (Ω; R d ) and > 0 there exist A, W ∈ L ∞ ((0, T ); R d ) and b ∈ L ∞ ((0, T ); R d ) such that the flow generated by (2.3), φ T (·; A, W, b), satisfies: Theorem 3 is a direct consequence of the following theorem: Setting T = 1, the norm the controls and their switches satisfy: The number of switches of the controls A, W, b is of the order of −2d(d−1) .
Proof. Let us first explain the main ideas of the proof.
In the NODE formulation of this result, each input x ∈ R d has a target f (x) ∈ R d . The actual value of f (x) depends on the characteristic set in which x lies. Ideally, we would like to find controls A, W and b such that for every x ∈ Ω we have that But, due to the discontinuity of f (·) and the continuous dependence of the solutions of the ODE with respect to the initial data, this task is impossible. Thus, we have to relax the problem so to formulate our goal in an approximate manner.
First, we will define for any h > 0 small enough, a partition of Ω in a way that Then, the control problem we seek to solve is the following. For all ν > 0, find A, W and b (that depend on h) satisfying that the associated controlled flow φ T (·, A, W, b) satisfies: with K independent of h and ν. This result can be achieved by means of a careful implementation and further development of the ideas of the proof of the simultaneous controllability result of Theorem 2. The proof is structured as follows. In Step 1 we construct the partition (Ω h , Ω c ). Furthermore, we will split Ω c in d-dimensional hyperrectangles, each of them having a single target associated in R d .
Step 2 concerns the compression of Ω c and of its hyperrectangles plus an adaptation of Theorem 1. The goal is to arrive to a configuration in which one will be able to apply the simultaneous controllability, Theorem 2. Finally, in Step 3, we apply the simultaneous controllability, Theorem 2.

Covering of boundaries of the characteristic sets.
Thanks to the compactness of the boundary Γ = M m=1 ∂Ω m it can be covered with a finite number of hypercubes of side h. Let us denote this cover by Γ h . The cover Γ h can be chosen in a way that each edge of each hypercube is oriented in the direction of an element of the canonical basis. Since the length of Γ is finite, and its dimension is d − 1, the number, N Γ , of needed hypercubes of size h in the cover Γ h is: Since each hypercube has volume h d , the volume of the cover Γ h is of the order of h. Note that, Γ h (represented by in white in Figure 5.2a) contains the jumps of the simple function f to be approximated. ∀x, y ∈ H l , f (x) = f (y), i.e. each hyperrectangle has at most a single associated value α m . For finding such partition we proceed as follows. We consider all hyperplanes that define the hypercubes in Γ h (see Figure 5.2b for a qualitative illustration). Altogether we have where N Γ,k is the number of hyperplanes orthogonal to the k-th vector of the canonical basis and let c k,N Γ,k +1 = +∞ . By construction one has that N Γ,k ≤ 2N Γ and therefore All these hyperplanes generate d-dimensional hyperrectangles {H l } N l=1 with the property that every H l has a associated target in α m(l) ∈ R d . Then, the number of final hyperrectangles N can be estimated by: Hence We consider sufficiently thin d-dimensional strip around each hyperplane {x (k) = c k,n } of side ζ = h d given by (5.9) S n,k := (R × R × · · · × [c n,k − ζ, c n,k + ζ] × · · · × R) ∩ Ω.
These strips, of measure of the order of h d , belong to Ω h together with Γ h (see Figure 5.3.). The total number of hyperplanes is of the order of h −(d−1) , (5.7), and therefore also the number of the thin strips. Hence, one can see that the measure of Ω h is of the order of h. After removing the d-dimensional strips, we are left with at most N connected sets.
Then, making an abuse of notation, we denote the number by N the number of connected components of Ω \ ∪ k,n S n,k and its connected components by {H l } N l=1 . Now the L 2 norm can be split in two: We see that it suffices to pick: to ensure that Step 2: Compression of Ω c . In this Step we apply a compression Lemma that allows us to apply the simultaneous controllability of Theorem 2. Given a set ω ⊂ R d , let us define diam(ω) := max x1,x2∈ω Furthermore, let us denote by (ω) (k) the k-th component of the set ω, i.e. ω (k) := {x ∈ R : ∃y ∈ ω such that y (k) = x}.
As done in Theorem 2, in step 1, since the backward NODE is a solution for the forward NODE for the appropiate controls, we can find an equivalent set of targets {α m } M m=1 such that α m if m = m . Without loss of generality, we will assume that the set of targets fulfills such contidition.
Let us summarize the properties we have built in the following (a) Consider hyperplanes {x (k) = c k,n } in the form of (5.6). (b) Let ζ > 0 be the width of the strips S n,k in (5.9). (c) Consider the sets {H l } N l=1 defined by the connected components of Ω \ ∪ n,k S n,k . (d) Assume that each set H l has a single associated target α m(l) .

Lemma 1 (Fundamental Compression Lemma).
Let Ω ⊂ R d be bounded, and consider hyperplanes, strips and sets as in (H). Assume that the set of targets For any T > 0 and any 0 < η < 1 /N small enough, there exist piecewise-constant controls A, W ∈ L ∞ ((0, T ); R d×d ) and b ∈ L ∞ ((0, T ); R d ) such that with K and C independent of N, N Γ and η. Furthermore, fixing T = 1, one has that: The number of discontinuities of the controls A, W, b is of the order of N.
Here, the subscript α denotes the dependence on the set {α m } M m=1 .
For the sake of readability, the proof is postponed to Appendix A. The proof uses a concatenation of flows as the ones described in Subsection 2.2 to achieve the desired result.
Note that, Lemma 1, aside of compressing, for η small enough, also contains a generalization of Theorem 1 for the hyperrectangles. 1 Step 3: Simultaneous Control. For η small enough, one can find a collection of M sets, Therefore we can apply Theorem 2. Let us denote the flow from Theorem 2 by ψ T . The controllability time of Theorem 2, T M does not depend on h and by Grönwall inequality to obtain that where C M is a constant that depends on M and on the set {α m } M m=1 but not on h. Hence ν ∼ M,α η. Furthermore, the extra number of switches will depend on M as in Theorem 2 but not on N .
We summarize the above discussion in a lemma that will be used in the next section.

Lemma 2.
Let Ω ⊂ R d be bounded, and consider hyperplanes, strips and sets as in (H). Let f be a simple function as in (5.1). For all ν > 0 small enough, there exist piecewise-constant controls A, W ∈ L ∞ ((0, T ); R d×d ) and b ∈ L ∞ ((0, T ); R d ) such that the associated flow of (2.3) satisfies Moreover, control norms are with K and C independent of N, N Γ and η. Furthermore, fixing T = 1, one has that: The number of discontinuities of the controls A, W, b is of the order of N.
Step 4: Estimates. Recalling (5.10a)(5.10b)(5.8) plus the requirement that η < N −1 applied to (5.12) we obtain that: The number of switches of the controls A, W, b is of the order of −2d(d−1) .

Remark 5.1 (General simple functions).
In the previous theorem we have assumed that the perimeter of the characteristic sets is finite. However, the same proof works for certain sets with infinite perimeter. Let D be the upper box-counting dimension of Γ (see [15,Chapter 2], see also [4, Chapter 1]). As in the proof above, let N Γ (h) be the number of hypercubes of side h needed to cover the boundary Γ, the upper box-counting dimension is defined as Hence, by definition one has that: Therefore, the estimates (5.3)(5.4) and (5.5) are replaced by The number of switches of A, W, b will be of the order of − 2dD d−D . Note that the estimates blow up when D = d, in such case, we cannot provide an estimate. Boundaries with dimension D = d exist, for instance, the Mandelbrot set [34]. Generally speaking, one can have a Jordan curve that has dimension equals to the ambient dimension [26].
Taking the box-counting dimension as a measure of complexity of the simple function we see that as the characteristic sets are more complex the cost of the control grows.

Control and simultaneous control of Neural transport equations
In this section we use the techniques of the proof of Theorem 4 to control a Neural transport equation (NTE) of the form where C c (R d , R + ) stands for the space of compactly supported continuous nonnegative functions. Our proof of the Universal Approximation Theorem can be interpreted in the context of the simultaneous control of these NTEs. This will also give rise to a transport formulation of the classification problem.
In this section the approximation will be measured in the sense of the Wasserstein-1 distance: is defined by as: where Lip(g) ≤ 1 stands for the class of Lipschitz functions with Lipschitz constant less or equal than 1.
Wasserstein distances play a central role in the theory of optimal transport (see [41,Chapter 5] for a general reference in this subject).
6.1. Control of the Neural Transport equation. . Since the vector field W σ(Ax+b) of the Neural Transport equation (6.1) is Lipschitz, it preserves the mass for each set H along the characteristic flow φ T (H) The control of the (6.1) can be done by controlling its characteristics plus a compression argument as done in Lemma 1 in Theorem 4.
We will consider target measures ρ * that are finite combinations of Dirac masses: Theorem 5. Let T > 0, d ≥ 2, σ be as (2.1) and ρ * as in (6.2) be a target positive measure satisfying: Then, for every > 0, there exist piecewise-constant control functions W, A ∈ L ∞ ((0, T ); R d×d ) and b ∈ L ∞ ((0, T ); R d ) such that the solution of (6.1) satisfies: where W 1 is the Wasserstein-1 distance. Furthermore, one has Proof. We proceed in several steps: (1) Assigning targets to the initial datum. Consider F : R → R be defined as: (1) .
Using the fact that ρ 0 ∈ C c (R d , R + ), there exist two real numbers, name them c 0 and c M , such that ∀s ≤ c 0 , F (s) = 0, ∀s ≥ c M , F (s) = 1.
Since ρ 0 ∈ C c (R d , R + ), F is continuous and monotone increasing, i.e. we can find M − 1 real numbers such that: Hence, we will divide ρ 0 with M − 1 hyperplanes of the form: Figure 6.1). The intention is to find a flow such that, by means of the characteristics, can transport all the mass β m contained between the hyperplanes {x (1) = c m } and {x (1) = c m−1 } to the target α m . As in the previous section, due to the continuity with respect of the initial data of the characteristics, this task will have to be carried out in an approximate manner. for some ζ > 0 in a way that: In this way, we ensure that the mass lying in the strips is of the order of . Using Lemma 2, we can ensure that the mass lying in the d-dimensional strips will lie in a bounded set independently of ζ and consequently independently of .
Let H ⊂ R d be a subset, and let ρ(t) be the solution of (6.1) for certain A, W, b, let φ T be the solution of (2.3) for the same controls. Due to the divergence structure of Equation (6.1) and that the field W σ(Ax + b) is Lipschitz, one has that the mass is preserved along the characteristics and . Now we may apply Lemma 2 to obtain the approximate controllability of all characteristics in every H m while having the mass allocated in the strips has been transported in a bounded set B(0, K) with K independent of ζ and ν. Moreover, choose ν small enough so that Now, the objective is to quantify, in terms of the Wasserstein distance, the proximity of ρ(T ) to ρ * . The conclusion will follow from a straightforward computation using the triangular inequality with a suitable ρ . We will find a ρ of the form Observe that the function g m : R + → [0, 1] defined as is an increasing function. Therefore, if the minimizer (6.4) is unique and if (6.4) satisfies ξ m < ν, then we have that g m (ξ m ) = β m . We, then, define By definition, β m > 0. Moreover, taking into consideration that the mass transported, through the application of Lemma 2, into B(α m , ν) is at least β m − 2 /M, one has that Now define the remainders of mass β m as follows (2) ...dx (d) dx (1) .
As in step 1 of this proof, there exist real numbers c 0 and c M such that: and as before, we find M − 1 real numbers such that Take a collection of {z m } M m=1 such that for every m ∈ {1, ..., M }, satisfies z m ∈ T m ∩ B(0, K). Then, choosing ν = the result follows by the triangular inequality.
(3) Estimates. The control cost, also depends on the initial density function ρ 0 . Define the map c : R + → R + as (1) .
Formally, the derivative of inverse map of c , which we denote by ζ c : R + → R + , at zero fulfills: Therefore, it suffices to require that ρ 0 (c m , x (2) , ..., x (d) )dx (2) ...dx (d) From Lemma 2 we obtain that: Furthermore, A L ∞ = 1, b L ∞ ≤ C and the number of discontinuities of the controls is bounded independently of but dependent on the target ρ * .
Remark 6.1 (More general target configuarations). We stated the theorem by setting targets that are Dirac masses, but, it is well known that, in particular, one can approximate any compactly supported probability measure by a finite number of Dirac masses.
Remark 6.2 (Influence of the initial condition on the control cost). The proof is valid for more general settings, namely for ρ 0 ∈ L 1 c (R d ; R + ), where by L 1 c (R d ; R + ) we understand the L 1 compactly supported nonnegative functions. However, the quantification of the control cost depends on the singularities of ρ 0 .
In (6.5) one can observe that as the density at a point c m increases, ζ should be smaller. If ρ 0 ∈ L 1 c (R d ) the cost will depend on the "strength" of the singularities of ρ 0 . The precise quantification of the control cost depending on the regularity of F (or singularity of ρ 0 ) will not be treated in this article. For more details and a fine analysis of singularities we refer to [20,Chapter I and II]. Remark 6.3. Due to the uniqueness of the solution of the underlying characteristic system, we cannot consider, in general, initial data that are Dirac deltas. More precisely, by the uniqueness of the characteristics, it is impossible to bring the initial datum ρ 0 = δ 0 to any configuration of the type ρ * = λδα 1 + (1 − λ)δα 2 for any λ ∈ (0, 1) and any α 1 , α 2 ∈ R d . Remark 6.4. In Theorems 4 and 5, to object of study was the input output map given by the flow φ T associated to (2.3). In Theorem 5 we use the same principles for another objective, to control in approximate manner from a given density function to another.
In Theorem 4 and Remark 5.1 we observed how the geometry of the supports plays a crucial role in the control cost. However, in Theorem 5, since we are controlling a single probability density, we do not observe such dependence. In contrast, in (6.5) we observe how high concentrations of mass can increase the cost of control.
Remark 6.5. The restriction of d ≥ 2 comes from the limitation pointed out in Remark 3.1. However, the transport equation with V as a control, can be controlled in the one-dimensional case approximately with respect to the Wasserstein-1 distance. It would be enough to design appropriate locations of attractors and repulsors depending on the mass distribution of ρ 0 . This would allow concentrating the mass in a finite number of points. Later, one needs to design a dynamical system that brings each mass approximately on the approximation of the target. Neural Transport Equations can achieve this in dimension d ≥ 2 just with the controls A, W and b. Remark 6.6. If the target Dirac masses have the same mass, one does not need the simultaneous controllability. Given M initial data {x m } M m=1 and M targets {α m } M m=1 , it is enough that for every > 0 one can find controls A, W and b such that Note that, one can also approximate any continuous compactly supported function by a finite number of Dirac masses whose support is disjoint and all of them have the same mass. Remark 6.7 (Other Wasserstein metrics). The definition we gave of the Wasserstein-1 metric relies on the Kantorovich-Rubinstein Theorem [40,Theorem 1.14]. In the Euclidean space, a Wasserstein-p distance is defined as: where Γ(µ 1 , µ 2 ) denotes the set of all measures γ that satisfy µ 1 = R d γ(x, y)dx and µ 2 = R d γ(x, y)dy. Theorem 5 states the result in the Wasserstein-1 metric, however, since the supports of the final datum ρ(T ) and the target are compact, one has equivalence of all W p distances [40, Chapter 7, Section 7.1.2]. Remark 6.8. For the ReLU case, and for any Lipschitz nonlinearity, fixing controls A, W and b, we see that the characteristics, at most, have an exponential decay. This is the impediment for Theorem 3 to be a universal approximation theorem for functions in L 2 (R d , R d ).
6.2. Simultaneous control of Neural transport equations and classification. .
The scalar transport equation does not allow distinguishing among labels. One should take a vectorial structure in order to have a transport formulation for classification.
Let us consider M classes, and M compactly supported probability densities ρ m for m = 1, ..., M . Moreover, we assume, as in Section 2, that for every x there exists, at most, a unique label y. This implies that the supports of the probability measures ρ m are disjoint: Theorem 4 is based on approximating a simple function. Therefore, in Theorem 4 we have already treated the problem of sending uniform disjoint compactly supported distributions to precisely located Dirac deltas (in an approximate manner). Combining the ideas of Theorem 4 and 5 one has the following: Let M ∈ N be a natural number and Ω ⊂ R d for d ≥ 2 and consider the system: Theorem 6. Let T > 0 and ρ * m be as in (6.2) the target measure for the the m-th equation. Assume that the initial conditions satisfy: .., M } where by Per we understand the perimeter. In addition, assume that the target functions satisfy: Then, for any > 0, there exist controls W, A ∈ L ∞ ((0, T ), R d×d ) and b ∈ L ∞ ((0, T ), R d ) such that the solution of (6.6) satisfies: where W 1 is the Wasserstein-1 distance.
The proof follows combining the arguments of Theorem 4 and 5. The only difference with respect to Theorem 5 is that all initial data have to be separated by several extra hyperrectangles as done in Step 1 in Theorem 4. This step was not needed for controlling just one transport equation, but it had to be done in the universal approximation Theorem 4. This implies that, in this case, the control cost will depend both, on the geometry of the supports of each density function, and on the concentration of mass of such density functions.
The classification interpretation is the following. Let S = {S m } M m=1 be a partition of R d and consider system (6.6). The classification problem consists of: for every > 0, being able to find controls (that depend on ) W, A ∈ L ∞ ((0, T ), R d×d ) and b ∈ L ∞ ((0, T ), R d ) such that the solution of (6.6) satisfies: Roughly speaking, the classification problem, in terms of transport equations, is equivalent to send most of the mass of each probability measure to prefixed sets and to control its expectation to the prefixed set.
Remark 6.9. Note that, since the supports are disjoint, if we consider ρ = N m=1 ρ m , the controls A, W and b of Theorem 6 control also the scalar equation (6.1). Remark 6.10 (Robustness with respect to large number finite dimensional samples). The controls obtained from Theorem 6 are robust with respect to large number finite dimensional samples that follow the probabilities distributions of the initial data of system (6.6). This implies that, given an error , the control cost of classification for finite dimensional samples as in Theorem 1 does not blow up when the cardinality of sample tends to infinity. Remark 6.11. Note that Remark 5.1 also applies in this case for having more irregular boundaries. However, if the boundaries of the supports of each distribution are at positive distance, one does not need to consider a fine meshing of the boundary.

Conclusions and Perspectives
We now present a number of conclusions and open problems.
(1) The role of the activation function. The key property that allows to prove all the results in this paper is that the activation functions are able to leave one halfspace invariant while moving the other half. As mentioned in the Introduction, the simultaneous controllability is a rare property for dynamical systems. In particular, it can never occur in linear systems when all trajectories solve the same system. Our methods exploit this nonlinear feature of activation functions to obtain simultaneous control properties of NODEs that lead to the desired results in the Machine Learning context.
(2) Algorithmic complexity. All proofs developed in this paper are algorithmic. The proofs in Theorem 1 and 2 require a number of steps of the order of O(N ), where N is the number of points to be classified (see Remarks 3.3 and 4.4). In the classification Theorem 1 each point to be classified is handled in an iterative manner. Thus, our estimate on the number of needed iterations corresponds to the worst-case scenario. In practice, the optimal number of required switches will be lower, depending on the structure of the data to be classified (the initial data of the Neural ODE) and also on the partition chosen (or target in the case of Theorem 2). In particular, when the initial data are structured into clusters our proof can take advantage of that to classify the data in fewer iterations, keeping them grouped in packages during the evolution.
Note that, even in the context of Theorem 1 and 2, where the number of switches is of order of O(N ), the control cost depends on the initial data and target destination. In particular, the cost of control increases when the distance between data corresponding to different labels decreases or their mixing increases.
The analysis is more intricate for the universal approximation theorem (Section 5) and for the control of Neural transport equations (Section 6). Even if the proofs are also algorithmic, the L ∞ -norm of the control depends on the precision required when approximating the target as well as on the complexity of initial and final data, see Remark 5.1.
(3) Optimal control strategies. As we have explained above, by time-scaling, all results are valid for an arbitrary T > 0. For the analysis of the complexity of the control dynamics and the cost of control it is convenient to normalize the time to T = 1. The complexity of the controls can be then estimated in terms of the number of switches and their T V -norm. The controls we build are not necessarily the optimal ones. In fact, as we have mentioned above, when the data to be classified present clustering phenomena, the complexity of the needed control diminishes.
Optimal controls could be defined setting the time horizon to be T = 1 and minimizing their L ∞ -norm. This would lead to an optimality system or Pontryagin maximum principle ( [30]) characterizing optimal controls that could be expected to present a bang-bang structure. The analysis of the complexity of optimal controls through such optimality system, which will strongly depend on the configuration and structure of the data to be classified, is an interesting and complex open problem. Of course such characterization does not lead to any explicit expression, it is of a purely implicit nature, making its posterior use rather complex ( [36]). Note that this issue is even more complex in the context of the ReLU activation function, because of its lack of regularity.
Other norms and cost criteria can also be used to define optimal controls. For instance, other than penalizing solely the norm of the control, one could also penalize the controlled trajectories, enhancing the stability and the turnpike properties of the control processes, as shown in [14] or in [43] where the L 1 penalization of the state gives rise to sparse bang-bang controls as the ones shown in this article.
In summary, the development of computational methodologies to derive control strategies of minimal complexity is an interesting and challenging topic.
(4) Optimal activation functions. In the previous point, we were proposing to find optimal controls with prescribed initial and final data. One general open question would be the following. How to construct a vector field such that the norm of the control or the number of switches is minimal? or more precisely, what kind of vector fields can classify with minimal cost for any target function? This question is analogous to the optimal observation and location of sensors in linear systems, addressed for instance in [31,32], but optimizing with respect to the nonlinearity. This question is not typically addressed in the control literature since, normally, in mechanics, we are not allowed to choose the dynamics that we want to control. On the other hand, the transport interpretation made in Section 6 allows for a more synthetic understanding of supervised learning. The time-continuous perspective can also be more adequate to analyze optimal control strategies.
Once the NODE setting is well understood, a posteriori, the time-discrete setting allows achieving similar results by means of fine discretizations of the Neural ODE, i.e., with sufficiently many deep layers.
On the other hand, as mentioned in Remark 3.1, when the time step is large enough, the discrete dynamics is richer in the sense that they are able to express features that the continuous NODEs cannot.
(7) On the topology approximation. Our proof of universal approximation guarantees approximation in L 2 . But it does not allow to obtain density results in spaces of higher regularity since it is fundamentally based on the approximation of simple functions. The approximation in W k,p , if possible, would require more sophisticated techniques.
(8) Other Neural ODEs. In our poofs we have considered controls W (t), A(t), b(t). But as we have seen, in each step of the control iteration only a few of the available control components were activated. This suggests that similar results can be achieved with fewer controls. For instance, • W = Id. We could consider the model on which the vector field σ is always pointing towards the first quadrant, which makes it impossible to achieve purposes such as the Universal Approximation (Section 5) or the approximate control of NTEs (Section 6). • We could also consider the varianṫ Our results could be extended for this type of Neural ODEs, using in an essential way the fact that a quadrant of the phase space remains invariant under the dynamics. However, this structure of the NODE introduces further rigidity in the dynamics and the proofs require more complex arguments. Another perspective is to generalize the results for second-order Neural ODEs such as continuous versions of Momentum Residual Neural Networks [33].
(9) Optimal Transport. Theorems 5 and 6 naturally establish a path towards optimal transport theory. Optimal transport can be formulated in a dynamic framework. Namely, given two positive measures ρ 0 and ρ 1 with the same mass, the goal is to find an optimal vector field V and a solution ρ of the continuity equation 3,7,27]). This can be interpreted as an optimal controllability problem, the vector field V being aimed to optimize the energy subject to (7.1). Here W 2 stands for the Wasserstein-2 distance. The control problem (7.1) is equivalent on finding geodesics in the measure space. Theorems 5 and 6 guarantee that Neural transport equations (6.1) and (6.6) allow to transport approximately ρ 0 into ρ 1 . Analyzing how close are the optimal flows of Neural transport equations from the actual optimal solution of the dynamic optimal transport is an interesting and challenging problem.
if m = m . For any T > 0 and any 0 < η < 1 /N small enough, there exist piecewise-constant controls A, W ∈ L ∞ ((0, T ); R d×d ) and b ∈ L ∞ ((0, T ); R d ) such that for every 1 ≤ l ≤ N, with K and C independent of N, N Γ and η. Furthermore, fixing T = 1, one has that: The number of discontinuities of the controls A, W, b is of the order of N.
Here, the subscript α denotes the dependence on the set {α m } M m=1 . Proof. In the whole proof we fix the L ∞ norm of W to be 1, and we study the time needed for each step to be completed. By rescaling, this will give us the L ∞ norm of W . 1. Compression of Ω c . In this step we compress Ω c . Notice that in all this process we will be transforming also Ω, being able to bound the transformation of Ω h whose diameter will remain very close to the diameter of Ω c . We will compress the sets H l so that: • All hyperrectangles, after applying the controls, will have a diameter of the order of η as small as we want. • The compressed sets will be separated by a fix constant independent of η. This will allow us to apply the simultaneous controllability after a refinement of Theorem 1 made in the next step of the proof.
There will be two phases, firstly we compress the x (k) -coordinates for k = 1, ..., d − 1. As we will see later on, with this compression the sets H l will have been transformed and moved in a particular configuration that will allow the the compression in the x (d) -coordinate. Moreover, while compressing the the x (d) -coordinate we will create a separation between the compressed sets.
For each k from 1 to d − 1 apply successively the following steps: 1.1 Compression of the x (k) -coordinates 1, ..., d − 1.
Here below we describe a procedure that will be recursively applied using successively all hyperplanes {x (k) = c k,n } starting from n = 1 until n = N Γ,k .
1.1.1 First Compression. We set the hyperplane {x (k) = c k,n } and the controls We choose W so that the field points towards the hyperplane, generating a contraction. We solve (2.3) for T large enough so that the hyperrectangles that have been affected by the nonzero vector field have a diam (k) lower than δ (see Figure A. By φ T (H l ) we understand the flow applied to the set H l until time T with the aforementioned controls. We make an abuse of notation, and we redefine H l := φ T (H l ). This will be done after each step in the proof to avoid the long notation of applying sequentially different flows with different controls. The time needed for making this compression is of the order of: We choose W so that the field points to −∞ in the x (d) -axis. We solve (2.3) for T large enough so that the sets that have been affected by the nonzero vectorfield, have lower values in the x (d) -coordinate than the sets that have not been affected by the vectorfield (in the region in which the field has been 0) (see Figure A.2), i.e. max x∈φ T (Hi) Given that the side of the removed strip is of the order of ζ, the field will be bounded by below by a quantity of the order of ζ. The space needed for the translation is a constannt depending on Ω. Therefore, the time needed for achieving the separation is of the order of (A.4) and W so that the field points towards the hyperplane, making it attractive (see Figure  A.3). We solve (2.3) for T large enough so that: The needed time is: As said before, this procedure is applied sequentially for every 1 ≤ k ≤ d − 1 and every 1 ≤ n ≤ N Γ,k − 1 obtaining that the resulting sets satisfy The set Ω has transformed, and diam (k) (φ T (Ω)) < δ for 1 ≤ k ≤ d − 1. Before continuing, we will "move the sets away from the targets". The reason of this particular action will be better understood in Step 2 of this proof.
Choose A and W so that the vector field, in its nonzero region in the phase space, follows the Cartesian direction of the first element of the canonical basis and makes the hyperplane (A.6) attractive. We solve (2.3) for T long enough so that: The overall required time of this step is the following where by ∼ α we denote the dependence on the targets. In the separation step, Step 1.1.2, the set deforms since the field is not homogeneous along the set. We must ensure that a hyperplane with normal vector equal to an element of the canonical basis can fit between the transformed sets. This is the main reason of the first compression step, Step 1.1.1. This step reduces the diam (k) of the set making that the field in the next step, Step 1.1.2 is almost constant on the hyperrectagle. Let us quantify the present discussion.
Before Step 1.1.2 the distance between two neighbor rectangles is ζ. After applying Step 1.1.2 (d − 1)N Γ times the separation is: so it will be enough if δ ≤ CN −1 Γ ζ for a certain constant C to guarantee that after the whole process the distance between two neighbor rectangles is of the order of ζ.
Notice that A has been chosen to be unitary, also b is bounded since we have chosen only among the original hyperplanes defined in Step 1. Furthermore, the number of switches we employed are explicit: The number of switches in the A control is d − 1.
(A.8) The number of switches in the b control and W control are of the order of N Γ .

1.2
Compression of the d-coordinate. Now, the task is to compress the x (d) -coordinate and separate enough each set while keeping a small diameter. After the transformations made above, the sets can be separated by N − 1 hyperplanes of the form For doing so, for l from 1 until N we proceed as follows: The diameter of φ T (H l ) is lower than diam(φ T (H l )) < √ dδ.

Separation.
We fix a parallel hyperplane to the one before. We set it slightly below at a distance of the order of the removed band, ζ, {x (d) = C d,l − ( 1 /4)ζ}, and we choose A ij = δ di δ dj , b i = (−C d,l + ( 1 /4)ζ)δ di and W such that the field points to the x (1) -coordinate. We solve (2.3) for T * large enough so that inf x∈φ(H l ) This would require T * = 8/(N ζ).
When separating the sets one should ensure that at the end of the process we stay in a bounded set independently of N . After applying this movement along the x (1) -coordinate, the diameter of the set H l has grown being: The time needed to perform this steps, the number of switches and the norm of b are estimated by the following: The number of switches of the A control has been 1. The number of switches of the b and W control is of the order of N .
We apply iteratively this process (see Figure A.6) for all hyperplanes {x (d) = C d,l } (obtaining at the end Figure A.7). Since the velocity field has different strength at different parts of the set, the distance between two H l can diminish.
In order to guarantee that after the whole process we have enough distance between the sets, we have to require δ to be small enough. When we apply the horizontal translation, the difference in the vector field inside each set is δ. Hence, we know that the distance in their x (1) -coordinate between the sets reduces as: If we impose that, after T * N units of time, the distance in the x (1) -coordinate is bigger than 1/N , this implies that the requirement is the following: (A.11) δ < ζ 8N The requirement (A.11) will be not the only one. Now the task is to guarantee that diam(φ T (Ω h )) is uniformly bounded for h small enough and that the maximum diameter of all the sets φ T (H l ) is as small as we desire.
Choosing any η such that (A.13) η < 1 N we would guarantee that φ T (Ω) remains bounded. With all the process described, we have moved the sets downwards in the x (d) -component, this translation depends on N . However, by putting an attractive hyperplane at the of the form x (d) = c at the desired place, we can move the set without having a dependence on N . Now the sets are at distances bigger than 1/N between them and each one of them has a bound in its diameter (A.12). We can choose δ as small as we want and we can guarantee that: Using the expressions (A.9) and (A), we find that the time horizon needed in Step 1 of this proof scales in the following manner: (A.14) T step 1 ∼ Ω,d,α N log The norm of the control b depends on N as seen in the treatment of the x (d) -coordinate (Step 1. Let φ T be the flow associated to the compression phase, and let Ω and H l be defined as in Step 1 of this proof. We redefine Ω and H l as Ω := φ T (Ω) and H l := φ T (H l ). Without loss of generality, we make a change of coordinates, for simplifying notation, in a way that: Remind that by (ω) (k) we denote the k-th component of the set ω, i.e. ω (k) := {x ∈ R : ∃y ∈ ω such that y (k) = x}.
Now we build a control that will control the first component of H l to their corresponding target while ensuring that after this control process, (A.18) still holds. This will allow us to proceed recursively. The procedure we are going to show is, in particular, an improvement of Theorem 1 for keeping track properly of Ω h while applying our controls.
x (2) = −1 x (1) = d l + 2η x (2) = −1 x (1) = d l + 2η x (2) = −1 x (1) = dl + 2η The time horizon will depend on the number of elements N . The mission is to make use of the explicit controlled flows to obtain sharper bounds and independent on the number of elements. Since the time horizon depends on N , the direct application of Grönwall's inequality is not of practical use for bounding Ω h independently of N . This is the reason why we examine the dynamics in more detail.
This Lemma that we are proving, leaves the system in a configuration in which the simultaneous controllability of Theorem 2 can be applied. In the application of Theorem 2, after applying Theorem 1, we only set controls that endow flows that follow directions in the orthogonal subspace of Span((1, 0, ..., 0)). Due to the Cartesian structure of our flows see that we only have to ensure that after the refinement of Theorem 1 the condition diam (k) (Ω) = O(η) for any k ≥ 2 still holds.
For every l ∈ {1, ..., N } consider the center of H l , d l = 1 |H l | H l xdx and we proceed sequentially.
(1) Observe that Ω is bounded in the second coordinate by: Fix the following bounds: x (1) = d (1) l + 2η (A.19a) x (1) = d Remark A.1 (The Bottleneck). The compression process done in the Lemma above is constituted by several stages, combining contractions and separations in different manner. However, there is a control stage which is giving the bound on the control cost which is the Ordering and Grouping step (Step 2 of the proof ). In producing certain flows, such as separations, NODEs require to have higher controls than when producing contractions.