# Interplay between depth and width for interpolation in neural ODEs

Antonio Álvarez-López<EMAIL_ADDRESS>Arselane Hadj Slimane
<EMAIL_ADDRESS>Enrique Zuazua
<EMAIL_ADDRESS>
###### Abstract

Neural ordinary differential equations (neural ODEs) have emerged as a natural
tool for supervised learning from a control perspective, yet a complete
understanding of their optimal architecture remains elusive. In this work, we
examine the interplay between their width $p$ and number of layer transitions
$L$ (effectively the depth $L+1$). Specifically, we assess the model
expressivity in terms of its capacity to interpolate either a finite dataset
$\mathcal{D}$ comprising $N$ pairs of points or two probability measures in
$\mathbb{R}^{d}$ within a Wasserstein error margin $\varepsilon>0$. Our
findings reveal a balancing trade-off between $p$ and $L$, with $L$ scaling as
$O(1+N/p)$ for dataset interpolation, and
$L=O\left(1+(p\varepsilon^{d})^{-1}\right)$ for measure interpolation.

In the autonomous case, where $L=0$, a separate study is required, which we
undertake focusing on dataset interpolation. We address the relaxed problem of
$\varepsilon$-approximate controllability and establish an error decay of
$\varepsilon\sim O(\log(p)p^{-1/d})$. This decay rate is a consequence of
applying a universal approximation theorem to a custom-built Lipschitz vector
field that interpolates $\mathcal{D}$. In the high-dimensional setting, we
further demonstrate that $p=O(N)$ neurons are likely sufficient to achieve
exact control.

###### keywords:

Neural ODEs , Depth , Width , Simultaneous controllability , Transport
equation , Wasserstein distance

###### MSC:

34H05 , 35Q49 , 68T07 , 93B05

[1]organization=Universidad Autónoma de Madrid, addressline=C. Francisco Tomás
y Valiente, 7, city=Madrid, postcode=28049, country=Spain [2]organization=ENS
Paris Saclay, addressline=Avenue des science, 4, city=Gif-sur-Yvette,
postcode=91190, country=France [3]organization=Chair for Dynamics, Control,
Machine Learning, and Numerics, Alexander von Humboldt-Professorship,
Department of Mathematics, Friedrich-Alexander-Universität Erlangen-Nürnberg,
addressline=Cauerstraße, 11, city=Erlangen, postcode=91058, country=Germany
[4]organization=Fundación Deusto, addressline=Av. de Universidades, 24,
city=Bilbao, postcode=48007, country=Spain

##  1 Introduction

Residual neural networks (ResNets) are formally defined as the family of
discrete systems

$\displaystyle\begin{cases}\mathbf{x}_{k+1}&=\mathbf{x}_{k}+W_{k}\boldsymbol{\sigma}\left(A_{k}\mathbf{x}_{k}+\mathbf{b}_{k}\right),\\\
\mathbf{x}_{0}&\in\mathbb{R}^{d},\end{cases}$ (1.1)

where $k=0,\ldots,L$, $W_{k}\in\mathbb{R}^{d\times p}$,
$A_{k}\in\mathbb{R}^{p\times d}$ and $\mathbf{b}_{k}\in\mathbb{R}^{p}$, for
some $d\geq\leavevmode\nobreak\ 1$, $L\geq 0$ and $p\geq 1$. Each time step
$k$ identifies a layer of the network. The number of layers $L+1$ is the
_depth_ of (1.1). The parameter $p$ is the _width_ of (1.1), identifying the
number of neurons per layer. The activation function
$\boldsymbol{\sigma}:\mathbb{R}^{p}\to\mathbb{R}^{p}$ is defined as the column
vector
$\boldsymbol{\sigma}(\mathbf{y})=\big{(}\sigma(y^{(1)}),\dots,\sigma(y^{(p)})\big{)}^{\top}$
from a chosen nonlinear function $\sigma:\leavevmode\nobreak\
\mathbb{R}\to\mathbb{R}$. We consider the Rectified Linear Unit (ReLU), given
by $\sigma(z)=\max\\{z,0\\},\;z\in\mathbb{R}.$

It has been noted [10, 17, 5, 4, 29] that (1.1) can be identified with the
forward Euler discretization scheme for the class of continuous models known
as neural ordinary differential equations (neural ODEs),

$\displaystyle\begin{cases}\dot{\mathbf{x}}(t)&=W(t)\boldsymbol{\sigma}\left(A(t)\,\mathbf{x}(t)+\mathbf{b}(t)\right),\\\
\mathbf{x}(0)&=\mathbf{x}_{0}\in\mathbb{R}^{d},\end{cases}$ (1.2)

where $(W,A,\mathbf{b})\in L^{\infty}\left((0,T),\mathbb{R}^{d\times
p}\times\mathbb{R}^{p\times d}\times\mathbb{R}^{p}\right)$ for some $T>0$.
Here, $t\in(0,T)$ parameterizes the evolution of the states through a
continuous range of layers. As discussed in [27, 2], it is common to assume
that $(W,A,\mathbf{b})$ is a step function over $(0,T)$, to align closer with
the dynamics of (1.1). Then, since $\sigma$ is Lipschitz, existence and
uniqueness of solutions hold for any $(W,A,\mathbf{b})$ and initial condition
$\mathbf{x}_{0}$. Equation (1.2) can be equivalently written as

$\dot{\mathbf{x}}=\sum_{i=1}^{p}\mathbf{w}_{i}(t)\sigma(\mathbf{a}_{i}(t)\cdot\mathbf{x}+b_{i}(t)),$
(1.3)

where $\mathbf{w}_{i}$ and $\mathbf{a}_{i}$ are respectively the $p$ columns
of $W$ and the $p$ rows of $A$, both seen as column vectors in
$\mathbb{R}^{d}$, while $b_{i}$ is the $i$-th coordinate of $\mathbf{b}$, for
$i=1,\dots,p$. In this work, we use formulation (1.3), although, for
simplicity, we represent $(\mathbf{w}_{i},\mathbf{a}_{i},b_{i})_{i=1}^{p}$ in
their matrix form $(W,A,\mathbf{b})$, which corresponds to the equivalent
system (1.2).

A Shallow ResNet

B Narrow ResNet

Figure 1: Qualitative representation of models (1.5) and (1.6) as discrete
systems. Blue circles represent the input $\mathbf{x}$; switches depict ReLU
functions; green circles indicate the result of
$W\boldsymbol{\sigma}(A\mathbf{x}+\mathbf{b})$; orange circles represent the
output after residual connections.

Equation (1.3) can be naturally extended to handle probability distributions,
rather than points in $\mathbb{R}^{d}$, by interpreting its right-hand side as
the advection field that drives the evolution of a measure $\rho$. This
extension gives rise to the neural transport equation [28]:

$\partial_{t}\rho+\operatorname{div}_{\mathbf{x}}\Big{(}\sum_{i=1}^{p}\mathbf{w}_{i}\,\sigma(\mathbf{a}_{i}\cdot\mathbf{x}+b_{i})\rho\Big{)}=0.$
(1.4)

Prior research indicates that control theory offers significant potential for
examining the properties of neural ODEs, for instance, via optimal control
[10, 13, 14] or geometric control techniques [31, 1, 30]. A fundamental
problem still open is to develop a comprehensive understanding of the roles
played by depth and width with respect to the expressive power of the model,
see [18, 23, 15]. This property is often evaluated by its capacity to
interpolate either a finite set of point pairs or two given probability
measures.

The first, commonly referred to as finite-sample expressivity [34], is
associated with the approximation power of the model, see [21]. It essentially
amounts to a simultaneous control problem, where the aim is to find a control
function $(W,A,\mathbf{b})$ such that the associated input-output map, given
by the flow of (1.3), maps $N$ specified points to their $N$ corresponding
target points in $\mathbb{R}^{d}$. Throughout this work, figures represent
each input point as a colored solid circle and its corresponding target point
as an empty circle of the same color.

In the second scenario, we aim to control the transport dynamics described by
(1.4) in order to transform a given initial density $\rho_{0}$ into another
density $\rho_{T}$. This task is highly relevant for probabilistic modeling or
the generation of synthetic data via normalizing flows [20, 25, 16]. We
approach it as an approximate control problem in the Wasserstein-$q$ metric
space for $q\geq 1$, extending prior work focused on $W_{2}$, see [11], or in
$W_{1}$ with $p=1$, see [27].

Our main objective is the development of a comprehensive theory of
interpolation for the family of models described by (1.3), linking the error
to the specific architecture given by $p$ and $L$. Both numerical and
theoretical studies [19, 12, 24] suggest that networks with greater depth
usually achieve better performance. This tendency is particularly noticeable
in training [35, 34], which reinforces the intuition that a deeper network
should possess greater expressivity, i.e., an enhanced ability to learn more
complex non-linear functions. Understanding the balance between width and
depth is thus vital for the optimal design of networks. We tackle this
significant question using the continuous framework of neural ODEs, where
depth is expressed as $L+1$, $L$ being the number of time discontinuities of
the control $(W,A,\mathbf{b})$. As we vary $L$ and $p$, two limiting models
emerge.
Shallow neural ODEs. $L=0$ is fixed, while the width $p$ can be as large as
required:

$\dot{\mathbf{x}}=\sum_{i=1}^{p}\mathbf{w}_{i}\sigma\left(\mathbf{a}_{i}\cdot\mathbf{x}+b_{i}\right),$
(1.5)

where
$\\{(\mathbf{w}_{i},\mathbf{a}_{i},b_{i})\\}_{i=1}^{p}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}$
are constant controls, making the equation autonomous. The field on the right-
hand side of (1.5) corresponds to a one hidden layer neural network with $d$
components. The approximation capacity of this class of functions has been
extensively studied (see [6, 26, 7]). The discrete version of (1.5) can be
identified with a one hidden layer ResNet (see fig. 1A).
Narrow neural ODEs. $p=1$ is fixed, while the depth $L+1$ can be as large as
required:

$\dot{\mathbf{x}}=\mathbf{w}(t)\sigma\left(\mathbf{a}(t)\cdot\mathbf{x}+b(t)\right),$
(1.6)

where $(\mathbf{w},\mathbf{a},b)\in
L^{\infty}\left((0,T),\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\right)$.
The ability of this model to interpolate data and approximate functions has
been explored in [27, 22]. It offers the advantage of easier construction of
explicit controls compared to (1.5), owing to its simplified dynamics, albeit
at the expense of increased depth, which scales with the cardinal $N$ of the
dataset. The discrete version of (1.6) corresponds to a deep ResNet with one
neuron per hidden layer, so it alternates layers of dimension 1 and $d$ (see
fig. 1A).

Developing a unified theory that bridges shallow and narrow neural ODEs would
combine the vast work done for (1.5) with the intuitive dynamics of (1.6).
Moreover, it would facilitate the optimal design of a neural ODE through the
strategic choice of depth and width. This entails optimizing the complexity
$\kappa$, defined as the total number of parameters in (1.3):

$\kappa\coloneqq(L+1)\times p\times(2d+1).$ (1.7)

Indeed, on each of the $L+1$ hidden layers, $p$ neurons of dimension $2d+1$
need to be determined.

###  1.1 Roadmap

In section 2, we present the main results of our work in two parts. First, in
section 2.1, we study the problem of interpolating a finite dataset in (1.3),
which is recast as the property of simultaneous control for neural ODEs.
Second, in section 2.2, we approach approximate controllability of probability
measures using the dynamics provided by the neural transport equation (1.4).
In section 3, we discuss the main implications of our work and pose some open
questions. In section 4, we prove the main results and provide the necessary
tools as lemmas.

###  1.2 Notation

  * 1.

We use subscripts to identify the particular elements from a dataset and
superscripts for the coordinates of a vector. In addition, (column) vectors
are denoted with bold letters and matrices with capital letters.

  * 2.

We denote by $\mathbf{x}\cdot\mathbf{y}$ the scalar product of
$\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$.

  * 3.

We denote by $\lceil z\rceil$ the lowest integer greater than or equal to
$z\in\mathbb{R}$, and by $\lfloor z\rfloor$ the highest integer lower than or
equal to $z$.

  * 4.

We denote by $\mathbb{S}^{d-1}$ the $(d-1)$-dimensional sphere in
$\mathbb{R}^{d}$.

  * 5.

We denote by $\operatorname{Lip}\left(\mathbb{R}^{d},\mathbb{R}^{d}\right)$
the space of Lipschitz-continuous vector fields in the usual norm, and by
$L_{V}$ the Lipschitz constant of each
$\mathbf{V}\in\operatorname{Lip}\left(\mathbb{R}^{d},\mathbb{R}^{d}\right)$.

  * 6.

We denote by $\mathcal{B}(\mathbb{R}^{d})$ the space of Borel measurable sets
in $\mathbb{R}^{d}$, and by $\mathcal{M}_{b}(\mathbb{R}^{d})$ the space of
Borel measures. Moreover, $\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ denotes the
set of compactly supported and absolutely continuous Borel probability
measures.

  * 7.

For any $\mu\in\mathcal{M}_{b}(\mathbb{R}^{d})$ and any measurable function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$,
$f_{\\#}\mu\in\mathcal{M}_{b}(\mathbb{R}^{d})$ denotes the pushforward measure
of $\mu$, defined for all $A\in\mathcal{B}(\mathbb{R}^{d})$ by

$f_{\\#}\mu(A)=\mu(f^{-1}(A))\in\mathcal{B}(\mathbb{R}^{d}).$

##  2 Main results

###  2.1 Simultaneous control

Let $N\geq 1$, $d\geq 2$, and consider a dataset

$\mathcal{D}=\\{(\mathbf{x}_{n},\mathbf{y}_{n})\\}_{n=1}^{N}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$
(2.1)

with $\mathbf{x}_{n}\neq\mathbf{x}_{m}$ and $\mathbf{y}_{n}\neq\mathbf{y}_{m}$
for all $n\neq m$. First, we study the finite-sample expressivity of the
general model (1.3), recast as a problem of simultaneous control.

######  Problem-Definition.

For any fixed time horizon $T>0$, find controls

$\left\\{(\mathbf{w}_{i},\mathbf{a}_{i},b_{i})\right\\}_{i=1}^{p}\subset
L^{\infty}\left((0,T);\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\right),$

for some $p\geq 1$, such that the flow $\Phi_{T}(\cdot;W,A,\mathbf{b})$ of
(1.2) _interpolates_ the dataset $\mathcal{D}$, i.e., it simultaneously drives
each data point from its initial position $\mathbf{x}_{n}$ to its target
$\mathbf{y}_{n}$. This is fulfilled when

$\Phi_{T}(\mathbf{x}_{n};W,A,\mathbf{b})=\mathbf{y}_{n}\qquad\text{ for all
}n=1,\dots,N,$

where $W,A,\mathbf{b}$ are respectively the matrix with columns
$\mathbf{w}_{i}$, the matrix with rows $\mathbf{a}_{i}$ and the vector with
components $b_{i}$, for $i=1,\dots,p.$

Our first result provides a relationship between $L$ and $p$ that ensures
interpolation of $\mathcal{D}$:

######  Theorem 1.

Let $N\geq 1$, $d\geq 2$ and $T>0$ be fixed. Consider the dataset
$\mathcal{D}$ as defined in (2.1). For any $p\geq 1$, there exists a piecewise
constant control

$\left(W,A,\mathbf{b}\right)\in L^{\infty}\left((0,T);\mathbb{R}^{p\times
d}\times\mathbb{R}^{p\times d}\times\mathbb{R}^{p}\right)$

such that the flow $\Phi_{T}(\cdot;W,A,\mathbf{b})$ generated by (1.3)
interpolates the dataset $\mathcal{D}$, i.e.,

$\Phi_{T}(\mathbf{x}_{n};W,A,\mathbf{b})=\mathbf{y}_{n},\qquad\text{for all
}n=1,\dots,N.$

Furthermore, the number of discontinuities of $(W,A,\mathbf{b})$ is

$L=2\left\lceil N/p\right\rceil-1.$ (2.2)

######  Remark 1.

If the target points $\\{\mathbf{y}_{n}\\}_{n=1}^{N}$ in (2.1) are not
distinct, interpolation is not achievable due to the uniqueness of solutions
in the system (1.3). In such cases, we relax the statement from exact to
approximate controllability by applying theorem 1 to an
$\varepsilon$-perturbation of the targets, for some $\varepsilon>0$.

Let us briefly describe the algorithm. First, we pivot around the
$x^{(1)}$-coordinate and controls the remaining $d-1$ coordinates.
Consequently, the trajectory of each data point $\mathbf{x}_{n}$ is confined
within the hyperplane defined by the equation $x_{n}^{(1)}=x^{(1)}_{n}$. Then,
we pivot using the controlled coordinates to adjust $x^{(1)}$. This algorithm
requires a depth of $2\lceil N/p\rceil$ layers, which is independent of the
dimension $d$, since a constant control suffices to simultaneously steer $d-1$
coordinates in the first step, assuming that $p>N$. This approach is broadly
applicable to any activation function, provided it meets the following three
conditions:

$\text{1. }\sigma\text{ is Lipschitz;}\quad\text{2. }\sigma(z)|_{z\leq
0}=0;\quad\text{3. }\sigma(z)|_{z>0}>0.$

Although this generalization holds theoretical interest, the ReLU function
appears to be superior due to its ability to work effectively with parameters
of smaller magnitude.

In (2.2) we can see that, as the width $p$ increases, the number of
discontinuities $L$ decreases with the same rate, meaning that width and depth
play a similar role in the steering. Nevertheless, a result on the optimal
design of our interpolating models can be derived:

######  Corollary 2.

For the family of controls defined in the proof of theorem 1 that ensure
interpolation of $\mathcal{D}$, the minimal value of complexity is

$\kappa_{\text{min}}=(4d+2)(N+1),$

obtained when $p=1$, i.e., when the neural ODE belongs to the narrow model
(1.6).

The complete transition from the narrow model (1.6) to the shallow model
(1.5), characterized by $L=0$ is not attained in (2.2). Due to the division
into two steps in the proposed algorithm, whenever $p>N$ the selected control
will exhibit a single switch ($L=1$), reaching a two-layer architecture,
rather than the autonomous ansatz (1.5). The restriction naturally raises the
question of whether simultaneous control is possible in shallow neural ODEs
(1.5). For this task, a reconsideration of the algorithm presented in [27]
becomes necessary. In the high dimensional setting, and more precisely, when
the dimension exceeds the number of data points ($d>N$), we can refine the
statement of theorem 1 to include the case $L=0$:

######  Corollary 3.

Let $N\geq 1$, $d\geq 2$ with $d>N$, and $T>0$ be fixed. Consider the dataset
$\mathcal{D}$ as defined in (2.1). For any $p\geq 1$, there exists a piecewise
constant control

$\left(W,A,\mathbf{b}\right)\in L^{\infty}\left((0,T);\mathbb{R}^{d\times
p}\times\mathbb{R}^{p\times d}\times\mathbb{R}^{p}\right)$

such that the flow $\Phi_{T}(\cdot;W,A,\mathbf{b})$ generated by (1.3)
interpolates the dataset $\mathcal{D}$. Furthermore, the number of
discontinuities of $(W,a,\mathbf{b})$ is

$L=2\big{(}\left\lceil N/p\right\rceil-1\big{)}.$

The key idea is that, when $d>N$, the first step in the proof of theorem 1 can
be suppressed. This is done by transforming the $x^{(1)}-$axis so that each
$\mathbf{x}_{n}$ shares the same first coordinate with $\mathbf{y}_{n}$, or,
equivalently, by optimally reorienting the hyperplanes represented in fig. 5.
For more insights on the proof, see fig. 6 in section 4.

######  Remark 2.

Corollary 3 suggests ideas similar to those in [13, Theorem 5.1]. In that
result, interpolation is established for $d\geq N$ in a simplified neural ODE,
when $\sup_{n=1,\dots,N}|\mathbf{x}_{n}-\mathbf{y}_{n}|<\varepsilon$ for a
sufficiently small $\varepsilon>0$, under a geometric assumption on the images
of the targets through $\boldsymbol{\sigma}$. Moreover, an estimation of the
control cost is obtained, which is linear with respect to $\varepsilon$. The
generation of new synthetic coordinates until $d\geq N$ is not typically a
problem, as discussed in [8], where the technique of embedding the dataset in
$\mathbb{R}^{d}\times\\{0,\dots,0\\}$ is proposed and its computational
advantages are studied.

In practice, $N$ tends to be larger than $d$. In that case, interpolation with
constant controls can be obtained for $p=N$ under a certain separability
hypothesis on $\mathcal{D}$:

######  Assumption 1.

Let
$\mathcal{D}=\\{(\mathbf{x}_{n},\mathbf{y}_{n})\\}_{n=1}^{N}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$
as defined in (2.1). There exist a vector $\mathbf{a}\in\mathbb{S}^{d-1}$, a
permutation $\tau$ of $N$ elements and a sequence
$-\infty<b_{N+1}<b_{N}<\cdots<b_{1}<\infty$ such that, for all
$n=1,\dots,N-1$,

$-b_{n}<\mathbf{a}\cdot\mathbf{x}_{\tau(n)}<-b_{n+1}\quad\text{and}\quad-
b_{n}<\mathbf{a}\cdot\mathbf{y}_{\tau(n)}<-b_{n+1}.$

1 claims that we can diagonally separate each pair
$(\mathbf{x}_{n},\mathbf{y}_{n})$ from the rest, in the sense that we can
define $N+1$ parallel hyperplanes
$H_{n}=\left\\{\mathbf{a}\cdot\mathbf{x}+b_{n}=0\right\\}$ such that the strip
$S_{n}$ bounded by $H_{n}$ and $H_{n+1}$ contains only the point
$\mathbf{x}_{n}$ and its target $\mathbf{y}_{n}$, for $n=1,\dots,N$ (see fig.
2). While the hypothesis might seem overly restrictive, it is noteworthy that
if the points are randomly sampled from a compact set, the probability that
the condition is fulfilled converges to 1 when the dimension grows:

######  Proposition 4.

Let $\mu\in\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ such that the random
variables $\pi_{i}X$ are independent and identically distributed (i.i.d.) for
$i=1,\ldots,d$, where $X\sim\mu$ and $\pi_{i}$ is the canonical projection on
the $i$-th coordinate. If every $\mathbf{x}_{n}$ and $\mathbf{y}_{n}$ in
$\mathcal{D}$ is sampled from $\mu$, and $N$ is sufficiently large, then the
probability $P$ that 1 is satisfied is bounded as

$1-\left[1-\frac{1}{\sqrt{2}}\left(\frac{e}{2N}\right)^{N}\right]^{d}\leq
P\leq 1.$

The hypothesis of proposition 4 are fulfilled by the uniform probability
measure in any hypercube, or by any isotropic Gaussian distribution.

Now, under 1, we can build a constant control such that the flow of (1.5),
taking a width $p=N$, interpolates $\mathcal{D}$. Our result is somehow a
dynamic version of [35, Theorem 1], under certain geometric conditions. A
representation of the trajectories can be seen in fig. 2.

######  Corollary 5.

Consider a dataset $\mathcal{D}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$ for
$d\geq 2$, under 1. For any fixed $T>0$, there exists a control
$(W,A,\mathbf{b})\in\mathbb{R}^{d\times N}\times\mathbb{R}^{N\times
d}\times\mathbb{R}^{N}$ such that the flow $\Phi_{T}$ generated by (1.5)
interpolates the dataset $\mathcal{D}$.

Figure 2: Left: separability condition in 1, for $\mathbf{a}=\mathbf{e}_{1}$.
Right: trajectories for exact control in the same example.

All in all, new strategies are required to study simultaneous control in the
autonomous model (1.5) under general conditions. A natural starting point to
assess the problem’s feasibility is to relax it by admitting an error
$\varepsilon>0$, and use density tools provided by universal approximation
theorems (UATs); see [6, 26]. In deep learning, UATs establish the density of
neural networks in function spaces over compact domains. The decay rate of the
error in relation to the number of parameters of the network has been
quantified for certain spaces [3, 7]. Specifically, these studies bound the
uniform error decay rate when the target function is Lipschitz continuous in a
compact domain.

We will approach the UAT, often interpreted in a static manner, from our
dynamic control perspective. In this regard, shallow neural ODEs provide a
vector field that transitions initial data to final data, with its flow at
time $T$ approximating the target function. First, we establish the existence
of a time-independent field whose integral curves guide each input point
$\mathbf{x}_{n}$ in $\mathcal{D}$ to its corresponding target $\mathbf{y}_{n}$
within a fixed time $T$. This field is constructed based on purely geometric
considerations (see fig. 3) and can be chosen to be Lipschitz continuous.

######  Proposition 6.

Let $N\geq 1$, $d\geq 2$ and $T>0$ be fixed. Consider the dataset
$\mathcal{D}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$ as defined in (2.1),
and any compact subset $\Omega\subset\mathbb{R}^{d}$ such that
$\operatorname{Int}(\Omega)$ is connected and
$\mathcal{D}\subset\operatorname{Int}(\Omega)\times\operatorname{Int}(\Omega).$
Then, there exists a vector field
$\mathbf{V}\in\operatorname{Lip}\left(\mathbb{R}^{d},\mathbb{R}^{d}\right)$
such that the flow $\Psi_{T,\mathbf{V}}$ of the equation

$\dot{\mathbf{x}}=\mathbf{V}(\mathbf{x})$

interpolates the dataset $\mathcal{D}$, and the $N$ curves given by

$\mathcal{C}_{n}\coloneqq\\{\Phi_{t,\mathbf{V}}(\mathbf{x}_{n}):t\in[0,T]\\},$

for $n=1,\dots,N$, are contained in $\operatorname{Int}(\Omega)$.

The subset $\Omega$, which will serve as our approximation domain, can always
be established as $\Omega=[-R,R]^{d}$ for a sufficiently large $R>0$.
Consequently, for any dataset
$\mathcal{D}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$, there exists a field
$\mathbf{V}\in\operatorname{Lip}(\mathbb{R}^{d},\mathbb{R}^{d})$ whose
integral curves $\mathcal{C}_{n}$ interpolate $\mathcal{D}$, i.e., the space

$\mathcal{V}_{\mathcal{D}}\coloneqq\big{\\{}\mathbf{V}\in\operatorname{Lip}\left(\mathbb{R}^{d},\mathbb{R}^{d}\right):\Psi_{T,\mathbf{V}}\text{
interpolates }\mathcal{D}\big{\\}}$

is non-empty. Moreover, we can define

$L_{0}\coloneqq\inf_{\mathbf{V}\in\mathcal{V}_{\mathcal{D}}}L_{V}.$

It suffices then to combine the UAT from [7] (see lemma 12 in section 4) with
classical results on the stability of ODEs to obtain the following theorem:

Figure 3: Construction of the Lipschitz field $\mathbf{V}$ in proposition 6
which interpolates $\mathcal{D}$ in a compact domain $\Omega$ that contains
all the points and curves.

######  Theorem 7.

Let $N\geq 1$, $d\geq 2$ and $T>0$ be fixed. Consider the dataset
$\mathcal{D}$ as defined in (2.1). For each $p\geq 1$, there exists a control
$(W,A,\mathbf{b})\in\mathbb{R}^{d\times p}\times\mathbb{R}^{p\times
d}\times\mathbb{R}^{p}$ such that the flow $\Phi_{T}$ generated by (1.5)
satisfies

$\sup_{i=1,\dots,N}|\mathbf{y}_{n}-\Phi_{T}(\mathbf{x}_{n};W,A,\mathbf{b})|\leq
C_{d,\mathcal{L},T}\,\frac{\log_{2}(\kappa)}{\kappa^{1/d}},$ (2.3)

where $\kappa=(d+2)dp$ is the complexity of the NODE, and

$C_{d,\mathcal{L},T}=C_{d,L_{0}}T\,\exp\big{(}\mathcal{L}\,T\big{)},$

being $\mathcal{L}=\min\\{L_{0},L_{NN}\\}$, $L_{NN}$ the Lispchitz constant of
the neural vector field, and $C_{d,L_{0}}>0$ a constant depending on $d$ and
$L_{0}$ but independent of $\kappa$.

######  Remark 3.

The argument employed in this theorem extends beyond neural networks. Since we
solely rely on a density result that provides a convergence rate, other dense
families of functions like polynomials, trigonometric, finite element methods
or wavelets could also be considered, with their corresponding convergence
rates.

######  Remark 4.

Given a domain of approximation $\Omega$, for the bound (2.3) to be optimal,
it is natural to pose the problem of finding the interpolating field
$\mathbf{V}\in\mathcal{V}_{\mathcal{D}}$ which has the smallest possible
Lipschitz constant $L_{V}$ within $\Omega$.

######  Remark 5.

When $d\geq 3$, the construction of a field
$\mathbf{V}\in\mathcal{V}_{\mathcal{D}}$ is generally very simple. Since, in
that case, two arbitrary curves are unlikely to intersect, we can generally
consider the $N$ segments that connect each pair
$(\mathbf{x}_{n},\mathbf{y}_{n})\in\mathcal{D}$ and build $\mathbf{V}$ as one
of the piecewise constant fields having these segments as integral curves.
Selecting the optimal field then becomes a combinatorial problem.

###  2.2 Transport control

So far, we have considered the system (1.3) with a finite set of points
$\mathcal{D}$ as initial data. A natural extension of this setting,
particularly pertinent when the data points are sampled from an underlying
distribution, is to consider as input a probability measure $\mu_{0}$ on
$\mathbb{R}^{d}$. The scenario where this distribution is a finite combination
of Dirac deltas aligns with the study previously conducted in section 2.1.

Consequently, our goal is to study the transformation of $\mu_{0}$ into a
fixed target probability measure $\mu_{*}$ through the push-forward map
generated by a neural ODE, i.e.,

$\Phi_{T}(\cdot;W,A,\mathbf{b})_{\\#}\mu_{0}=\mu_{*}.$

This question can be recast as the control problem of a transport equation.
Consider the family of probability measures $\mu(t)=\Phi_{t\\#}\mu_{0}$, where
$\Phi_{t}$ is the flow in time $t$ generated by (1.3). Since the neural vector
field is Lipschitz, if $\mu_{0}$ is absolutely continuous, then so is the
measure $\mu(t)$ for each $t>0$, and the evolving density $\rho(t)$ satisfies
the neural transport equation

$\displaystyle\begin{cases}\partial_{t}\rho+\operatorname{div}_{\mathbf{x}}\left(\rho\sum_{i=1}^{p}\mathbf{w}_{i}\,\sigma(\mathbf{a}_{i}\cdot\mathbf{x}+b_{i})\right)=0\\\
\rho(0)=\rho_{0}.\end{cases}$ (2.4)

Here, we have assumed that $\mu_{0}$ has density $\rho_{0}$, and

$(\mathbf{w}_{i},\mathbf{a}_{i},b_{i})_{i=1}^{p}\subset
L^{\infty}\left((0,T);\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\right)$

serve again as control functions. The projected characteristics of (2.4) solve
the neural ODE (1.3) in $(0,T)\times\mathbb{R}^{d}$. If the controls are step
functions, and since the ReLU function is Lipschitz, the continuity equation
(2.4) is well-posed and the total mass is conserved. Therefore, we aim to find
some controls such that the corresponding solution of (2.4) with initial
condition $\rho_{0}$ satisfies

$\rho(T)=\rho_{*}.$

This task, however, can be very hard to achieve, so we consider a relaxation
of the problem to approximate control of (2.4). For this purpose, first we
must choose a function to quantify the difference between any two measures.

######  Definition 6.

For any $q\geq 1$, the _Wasserstein-q distance_ between
$\mu,\nu\in\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ is defined as

$W_{q}(\mu,\nu)\coloneqq\Big{(}\min_{\gamma\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}|\mathbf{x}-\mathbf{y}|^{q}d\gamma(x,y)\Big{)}^{1/q},$
(2.5)

where $\Pi(\mu,\nu)$ denotes the set of measures $\gamma$ on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ that couple $\mu$ and $\nu$ in the sense
that $\gamma(\cdot\times\mathbb{R}^{d})=\mu(\cdot)$ and
$\gamma(\mathbb{R}^{d}\times\cdot)=\nu(\cdot)$. Note that
$\mu\in\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ has finite $q$-th momentum for
every $q\geq 1$, hence the Wasserstein$-q$ distance is well-defined in this
space.

Moreover, recalling the Monge formulation of optimal transport, when $\mu$ and
$\nu$ belong to $\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ we have:

$W_{q}(\mu,\nu)=\Big{(}\min_{T}\Big{\\{}\int_{\mathbb{R}^{d}}|\mathbf{x}-T(\mathbf{x})|^{q}d\mu:T_{\\#}\mu=\nu\Big{\\}}\Big{)}^{1/q},$
(2.6)

where $T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ measurable, see [33].

######  Problem-Definition.

Let $\mu_{0}$ and $\mu_{*}$ be two compactly supported, absolutely continuous
probability measures with respective densities $\rho_{0}$ and $\rho_{*}$. For
any fixed time horizon $T>0$ and $\varepsilon>0$, find controls

$\left\\{(\mathbf{w}_{i},\mathbf{a}_{i},b_{i})\right\\}_{i=1}^{p}\subset
L^{\infty}\left((0,T);\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}\right),$

for some $p\geq 1$, such that the solution of (2.4) in time $T$ _approximately
interpolates_ the initial condition $\rho_{0}$ to the target density
$\rho_{*}$. This is achieved when the $W_{q}$-error of the corresponding
measures (for some $q\geq 1$) satisfies:

$W_{q}(\mu(T),\mu_{*})<\varepsilon.$

The following theorem offers a partial solution to this problem. It assumes
that $1\leq q<\frac{d}{d-1}$ and targets the uniform measure in $[0,1]^{d}$.
While this bears resemblance to achieving null controllability, the nonlinear
nature of the problem prevents from directly extending this result to
arbitrary targets.

######  Theorem 8.

Let $d\geq 1$, $\mu_{0}\in\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ with density
$\rho_{0}$, $\mu_{*}$ the uniform measure in $[0,1]^{d}$, and $T>0$ be fixed.
For any $\varepsilon>0$, $q\in[1,\frac{d}{d-1})$ and $p\geq 1$, there exists a
piecewise constant control

$\left(W,A,\mathbf{b}\right)\in L^{\infty}\left((0,T);\mathbb{R}^{d\times
p}\times\mathbb{R}^{p\times d}\times\mathbb{R}^{p}\right)$

such that the solution $\rho(t)$ of the transport equation (1.4), taking
$\rho_{0}$ as initial condition, satisfies

$W_{q}(\mu(T),\mu_{*})<\varepsilon.$

Furthermore, the number of discontinuities of $(W,A,\mathbf{b})$ is

$L=\left\lceil 2d/p\right\rceil+\max\\{\lceil n/p_{1}\rceil,\dots,\lceil
n^{d}/p_{d}\rceil\\}-1,$

for any $p_{1},\dots,p_{d}\geq 1$ such that $p_{1}+\cdots+p_{d}=p$, and

$n\coloneqq\left(\frac{3d^{1/2+1/q}}{\varepsilon}\right)^{\frac{1}{1+d/q-d}}.$

For a given $\varepsilon>0$, the behavior of $L$ resembles that described in
theorem 1, as it decreases with an increase in $p$, reaching $L=1$ when $p$ is
large enough. Our proof is similar to a strategy from [9], and based on the
specific movements that the neural ode (1.3) allows us to generate. We
compress the support of $\mu_{0}$ to $[0,1]^{d}$ and divide it into
hyperrectangles, each with a mass of $O(\varepsilon^{d})$. These subsets are
then transformed to match a similar partition of $[0,1]^{d}$ corresponding to
the uniform measure $\mu_{*}$.

######  Remark 7.

If $\varepsilon>0$ is sufficiently small, and we choose
$p_{1}=\cdots=p_{d-1}=1$, $p_{d}=p-d+1$, it follows that

$L=\lceil
2d/p\rceil+\left\lceil\frac{1}{p-d+1}\left(\frac{3^{1+d/q}\sqrt{d}}{\varepsilon}\right)^{\frac{d}{1+d/q-d}}\right\rceil-1.$

For $q=1$, this expression simplifies to:

$L=\left\lceil
2d/p\right\rceil+\left\lceil\frac{1}{p-d+1}\left(\frac{3^{1+d}\sqrt{d}}{\varepsilon}\right)^{d}\right\rceil-1.$

##  3 Discussion

###  3.1 Conclusions

We have established several results on the capacity of neural ODEs for
interpolation and its relationship with the chosen architecture, determined by
the depth $p$ and width $L$. More precisely, we have provided explicit
dependencies between these two parameters that are sufficient to (exactly or
approximately) interpolate either two sets of $N$ different points in
$\mathbb{R}^{d}$ or any compactly supported, absolutely continuous probability
measure in $\mathbb{R}^{d}$ with the uniform measure in $[0,1]^{d}$. Our work
reveals that $p$ and $L$ can play similar roles in the algorithms, thereby
exhibiting a degree of exchangeability in the network’s structure.

Specifically, theorem 1 proves that a neural ODE with $p$ neurons can
interpolate any dataset of $N$ pairs of points using a piecewise constant
control with $L=2\left\lceil N/p\right\rceil-1$ discontinuities. Although
increasing $p$ reduces the number of discontinuities, we find a limiting case
of a 2-hidden layer neural ODE when $p\geq N$. Explicit controls for
interpolation with a shallow neural ODE $(L=0)$ are obtained in corollary 3
when $d>N$; or in corollary 5, with $p=N$ under 1. More generally, theorem 7
provides an error decay rate with respect to the number of parameters for
shallow neural ODEs. Finally, theorem 8 explores the Wasserstein-$q$
approximate control of the neural transport equation to a uniform distribution
on $[0,1]^{d}$ using piecewise constant controls. As in theorem 1, the number
of discontinuities diminishes as $p$ increases.

###  3.2 Open questions

Some new objectives can be derived from our work:
1\. Approaching the autonomous regime. As we have discussed, both $d>N$ and 1
are only special cases where we can find controls to interpolate in the
autonomous regime of shallow neural ODE. The question of finding such a
construction for any $d\geq 1$, or at least under a less restrictive
hypothesis on the dataset than 1, is still open. A first step could involve
assuming the relaxed condition that the projections of points onto a line with
direction $\mathbf{a}\in\mathbb{S}^{d-1}$ are ordered as
$\mathbf{a}\cdot\mathbf{x}_{\tau(1)}<\cdots<\mathbf{a}\cdot\mathbf{x}_{\tau(N)}$
and
$\mathbf{a}\cdot\mathbf{y}_{\tau(1)}<\cdots<\mathbf{a}\cdot\mathbf{y}_{\tau(N)}$,
for a certain permutation $\tau$ of $N$ elements. The strategy to be adopted
is clear. It entails a combination of the one-dimensional control delineated
in lemma 13 with the transversal control required to prove corollary 5,
followed by a case-by-case analysis.
2\. Universal approximation. In [27], the authors demonstrate that a narrow
neural ODE can approximate any simple function with compact support. The proof
hinges on three key aspects: interpolation capacity, the compressive nature of
neural ODEs, and a control strategy that bounds the support of the target
function irrespective of required approximation accuracy. However, minimizing
the number of time discontinuities required for this control, particularly by
increasing $p$, is a non-trivial task that may require an entirely different
approach.
3\. Neural transport equation. The problem of controlling the neural transport
equation (2.4) using constant controls is as yet unresolved. One possible
approach involves approximating both the initial and target measures with
atomic measures of the form
$\rho_{N}=\frac{1}{N}\sum_{n=1}^{N}\beta_{n}\delta_{\alpha_{n}}$, with
$\beta_{n}>0$ and $\alpha_{n}\in\mathbb{R}^{d}$, and then interpolating those
Dirac deltas by controlling the characteristic curves. However, a potential
issue arises as $N\to\infty$: the distance between $\rho(T)$ and
$\rho^{N}(T)$, the solutions to the transport equation with initial conditions
$\rho_{0}$ and $\rho_{0}^{N}$, respectively, may diverge significantly. This
error can be quantified using the Grönwall inequality, which suggests that the
Lipschitz constant could increase unboundedly with $N$, especially as the
number of controlled points grows.
4\. Minimizing the number of time jumps. Another interesting question is how
to frame the reduction of discontinuities as an optimal control problem. For
instance, one could penalize the frequency of time jumps by targeting the
total variation seminorm. However, this seminorm lacks regularity, and
moreover the class of piecewise constant functions is not a closed set of
admissible controls.
5\. Switching dimensions. In our simplified ResNet (1.1), the dimension
remains constant across layers. However, strategically varying the hidden
dimension by defining $p=p(t)$ could offer advantages, either by reducing
complexity through dimension shrinkage or by creating space through dimension
increase. Exploring effective methods to implement these transitions, whether
by employing projections or by applying nonlinear transformations to the data,
constitutes a compelling area for research.

##  4 Proofs

###  4.1 Basic dynamics

Figure 4: Left to right: Compression, parallel motion, expansion.

We describe the simplest dynamics that we can generate via (1.6) by
conveniently choosing $(\mathbf{w},\mathbf{a},b)$:

1\. For each $t>0$, the term $\mathbf{a}(t)\cdot\mathbf{x}+b(t)$ identifies a
hyperplane in $\mathbb{R}^{d}$. For instance, taking
$\mathbf{a}=\mathbf{e}_{k}$ and $b=-c$, we fix the hyperplane $h$ with the
equation $x^{(k)}-c=0$.

2\. The application of $\sigma$ and the product with the vector
$\mathbf{w}(t)$ yields the field $\mathbf{w}(t)\max\\{x^{(k)}-c,0\\}$, which
exhibits distinct dynamics in two complementary half-spaces:
$H^{+}\equiv\\{x^{(k)}-c>0\\}$, where the field equals
$\mathbf{w}(t)(x^{(k)}-c)$, and $H^{-}\equiv\\{x^{(k)}-c\leq 0\\}$, where the
field is zero, meaning this set remains stationary under the flow.

3\. The choice of $\mathbf{w}(t)$ specifies the orientation and magnitude of
the field. For example, $\mathbf{w}(t)=\pm\mathbf{e}_{k}$ yields
$\dot{x}^{(k)}(t)=\pm\max\\{(x^{(k)}-c),0\\}$, so the points in $H^{+}$ can
either be attracted to or repelled from $h$, enabling compression or expansion
along the $k$-th coordinate. Conversely, $\mathbf{w}(t)=\mathbf{e}_{i}$, with
$i\not=k$, results in $\dot{x}^{(i)}(t)=\max\\{(x^{(k)}-c),0\\}$. In this
case, we generate in $H^{+}$ a movement that is parallel to $h$, i.e., along
the coordinate $i$.

The three basic operations of compression, expansion and movement in parallel
with the hyperplane (represented in fig. 4) constitute our toolbox for many
subsequent proofs.

###  4.2 Proof of theorem 1.

We will employ the following lemma, whose proof we postpone to the end of this
subsection:

######  Lemma 9.

Let $N\geq 1$, $d\geq 2$ and consider the dataset
$\mathcal{D}=\\{(\mathbf{x}_{n},\mathbf{y}_{n})\\}_{n=1}^{N}$ as defined in
(2.1). There exists a change of coordinates in $\mathbb{R}^{d}$ such that

$x_{n}^{(1)}\not=x_{m}^{(1)}\qquad\text{and}\qquad
y_{n}^{(2)}\not=y_{m}^{(2)},\qquad\text{if }n\not=m.$ (4.1)

Under the separability condition (4.1), we achieve the exact control by
building on the methods developed in [27]. Let $p\geq 1$ be fixed. We divide
the proof in two steps, illustrated in fig. 5 and fig. 5.
Step 1: Control of $d-1$ coordinates. By (4.1), we can relabel the data
$\\{\mathbf{x}_{n}\\}_{n=1}^{N}$ to impose the ordering

$x_{1}^{(1)}<\cdots<x_{N}^{(1)}.$

We define a partition of $\\{\mathbf{x}_{n}\\}_{n=1}^{N}$ in $\lceil
N/p\rceil$ subsets by increasing order of the $x^{(1)}$-coordinate. The $j$-th
subset is

$X_{j}\coloneqq\\{\mathbf{x}_{(j-1)\cdot p+1},\dots,\mathbf{x}_{j\cdot
p}\\},\quad\text{for }j=1,\dots,\lceil N/p\rceil-1,$

and $X_{\lceil N/p\rceil}$ contains the remaining $N-p\lfloor N/p\rfloor$
points. We describe the control of the first subset $X_{1}$. We take controls
$\mathbf{a}_{i}=\mathbf{e}_{1}$ and $b_{i}\in\mathbb{R}$, for $i=1,\dots,p$,
satisfying

$-b_{1}<x_{1}^{(1)}<-b_{2}<x_{2}^{(1)}<\cdots<-b_{p}<x_{p}^{(1)}.$

These controls define a family of parallel hyperplanes, given by
$\mathbf{a}_{i}\cdot\mathbf{x}+b_{i}=x^{(1)}+b_{i}=0$, which separate the
points of $X$. In this way, $i-1$ terms of the sum in (1.3) cancel inside the
strip $-b_{i}<x^{(1)}<-b_{i+1}$ for each $i=1,\dots,p$, so (1.3) simplifies to

$\dot{\mathbf{x}}=\sum_{l=1}^{i}\big{\\{}\mathbf{w}_{l}x^{(1)}+\mathbf{w}_{l}b_{l}\big{\\}}.$
(4.2)

Figure 5: Left: Step 1. Fix $x^{(1)}$ and control $x^{(2)},\dots,x^{(d)}$.
Right: Step 2. Control $x^{(1)}$ while $x^{(2)},\dots,x^{(d)}$ are fixed.

We consider velocities of the form
$\mathbf{w}_{i}=(0,w_{i}^{(2)},\dots,w_{i}^{(d)}),$ where the components
$w_{i}^{(k)}\in\mathbb{R}$ have to be defined in order to achieve the exact
control in time $T=1$. The first point, $\mathbf{x}_{1}$, is subject only to
one velocity, $\mathbf{w}_{1}$, so

$x_{1}^{(k)}(t)=w_{1}^{(k)}(x_{1}^{(1)}+b_{1})t+x_{1}^{(k)},$

while $x_{1}^{(1)}$ remains fixed. Therefore, it is enough to take

$w_{1}^{(k)}=\frac{y_{1}^{(k)}-x_{1}^{(k)}}{x_{1}^{(1)}+b_{1}}.$

Similarly, for $i=2,\dots,p$, having fixed
$\mathbf{w}_{1},\dots,\mathbf{w}_{i-1}$ it is enough to take

$w_{i}^{(k)}=\frac{y_{i}^{(k)}-x_{i}^{(k)}-\sum_{l=1}^{i-1}w_{l}^{(k)}(x_{i}^{(1)}+b_{l})}{x_{i}^{(1)}+b_{i}},$

for $k=1,\dots,d$. The described procedure can be simultaneously done for each
$X_{j}$, with $j=2,\dots\lceil N/p\rceil$, taking into account that the fields
used to control $X_{1},\dots,X_{j-1}$ (all of them orthogonal to
$\mathbf{e}_{1}$) will be added as new terms in (4.2). In the end, we will
have, for every $n=1,\dots,N$:

$\Phi_{1}(\mathbf{x}_{n};W,A,\mathbf{b})^{(k)}=y_{n}^{(k)},\qquad\text{for
}k=2,\dots,d.$

The total number of iterations employed in this step is $\left\lceil
N/p\right\rceil$, which corresponds to $\left\lceil N/p\right\rceil-1$
switches.
Step 2: Control of the remaining coordinate. In a slight abuse of notation, we
redefine $\mathbf{x}_{n}\coloneqq\Phi_{1}(\mathbf{x}_{n})$, where $\Phi_{1}$
is the flow resulting from step 1. Once again, we can relabel the data, now
assuming

$x_{1}^{(2)}<\dots<x_{N}^{(2)}.$

Following the increasing order of the $x^{(2)}$-coordinate, we define
$X_{1},\dots,X_{\lceil N/p\rceil-1}$, each being a subset of
$\\{\mathbf{x}_{n}\\}_{n=1}^{N}$ with $p$ points, and $X_{\lceil N/p\rceil}$,
which contains the remaining $N-p\lfloor N/p\rfloor$ points.

We follow an analogous methodology to step 1. For each $j$, we define controls
$\mathbf{a}_{i}=\mathbf{e}_{2}$ and $b_{i}$ (for $i=1,\dots,p$) that separate
the points of $X_{j}$ using $p$ parallel hyperplanes, each described by the
equation $x^{(2)}=b_{i}$. Now, we consider velocities of the form
$\mathbf{w}_{i}=w_{i}\mathbf{e}_{1}$, where the values $w_{i}$ are determined,
as in step 1, to ensure

$\Phi_{1}(\mathbf{x}_{n};W,A,\mathbf{b})^{(1)}=y_{n}^{(1)},\qquad\text{for
}n=1,\dots,N.$

The number of switches employed in step 2 is $\left\lceil N/p\right\rceil-1.$
Then, by adding one more to transition between steps, the whole control
requires $L=2\left\lceil N/p\right\rceil-1$ switches, hence proving theorem 1.

###### Proof of lemma 9.

The set of vectors in $\mathbb{R}^{d}$ that are orthogonal to any point
$\mathbf{x}_{n}$ from $\mathcal{D}$, is a finite union of hyperplanes in
$\mathbb{R}^{d}$. Therefore, one can always choose
$\mathbf{u}_{1}\in\mathbb{S}^{d-1}$ inside the complement of this set in
$\mathbb{R}^{d}$.

Consider the orthogonal subspace
$S=\langle\mathbf{u}_{1}\rangle^{\perp}\subset\mathbb{R}^{d}$, of dimension
$d-1$. With a similar argument, one can select a vector
$\mathbf{u}_{2}\in\mathbb{S}^{d-1}\cap S$ that has a non-zero scalar product
with all of the vectors $\mathbf{y}_{n}$ $(n=1,\dots,N)$. Completing the pair
$(\mathbf{u}_{1},\mathbf{u}_{2})$ to form a orthonormal basis
$\mathcal{B}=\\{\mathbf{u}_{1},\dots,\mathbf{u}_{d}\\}$ of $\mathbb{R}^{d}$,
the dataset $\mathcal{D}$ will satisfy the separability condition eq. 4.1 when
expressed in $\mathcal{B}$. ∎

###  4.3 Proof of corollary 3

We aim to eliminate the initial step in the algorithm defined in the proof of
theorem 1. To achieve this, we seek a new vector basis in $\mathbb{R}^{d}$
where the input-target pairs inherently share the first coordinate. When
$d\leq N$, this condition usually cannot be met. However, when $d>N$, there
exists an orthonormal vector basis $\mathcal{B}\subset\mathbb{R}^{d}$ such
that for all $n=1,\ldots,N$, the first coordinates of the $N$ pairs with
respect to $\mathcal{B}$ satisfy $x_{n}^{(1)}=y_{n}^{(1)}$.

To construct such a vector system, without loss of generality we can assume
that $d=N+1$. Let
$(\mathbf{x},\mathbf{y})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ with
$\mathbf{x}\neq\mathbf{y}$. We seek a vector $\mathbf{u}\in\mathbb{S}^{d-1}$
such that $\mathbf{u}\cdot\mathbf{x}=\mathbf{u}\cdot\mathbf{y}$. This
condition is equivalent to $\mathbf{u}\cdot(\mathbf{x}-\mathbf{y})=0$, which
is satisfied by any unit vector $\mathbf{u}$ contained in the linear
hyperplane orthogonal to $\mathbf{x}-\mathbf{y}$.

For $d-1$ input-target pairs of points, we consider the corresponding
hyperplanes $\\{H_{n}\\}_{n=1}^{d-1}$. Note that some of these hyperplanes can
be repeated. So, the intersection $\bigcap_{i=1}^{d-1}H_{i}$ yields a linear
subspace of dimension at most $d-1$. We choose any unit vector
$\mathbf{e}_{1}^{\prime}$ contained in that subspace. Any completion to an
orthonormal basis
$\mathcal{B}=\\{\mathbf{e}_{1}^{\prime},\dots,\mathbf{e}_{d}\\}\subset\mathbb{R}^{d}$
will satisfy the desired condition for $d-1$ points. The procedure is
illustrated in fig. 6.

Figure 6: For $N=2$, $d=3$, construction of a new basis of $\mathbb{R}^{d}$ in
which the first coordinates of the pairs $(\mathbf{x}_{n},\mathbf{y}_{n})$ are
matched. theorem 1 is applied afterwards.

###  4.4 Proof of corollary 5

###### Proof of proposition 4.

We carry out a similar analysis to the one in [2], where the probability of
requiring $k$ hyperplanes to separate two sets of $N$ points for binary
classification was estimated. We introduce the random variable
$Z^{p}_{d,2N}(\mathcal{D})$, which assigns to each possible dataset

$\mathcal{D}=\\{(\mathbf{x}_{n},\mathbf{y}_{n})\\}_{n=1}^{N}\subset\operatorname{supp}(\mu)\times\operatorname{supp}(\mu)$

the minimum number of parallel hyperplanes needed to separate in
$\mathbb{R}^{d}$ every pair $(\mathbf{x}_{n},\mathbf{y}_{n})$ from the others.
Note that $\min_{\mathcal{D}}Z^{p}_{d,2N}(\mathcal{D})=N-1$. We estimate the
probability $P_{d,N}\coloneqq P(Z_{d,2N}^{p}(\mathcal{D})=N-1)$ for any $d$
and $N$.

First, we consider the one-dimensional case. Since all the points are sampled
from the same distribution, every possible configuration of the $2N$ points in
the real line will have the same probability, i.e., their distribution is
uniform on the finite space of all possible orderings. Therefore, we can
compute:

$P_{1,N}=\frac{\text{favorable configurations}}{\text{total
configurations}}=\frac{N!\,2^{N}}{(2N)!}.$

Now, we apply Stirling’s formula $\frac{n^{n}\sqrt{2\pi
n}}{n!e^{n}}\xlongrightarrow{n\to\infty}1$ to approximate, for sufficiently
large $N$:

$P_{1,N}\approx\frac{N^{N}\sqrt{2\pi N}(2/e)^{N}}{(2N)^{2N}\sqrt{4\pi
N}/e^{2N}}\\\ =\frac{1}{\sqrt{2}}\left(\frac{e}{2N}\right)^{N}.$

Let $Z_{d,2N}^{p,c}$ be similarly defined to $Z_{d,2N}^{p}$ but restricting
the hyperplanes to be orthogonal to one of the $d$ canonical axes. Then, for
any $d\geq 1$, we can bound:

$P(Z_{d,2N}^{p,c}=N-1)\leq P(Z_{d,2N}^{p}=N-1).$ (4.3)

By hypothesis, the $d$ random variables defined as $Z_{1,2N}^{p}$ over the
projection of $\mathcal{D}$ on each canonical axis are i.i.d. to
$Z_{1,2N}^{p}$, so we can write:

$\displaystyle P(Z_{d,2N}^{p,c}>N-1)$
$\displaystyle=\left[1-P(Z_{1,2N}=N-1)\right]^{d}$
$\displaystyle=\frac{N!\,2^{N}}{(2N)!}\approx\left[1-\frac{1}{\sqrt{2}}\left(\frac{e}{2N}\right)^{N}\right]^{d},$

if $N\gg 1$. By (4.3), the complementary provides the desired lower bound for
$P_{d,N}$. ∎

###### Proof of corollary 5.

Let $\mathbf{a}\in\mathbb{S}^{d-1}$,
$\\{b_{n}\\}_{n=1}^{N+1}\subset\mathbb{R}$ and $\tau$ be given by 1. With no
loss of generality, we can assume that $\mathbf{a}=\mathbf{e}_{1}$ and $\tau$
is the identity permutation. The argument that we will use is similar to the
one employed in the proof of theorem 1, but now the motion must be
longitudinal as well as transverse. It also hinges on the fact that, , inside
the $n$-th strip

$S_{n}\coloneqq\\{\mathbf{x}\in\mathbb{R}^{d}:-b_{n}<\mathbf{a}\cdot\mathbf{x}<-b_{n+1}\\},$

the equation (1.5) simplifies to (4.2). The simultaneous control of the data
points is achieved inductively, in increasing order of the first coordinates,
by appropriately defining the field $\mathbf{w}_{n}$ associated with each
hyperplane $H_{n}$. Both the base case and the inductive step are established
by the following two lemmas, which will be proven later.

######  Lemma 10.

Consider two points $\mathbf{x}_{1},\mathbf{y}_{1}\in\mathbb{R}^{d}$ with
$\mathbf{x}_{1}\neq\mathbf{y}_{1}$. For any $T>0$ and $b\in\mathbb{R}$
satisfying $x_{1}^{(1)}+b>0$ and $y_{1}^{(1)}+b>0$, there exists a unique
$\mathbf{w}\in\mathbb{R}^{d}$ such that the solution of

$\displaystyle\begin{cases}\dot{\mathbf{x}}&=\mathbf{w}\sigma(x^{(1)}+b),\\\
\mathbf{x}(0)&=\mathbf{x}_{1}\in\mathbb{R}^{d}\end{cases}$ (4.4)

reaches $\mathbf{x}(T)=\mathbf{y}_{1}$.

Having controlled $\mathbf{x}_{1},\dots,\mathbf{x}_{n-1}$ to
$\mathbf{y}_{1},\dots,\mathbf{y}_{n-1}$ in a time horizon $T>0$, and using
parameters $\\{\mathbf{w}_{i}\\}_{i=1}^{n-1}\subset\mathbb{R}^{d}$ and
$\\{b_{i}\\}_{i=1}^{n-1}\subset\mathbb{R}$ such that $b_{n-1}<\cdots<b_{1}$
and

$x_{i}^{(1)}+b_{i}>0\quad\text{and}\quad y_{i}^{(1)}+b_{i}>0\quad\text{for
}i=1,\dots,N,$

steering $\mathbf{x}_{n}$ involves overcoming an autonomous drift field

$\mathbf{d}(\mathbf{x})\coloneqq\sum_{i=1}^{n-1}\mathbf{w}_{i}\sigma(\mathbf{a}\cdot\mathbf{x}+b_{i})=\sum_{i=1}^{n-1}\mathbf{w}_{i}\sigma(x^{(1)}+b_{i}).$
(4.5)

The drift field $\mathbf{d}$ becomes more intense as the first coordinate
increases, owing to the characteristics of the ReLU function. However, the
following lemma shows that the control is possible:

######  Lemma 11.

Consider $\mathbf{x}_{n},\mathbf{y}_{n}\in\mathbb{R}^{d}$ with
$\mathbf{x}_{n}\neq\mathbf{y}_{n}$. With the above notation, for any $T>0$
there exists a unique $\mathbf{w}_{n}\in\mathbb{R}^{d}$ and some
$b_{n}\in\mathbb{R}$ satisfying

$b_{n}<b_{n-1},\quad x_{n}^{(1)}+b_{n}>0\quad\text{and}\quad
y_{n}^{(1)}+b_{n}>0$

such that the solution of the Cauchy problem

$\displaystyle\begin{cases}\dot{\mathbf{x}}&=\mathbf{d}(\mathbf{x})+\mathbf{w}_{n}\sigma(x^{(1)}+b_{n}),\\\
\mathbf{x}(0)&=\mathbf{x}_{n},\end{cases}$ (4.6)

where $\mathbf{d}$ is given by (4.5), reaches $\mathbf{x}(T)=\mathbf{y}_{n}$.

With lemmas 10 and 11, the inductive argument is almost complete. It is left
to show that the trajectory of each initial datum $\mathbf{x}_{n}$ will remain
in its corresponding strip $S_{n}$ for all $t\in(0,T)$.

On one hand, lemma 11 guarantees that the trajectory $\mathbf{x}(t)$
originating from $\mathbf{x}_{n}$ will reach the endpoint $\mathbf{y}_{n}$. On
the other hand, taking into account that the field is autonomous and invariant
along each hyperplane $x^{(1)}=\text{const}$, then $\dot{x}^{(1)}(t)$ cannot
change sign at any time, that is to say, the $x^{(1)}(t)$ does not change its
direction. Consequently, the entire trajectory will be contained within the
strip bounded by $x=x_{n}^{(1)}$ and $x=y_{n}^{(1)}$, which in turn is
contained in $S_{n}$. ∎

###### Proof of lemma 10.

In the half-space $\\{x^{(1)}+b>0\\}$, the equation (4.4) is written as

$\dot{\mathbf{x}}=\mathbf{w}\sigma(x^{(1)}+b)=\mathbf{w}x^{(1)}+\mathbf{w}\,b.$
(4.7)

We can assume that $b=0$, so the solution of (4.7) is

$\mathbf{x}(t)=\frac{x_{1}^{(1)}}{w^{(1)}}\mathbf{w}\left(e^{w^{(1)}t}-1\right)+\mathbf{x}_{1},$

which can be driven to $\mathbf{x}(T)=\mathbf{y}_{1}$ by taking

$w^{(1)}=\frac{1}{T}\ln\left(\frac{y_{1}^{(1)}}{x_{1}^{(1)}}\right)\quad\text{and}\quad
w^{(k)}=\frac{y_{1}^{(k)}-x_{1}^{(k)}}{y_{1}^{(1)}-x_{1}^{(1)}}\,w^{(1)},$

or $w^{(k)}=y_{1}^{(k)}-x_{1}^{(k)}$ if $x_{1}^{(1)}=y_{1}^{(1)}$, for
$k=2,\dots,d$. Moreover, $\mathbf{x}(t)$ stays in the half-space $x^{(1)}>0$
for $t\in[0,T]$ because $x^{(1)}(t)$ is monotone in that interval. ∎

###### Proof of lemma 11.

First, take any $b_{n}\in\mathbb{R}$ satisfying

$b_{n}<b_{n-1},\quad x_{n}^{(1)}+b_{n}>0\quad\text{and}\quad
y_{n}^{(1)}+b_{n}>0.$

For simplicity, we rewrite (4.4) as

$\dot{\mathbf{x}}=\left(\mathbf{s}_{n-1}+\mathbf{w}_{n}\right)x^{(1)}+\mathbf{w}_{n}b_{n}+\mathbf{c}_{n-1},$

where $\mathbf{s}_{n-1}=\sum_{i=1}^{n-1}\mathbf{w}_{i}$,
$\mathbf{c}_{n-1}=\sum_{i=1}^{n-1}\mathbf{w}_{i}b_{i}$. If we restrict to the
first coordinate, we have:

$\displaystyle\begin{cases}\dot{x}^{(1)}&=\left(s_{n-1}^{(1)}+w_{n}^{(1)}\right)x^{(1)}+w_{n}^{(1)}b_{n}+c_{n-1}^{(1)},\\\\[4.0pt]
x^{(1)}(0)&=x_{n}^{(1)},\end{cases}$

which has solution

$x^{(1)}(t)=\frac{w_{n}^{(1)}b_{n}+c_{n-1}^{(1)}}{s_{n-1}^{(1)}+w_{n}^{(1)}}\left[e^{\left(s_{n-1}^{(1)}+w_{n}^{(1)}\right)t}-1\right]\\\
+x_{n}^{(1)}e^{\left(s_{n-1}^{(1)}+w_{n}^{(1)}\right)t}.$

First, we want to see if there exists $\hat{w}_{n}^{(1)}\in\mathbb{R}$ such
that $x^{(1)}(T)=y_{n}^{(1)}$, or, equivalently, if the function

$f(z)=\frac{z\,b_{n}+c_{n-1}^{(1)}}{s_{n-1}^{(1)}+z}\left[e^{\left(s_{n-1}^{(1)}+z\right)T}-1\right]\\\
+x_{n}^{(1)}e^{\left(s_{n-1}^{(1)}+z\right)T}-y_{n}^{(1)}$ (4.8)

has a real root. For that task, we compute:

$\displaystyle\lim_{z\to\infty}f(z)$
$\displaystyle=\operatorname{sign}\left(b_{n}+x_{n}^{(1)}\right)\cdot\infty=+\infty,$
$\displaystyle\lim_{z\to-\infty}f(z)$ $\displaystyle=-b_{n}-y_{n}^{(1)}<0.$

On the other hand,

$\displaystyle\lim_{z\to{-s_{n-1}^{(1)}}^{+}}f(z)$
$\displaystyle=\big{(}c_{n-1}^{(1)}-b_{n}s^{(1)}_{n-1}\big{)}T+x_{n}^{(1)}-y_{n}^{(1)}$
$\displaystyle=\lim_{z\to{-s_{n-1}^{(1)}}^{-}}f(z),$

so $f$ is continuous. Hence, we can assure that there exists
$\hat{w}_{n}^{(1)}\in\mathbb{R}$ such that
$f\left(\hat{w}_{n}^{(1)}\right)=0$. Now, we denote

$s_{n}^{(1)}=\hat{w}_{n}^{(1)}+s_{n-1}^{(1)}\quad\text{and}\quad
c_{n}^{(1)}=\hat{w}_{n}^{(1)}b_{n}+c_{n-1}^{(1)}.$

For each component $j\in\\{2,\dots,d\\}$,

$\displaystyle\begin{cases}\dot{x}^{(j)}&=c_{n}^{(1)}\frac{s_{n-1}^{(j)}+w_{n}^{(j)}}{s_{n}^{(1)}}\left[e^{s_{n}^{(1)}t}-1\right]\\\
&\;+x_{n}^{(1)}\left(s_{n-1}^{(j)}+w_{n}^{(j)}\right)e^{s_{n}^{(1)}t}+w_{n}^{(j)}b_{n}+c_{n-1}^{(j)},\\\\[8.0pt]
x^{(j)}(0)&=x_{n}^{(j)},\end{cases}$

which has solution

$x^{(j)}(t)=x_{n}^{(j)}+\left[w_{n}^{(j)}b_{n}+c_{n-1}^{(j)}-c_{n}^{(1)}\frac{s_{n-1}^{(j)}+w_{n}^{(j)}}{s_{n}^{(1)}}\right]t\\\
+\left[c_{n}^{(1)}\frac{s_{n-1}^{(j)}+w_{n}^{(j)}}{s_{n}^{(1)^{2}}}+x_{n}^{(1)}\frac{s_{n-1}^{(j)}+w_{n}^{(j)}}{s_{n}^{(1)}}\right]e^{s_{n}^{(1)}t}.$

Now, we want to find a solution of $g(z)=0$ for

$g(z)=x_{n}^{(j)}+\left[zb_{n}+c_{n-1}^{(j)}-c_{n}^{(1)}\frac{s_{n-1}^{(j)}+z}{s_{n}^{(1)}}\right]T\\\
+\left[c_{n}^{(1)}\frac{s_{n-1}^{(j)}+z}{s_{n}^{(1)^{2}}}+x_{n}^{(1)}\frac{s_{n-1}^{(j)}+z}{s_{n}^{(1)}}\right]e^{s_{n}^{(1)}T}-y_{n}^{(j)}.$

This is an affine function in $z$, so $g(z)=0$ has a unique solution
$\hat{w}_{n}^{(j)}\in\mathbb{R}$ if and only if the slope is non-zero. Suppose
that the $b_{n}$ we have chosen yields a zero slope, i.e.,

$\left[\frac{c_{n}^{(1)}}{s_{n}^{(1)}}+x_{n}^{(1)}\right]e^{s_{n}^{(1)}T}+(b_{n}s_{n-1}^{(1)}-c_{n-1}^{(1)})T=0.$

Recalling that $f(\hat{w}_{n}^{(1)})=0$, we can write this equation as

$y_{n}^{(1)}+\frac{c_{n}^{(1)}}{s_{n}^{(1)}}+(b_{n}s_{n-1}^{(1)}-c_{n-1}^{(1)})T=0,$

so

$\frac{c_{n}^{(1)}}{s_{n}^{(1)}}=(c_{n-1}^{(1)}-b_{n}s_{n-1}^{(1)})T-y_{n}^{(1)}$
(4.9)

and we can solve for $\hat{w}_{n}$ in this equation as

$\hat{w}_{n}^{(1)}=\frac{\left[(c_{n-1}^{(1)}-b_{n}s_{n-1}^{(1)})T-y_{n}^{(1)}\right]s_{n-1}^{(1)}-c_{n-1}^{(1)}}{b_{n}-\left(c_{n-1}^{(1)}-b_{n}s_{n-1}^{(1)}\right)T-y_{n}^{(1)}}.$

Meanwhile, substituting (4.9) in $f(\hat{w}_{n}^{(1)})=0$, we get:

$\left[x_{n}^{(1)}-y_{n}^{(1)}\right]e^{s_{n}^{(1)}T}-T(c_{n-1}^{(1)}+b_{n}s_{n-1}^{(1)})[e^{s_{n}^{(1)}T}-1]=0,$

which can also be solved for $\hat{w}_{n}^{(1)}$ as

$\hat{w}_{n}^{(1)}=\frac{1}{T}\ln\left(\frac{(c_{n-1}^{(1)}-b_{n}s_{n-1}^{(1)})T}{(c_{n-1}^{(1)}-b_{n}s_{n-1}^{(1)})T-y_{n}^{(1)}+x_{n}^{(1)}}\right)-s_{n-1}^{(1)}.$

Equalizing both expressions of $\hat{w}_{n}^{(1)}$, we obtain an equation with
different analytic functions of $b_{n}$ and not involving $\hat{w}_{n}^{(1)}$.
Therefore, changing slightly the chosen value of $b_{n}$, these expressions
become different, independently of the corresponding value of
$\hat{w}_{n}^{(1)}$, so there is a perturbation of $b_{n}$ that satisfies the
statement of the lemma while $g(z)=0$ has a unique solution.

∎

###  4.5 Proof of theorem 7

First, we prove proposition 6 with an inductive argument of topological
nature.

###### Proof of proposition 6.

First, we aim to build a family of $N$ disjoint $C^{\infty}$ curves contained
in $\operatorname{Int}(\Omega)$, each connecting the two points of a
corresponding pair $(\mathbf{x}_{n},\mathbf{y}_{n})\in\mathcal{D}$. When
$d\geq 2$, any connected open set in $\mathbb{R}^{d}$ is path-connected.
Therefore, we can take a continuous path $C$ that connects any two given
points $(\mathbf{x},\mathbf{y})$ inside
$U\coloneqq\operatorname{Int}(\Omega)\setminus K$, where $K$ represents any
finite union of disjoint curves contained in $\operatorname{Int}(\Omega)$.
Moreover, by a well-known approximation argument this path can be chosen to be
$C^{\infty}$.

Now, we have $N$ disjoint $C^{\infty}$ curves $\\{C_{n}\\}_{n=1}^{N}$
contained in $\operatorname{Int}(\Omega)$, each connecting a corresponding
pair of points $(\mathbf{x}_{n},\mathbf{y}_{n})$. The tangent velocity field
of each curve is also $C^{\infty}$, so the assembled field
$\mathbf{V}^{\prime}$, defined in $\bigcup_{n=1}^{N}C_{n}$, is smooth too.
Since its domain is compact, it is also Lipschitz-continuous. The required
vector field $\mathbf{V}$ is provided by Kirzsbraun’s Theorem (see [32]),
which ensures the existence of a Lipschitz-continuous map
$\mathbf{V}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ that extends
$\mathbf{V}^{\prime}$ sharing the same Lipschitz constant. ∎

The proof of theorem 7 employs the following lemma from [7, Section 7.2.2]:

######  Lemma 12 (Approximation rate for Lipschitz functions).

Let $K$ be the unit ball in $\operatorname{Lip}(\Omega,\mathbb{R})$, where
$\Omega=[-R,R]^{d}$. We have

$\frac{1}{[\kappa\log_{2}\kappa]^{1/d}}\,C_{d,R}\leq
E(K,\Sigma_{\kappa})_{C(\Omega)}\leq
C_{d,R}\,\frac{\log_{2}\kappa}{\kappa^{1/d}},$

where:

  * 1.

$\Sigma_{\kappa}\coloneqq\left\\{S_{\kappa}:\mathbb{R}^{d}\to\mathbb{R}\right\\}$
is the space of shallow neural networks with $\kappa=(d+2)p$ parameters, $p$
being the number of neurons in the hidden layer;

  * 2.

$E(K,\Sigma_{\kappa})_{C(\Omega)}=\sup_{f\in
K}\inf_{S\in\Sigma_{\kappa}}\|f-S\|_{C(\Omega)}$ measures the capacity of
$\Sigma_{\kappa}$ to approximate any function in $K$.

This result was in turn derived from [3, Proposition 6], where the upper bound
is obtained for $K_{L}$, the ball of radius $L$ in
$\operatorname{Lip}(\Omega,\mathbb{R})$:

$E(K_{L},\Sigma_{\kappa})_{C(\Omega)}\leq
C_{d,R}\,L\;\frac{\log_{2}\kappa}{\kappa^{1/d}}.$ (4.10)

###### Proof of theorem 7.

Proposition 6 ensures that we can find a Lipschitz-continuous field
$\mathbf{V}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ such that the flow $\Psi_{T}$ of
the ODE

$\dot{\mathbf{x}}=\mathbf{V}(\mathbf{x})$ (4.11)

interpolates the dataset. The classical UAT result in [6] guarantees that we
can uniformly approximate in $\Omega$ with precision $\varepsilon/\sqrt{d}$
each of its components $V^{(i)}$ using a corresponding shallow neural network
$S^{i}_{\kappa_{i}}:\Omega\rightarrow\mathbb{R}$, with $\kappa_{i}\geq 1$ for
$i=1,\dots,d$. Moreover, lemma 12 quantifies the dependence of the error with
respect to the number of parameters of each $S^{i}_{\kappa_{i}}$, ensuring
that

$\kappa_{1}=\kappa_{2}=\cdots=\kappa_{d}=(d+2)p$

parameters suffice to ensure (4.10) on each component. The assembled field

$\mathbf{V}_{NN}=\left(S_{\kappa_{1}}^{1},\dots,S_{\kappa_{d}}^{d}\right):\mathbb{R}^{d}\to\mathbb{R}^{d}$

has $\kappa=(d+2)pd$ parameters and governs the equation

$\dot{\mathbf{x}}=\mathbf{V}_{NN}(\mathbf{x})$ (4.12)

and satisfies
$\sup_{\mathbf{x}\in\Omega}|\mathbf{V}(\mathbf{x})-\mathbf{V}_{NN}(\mathbf{x})|<\varepsilon.$
Suppose that the Lipschitz constant of $\mathbf{V}_{NN}$ is
$L_{NN}\leavevmode\nobreak\ >\leavevmode\nobreak\ 0$. Let
$\mathbf{X}_{V}(t;\mathbf{x}_{0})$ and $\mathbf{X}_{NN}(t;\mathbf{x}_{0})$ be
the respective trajectories in time $t>0$ that a point
$\mathbf{x}_{0}\leavevmode\nobreak\ \in\leavevmode\nobreak\ \Omega$ will
follow under the dynamics provided by (4.11) and (4.12). The deviation

$\mathbf{z}(t)=|\mathbf{X}_{V}(t;\mathbf{x}_{0})-\mathbf{X}_{NN}(t;\mathbf{x}_{0})|$

is bounded as

$\displaystyle\mathbf{z}(t)$
$\displaystyle\leq\int_{0}^{t}\left|\mathbf{V}(\mathbf{X}_{V}(s;\mathbf{x}_{0}))-\mathbf{V}_{NN}(\mathbf{X}_{NN}(s;\mathbf{x}_{0}))\right|\mathrm{ds}$
$\displaystyle\leq\int_{0}^{t}\big{\\{}\left|\mathbf{V}(\mathbf{X}_{V}(s;\mathbf{x}_{0}))-\mathbf{V}(\mathbf{X}_{NN}(s;\mathbf{x}_{0}))\right|$
$\displaystyle\;+\left|\mathbf{V}(\mathbf{X}_{NN}(s;\mathbf{x}_{0}))-\mathbf{V}_{NN}(\mathbf{X}_{NN}(s;\mathbf{x}_{0}))\right|\big{\\}}\mathrm{ds}$
$\displaystyle\leq
L_{V}\int_{0}^{t}\left|\mathbf{X}_{V}(s;\mathbf{x}_{0})-\mathbf{X}_{NN}(s;\mathbf{x}_{0})\right|\mathrm{ds}+\varepsilon\;t$
$\displaystyle=L_{V}\int_{0}^{t}\mathbf{z}(s)\mathrm{ds}+\varepsilon\;t.$

By Grönwall’s inequality, it follows that

$\mathbf{z}(t)\leq\varepsilon\;t\;\exp\left\\{L_{V}\;t\right\\}.$

On the other hand, in the second line, we could alternatively add and subtract
$\mathbf{V}_{NN}(\mathbf{X}_{V}(s;\mathbf{x}_{0}))$, so:

$\displaystyle\mathbf{z}(t)$
$\displaystyle\leq\int_{0}^{t}\big{\\{}\left|\mathbf{V}(\mathbf{X}_{V}(s;\mathbf{x}_{0}))-\mathbf{V}_{NN}(\mathbf{X}_{V}(s;\mathbf{x}_{0}))\right|$
$\displaystyle\;+\left|\mathbf{V}_{NN}(\mathbf{X}_{V}(s;\mathbf{x}_{0}))-\mathbf{V}_{NN}(\mathbf{X}_{NN}(s;\mathbf{x}_{0}))\right|\big{\\}}\mathrm{ds}$
$\displaystyle\leq\varepsilon\;t+L_{NN}\int_{0}^{t}\left|\mathbf{X}_{NN}(s;\mathbf{x}_{0})-\mathbf{X}_{V}(s;\mathbf{x}_{0})\right|\mathrm{ds}$
$\displaystyle=\varepsilon\;t+L_{NN}\int_{0}^{t}\mathbf{z}(s)\mathrm{ds}.$

By Grönwall’s inequality, it follows that

$\mathbf{z}(t)\leq\varepsilon\;t\;\exp\left\\{L_{NN}\;t\right\\}.$

Taking $\mathbf{x}_{0}=\mathbf{x}_{n}$ and $t=T$, the two bounds for
$\mathbf{z}(t)$ give:

$|\mathbf{y}_{n}-\Phi_{T}(\mathbf{x}_{n})|\leq\varepsilon\;T\;\exp\big{\\{}\min\\{L_{V},L_{NN}\\}\,T\big{\\}}.$

The approximation rate in $\varepsilon$ with respect to the number of
parameters is obtained by direct application of (4.10). ∎

###  4.6 Proof of theorem 8

We seek to find $(W,A,\mathbf{b})$ such that the generated vector field moves,
compresses and stretches the mass distributed following $\rho_{0}$, to drive
it approximately to the target density, given by $\rho_{*}$. This is achieved
in four steps, illustrated in figs. 7, 8, 9 and 10 for the case $d=2$.

Figure 7: Illustration of the division of $[0,1]^{2}$ into rectangles, each
containing a mass of $1/n^{2}$ following the distribution given by $\rho_{0}$.

1\. Preparation. We compress $\operatorname{supp}(\mu_{0})$ into $[0,1]^{d}$.
To do this, we aim to find a control

$\left(W,A,\mathbf{b}\right)\in L^{\infty}\left((0,T);\mathbb{R}^{d\times
p}\times\mathbb{R}^{p\times d}\times\mathbb{R}^{p}\right)$

such that, in a time $T_{1}>0$, the flow $\Phi_{T_{1}}$ of (1.3) satisfies

$\Phi_{T_{1}}(\operatorname{supp}(\mu_{0}))\subset[0,1]^{d}$

For $k=1,\dots,d$, we fix the hyperplane $x^{(k)}=0$ and a compressive
velocity field by taking
$(\mathbf{w},\mathbf{a},b)=(-\mathbf{e}_{k},\mathbf{e}_{k},0)$. We choose
$T_{1,k}>0$ sufficiently large to ensure

$\Phi_{T_{1,k}}\left(\operatorname{supp}(\mu_{0})\cap\\{x^{(k)}\geq
0\\}\right)\subset\\{0\leq x^{(k)}\leq 1\\}.$

We repeat the operation with the hyperplanes $x^{(k)}=1$ for $k=1,\dots,d$,
taking $(\mathbf{w},\mathbf{a},b)=(\mathbf{e}_{k},-\mathbf{e}_{k},1)$ and
$T_{1,k}^{\prime}>0$ such that

$\Phi_{T_{1,k}^{\prime}}\circ\Phi_{T_{1,k}}\left(\operatorname{supp}(\mu_{0})\cap\\{x^{(k)}\leq
0\\}\right)\subset\\{0\leq x^{(k)}\leq 1\\}.$

Both operations are possible in a finite time because $\mu_{0}$ has compact
support. In the end, we will have built piecewise constant controls
$(\mathbf{w},\mathbf{a},b)$ that take $2d$ values, such that

$\Phi_{T_{1}}(\operatorname{supp}(\mu_{0}))\subset[0,1]^{d},\quad\text{for
}T_{1}\coloneqq\sum_{k=1}^{d}\left(T_{1,k}+T_{1,k}^{\prime}\right).$

Using $p$ neurons, we can simultaneously apply $p$ controls, because the
characteristic curves of (2.4) when
$(\mathbf{w},\mathbf{a},b)=(\pm\mathbf{e}_{k},\pm\mathbf{e}_{k},1)$ are
parallel to the hyperplanes $\\{x^{(l)}=0\\}$ for every $l\neq k$. So, the
total number of values taken by $(W,A,\mathbf{b})$ is $\lceil 2d/p\rceil$.
2\. Partition. We aim to divide $[0,1]^{d}$ into a collection of $n^{d}$
hyperrectangles, each containing a mass of $1/n^{d}$, as distributed by
$\mu_{0}$. The process can be visualized in fig. 7 for $d=2$. For simplicity,
we redefine $\mu_{0}:=\Phi_{T_{1}\\#}\mu_{0}$ with density $\rho_{0}$, now
satisfying $\operatorname{supp}(\mu_{0})\subset[0,1]^{d}$. Let $n\geq 1$ and
consider the function

$t\longmapsto\int_{[0,t)\times[0,1]^{d-1}}d\mu_{0}=\int_{[0,t)\times[0,1]^{d-1}}\rho_{0}.$

This function is continuous, strictly increasing (by absolute continuity),
equal to $0$ at $t=0$ and equal to $1$ at $t=1$. Therefore, we can choose
$n+1$ numbers

$c_{0}=0<c_{1}<\dots<1=c_{n}$

such that, for $i_{1}=0,\dots,n-1,$

$\int_{[c_{i_{1}},c_{i_{1}+1}]\times[0,1]^{d-1}}\rho_{0}=\frac{1}{n}.$

Similarly, for each $i_{1}=0,\dots,n-1$ we choose $n+1$ numbers

$c_{i_{1},0}=0<c_{i_{1},1}<\dots<1=c_{i_{1},n}$

such that, for $i_{2}=0,\dots,n-1,$

$\int_{[c_{i_{1}},c_{i_{1}+1}]\times[c_{i_{1},i_{2}},c_{i_{1},i_{2}+1}]\times[0,1]^{d-2}}\rho_{0}=\frac{1}{n^{2}}.$

Repeating this operation recursively for each coordinate, we end up with
$n^{d}$ hyperrectangles

$C_{i_{1},\dots,i_{d}}^{0}\coloneqq\left[c_{i_{1}},c_{i_{1}+1}\right]\times\cdots\times\left[c_{i_{1},\dots,i_{d}},c_{i_{1},\dots,i_{d}+1}\right]\subset\mathbb{R}^{d},$

with $i_{k}\in\\{0,\dots,n-1\\}$ for every $k=1,\dots,d$, such that

$\int_{C_{i_{1},\dots,i_{d}}^{0}}\rho_{0}=\frac{1}{n^{d}}.$

The analogous partition for the uniform measure $\mu_{*}$ is

$G_{i_{1},\dots,i_{d}}\coloneqq\left[\frac{i_{1}}{n},\frac{i_{1}+1}{n}\right]\times\dots\times\left[\frac{i_{d}}{n},\frac{i_{d}+1}{n}\right].$

For the sake of readability, we will denote each multi-index by
$I=(i_{1},\dots,i_{d})\in\\{0,\dots,n-1\\}^{d}$ and write $C_{I}^{0}$ and
$G_{I}$.
3\. Control. We aim to define the controls that expand and compress the mass
until the hyperrectangles $C_{I}^{0}$ approximate a corresponding collection
of $n^{d}$ hypercubes, each of them containing the same mass $1/n^{d}$, as
distributed by $\mu_{*}$. Ideally, we would build
$\left(W,A,\mathbf{b}\right)$ such that the flow of the ODE (1.3) satisfied

$\Phi_{T}(C_{I}^{0};W,A,\mathbf{b})=G_{I}$ (4.13)

for each $I\in\\{0,\dots,n-1\\}^{d}$. This would be done by transforming each
hyperplane $\\{x^{(k)}=c_{i_{1},...,i_{k}}\\}$ into a target hyperplane
$\\{x^{(k)}=i_{k}/n\\}$.

Figure 8: Representation of the ideal transformation of the rectangles $C_{I}$
into the corresponding ones $G_{I}$.

However, this task is not possible in general, since any
$\\{x^{(k)}=i_{k}/n\\}$ can be a target for multiple distinct hyperplanes
$\\{x^{(k)}=c_{i_{1},...,i_{k}}\\}$. We therefore relax the problem to
$\delta$-approximate control by considering $\delta$-displacements of the
target hyperplanes $\\{x^{(k)}=i_{k}/n\\}$. Now we aim to control each
$C_{I}^{0}$ to a corresponding target $G_{I}^{\delta}$ that is $\delta$-close
to $G_{I}$, for a sufficiently small $\delta>0$.

Figure 9: Construction of the partition in rectangles $G_{I}^{\delta}$.

Let us build the new hyperrectangles $G_{I}^{\delta}$, a process shown in fig.
9. For each $k\in\\{2,\dots,d\\}$,
$(i_{1},\dots,i_{k-1})\in\\{0,\dots,n-1\\}^{k-1}$ and
$i_{k}\in\\{0,\dots,n\\}$, we define

$g_{i_{1},\dots,i_{k}}^{\delta}\coloneqq
i_{k}/n+\delta(c_{i_{1},\dots,i_{k}}-\tilde{c}_{i_{k}}),$

where

$\tilde{c}_{i_{k}}\coloneqq\min\\{c_{i_{1}^{\prime},\dots,i_{k-1}^{\prime},i_{k}}:(i_{1}^{\prime},\dots,i_{k-1}^{\prime})\in\\{0,\dots,n-1\\}^{k-1}\\}.$

Note that $g_{i_{1},\dots,i_{k}}^{\delta}=1$ whenever $i_{k}=n$. By
construction,

$g_{i_{1},\dots,i_{k-1},i_{k}}^{\delta}<g_{i_{1}^{\prime},\dots,i_{k-1}^{\prime},i_{k}}^{\delta}\iff
c_{i_{1},\dots,i_{k-1},i_{k}}<c_{i_{1}^{\prime},\dots,i_{k-1}^{\prime},i_{k}}$

and

$g_{i_{1},\dots,i_{k-1},i_{k}}^{\delta}=g_{i_{1}^{\prime},\dots,i_{k-1}^{\prime},i_{k}}^{\delta}\iff
c_{i_{1},\dots,i_{k-1},i_{k}}=c_{i_{1}^{\prime},\dots,i_{k-1}^{\prime},i_{k}}.$

By recursion, we define a new partition of $[0,1]^{d}$ into a collection of
rectangles $G_{I}^{\delta}$ with $I\in\\{0,\dots,n-1\\}^{d}$, where

$G_{I}^{\delta}\subset G_{I}+\\{0\\}\times[-\delta,\delta]^{d-1}.$

Moreover, this partition mimics the structure of the partition defined for
$\mu_{0}$, in the sense that there is the same number of distinct target
hyperplanes as initial ones to be controlled. To sum up, if we take
$\delta<1/n$, we end up with:

$\displaystyle\left\\{C_{I}^{0}:I\in\\{0,\dots,n-1\\}^{d}\right\\}\text{ s.t.
}$ $\displaystyle\int_{C_{I}^{0}}d\mu_{0}=\frac{1}{n^{d}},$ (4.14)
$\displaystyle\left\\{G_{I}^{\delta}:I\in\\{0,\dots,n-1\\}^{d}\right\\}\text{
s.t. }$ $\displaystyle\int_{G^{\delta}_{I}}d\mu_{*}\leq\frac{3^{d}}{n^{d}},$
and $\displaystyle\operatorname{diam}(G^{\delta}_{I})\leq
3\sqrt{d}\frac{1}{n}.$

It is left to map $C_{I}^{0}$ to $G_{I}^{\delta}$ for each $I$. This is based
on the following lemma, whose proof we postpone for readability.

######  Lemma 13.

Let $d\geq 2$, $\mu_{0}\in\mathcal{P}_{ac}^{c}(\mathbb{R}^{d})$ with density
$\rho_{0}$, $\rho_{*}$ the uniform density in $[0,1]^{d}$, and $T>0$ be fixed.
Let $n\geq 1$ and consider a family of hyperrectangles such as (4.14). For any
$p_{1},\dots,p_{d}\geq 1$, there exists a piecewise constant control

$(W,A,\mathbf{b})\in L^{\infty}((0,T);\mathbb{R}^{p\times
d}\times\mathbb{R}^{d\times p}\times\mathbb{R}^{p})$

with $p=p_{1}+\dots+p_{d}$ such that, for each
$I=(i_{1},\dots,i_{d})\in\\{0,\dots,n-1\\}^{d},$ the flow $\Phi_{T}$ generated
by (1.3) satisfies

$\Phi_{T}(C^{0}_{I};W,A,\mathbf{b})=G^{\delta}_{I}$

Furthermore, the number of discontinuities of $(W,A,\mathbf{b})$ is

$L=\max\\{\lceil n/p_{1}\rceil,\dots,\lceil n^{d}/p_{d}\rceil\\}-1.$

4\. Estimates. We compute the $W_{q}$-distance between both measures to verify
the approximate control. Let $\Phi_{T}$ be the flow given by lemma 13,
satisfying

$\Phi_{T}(C_{I})=G_{I}^{\delta},\;\text{ for
}I=(i_{1},\dots,i_{d})\in\\{0,\dots,n-1\\}^{d}.$

Let us quantify, in the Wasserstein-$q$ distance, the proximity of
$\Phi_{T\\#}\mu_{0}$ to $\mu_{*}$. We have

$W_{q}(\mu(T),\mu_{*})=W_{q}(\Phi_{T\\#}\mu_{0},\mu_{*})\\\
\leq\sum_{I\in\\{0,\dots,n-1\\}^{d}}W_{q}(\Phi_{T\\#}\mu_{0}|_{G_{I}^{\delta}},\mu_{*}|_{G_{I}^{\delta}}),$
(4.15)

see [33] for the inequality. For each $I\in\\{0,\dots,n-1\\}^{d}$, let
$\gamma_{I}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ be the measurable
function that satisfies

$\gamma_{I\\#}(\Phi_{T\\#}\mu_{0}|_{G_{I}^{\delta}})=\mu_{*}|_{G_{I}^{\delta}},$

attaining the minimum in Monge formulation (2.6) for $W_{q}$. In particular,
$\gamma_{I}$ only redistributes the mass inside $G_{I}^{\delta}$, so

$\displaystyle\int_{\mathbb{R}^{d}}|x-\gamma_{I}(x)|^{q}d\mu_{*}|_{G_{I}^{\delta}}$
$\displaystyle=\int_{G_{I}^{\delta}}|x-\gamma_{I}(x)|^{q}d\mu_{*}$
$\displaystyle\leq\operatorname{diam}(G_{I}^{\delta})^{q}\int_{G_{I}^{\delta}}d\mu_{*}\leq\frac{3^{q+d}d^{q/2}}{n^{q+d}}.$

Plugging this bound into (4.15), we conclude that

$\displaystyle W_{q}(\mu(T),\mu_{*})$ $\displaystyle\leq
3^{1+d/q}\sqrt{d}\,n^{d}\left(\frac{1}{n^{q+d}}\right)^{1/q}$
$\displaystyle=3^{1+d/q}\sqrt{d}\,n^{-(1+d/q-d)}.$

It follows that

$W_{q}(\mu(T),\mu_{*})\xlongrightarrow{n\to\infty}0\;\iff\;q<\frac{d}{d-1},$

and, in that case, $W_{q}(\mu(T),\mu_{*})<\varepsilon$ is obtained for

$n>\left(\frac{3^{1+d/q}\sqrt{d}}{\varepsilon}\right)^{\frac{1}{1+d/q-d}},$

hence proving theorem 8.

###### Proof of lemma 13.

We control the rectangles by mapping each hyperplane to its corresponding
target hyperplane. This strategy is illustrated in fig. 10.

The proof is divided into three steps. First, we achieve simultaneous control
for any $p$ hyperplanes orthogonal to a fixed direction. Second, for $N\geq p$
hyperplanes, we iteratively control $p$-subsets. Third, we show that this
approach is applicable to all $d$ canonical directions simultaneously.

Figure 10: From left to right, by rows: We transform the rectangles $C_{I}$ to
approximate the corresponding ones $G_{I}^{\delta}$ (in green)

Step 1. Let $k\in\\{1,\dots,d\\}$, $p\geq 1$ and

$\displaystyle-\infty<c_{1}<\cdots<c_{p}<\infty,$
$\displaystyle-\infty<g_{1}<\cdots<g_{p}<\infty.$

We aim to find constant controls
$\\{(\mathbf{w}_{i},\mathbf{a}_{i},b_{i})\\}_{i=1}^{p}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}$
such that the flow generated by (1.3) satisfies

$\displaystyle\Phi_{T}(\\{x^{(k)}=c_{i}\\})=\\{x^{(k)}=g_{i}\\},\quad
i=1,\dots,p.$ (4.16)

In particular, we will take $\mathbf{a}_{i}=\mathbf{e}_{k}$ for all $i$, and
$b_{i}$ such that $-b_{i}<\min\\{c_{i},g_{i}\\}$. Since the field

$\sum_{i=1}^{p}\mathbf{w}_{i}\sigma(\mathbf{a}_{i}\cdot\mathbf{x}+b_{i})=\sum_{i=1}^{p}\mathbf{w}_{i}\sigma(x^{(k)}+b_{i})$

only depends on the $x^{(k)}$-coordinate, it is projectable onto the
$x^{(k)}$-axis, i.e., the forward evolution of a hyperplane orthogonal to any
coordinate-axis is a hyperplane orthogonal to the same coordinate-axis, for
any time. Therefore, we can identify each $\\{x^{(k)}=c_{i}\\}$ and
$\\{x^{(k)}=g_{i}\\}$ with the point $c_{i}\in\mathbb{R}$ or
$g_{i}\in\mathbb{R}$ and study their evolution in the real line, so the
problem becomes one-dimensional. Thus, we identify $x^{(k)}\equiv x$ and fix
$\mathbf{w}_{i}=w_{i}\mathbf{e}_{k}$ with $w_{i}\in\mathrm{R}$, so we aim to
find $(w_{i})_{i=1}^{p}\subset\mathbb{R}$ and
$(b_{i})_{i=1}^{p}\subset\mathrm{R}$ such that $\Phi_{T}(c_{i})=g_{i}$ for
$i=1,\dots,p$. We proceed by induction on $p$.
First, we consider $p=1$ and let $c_{1},g_{1}\in\mathbb{R}$. Take any
$-b<\min\\{c_{1},g_{1}\\}$, so $\\{c_{1},g_{1}\\}\subset\\{x+b>0\\}$, and
$w=\operatorname{sign}(g_{1}-c_{1}),$ pointing from $c_{1}$ to $g_{1}$. Then,

$\Phi_{\tilde{T}}\big{(}c_{1}\big{)}=g_{1},\quad\text{for
}\tilde{T}=\left|\log\left(\frac{g_{1}-b}{c_{1}-b}\right)\right|.$

Using a time rescaling argument identical to the one in the proof of 1, this
control can be done in any time $T>0$. Also note that $\\{x+b\leq 0\\}$ is
fixed by the flow $\Phi_{t}$.
In the inductive step, we assume that the statement is true for some $p\geq
1$, and consider

$\displaystyle-\infty<c_{1}<\cdots<c_{p+1}<\infty,$
$\displaystyle-\infty<g_{1}<\cdots<g_{p+1}<\infty.$

Let $\\{(w_{i},b_{i})\\}_{i=1}^{p}\subset\mathbb{R}\times\mathbb{R}$ with
$-b_{i}<\min\\{c_{i},g_{i}\\}$ for all $i$, such that the flow of the one-
dimensional problem satisfies

$\Phi_{T}(c_{i})=g_{i},\quad i=1,\dots,p.$

We want to add a new pair $(w_{p+1},b_{p+1})\in\mathbb{R}\times\mathbb{R}$
such that

$\Phi_{T}(c_{i})=g_{i},\quad i=1,\dots,p+1.$

_Case 1: $c_{p}<g_{p}$._ Take $b_{p+1}=-g_{p}$, so (1.3) becomes

$\dot{x}=\sum_{i=1}^{p}w_{i}(x+b_{i})\mathbbm{1}_{c_{i}<x}(x)+w_{p+1}(x-g_{p})\mathbbm{1}_{g_{p}<x}(x).$

The added velocity $w_{p+1}(x-g_{p})$ only acts on the half-space $\\{x\geq
g_{p}\\}$, so the points $\\{c_{i}\\}_{i=1}^{p}$ are only subject to the drift
field

$d(x)\coloneqq\sum_{i=1}^{p}w_{i}(x+b_{i})\mathbbm{1}_{c_{i}<x}(x).$

Therefore, when adding $(w_{p+1},b_{p+1})$, we still have

$\Phi_{T}(c_{i})=g_{i},\quad i=1,\dots,p.$

If $c_{p+1}\leq g_{p}$, there is $0\leq s<T$ such that
$\Phi_{s}(c_{p+1})=g_{p}$. Otherwise, if $c_{p+1}>g_{p}$, consider $s=0$. Note
that $s$ only depends on the $p$ first neurons, on $c_{p+1}$ and on $g_{p}$,
and it is thus independent of $(w_{p+1},b_{p+1})$. Therefore, $c_{p+1}$ is
only subject to $d(x)$ for $t\in(0,s)$, and to $d(x)+w_{p+1}(x-g_{p})$ for
$t\in(s,T)$. More precisely:

$\frac{d}{dt}\Phi_{t}(c_{p+1})=d(\Phi_{t}(c_{p+1}))\\\
+w_{p+1}(\Phi_{t}(c_{p+1})-g_{p})\mathbbm{1}_{s\leq t<T}(t).$

Then, with a similar computation to that of lemma 11:

$\displaystyle\Phi_{T}(c_{p+1})=$
$\displaystyle\Phi_{T-s}\circ\Phi_{s}(c_{p+1})$ $\displaystyle=$
$\displaystyle\Big{(}g_{p}+\frac{\sum_{i=1}^{p}w_{i}b_{i}-w_{p+1}g_{p}}{\sum_{i=1}^{p}w_{i}+w_{p+1}}\Big{)}e^{(T-s)\sum_{i=1}^{p+1}w_{i}}$
$\displaystyle-\frac{\sum_{i=1}^{p}w_{i}b_{i}-w_{p+1}g_{p}}{\sum_{i=1}^{p}w_{i}+w_{p+1}}.$

The expression converges to $g_{p}$ when $w_{p+1}\to-\infty$, and diverges to
$\infty$ when $w_{p+1}\to\infty$. Therefore, by continuity, and since
$g_{p}<g_{p+1}$, there exists $w_{p+1}\in\mathbb{R}$ such that
$\Phi_{T}(c_{p+1})=g_{p+1}$.
_Case 2: $g_{p}\leq c_{p}<g_{p+1}$._ Take $b_{p+1}=-c_{p}$, so (1.3) becomes

$\dot{x}=\sum_{i=1}^{p}w_{i}(x+b_{i})\mathbbm{1}_{c_{i}<x}(x)+w_{p+1}(x-c_{p})\mathbbm{1}_{c_{p}<x}(x),$

Again, the points $\\{c_{i}\\}_{i=1}^{p}$ are only subject to $d(x)$ because
the added velocity only acts on $\\{x\geq c_{p}\\}$, so

$\Phi_{T}(c_{i})=g_{i},\quad i=1,\dots,p.$

The point $c_{p+1}$ is subject to the total velocity at $t=0$. A similar
computation to Case 1 leads to

$\displaystyle\Phi_{T}(c_{p+1})=$
$\displaystyle\left(c_{p}+\frac{\sum_{i=1}^{p}w_{i}b_{i}-w_{p+1}c_{p}}{\sum_{i=1}^{p}w_{i}+w_{p+1}}\right)e^{T\sum_{i=1}^{p+1}w_{i}}$
$\displaystyle-\frac{\sum_{i=1}^{p}w_{i}b_{i}-w_{p+1}c_{p}}{\sum_{i=1}^{p}w_{i}+w_{p+1}}.$

$\Phi_{T}(c_{p+1})$ tends to $c_{p}$ when $w_{p+1}\to-\infty$, and diverges to
$\infty$ when $w_{p+1}\to\infty$. By continuity, and since $c_{p}<g_{p+1}$,
there exists $w_{p+1}\in\mathbb{R}$ such that $\Phi_{T}(c_{p+1})=g_{p+1}$.
_Case 3: $g_{p+1}\leq c_{p}$._ Take $b_{p+1}=-c_{p}$, so (1.3) becomes

$\dot{x}=\sum_{i=1}^{p}w_{i}(x+b_{i})\mathbbm{1}_{c_{i}<x}(x)+w_{p+1}(x-c_{p})\mathbbm{1}_{c_{p}<x},$

and

$\Phi_{T}(c_{i})=g_{i},\quad i=1,\dots,p.$

Let $0<s\leq T$ be the first time such that $\Phi_{s}(c_{p+1})=c_{p}$. An
analogous computation to the previous cases gives

$\displaystyle\Phi_{s}(c_{p+1})=$
$\displaystyle\left(c_{p+1}+\frac{\sum_{i=1}^{p}w_{i}b_{i}-w_{p+1}c_{p}}{\sum_{i=1}^{p+1}w_{i}}\right)e^{s\sum_{i=1}^{p+1}w_{i}}$
$\displaystyle-\frac{\sum_{i=1}^{p}w_{i}b_{i}-w_{p+1}c_{p}}{\sum_{i=1}^{p+1}w_{i}},$

so

$s=\frac{1}{\sum\limits_{i=1}^{p+1}w_{i}}\log\left(\frac{\sum\limits_{i=1}^{p}w_{i}(c_{p}-b_{i})+2w_{p+1}c_{p}}{\sum\limits_{i=1}^{p}w_{i}(c_{p+1}+b_{i})+w_{p+1}(c_{p+1}-c_{p})}\right).$

By varying $w_{p+1}$ in $(-\sum_{i=1}^{p}w_{i},+\infty)$, we can ensure that
$s$ can take any value in $(0,T)$. Then,

$\displaystyle\Phi_{T}(c_{p+1})=$
$\displaystyle\Phi_{T-s}\circ\Phi_{s}(c_{p+1})$ $\displaystyle=$
$\displaystyle\left(c_{p}+\frac{\sum_{i=1}^{p}w_{i}b_{i}}{\sum_{i=1}^{p}w_{i}}\right)e^{(T-s)\sum_{i=1}^{p}w_{i}}$
$\displaystyle-\frac{\sum_{i=1}^{p}w_{i}b_{i}}{\sum_{i=1}^{p}w_{i}}.$

For $s=T$, the expression is equal to $c_{p}$, and for $s=0$ it equals
$\Phi_{T}(c_{p})=g_{p}$. By continuity, and also because $g_{p}<g_{p+1}\leq
c_{p}$, an argument like in case 1 ensures the existence of
$w_{p+1}\in\mathbb{R}$ such that $\Phi_{T}(c_{p+1})=g_{p+1}.$
Step 2. Let $k\in\\{1,\dots,d\\}$, $p\geq 1,$ and $N\geq p$ and

$-\infty<c_{1}<\dots<c_{N}<\infty,\quad-\infty<g_{1}<\dots<g_{N}<\infty.$

We will show that there exist piecewise constant controls
$\left(w_{i},a_{i},b_{i}\right)_{i=1}^{p}$ such that the flow of (1.3)
satisfies

$\Phi_{T}(c_{i})=g_{i},\quad i=1,\dots,N,$

and the number of discontinuities is $\lceil N/p\rceil-1$. We use a similar
argument to the one in the proof of theorem 1. We divide
$\\{c_{i}\\}_{i=1}^{N}$ and $\\{g_{i}\\}_{i=1}^{N}$ into subsets of $p$ points

$C_{j}\coloneqq\\{c_{(j-1)\cdot p+1},\dots,c_{j\cdot
p}\\},\;G_{j}\coloneqq\\{g_{(j-1)\cdot p+1},\dots,g_{j\cdot p}\\},$

for $j=1,\dots,\lceil N/p\rceil-1$, and $C_{\lceil N/p\rceil}$, $G_{\lceil
N/p\rceil}$ with the remaining $N-p\lfloor N/p\rfloor$ points.

The piecewise constant controls are obtained by induction on $j$. In each
iteration, we apply step 1 to define $p$ constant controls
$(w_{i}^{j},a_{i}^{j},b_{i}^{j})_{i=1}^{p}$ that map the $p$ points of $C_{j}$
to the corresponding ones in $G_{j}$ in time $\frac{T}{\lceil N/p\rceil}$.
Note that the initialization of induction in step 1 ensures that the
previously controlled subsets $C_{1},\dots,C_{j-1}$ can remain fixed during
the subsequent iterations, which trivializes the induction. Finally, we have
the piecewise constant controls

$(w_{i},a_{i},b_{i})_{i=1}^{p}=\sum_{j=1}^{\lceil
N/p\rceil}(w_{i}^{j},a_{i}^{j},b_{i}^{j})_{i=1}^{p}\mathbbm{1}_{\left(\frac{(j-1)T}{\lceil
N/p\rceil},\frac{jT}{\lceil N/p\rceil}\right)}(t),$

which achieve the desired objective.
Step 3. For every $k=1,\dots,d$, let $p_{k},N_{k}\geq 1$, and

$-\infty<c_{1}^{k}<\cdots<c_{N_{k}}^{k}<\infty,\quad-\infty<g_{1}^{k}<\cdots<g_{N_{k}}^{k}<\infty.$

For each fixed direction $k\in\\{1,\dots,d\\}$, step 2 is used to build
piecewise constant controls

$(\mathbf{w}_{j,k},\mathbf{a}_{j,k},b_{j,k})_{j=1}^{p_{k}}=(w_{j,k}\,\mathbf{e}_{k},\mathbf{e}_{k},b_{j,k})_{j=1}^{p_{k}}$

with $\lceil N_{k}/p_{k}\rceil-1$ discontinuities, such that

$\Phi_{T}(\\{x^{(k)}=c_{i}^{k}\\})=\\{x^{(k)}=g_{i}^{k}\\},\quad i=1,\dots,N.$

Moreover, only the $k$-th coordinate is varying on each flow, as argued in
steps 1 and 2 when we simplified the problem to one dimension. We define each
tern of the assembled control
$(\mathbf{w}_{j},\mathbf{a}_{j},b_{j})_{j=1}^{p}$, with
$p=\sum_{k=1}^{d}p_{k}$, by

$(\mathbf{w}_{j},\mathbf{a}_{j},b_{j})=\big{(}\mathbf{w}_{j-\sum_{i=1}^{k-1}p_{i},p_{k}},\mathbf{a}_{j-\sum_{i=1}^{k-1}p_{i},p_{k}},b_{j-\sum_{i=1}^{k-1}p_{i},p_{k}}\big{)},$

for $\sum_{i=1}^{k-1}p_{i}+1\leq j\leq\sum_{i=1}^{k}p_{i}$. Therefore, the
resulting neural ODE (1.3) on each coordinate writes

$\dot{x}^{(k)}=w_{1,k}^{(k)}(x^{(k)}+b_{1,k})+\dots+w_{p_{k},k}^{(k)}(x^{(k)}+b_{p_{k},k})$

All the equations of the system are independent, so each movement does not
interfere with the other $d-1$ movements. Therefore, the corresponding flow of
(1.3) satisfies

$\Phi_{T}(\\{x^{(k)}=c_{i}^{k}\\})=\\{x^{(k)}=g_{i}^{k}\\}$

for $k=1,\dots,d$ and $i=1,\dots,N_{k}$, and moreover, the number of
discontinuities in the controls is

$\displaystyle L$ $\displaystyle=\max\\{\lceil
N_{1}/p_{1}\rceil-1,\dots,\lceil N_{d}/p_{d}\rceil-1\\}$
$\displaystyle=\max_{k=1,\dots,d}\lceil N_{k}/p_{k}\rceil-1.$

Recalling that $N_{k}=n^{k}$ for $k=1,\dots,d$ (by construction of the
rectangles $C_{I}^{0}$ and $G_{I}^{\delta}$), it follows the desired result.

∎

## References

  * Agrachev and Sarychev [2021] Agrachev, A.A., Sarychev, A.V., 2021\.  Control on the manifolds of mappings with a view to the deep learning.  Journal of Dynamical and Control Systems 28, 989–1008.
  * Alvarez-Lopez et al. [2023] Alvarez-Lopez, A., Orive-Illera, R., Zuazua, E., 2023.  Optimized classification with neural odes via separability.  arXiv:2312.13807.
  * Bach [2017] Bach, F., 2017.  Breaking the curse of dimensionality with convex neural networks.  Journal of Machine Learning Research 18, 1–53.
  * Chang et al. [2018] Chang, B., Meng, L., Haber, E., Tung, F., Begert, D., 2018.  Multi-level residual networks from dynamical systems view.  arXiv:1710.10348.
  * Chen et al. [2018] Chen, R.T.Q., Rubanova, Y., Bettencourt, J., Duvenaud, D., 2018\.  Neural ordinary differential equations, in: Proceedings of the 32nd International Conference on Neural Information Processing Systems, Curran Associates Inc.. p. 6572–6583.
  * Cybenko [1989] Cybenko, G., 1989.  Approximation by superpositions of a sigmoidal function.  Mathematics of Control, Signals and Systems 2, 303–314.
  * DeVore et al. [2021] DeVore, R., Hanin, B., Petrova, G., 2021.  Neural network approximation.  Acta Numerica 30, 327–444.
  * Dupont et al. [2019] Dupont, E., Doucet, A., Teh, Y.W., 2019.  Augmented neural odes, in: Proceedings of the 33rd International Conference on Neural Information Processing Systems, Curran Associates Inc.. p. 3140–3150.
  * Duprez et al. [2017] Duprez, M., Morancey, M., Rossi, F., 2017.  Controllability and optimal control of the transport equation with a localized vector field, in: 2017 25th Mediterranean Conference on Control and Automation (MED), pp. 74–79.
  * E [2017] E, W., 2017.  A proposal on machine learning via dynamical systems.  Communications in Mathematics and Statistics 5, 1–11.
  * Elamvazhuthi et al. [2022] Elamvazhuthi, K., Gharesifard, B., Bertozzi, A.L., Osher, S., 2022\.  Neural ode control for trajectory approximation of continuity equation.  IEEE Control Systems Letters 6, 3152–3157.
  * Eldan and Shamir [2015] Eldan, R., Shamir, O., 2015\.  The Power of Depth for Feedforward Neural Networks.  JMLR: Workshop and Conference Proceedings 49, 1–34.
  * Esteve et al. [2021] Esteve, C., Geshkovski, B., Pighin, D., Zuazua, E., 2021\.  Large-time asymptotics in deep learning.  arXiv:2008.02491.
  * Esteve-Yagüe and Geshkovski [2023] Esteve-Yagüe, C., Geshkovski, B., 2023\.  Sparsity in long-time control of neural ODEs.  Systems Control Lett. 172, Paper No. 105452, 14.
  * Fan et al. [2020] Fan, F., Lai, R., Wang, G., 2020.  Quasi-equivalence of width and depth of neural networks.
  * Grathwohl et al. [2019] Grathwohl, W., Chen, R.T.Q., Bettencourt, J., Sutskever, I., Duvenaud, D., 2019.  FFJORD: free-form continuous dynamics for scalable reversible generative models, in: 7th International Conference on Learning Representations, ICLR 2019, New Orleans.
  * Haber and Ruthotto [2017] Haber, E., Ruthotto, L., 2017\.  Stable architectures for deep neural networks.  Inverse Problems 34, 014004\.
  * Hardt and Ma [2017] Hardt, M., Ma, T., 2017.  Identity matters in deep learning, in: International Conference on Learning Representations.
  * Huang [2003] Huang, G.B., 2003.  Learning capability and storage capacity of two-hidden-layer feedforward networks.  IEEE Transactions on Neural Networks 14, 274–281.
  * Kobyzev et al. [2021] Kobyzev, I., Prince, S.J., Brubaker, M.A., 2021.  Normalizing flows: An introduction and review of current methods.  IEEE Transactions on Pattern Analysis and Machine Intelligence 43, 3964–3979.
  * Li et al. [2022] Li, Q., Lin, T., Shen, Z., 2022.  Deep learning via dynamical systems: An approximation perspective.  Journal of the European Mathematical Society 25, 1671–1709.
  * Lin and Jegelka [2018] Lin, H., Jegelka, S., 2018.  Resnet with one-neuron hidden layers is a universal approximator, in: Proceedings of the 32nd International Conference on Neural Information Processing Systems, p. 6172–6181.
  * Lu et al. [2017] Lu, Z., Pu, H., Wang, F., Hu, Z., Wang, L., 2017.  The expressive power of neural networks: A view from the width, in: Proceedings of the 31st International Conference on Neural Information Processing Systems, Curran Associates Inc., Red Hook, NY, USA. p. 6232–6240.
  * Mhaskar et al. [2017] Mhaskar, H., Liao, Q., Poggio, T., 2017.  When and why are deep networks better than shallow ones?, in: Proceedings of the AAAI conference on artificial intelligence.
  * Papamakarios et al. [2021] Papamakarios, G., Nalisnick, E., Rezende, D.J., Mohamed, S., Lakshminarayanan, B., 2021.  Normalizing flows for probabilistic modeling and inference.  J. Mach. Learn. Res. 22.
  * Pinkus [1999] Pinkus, A., 1999.  Approximation theory of the mlp model in neural networks.  Acta Numerica 8, 143–195.
  * Ruiz-Balet and Zuazua [2023] Ruiz-Balet, D., Zuazua, E., 2023\.  Neural ODE control for classification, approximation, and transport.  SIAM Rev. 65, 735–773.
  * Ruiz-Balet and Zuazua [2024] Ruiz-Balet, D., Zuazua, E., 2024\.  Control of neural transport for normalising flows.  Journal de Mathématiques Pures et Appliquées 181, 58–90.
  * Ruthotto [2024] Ruthotto, L., 2024.  Differential equations for continuous-time deep learning.  arXiv:2401.03965.
  * Scagliotti [2023] Scagliotti, A., 2023.  Deep learning approximation of diffeomorphisms via linear-control systems.  Mathematical Control and Related Fields 13, 1226–1257.
  * Tabuada and Gharesifard [2023] Tabuada, P., Gharesifard, B., 2023\.  Universal approximation power of deep residual neural networks through the lens of control.  IEEE Transactions on Automatic Control 68, 2715–2728.
  * Valentine [1945] Valentine, F.A., 1945.  A lipschitz condition preserving extension for a vector function.  American Journal of Mathematics 67, 83–93.
  * Villani [2008] Villani, C., 2008.  Optimal transport – Old and new. Springer Berlin, Heidelberg. volume 338\.  pp. xxii+973.
  * Yun et al. [2019] Yun, C., Sra, S., Jadbabaie, A., 2019.  Small relu networks are powerful memorizers: A tight analysis of memorization capacity, in: Proceedings of the 33rd International Conference on Neural Information Processing Systems, Curran Associates Inc.. p. 15558–15569.
  * Zhang et al. [2016] Zhang, C., Bengio, S., Hardt, M., Recht, B., Vinyals, O., 2016.  Understanding deep learning requires rethinking generalization.  Communications of the ACM 64.