# Dynamical systems’ based neural networks

Elena Celledoni Department of Mathematical Sciences, NTNU, N-7491 Trondheim,
Norway<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Davide Murari 11footnotemark: 1 Brynjulf
Owren11footnotemark: 1 Carola-Bibiane Schönlieb Department of Applied
Mathematics and Theoretical Physics, University of Cambridge, Wilberforce
Road, Cambridge CB3 0WA, UK<EMAIL_ADDRESS><EMAIL_ADDRESS>Ferdia
Sherry22footnotemark: 2

###### Abstract

Neural networks have gained much interest because of their effectiveness in
many applications. However, their mathematical properties are generally not
well understood. If there is some underlying geometric structure inherent to
the data or to the function to approximate, it is often desirable to take this
into account in the design of the neural network. In this work, we start with
a non-autonomous ODE and build neural networks using a suitable, structure-
preserving, numerical time-discretisation. The structure of the neural network
is then inferred from the properties of the ODE vector field. Besides
injecting more structure into the network architectures, this modelling
procedure allows a better theoretical understanding of their behaviour. We
present two universal approximation results and demonstrate how to impose some
particular properties on the neural networks. A particular focus is on
1-Lipschitz architectures including layers that are not 1-Lipschitz. These
networks are expressive and robust against adversarial attacks, as shown for
the CIFAR-10 dataset.

###### keywords:

Neural networks, dynamical systems, Lipschitz networks, Structure preserving
deep learning, Universal approximation theorem.

65L05, 65L06, 37M15

##  1 Introduction

Neural networks have been employed to accurately solve many different tasks
(see, e.g., [4, 11, 48, 31]). Indeed, because of their excellent approximation
properties, ability to generalise to unseen data, and efficiency, neural
networks are one of the preferred techniques for the approximation of
functions in high-dimensional spaces.

In spite of this popularity, a substantial part of results and success stories
in deep learning still relies on empirical evidence and more theoretical
insight is needed. Recently, a number of scientific papers on the mathematical
foundations of neural networks have appeared in the literature, [8, 64, 53,
54, 59, 29]. In a similar spirit, many authors consider the design of deep
learning architectures taking into account specific mathematical properties
such as stability, symmetries, or constraints on the Lipschitz constant, [32,
27, 23, 56, 19, 24, 60, 30, 62, 66]. Even so, the imposition of structure on
neural networks is often done in an ad hoc manner, making the resulting input
to output mapping $F:\mathcal{X}\rightarrow\mathcal{Y}$ hard to analyse. In
this paper, we describe a general and systematic way to impose desired
mathematical structure on neural networks leading to an easier approach to
their analysis.

There have been multiple attempts to formulate unifying principles for the
design of neural networks. We hereby mention Geometric Deep Learning (see e.g.
[12, 11]), Neural ODEs (see e.g. [18, 41, 50, 65]), the continuous-in-time
interpretation of Recurrent Neural Networks (see e.g. [51, 17]) and of
Residual Neural Networks (see e.g. [64, 39, 15, 52, 1]). In this work, we
focus on Residual Neural Networks (ResNets) and build upon their continuous
interpretation.

Neural networks are compositions of parametric maps, i.e. we can characterise
a neural network as a map $\mathcal{N}=f_{\theta_{k}}\circ\ldots\circ
f_{\theta_{1}}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, with
$f_{\theta_{i}}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}^{n_{i+1}}$ being the
network layers. For ResNets the most important parametric maps are of the form

(1) $x\mapsto f_{\theta_{i}}(x)=x+h\Lambda(\theta_{i},x).$

The continuous-in-time interpretation of ResNets arises from the observation
that if $n_{i}=n_{i+1}$, $f_{\theta_{i}}$ coincides with one $h-$step of the
explicit Euler method applied to the non-autonomous ODE
$\dot{x}(t)=\Lambda(\theta(t),x(t))$. In this work, we consider piecewise-
autonomous systems, i.e. we focus on time-switching systems of the form

(2)
$\dot{x}(t)=f_{s(t)}(x(t)),\,\,s:[0,T]\rightarrow\\{1,\ldots,N\\},\,f_{i}\in\mathcal{F},$

with $s$ being piecewise constant (see e.g. [50, 36]), and $\mathcal{F}$ a
family of parametric vector functions. This simplification is not restrictive
and can help analyse and design neural networks, as we will clarify throughout
the paper.

This interpretation of ResNets gives the skeleton of our reasoning. Indeed, we
replace the explicit Euler method in Eq. 1 with suitable (structure-
preserving) numerical flows of appropriate vector fields. We call “dynamical
blocks” the groups of layers obtained with these numerical flows. The choice
of the vector field is closely related to the structure to impose. For
example, to derive symplectic neural networks we would apply symplectic time
integrators to Hamiltonian vector fields. This approach enables us to derive
new structured networks systematically and collocate other existing
architectures into a more general setting, making their analysis easier. For
instance, Section 2 presents a strategy to study the approximation
capabilities of some structured networks. Finally, we highlight the
flexibility and the benefits of this framework in Section 3, where we show
that to obtain expressive and robust neural networks, one can also include
layers that are not 1-Lipschitz.

There are multiple situations where one could be interested in networks with
some prescribed property. We report three of them here, where we refer to $F$
as the function to approximate:

  1. 1.

When $F$ has some known characterising property, e.g. $F$ is known to be
symplectic, compare Section 4.

  2. 2.

When the data we process has a particular structure, e.g. vectors whose
entries sum to one, as we present in Section 4.

  3. 3.

When we can approximate $F$ to sufficiently high accuracy with functions in
$\mathcal{G}$, a space that is well suited to model the layers of a network.
An example is using the space $\mathcal{G}$ of 1-Lipschitz functions to define
a classifier robust to adversarial attacks; see Section 3.

Thus, there are various applications where having neural networks structured
in a particular way is desirable. We will delve deeper into some of them in
the following sections. To be precise, we remark that all the properties we
focus on are preserved under composition, such as being 1-Lipschitz or
symplectic.

The paper is structured in five sections. First, in Section 2 we investigate
the universal approximation capabilities of some neural networks, thanks to
vector field decompositions, splitting methods and an embedding of the
dynamics into larger dimensional spaces. We then move, in Section 3, to a
neural network that has the property of being 1-Lipschitz. After the
mathematical derivation of the architecture, we present some numerical
experiments on adversarial robustness for the CIFAR-10 image classification
problem. Next, in Section 4, we introduce other neural networks with specific
designs. This last section aims to present in a systematic way how one can
impose certain properties on the architecture. We finally conclude the paper
in Section 5, mentioning some promising directions for further work.

Before moving on, we now report a numerical experiment that motivates our
investigation of structured neural networks. The results highlight how
imposing a structure does not have to degrade the approximation’s quality
considerably. Furthermore, this experiment suggests that not all the
constraining strategies perform similarly. Thus, a systematic process to
impose structure is essential since it allows changing the architecture in a
guided manner while preserving the property of interest.

###  1.1 Classification of points in the plane

We present a numerical experiment for the classification problem of the
dataset in Fig. 1(b). We consider neural networks that are 1-Lipschitz, as in
Section 3. We define the network layers alternatingly as contractive flow
maps, whose vector fields belong to
$\mathcal{F}_{c}=\\{-A^{T}\Sigma(Ax+b):\,A^{T}A=I\\}$, and as flows of other
Lipschitz continuous vector fields in
$\mathcal{F}=\\{\Sigma(Ax+b):\,A^{T}A=I\\}$, with
$\Sigma(z)=[\sigma(z_{1}),\ldots,\sigma(z_{n})]$ and
$\sigma(s)=\max\\{s,s/2\\}$111To impose the weight orthogonality, we set
$A=\mathrm{expm}(W-W^{T})$ with $\mathrm{expm}$ being the matrix exponential
and $W$ a trainable matrix.. The time steps for each vector field are network
parameters, and we constrain them to get a 1-Lipschitz network, see Section 3.
We report the results in Fig. 1(a) and Table 1.

The average classification test accuracy and final integration time, in
combination, get better by combining $\mathcal{F}_{c}$ with $\mathcal{F}$
instead of considering $\mathcal{F}_{c}$ alone. In particular, we see that the
final integration time $T$ with $\mathcal{F}_{c}\cup\mathcal{F}$ is the
smallest without significantly compromising the accuracy. The parameter $T$
measures how complex the network is; hence we get a more natural and efficient
solution by alternating the vector fields. In Section 2 we reinforce this
empirical result, proving results about theoretical approximation guarantees.
This renders the possibility of obtaining neural networks with prescribed
properties without compromising their approximation capabilities.

Adopted family of vector fields | Median accuracy | Median of $T$
---|---|---
$\mathcal{F}\cup\mathcal{F}_{c}$ | 98% | 1.84
$\mathcal{F}$ | 99% | 7.53
$\mathcal{F}_{c}$ | 97.33% | 19.82
Table 1: We perform 100 experiments alternating vector fields in
$\mathcal{F}_{c}$ with those in $\mathcal{F}$, 100 using just vector fields in
$\mathcal{F}_{c}$, and 100 with only those in $\mathcal{F}$. We work with
networks with ten residual layers throughout the experiments. In the table we
report the median final time $T$ and test accuracy for the three set of
experiments analysed

(a) Mean length of the time intervals along which the two dynamical regimes
are active.

(b) Dataset studied for the classification problem

Figure 1: Results from the experiments alternating the vector fields of
$\mathcal{F}$ and those of $\mathcal{F}_{c}$, together with the dataset of
interest

##  2 Universal approximation properties

As introduced before, neural network layers can be modelled by discretising
ordinary differential equations. In particular, this ODE-based approach can
also be beneficial for imposing some structure on neural networks and
providing a better theoretical understanding of their properties. In this
section, we follow this principle and prove the universal approximation
capabilities of two neural network architectures. Starting with the
continuous-in-time interpretation of neural networks, many approaches are
possible to prove such approximation properties, often based on merging the
broad range of results from dynamical systems theory and numerical analysis.
One can proceed, for example, in a constructive way as done in [50], where the
authors investigate the dynamics of some neural networks and explicitly
construct solutions to the problem of approximating a target function. Another
possibility is to study the approximation capabilities of compositions of flow
maps, as done in [35]. In this section, we focus on two solutions that, to the
best of our knowledge, are new. The first result is based on a vector field
decomposition, while the second is based on embedding vector fields into
larger dimensional spaces.

The approximation results that we cover rely on approximating vector fields
arbitrarily well. Consequently, this allows us to approximate their flow maps
accurately. This is based on the fact that for a sufficiently regular vector
field $X\in\operatorname{Lip}(\mathbb{R}^{n},\mathbb{R}^{n})$, if
$\tilde{X}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is such that for every
$x\in\mathbb{R}^{n}$

(3) $\|X(x)-\tilde{X}(x)\|<\varepsilon,$

then also their flow maps are close to one another for finite time intervals.
We formalise this reasoning in Proposition 2.1. In the proposition and
throughout the paper, we denote with $\Phi^{t}_{X}(z)$ the time-$t$ flow of
the vector field $X$, applied to $z$.

######  Proposition 2.1.

Let $X\in\operatorname{Lip}(\mathbb{R}^{n},\mathbb{R}^{n})$ and
$\tilde{X}\in\operatorname{Lip}(\mathbb{R}^{n},\mathbb{R}^{n})$ be as in Eq.
3. Then $\|\Phi_{X}^{t}(x)-\Phi_{\tilde{X}}^{t}(x)\|\leq\varepsilon
t\exp{\left(\operatorname{Lip}(X)t\right)},$ where $\operatorname{Lip}(X)$ is
the Lipschitz constant of $X$.

######  Proof 2.2.

We consider the integral equations associated to the ODEs $\dot{x}(t)=X(x(t))$
and $\dot{\tilde{x}}(t)=\tilde{X}(\tilde{x}(t))$ and study the difference of
their solutions both with the same initial condition $x\in\mathbb{R}^{n}$

$\displaystyle\|\Phi^{t}_{X}(x)-\Phi^{t}_{\tilde{X}}(x)\|$
$\displaystyle=\left\|x+\int_{0}^{t}X(\Phi_{X}^{s}(x))\,\mathrm{d}s-x-\int_{0}^{t}\tilde{X}(\Phi_{\tilde{X}}^{s}(x))\,\mathrm{d}s\right\|$
$\displaystyle\leq\int_{0}^{t}\left\|X(\Phi_{X}^{s}(x))-\tilde{X}(\Phi_{\tilde{X}}^{s}(x))\right\|\,\mathrm{d}s$
$\displaystyle=\int_{0}^{t}\left\|X(\Phi_{X}^{s}(x))-X(\Phi_{\tilde{X}}^{s}(x))+X(\Phi_{\tilde{X}}^{s}(x))-\tilde{X}(\Phi_{\tilde{X}}^{s}(x))\right\|\,\mathrm{d}s$
$\displaystyle\leq\operatorname{Lip}(X)\int_{0}^{t}\|\Phi_{X}^{s}(x)-\Phi_{\tilde{X}}^{s}(x)\|\,\mathrm{d}s+\varepsilon
t.$

Then we conclude that

$\|\Phi^{t}_{X}(x)-\Phi^{t}_{\tilde{X}}(x)\|\leq\varepsilon
t\exp{\left(\operatorname{Lip}(X)t\right)}.$

applying Gronwall’s inequality.

A particular consequence of this proposition is that if for every
$\varepsilon>0$ there is an $\tilde{X}\in\mathcal{F}$ making Eq. 3 true, then
we can approximate the $T-$flow map of $X$ arbitrarily well using elements of
$\mathcal{F}$:

$\|\Phi^{T}_{X}(x)-\Phi^{T}_{\tilde{X}}(x)\|\leq\varepsilon
T\exp(\operatorname{Lip}(X)T)=c\varepsilon.$

Because of this result, we now derive two approximation results for neural
networks working at the level of modelling vector fields.

###  2.1 Approximation based on a vector field decomposition

We now aim to show that, for a particularly designed neural network, we can
approximate arbitrarily well any continuous function in the $L^{p}$ norm and
any differentiable invertible function in the supremum norm on compact sets.
We also mention how to extend this last result to generic continuous
functions.

######  Theorem 2.3.

Let $F:\Omega\subset\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ be a continuous
function, with $\Omega\subset\mathbb{R}^{n}$ a compact set. Suppose that it
can be approximated, with respect to some norm $\|\cdot\|$, by a composition
of flow maps of $\mathcal{C}^{1}(\Omega,\mathbb{R}^{n})$ vector fields, i.e.
for any $\varepsilon>0$, $\exists
f_{1},\ldots,f_{k}\in\mathcal{C}^{1}(\Omega,\mathbb{R}^{n})$, such that

(4)
$\|F-\Phi_{f_{k}}^{h_{k}}\circ\ldots\circ\Phi_{f_{1}}^{h_{1}}\|<\varepsilon.$

Then, $F$ can be approximated arbitrarily well by composing flow maps of
gradient and sphere preserving vector fields, i.e. $\|F-\Phi_{\nabla
U^{k}}^{h_{k}}\circ\Phi_{X_{S}^{k}}^{h_{k}}\circ\ldots\circ\Phi_{\nabla
U^{1}}^{h_{1}}\circ\Phi_{X_{S}^{1}}^{h_{1}}\|<\varepsilon$.

By sphere preserving vector field, we mean a vector field $X_{S}$ having
$z^{T}z$ as a first integral, i.e. such that $z^{T}X_{S}(z)=0$ for any
$z\in\mathbb{R}^{n}$.

The norm $\|\cdot\|$ in Eq. 4 can be any norm that is well defined for
functions in $\mathcal{C}^{1}(\Omega,\mathbb{R}^{n})$. Two typical choices in
the literature are $L^{p}$ norms and the supremum norm

(5)
$\|F-\Phi_{f_{k}}^{h_{k}}\circ\ldots\circ\Phi_{f_{1}}^{h_{1}}\|:=\sup_{x\in\Omega}\|F(x)-\Phi_{f_{k}}^{h_{k}}\circ\ldots\circ\Phi_{f_{1}}^{h_{1}}(x)\|.$

Various works, like [10] and [35], have already proven the existence of vector
fields $f_{1},\ldots,f_{k}$ making Eq. 4 true when $\|\cdot\|$ is the $L^{p}$
norm and $F$ is a continuous function. Regarding the validity of hypothesis
Eq. 4 with the norm defined in Eq. 5, we mention [58] where the authors have
proven that if $F$ is a smooth invertible map with smooth inverse, then the
existence of $f_{1},\ldots,f_{k}$ can be guaranteed.

Theorem 2.3 is a consequence of the Presnov decomposition of vector fields,
introduced in [47], applied to the $k$ vector fields
$f_{1},\cdots,f_{k}\in\mathcal{C}^{1}(\Omega,\mathbb{R}^{n})$ in Eq. 4. The
Presnov decomposition is indeed a global decomposition of
$\mathcal{C}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector fields into the sum of
a gradient and a sphere preserving vector field. We now prove Theorem Theorem
2.3, and we specialise it to the subfamilies of vector fields we implement to
define neural networks.

######  Proof 2.4.

The vector fields $f_{1},\ldots,f_{k}$ are supposed to be continuously
differentiable. Thus, they all admit a unique Presnov decomposition, i.e. they
can be written as

$f_{i}(z)=\nabla U^{i}(z)+X^{i}_{S}(z),$

for a scalar function $U_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$, with
$U_{i}(0)=0$, and a sphere preserving vector field $X^{i}_{S}$. In general the
two vector fields $\nabla U^{i}(z)$ and $X_{S}^{i}(z)$ do not commute, i.e.
the Jacobi-Lie bracket $[\nabla U^{i},X_{S}^{i}]$ is not identically zero.
However, because of the Baker-Campbell-Hausdorff formula (see e.g. [26,
Section III.4.2]), as in splitting methods (see e.g. [42]) we can proceed with
an approximation of the form

$\Phi^{h}_{f_{i}}=\Phi_{\frac{h}{2}}^{\nabla
U^{i}}\circ\Phi_{h}^{X_{S}^{i}}\circ\Phi_{\frac{h}{2}}^{\nabla
U^{i}}+\mathcal{O}(h^{3}).$

We recall that
$\Phi_{f_{i}}^{h}=\Phi_{f_{i}}^{\alpha_{n}h}\circ\ldots\Phi_{f_{i}}^{\alpha_{1}h}$
for any tuple $\alpha=(\alpha_{1},\ldots,\alpha_{n})$ such that
$\sum_{i=1}^{n}\alpha_{i}=1$. Thus, up to choosing the $\alpha_{i}$s small
enough, we can approximate as accurately as desired $\Phi_{f_{i}}^{h}$ with
composition of flow maps of sphere preserving and gradient vector fields. This
concludes the proof.

Similar reasoning can be extended to other vector field decompositions, e.g.
the Helmholtz decomposition, as long as $f_{1},\ldots,f_{k}$ admit such a
decomposition. In Section 3, we adopt gradient vector fields whose flow maps
expand and contract distances to obtain 1-Lipschitz neural networks. We now
specialise Theorem 2.3 to the vector fields we use to model such neural
networks.

######  Corollary 2.5.

Consider the same assumptions of Theorem 2.3, and in particular the inequality
Eq. 4. Then, we can approximate $F$ arbitrarily well by composing flow maps of
expansive, contractive and sphere-preserving vector fields.

We first remark that with an expansive vector field we mean a vector field $X$
such that $\|\Phi^{t}_{X}(x)-\Phi^{t}_{X}(y)\|>\|x-y\|$ for any $t>0$, while
by contractive we mean that $\|\Phi^{t}_{X}(x)-\Phi^{t}_{X}(y)\|<\|x-y\|$. To
prove the corollary, we rely on a classical universal approximation theorem
with non-polynomial activation functions (see e.g. [46]). For clarity, we
report it here.

######  Theorem 2.6 (Universal approximation, [46]).

Let $\Omega\subset\mathbb{R}^{n}$ be a compact set and
$U\in\mathcal{C}^{1}(\mathbb{R}^{n})$. Assume
$\gamma\in\mathcal{C}^{1}(\mathbb{R})$ and $\gamma$ is not a polynomial. Then
for every $\varepsilon>0$ there is

$\tilde{U}(x)=\boldsymbol{\alpha}^{T}\Gamma(Ax+b),\quad\Gamma(z)=[\gamma(z_{1}),\ldots,\gamma(z_{n})],$

such that $\sup_{x\in\Omega}\|\tilde{U}(x)-U(x)\|<\varepsilon$, and
$\sup_{x\in\Omega}\|\nabla\tilde{U}(x)-\nabla U(x)\|<\varepsilon.$

We now prove Corollary 2.5.

######  Proof 2.7.

The proof follows the same reasoning of the one of Theorem 2.3. Indeed, we
first decompose each of the $f_{1},\ldots,f_{k}$ of equation Eq. 4 via the
Presnov decomposition as $f_{i}(z)=\nabla U^{i}(z)+X_{S}^{i}(z)$. Then, we
approximate each of the $U^{i}$ functions thanks to Theorem 2.6. To ease the
notation, we focus on one of the $f_{i}$ and denote it with $f$ from now on in
the proof.

Let $U:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and $X_{S}$ be so that
$f(z)=\nabla U(z)+X_{S}(z)$. Choose then $\sigma(x)=\max\\{ax,x\\}$,
$a\in(0,1)$, and $\gamma(x)=\int_{0}^{x}\sigma(s)\,\mathrm{d}s$. Since
$\gamma$ is not a polynomial and it is continuously differentiable, Theorem
2.6 for any $\varepsilon>0$ ensures the existence of a function

$\tilde{U}(z)=\alpha^{T}\Gamma(Ax+b),$

that satisfies $\sup_{x\in\Omega}\|U(x)-\tilde{U}(z)\|<\varepsilon$ and
$\sup_{x\in\Omega}\|\nabla U(x)-\nabla\tilde{U}(z)\|<\varepsilon$. We now
split $\nabla\tilde{U}(z)=A^{T}\operatorname{diag}(\alpha)\Sigma(Az+b)$ into a
contractive and an expansive part, exploiting the two following properties of
$\sigma$ and $\gamma$:

  1. 1.

$\sigma$ is positively homogeneous, i.e. $\sigma(\lambda s)=\lambda\sigma(s)$
for $\lambda,s\in\mathbb{R}$, $\lambda\geq 0$,

  2. 2.

$\gamma$ is strongly convex.

We decompose $\alpha$ as $\alpha^{+}-\alpha^{-}$, where
$(\alpha^{+})_{k}=\max\\{0,\alpha_{k}\\}$,
$(\alpha^{-})=-\min\\{0,\alpha_{k}\\}$ with $k=1,\ldots,n$. Because of the
positive homogeneity, $\nabla\tilde{U}(z)$ can be rewritten as

$\nabla\tilde{U}(z)=A_{1}^{T}\Sigma(A_{1}z+b_{1})-A_{2}^{T}\Sigma(A_{2}z+b_{2})=X_{E}(z)+X_{C}(z)$

where

$A_{1}=\operatorname{diag}(\alpha^{+})^{\frac{1}{2}}A,\,A_{2}=\operatorname{diag}(\alpha^{-})^{\frac{1}{2}}A,\,b_{1}=\operatorname{diag}(\alpha^{+})^{\frac{1}{2}}b,\,b_{2}=\operatorname{diag}(\alpha^{-})^{\frac{1}{2}}b.$

Because of the strong convexity of $\gamma$, we have

$\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left\|z(t)-y(t)\right\|^{2}=\langle
X_{E}(z(t))-X_{E}(y(t)),z(t)-y(t)\rangle>0$

with $z(t)=\Phi_{X_{E}}^{t}(z_{0})$ and $y(t)=\Phi_{X_{E}}^{t}(y_{0})$. This
means that the flow of $X_{E}$ is an expansive map. A similar reasoning shows
that $X_{C}$ has a contractive flow map. We can now conclude as in Theorem 2.3
since we have shown that every $f_{i}$ in Eq. 4 can be approximated
arbitrarily well as $f_{i}(z)\approx X_{E}^{i}(z)+X_{C}^{i}(z)+X_{S}^{i}(z)$.

As for the expansive and the contractive vector fields, to define neural
networks based on Corollary 2.5 one needs to parameterise the vector field
$X_{S}^{i}(z)$ that preserves spheres. Many possibilities are available, and
we report a couple of them. The first is

$\tilde{X}_{S}(z)=P(z)B^{T}\Sigma(Cz+d),\,\,B,C\in\mathbb{R}^{m\times
n},\,d\in\mathbb{R}^{m},\,P(z)=I_{n}-\frac{zz^{T}}{\|z\|^{2}},$

where $P(z):T_{z}\mathbb{R}^{n}\rightarrow T_{z}S^{2}_{\|z\|}$ is the
orthogonal projection on the space $\langle z\rangle^{\perp}$ and
$I_{n}\in\mathbb{R}^{n\times n}$ is the identity matrix. Another option is

$\tilde{X}_{S}(z)=\Lambda(z,\theta)z$

where $\Lambda(z,\theta)=A(z,\theta)-A(z,\theta)^{T}\in\mathbb{R}^{n\times n}$
with $A$ being a strictly upper triangular matrix whose entries are modelled
by $B\Sigma(Cx+b)\in\mathbb{R}^{N}$, $B\in\mathbb{R}^{N\times m}$,
$C\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$, $N=\frac{n(n-1)}{2}$.
These two possibilities allow us to approximate any sphere preserving vector
field arbitrarily because of classical universal approximation results, like
the one mentioned in Theorem 2.6. We prefer, for practical reasons, the second
one in the experiments reported in Appendix A. On the other hand, in the
experiment reported in Section 3 we replace the part
$X_{S}^{i}(z)+X_{E}^{i}(z)$ of the vector field decomposition with vector
fields of the form $\Sigma(Az+b)$. This choice gives good experimental
performances. The primary motivation behind this replacement is that, in those
experiments, the linear transformations of the inputs are based on convolution
operations, and the extension of sphere preserving vector fields to that
setting brings to involved architectures.

We now summarise the results presented in the context of neural networks.
Suppose that
$\|F-\Phi_{f_{k}}^{h_{k}}\circ\ldots\circ\Phi_{f_{1}}^{h_{1}}\|<\varepsilon$
and that $f_{i}\approx\tilde{f}_{i}=X^{i}_{C}+X^{i}_{E}+\tilde{X}_{S}^{i}$ for
$i=1,\ldots,k$. In Theorem 2.3 we have worked with the exact flows of the
vector fields. However, most of the times these are not available and hence a
numerical approximation is needed. This is exactly equivalent to applying a
splitting numerical integrator (see e.g. [26, Chapter 2] or [42]) to
approximate the $h_{i}$-flow map of $\tilde{f}_{i}$ (and hence also of
$f_{i}$) and get

(6)
$F(x)\approx\mathcal{N}(x)=\Psi^{h_{k}}_{X^{k}_{C}}\circ\Psi^{h_{k}}_{X^{k}_{E}}\circ\Psi^{h_{k}}_{g_{S}^{k}}\circ\ldots\circ\Psi^{h_{1}}_{X^{1}_{C}}\circ\Psi^{h_{1}}_{X^{1}_{E}}\circ\Psi^{h_{1}}_{g_{S}^{1}}(x).$

Here we denote with $\Psi^{h}_{f}$ a discrete approximation of the exact flow
$\Phi^{h}_{f}$ and $\mathcal{N}$ is the neural network that approximates the
target function $F$. Because of Corollary 2.5 and basic theory of numerical
methods for ODEs, $\mathcal{N}$ hence can approximate arbitrarily well $F$ in
the norm $\|\cdot\|$.

We remark that the neural network $\mathcal{N}$ defined in Eq. 6 does not
change the dimensionality of the input point, i.e. it is a map from
$\mathbb{R}^{n}$ to itself. However, usually Residual Neural Networks allow
for dimensionality changes thanks to linear lifting and projection layers. One
can extend all the results presented in this section to the dimensionality
changes, where instead of defining the whole network as the composition of
discrete flow maps, just the “dynamical blocks” are characterised in that way,
as represented in Fig. 2.

Figure 2: Representation of a ResNet made of two dynamical blocks, two lifting
layers $\mathcal{L}_{1}$, $\mathcal{L}_{2}$ and a final projection layer
$\mathcal{P}$.

Consequently, one can extend the results presented in [67]. In particular, one
can show that by composing flow maps of sphere-preserving and gradient vector
fields, generic continuous functions can also be approximated in the sense of
Eq. 5.

In Appendix A we show some numerical experiments where some unknown dynamical
systems and some target functions are approximated starting from the above
results. We now introduce another way to get expressivity results starting
from the continuous-in-time interpretation of neural networks.

###  2.2 Approximation based on Hamiltonian vector fields

Augmenting the dimensionality of the space where the dynamics is defined is a
typical technique for designing deep neural networks, see Fig. 2. Based on
this idea, we now study the approximation properties of networks obtained by
composing flow maps of Hamiltonian systems. For an introductory presentation
of Hamiltonian systems, see [33].

We now show that for any function $F$ for which hypothesis Eq. 4 holds, one
can approximate $F$ arbitrarily well, in the same function norm, by composing
flow maps of Hamiltonian systems. Consequently, symplectomorphisms like those
defined by SympNets ([30]) also approximate $F$ arbitrarily well.

This result relies on the embedding of a vector field
$f\in\mathcal{C}^{1}(\mathbb{R}^{n},\mathbb{R}^{n})$ into a Hamiltonian vector
field on $\mathbb{R}^{2n}$. To do so, we first define the linear map
$L:\mathbb{R}^{n}\rightarrow\mathbb{R}^{2n}$, as $z\mapsto L(z)=(z,0).$ We
then introduce the function $H_{f}(z,p)=p^{T}f(z)$, where $p$ is the conjugate
momentum of $z$. The gradient of such a function is

$\nabla
H_{f}(z,p)=\begin{bmatrix}\frac{\partial\left[p^{T}f(z)\right]}{\partial
z}&f(z)\end{bmatrix}.$

This implies that the Hamiltonian ODEs associated to $H_{f}$ are

$\begin{bmatrix}\dot{z}\\\
\dot{p}\end{bmatrix}=X_{H_{f}}(z,p)=\mathbb{J}\nabla
H_{f}(z,p)=\begin{bmatrix}f(z)\\\
-\frac{\partial\left[p^{T}f(z)\right]}{\partial z}\end{bmatrix}.$

Hence, we have $\Phi_{f}^{h}=P\circ\Phi_{X_{H_{f}}}^{h}\circ L$ where $P$ is
the projection on the first component $\mathbb{R}^{2n}\ni(z,p)\mapsto
z\in\mathbb{R}^{n}$. This construction, together with the fact that $P\circ
L=\mathrm{id}_{\mathbb{R}^{n}}$, implies that

$\|F-\Phi_{f_{k}}^{h_{k}}\circ\ldots\circ\Phi_{f_{1}}^{h_{1}}\|<\varepsilon\implies\|F-P\circ\Phi_{X_{H_{f_{k}}}}^{h_{k}}\circ\ldots\circ\Phi_{X_{H_{f_{1}}}}^{h_{1}}\circ
L\|<\varepsilon.$

One could be interested in such a lifting procedure and hence work with
Hamiltonian systems because discretising their flow maps with symplectic
methods might generate more stable networks or, as highlighted in [23], can
prevent vanishing and exploding gradient problems.

##  3 Adversarial robustness and Lipschitz neural networks

In this section, we consider the problem of classifying points of a set
$\mathcal{X}\subset\mathbb{R}^{n}$. More precisely, given a set
$\mathcal{X}=\cup_{i=1}^{C}\mathcal{X}_{i}$ defined by the disjoint union of
$C$ subsets $\mathcal{X}_{1},\ldots,\mathcal{X}_{C}$, we aim to approximate
the function $\ell:\mathcal{X}\rightarrow\\{1,\ldots,C\\}$ that assigns all
the points of $\mathcal{X}$ to their correct class, i.e. $\ell(x)=i$ for all
$x\in\mathcal{X}_{i}$ and all $i=1,\ldots,C$. Because of their approximation
properties, one can often choose neural networks to solve classification
problems, i.e. as models that approximate the labelling function $\ell$. On
the other hand, there is increasing evidence that trained neural networks are
sensitive to well-chosen input perturbations called adversarial attacks. The
first work that points this out is [57] and, since then, numerous others (see
e.g. [40, 14, 25]) have introduced both new ways to perturb the inputs
(attacks) and to reduce the sensitivity of the networks (defences). We first
formalise the problem of adversarial robustness from the mathematical point of
view and then derive a network architecture with inherent stability
properties.

Let $\mathcal{N}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{C}$ be a neural network
trained so that the true labelling map $\ell$ is well approximated by
$\hat{\ell}(x)=\arg\max_{i=1,\ldots,C}\mathcal{N}(x)_{i}$ for points
$x\in\mathcal{X}$. Furthermore, let us assume

$\|x-y\|\geq 2\varepsilon\quad\forall
x,y\in\mathcal{X},\,\,\ell(x)\neq\ell(y)$

for some norm $\|\cdot\|:\mathbb{R}^{n}\rightarrow\mathbb{R}^{+}$ defined on
the ambient space. With this setup, we say the network $\mathcal{N}$ is
$\varepsilon-$robust if

$\ell(x)=\hat{\ell}(x)=\hat{\ell}(x+\delta),\quad\forall\delta\in\mathbb{R}^{n},\,\,\|\delta\|<\varepsilon.$

In order to quantify the robustness of $\mathcal{N}$ we, first of all,
consider its Lipschitz constant $\operatorname{Lip}(\mathcal{N})$, i.e. the
smallest scalar value such that

$\|\mathcal{N}(x)-\mathcal{N}(y)\|_{2}\leq\operatorname{Lip}(\mathcal{N})\|x-y\|,\quad\forall
x,y\in\mathbb{R}^{n},$

where $\|\cdot\|_{2}:\mathbb{R}^{C}\rightarrow\mathbb{R}^{+}$ is the
$\ell^{2}$ norm. We also need a way to quantify how certain the network
predictions are. A measure of this certainty level is called margin in the
literature (see e.g. [3, 7, 61]) and it is defined as

$\mathcal{M}_{\mathcal{N}}(x)=\mathcal{N}(x)^{T}e_{\ell(x)}-\max_{j\neq\ell(x)}\mathcal{N}(x)^{T}e_{j},$

where $e_{i}$ is the $i-$th vector of the canonical basis of $\mathbb{R}^{C}$.
Combining these two quantities, in [61] the authors show that if the norm
$\|\cdot\|$ considered for $\mathcal{X}$ is the $\ell^{2}$ norm of the ambient
space $\mathbb{R}^{n}$, then

(7)
$\mathcal{M}_{\mathcal{N}}(x)\geq\sqrt{2}\varepsilon\operatorname{Lip}(\mathcal{N})\implies\mathcal{M}_{\mathcal{N}}(x+\delta)\geq
0\quad\forall\delta\in\mathbb{R}^{n},\,\|\delta\|\leq\varepsilon.$

Hence, for the points in $\mathcal{X}$ where Eq. 7 holds, the network is
robust to perturbations with a magnitude not greater than $\varepsilon$. This
result can be extended to generic $\ell^{p}$ metrics, but, in this section, we
focus on the case where $\|\cdot\|$ is the $\ell^{2}$ metric of
$\mathbb{R}^{n}$ and, from now on, we keep denoting it as $\|\cdot\|$.

Motivated by inequality Eq. 7, we present a strategy to constrain the
Lipschitz constant of Residual Neural Networks to the value of $1$.
Differently from [15, 66, 44], we impose such a property on the network
without relying only on layers that are 1-Lipschitz. This strategy relies on
the ODE based approach that we are presenting and is motivated by the interest
of getting networks that also have good expressivity capabilities. Indeed, we
remark that in Section 2 we studied the approximation properties of networks
similar to those we consider in this section. We conclude the section with
numerical experiments for the adversarial robustness with the CIFAR-10 dataset
to test the proposed network architecture.

###  3.1 Non–expansive dynamical blocks

Consider a scalar differentiable function
$V:\mathbb{R}^{n}\rightarrow\mathbb{R}$ that is also $\mu-$strongly convex,
for some $\mu>0$, i.e.

(8) $\langle\nabla V(x)-\nabla V(y),x-y\rangle\geq\mu\|x-y\|^{2},$

see e.g. [28, Chapter 6]. This said, it follows that the dynamics defined by
the ODE

(9) $\dot{x}(t)=-\nabla V(x(t))=X(x(t))$

is contractive, since

(10)
$\begin{split}\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\|x(t)-y(t)\|^{2}&=-\langle
x(t)-y(t),\nabla V(x(t))-\nabla V(y(t))\rangle\leq-\mu\|x(t)-y(t)\|^{2}\\\
\implies\|x(t)-y(t)\|&\leq e^{-\mu
t}\|x_{0}-y_{0}\|<\|x_{0}-y_{0}\|\,\,\forall\,t\geq 0,\end{split}$

where $x(t)=\Phi^{t}_{X}(x_{0})$ and $y(t)=\Phi^{t}_{X}(y_{0})$. A choice for
$V$ is $V(x)=\boldsymbol{1}^{T}\Gamma(Px+p)$, where
$\Gamma(x)=[\gamma(x_{1}),\ldots,\gamma(x_{n})]$,
$\gamma:\mathbb{R}\rightarrow\mathbb{R}$ is a strongly convex differentiable
function, and $\boldsymbol{1}=[1,\ldots,1]\in\mathbb{R}^{n}$. In this way, the
network we generate by concatenating explicit Euler steps applied to such
vector fields has layers of the type

$x\mapsto\Psi_{X}^{h}(x)=x-hP^{T}\Sigma(Px+p)$

where $\Sigma(x)=[\sigma(x_{1}),\ldots,\sigma(x_{n})]$ and
$\sigma(s)=\gamma^{\prime}(s)$.

If we discretise the ODE introduced above reproducing the non-expansive
behaviour at a discrete level, as presented for example in [20, 15, 44], we
get that the numerical flow $\Psi^{h}_{X}$ is non-expansive too. Consequently,
we can obtain 1-Lipschitz neural networks composing these non-expansive
discrete flow maps. A simple way to discretise Eq. 9 while preserving non-
expansivity is to use explicit-Euler steps with a small enough step size.
Indeed, assuming $\operatorname{Lip}(\sigma)\leq 1$, a layer of the form

$x\mapsto x-hP^{T}\Sigma(Px+p),\quad
h\leq\frac{2}{\|P\|^{2}},\,\,\|P\|=\sup_{\begin{subarray}{c}x\in\mathbb{R}^{n}\\\
\|x\|=1\end{subarray}}\|Px\|,$

is guaranteed to be 1-Lipschitz. We remark that, as highlighted in [15, 44],
it is not necessary to require strong convexity for $\gamma$ in order to make
$\Phi^{t}_{X}$ 1-Lipschitz. Indeed, it is enough to take $\gamma$ convex.
However, the strong convexity assumption allows us to include other layers
that are not 1-Lipschitz thanks to inequality Eq. 10.

We now shortly present the reasoning behind this statement. Consider another
ODE $\dot{x}(t)=Y(x(t))$ where $Y$ is again a vector field on $\mathbb{R}^{n}$
and suppose that $Y$ is $L$-Lipschitz. Then, we have that

$\|\Phi^{\bar{t}}_{Y}(x_{0})-\Phi^{\bar{t}}_{Y}(y_{0})\|\leq\exp{(L{\bar{t}})}\|x_{0}-y_{0}\|.$

This implies that, given $X$ as in Eq. 9, the map
$\Phi_{X}^{t}\circ\Phi_{Y}^{\bar{t}}=:C_{{\bar{t}},t}$ satisfies

(11)
$\|\Phi^{t}_{X}(\Phi^{\bar{t}}_{Y}(x_{0}))-\Phi^{t}_{X}(\Phi^{\bar{t}}_{Y}(y_{0}))\|\leq\exp\left(-\mu
t+L{\bar{t}}\right)\|x_{0}-y_{0}\|,$

so $C_{{\bar{t}},t}$ is Lipschitz continuous and will be 1-Lipschitz if
$\exp{(-\mu t+L{\bar{t}})}\leq 1$. This amounts to impose $L{\bar{t}}\leq\mu
t$ for the considered vector fields and time intervals on which corresponding
flow maps are active. The map $C_{{\bar{t}},t}$ can be seen as the exact
$(t+{\bar{t}})-$flow map of the switching system having a piecewise constant
(in time) autonomous dynamics. In particular, such a system coincides with $Y$
for the first time interval $[0,{\bar{t}})$ and with $X$ for the time interval
$[{\bar{t}},{\bar{t}}+t)$.

We could choose $Y$ to be the gradient vector field of a $L$-smooth scalar
potential. In other words, we ask for its gradient to be $L$-Lipschitz. Thus,
one possible way of building a dynamical block of layers that is 1-Lipschitz
is through a consistent discretisation of the switching system

(12) $\dot{x}(t)=(-1)^{s(t)}R_{s(t)}^{T}\Sigma(R_{s(t)}x(t)+r_{s(t)}),$

where $t\mapsto s(t)$ is a piecewise constant time-switching signal that,
following Eq. 11, balances the expansive and contractive regimes. In Eq. 12,
we are assuming that $\Sigma(z)=[\sigma(z_{1}),\ldots,\sigma(z_{n})]$ with
$\sigma$ Lipschitz and that $\gamma(s)=\int_{0}^{s}\sigma(t)dt$ is strongly
convex. In the numerical experiment we report at the end of the section we
design $s(t)$ giving an alternation between contractive and possibly non-
contractive behaviours. In the following subsection, we present two possible
approaches to discretise numerically the system in Eq. 12, mentioning how this
extends to more general families of vector fields.

In this subsection, we have worked to obtain dynamical blocks that are non-
expansive for the Euclidean metric. In Appendix B we show a way to extend this
reasoning to more general metrics defined in the input space.

###  3.2 Non–expansive numerical discretisation

As presented in the introductory section, it is not enough to have a
continuous model that satisfies some property of interest to get it at the
network level. Indeed, discretising the solutions to such an ODE must also be
done while preserving such a property. One approach that always works is to
restrict the step sizes to be sufficiently small so that the behaviour of the
discrete solution resembles the one of the exact solution. This strategy can
lead to expensive network training because of the high number of time steps.
On the other hand, this strategy allows weaker weight restrictions and better
performances. We remark how this translates for the dynamical system
introduced in Eq. 12, with $\sigma(x)=\max\\{x,ax\\}$. For that ODE, the one-
sided Lipschitz constant of contractive layers is
$\mu=a\lambda_{\min}(P^{T}P)$, $\lambda_{\min}$ being the smallest eigenvalue.
Thus, if $P$ is orthogonal, we get $\mu=a$. Under the same orthogonality
assumption, the expansive layers in Eq. 12 have Lipschitz constant $L=1$ and
this allows to specialise the non-expansivity condition Eq. 11 to
${\bar{t}}\leq at$. Thus, if we impose such a relationship and do sufficiently
small time steps, also the numerical solutions will be non expansive.

However, frequently smarter choices of discrete dynamical systems can lead to
leaner architectures and faster training procedures. We focus on the explicit
Euler method for this construction; however, one can work with other numerical
methods, like generic Runge–Kutta methods, as long as the conditions we derive
are adjusted accordingly. We concentrate on two time steps applied to equation
Eq. 12, but then the reasoning extends to every pair of composed discrete
flows and to other families of vector fields. Let

(13)
$\begin{split}\tilde{\Psi}^{h_{1}}(x)&=x-h_{1}P^{T}\Sigma(Px+p)=:x-h_{1}X(P,p,x)\\\
\Psi^{h_{2}}(x)&=x+h_{2}Q^{T}\Sigma(Qx+q)=:x+h_{2}X(Q,q,x).\end{split}$

The condition we want to have is that the map
$F_{h}=\tilde{\Psi}^{h_{1}}\circ\Psi^{h_{2}}$ is 1-Lipschitz, or at least that
this is true when $P$, $Q$ and $\Sigma$ satisfy some well-specified
properties. We first study the Lipschitz constant of both the discrete flow
maps and then upper bound the one of $F_{h}$ with their products. We take two
points $x,y\in\mathbb{R}^{n}$, define $\delta X(P,p,x,y)=X(P,p,y)-X(P,p,x)$,
and proceed as follows

$\displaystyle\|\tilde{\Psi}^{h_{1}}(y)-\tilde{\Psi}^{h_{1}}(x)\|^{2}$
$\displaystyle=\|y-x\|^{2}-2h_{1}\langle y-x,\delta
X(P,p,x,y)\rangle+h_{1}^{2}\|\delta X(P,p,x,y)\|^{2}$
$\displaystyle\leq_{(*)}\|y-x\|^{2}-2h_{1}\lambda_{\min}(P^{T}P)a\|y-x\|^{2}+h_{1}^{2}\|P\|^{4}\|y-x\|^{2},$

where $(*)$ is because we consider $\sigma=\max\\{ax,x\\}$. In the experiments
we present at the end of the section, we assume the weight matrices to be
orthogonal, hence $\lambda_{\min}(P^{T}P)=1$. Multiple works support this
weight constraint as a way to improve the generalisation capabilities, the
robustness to adversarial attacks and the weight efficiency (see e.g. [62,
60]). We will detail how we get such a constraint on convolutional layers in
the numerical experiments.

(a) Case $S=1$.

(b) Case $S=2$.

Figure 3: Representation of the non-expansivity region Eq. 15 for the choice
$a=0.5$.

The orthogonality of $P$ also implies $\|P\|=1$ and hence we get that a
Lipschitz constant of $\tilde{\Psi}^{h_{1}}$ is
$L_{1}=\sqrt{1-2h_{1}a+h_{1}^{2}}$. For the expansive flow map $\Psi^{h_{2}}$,
we have

(14)
$\begin{split}\|\Psi^{h_{2}}(y)-\Psi^{h_{1}}(x)\|&\leq\|x-y\|+h_{2}\operatorname{Lip}\left(Q^{T}\Sigma(Qz+q)\right)\|x-y\|\\\
&\leq(1+h_{2})\|x-y\|\end{split}$

under the orthogonality assumption for $Q$. The same result holds also if we
just have $\|Q\|\leq 1$. This leads to a region in the $(h_{1},h_{2})-$plane
where $L_{1}\cdot L_{2}\leq 1$ that can be characterised as follows

(15)
$\mathcal{R}=\\{(h_{1},h_{2})\in[0,1]^{2}:\;(1+h_{2})\sqrt{1-2h_{1}a+h_{1}^{2}}\leq
1\\}.$

This is represented in Fig. 3(a) for the case $a=0.5$, that is the one used
also in the numerical experiment for adversarial robustness. Thus, we now have
obtained a way to impose the 1-Lipschitz property on the network coming from
the discretisation of the ODE Eq. 12. It is clear that the result presented
here easily extends to different time-switching rules (i.e. a different choice
of $s(t)$), as long as there is the possibility of balancing expansive vector
fields with contractive ones. Furthermore, to enlarge the area in the
$(h_{1},h_{2})-$plane where non-expansivity can be obtained, one can decide to
divide into sub-intervals the time intervals $[0,h_{1}]$ and $[0,h_{2}]$.
Doing smaller steps, the allowed area increases. Indeed, instead of doing two
single time steps of length $h_{1}$ and $h_{2}$, one can perform $S$ time-
steps all of step-length $h_{1}/S$ or $h_{2}/S$. Thus, replacing
$\bar{h}_{1}=h_{1}/S$ and $\bar{h}_{2}=h_{2}/S$ into Eq. 15 it is immediate to
see that $h_{1}$ and $h_{2}$ are allowed to be larger than with the case
$S=1$. For example, if we again fix $a=0.5$ and set $S=2$, we get the area
represented in Fig. 3(b). The choice of $a=0.5$ and $S=2$ is the one we adopt
in the experiments reported in this section. We now conclude the section
showing how the derived architecture allows to improve the robustness against
adversarial attacks for the problem of image classification.

###  3.3 Numerical experiments with adversarial robustness

We now apply the reasoning presented above to the problem of classifying
images of the CIFAR-10 dataset. The implementation is done with PyTorch and is
available at the GitHub repository associated to the
paper222https://github.com/davidemurari/StructuredNeuralNetworks. We work with
convolutional neural networks, and with the activation function
$\sigma(x)=\max\left\\{x,\frac{x}{2}\right\\}$. The residual layers are
defined by dynamical blocks based on the discrete flow maps

$\displaystyle\tilde{\Psi}^{h_{1}}(x)$
$\displaystyle=x-h_{1}P^{T}\Sigma(Px+p)=:x-h_{1}X(P,p,x),\,\,P^{T}P=I$
$\displaystyle\Psi^{h_{2}}(x)$
$\displaystyle=x+h_{2}\Sigma(Qx+q)=:x+h_{2}Y(Q,q,x),\,\,Q^{T}Q=I$ (16)
$\displaystyle x$
$\displaystyle\mapsto\Psi^{h_{2}/2}\circ\tilde{\Psi}^{h_{1}/2}\circ\Psi^{h_{2}/2}\circ\tilde{\Psi}^{h_{1}/2}(x)$

We prefer this alternation regime rather than the one introduced in Eq. 13
because experimentally it gives better performances and it allows to generate
networks that are easier to train. We remark that all the calculations and
step restrictions presented before, apply with no changes to the alternation
regime in Section 3.3333This highlights the flexibility of the ODE based
approach to impose structure on the network. Indeed, this strategy just a
guided procedure to impose structure, but usually it does not limit
considerably the options one has available.. Indeed, what should change is
just the estimate in equation Eq. 14, but the Lipschitz constant of
$\Sigma(Qx+q)$ can be bounded by $1$ as in that equation as long as $\|Q\|\leq
1$. The step restriction is imposed after every training iteration, projecting
back the pairs $(h_{1},h_{2})$ in the region represented in Fig. 3(b) if
needed.

We implement an architecture that takes as inputs tensors of order three and
shape $3\times 32\times 32$. The first dimensionality of the tensor increases
to $32-64-128$ feature maps throughout the network via convolutional layers.
For each fixed number of filters, we have four layers of the form specified in
Section 3.3. To be precise, convolutional layers replace the matrix-vector
products. As a comparison, we also consider a baseline ResNet with the same
layer sequence and lifting pattern. The residual layers are of the form
$x\mapsto x+A^{T}\Sigma(Bx+b)$. This ResNet has a number of parameters that is
roughly the double of those on which we constrain the Lipschitz constant, i.e.
based on Section 3.3.

The network architecture based on Eq. 12 gets close to 90% test accuracy when
trained with cross-entropy loss and without weight constraints. However, as
presented at the beginning of this section, one could consider its Lipschitz
constant and its margin at any input point of interest to get robustness
guarantees for a network. For this reason, we now focus on constraining the
Lipschitz constant of the architecture, and we introduce the loss function we
adopt to promote higher margins. As in [3], we train the network architecture
with the multi-class hinge loss function defined as

$\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\sum\limits_{\begin{subarray}{c}j=1\\\
j\neq\ell(x_{i})\end{subarray}}^{10}\max\left\\{0,\text{margin}-\left(\mathcal{N}(x_{i})^{T}e_{\ell(x_{i})}-\mathcal{N}(x_{i})^{T}e_{j}\right)\right\\},$

where margin is a parameter to tune. In Fig. 4, we report how the robustness
of the architecture changes when we change such a margin value. We also train
the baseline ResNet with this loss function for a fair comparison. Having
predictions with higher margins allows us to get more robust architectures if
we fix the Lipschitz constant. Still, too high margins can lead to poor
accuracy. In Fig. 4, we report the results obtained constraining all the
dynamical blocks and the final projection layer. However, we do not constrain
the lifting layers. In this way, we can still control the full network’s
Lipschitz constant, just considering the norms of those lifting layers. On the
other hand, we leave some flexibility to the network, which can train better
also when we increase the hinge-loss margin. We notice that the dynamical
blocks usually get a small Lipschitz constant. Thus, even when we do not
constrain all the layers, the network will still be 1-Lipschitz.

To get orthogonal convolutional filters, we apply the regularisation strategy
proposed in [62]. This strategy is not the only one possible. Still, for this
experiment, we preferred it to more involved ones since the main focus has
been on the architecture, and the obtained results are satisfactory. Various
works, e.g. [34, 45, 22], highlight how one can directly constrain the
optimisation steps without having to project the weights on the right space or
to add regularisation terms. We have not experimented with these kinds of
strategies, and we leave them for future study. We also work with an
orthogonal initialisation for the convolutional layers. The lifting layers are
modelled as $x\mapsto\alpha Wx$ for a convolutional filter $W$ with
$\|W\|_{2}\leq 1$. To constrain the norm to 1, we add a projection step after
the stochastic gradient descent (SGD) method updates the weights, i.e. we
normalise the weights as $W\mapsto W/\max\\{1,\|W\|_{2}\\}$. Here, the
$2-$norm of the convolutional filters is computed with the power method as
described, for example, in [44]. Furthermore, we work with the SGD having a
learning rate scheduler that divides the learning rate by $10$ after a fixed
number of epochs. Finally, we generate the adversarial examples with the
library “Foolbox” introduced in [49]. We focus on the $L^{2}-$PGD attack and
perform ten steps of it. We test different magnitudes of the adversarial
perturbations.

In Fig. 4 we see that the robustness of the constrained neural networks
improves compared to the baseline ResNet, whatever the margin we choose.
Furthermore, the flexibility in the Lipschitz constant induced by the
unconstrained lifting layers makes the training procedure less sensitive to
the changes in the $\mathrm{margin}$ parameter. This experimental evidence
supports the trade-off between these two scalar values reported in equation
Eq. 7.

Figure 4: Plot of the accuracy against the magnitude of the adversarial PGD
perturbation of $1024$ test clean images, comparing various values for
$\mathrm{margin}$. All the networks have been trained with hinge-loss
function, for three different values of the margin parameter.

##  4 Imposition of other structure

Depending on the problem and the application where a neural network is
adopted, the properties that the architecture should satisfy may be very
different. We have seen in the previous section a strategy to impose Lipschitz
constraints on the network to get some guarantees in terms of adversarial
robustness. In that context, the property is of interest because it is
possible to see that even when imposing it, we can get sufficiently accurate
predictions. Moreover, this strategy allows controlling the network’s
sensitivity to input perturbations. As mentioned in the introduction, there
are at least two other situations where structural constraints might be
desirable. The first one is when one knows that the function to approximate
has some particular property. The second is when we know that the data we
process induces some structure on the function we want to approximate.

This section supports the claim that combining ODE models with suitable
numerical integrators allows us to define structured networks. More precisely,
we derive multiple architectures putting together these two elements. Some of
these have already been presented in other works, and others are new. The
properties that we investigate are symplecticity, volume preservation and mass
preservation. For the first two, we describe how to constrain the dynamical
blocks. For the third, we propose how to structure also the linear lifting and
projection layers. Moreover, for this latter example, we also report some
numerical experiments.

###  4.1 Symplectic dynamical blocks

A function $F:\mathbb{R}^{2n}\rightarrow\mathbb{R}^{2n}$ is said to be
symplectic if it satisfies the identity

$\frac{\partial F(x)}{\partial x}^{T}\mathbb{J}\frac{\partial F(x)}{\partial
x}=\mathbb{J}\quad\forall
x\in\mathbb{R}^{2n},\,\,\mathbb{J}=\begin{bmatrix}0_{n}&I_{n}\\\
-I_{n}&0_{n}\end{bmatrix}\in\mathbb{R}^{2n\times 2n}$

with $0_{n},I_{n}\in\mathbb{R}^{n\times n}$ being the zero and the identity
matrices. Symplectic maps are particularly important in classical mechanics,
because the flow map $\Phi^{t}$ of a Hamiltonian system
$\dot{x}(t)=\mathbb{J}\nabla H(x(t))$ is symplectic, see e.g. [33]. This fact
implies that if one is interested in approximating such a flow map with a
neural network, then structuring it to be symplectic might be desirable. In
this direction there are a considerable number of works (see e.g. [30, 19,
68]). We mention in particular [30] where the authors construct layers of a
network to ensure the symplectic property is satisfied. On the other hand, in
[19] the authors consider a neural network as the Hamiltonian function of a
system and approximate its flow map with a symplectic numerical
integrator444We remark that a one-step numerical method $\Psi^{h}$ is
symplectic if and only if, when applied to any Hamiltonian system, it is a
symplectic map.. The simplest symplectic integrator is symplectic-Euler, that
applied to $H(q,p)=V(q)+K(p)$ writes

$q_{n+1}=q_{n}+h\partial_{p}K(p_{n}),\quad
p_{n+1}=p_{n}-h\partial_{q}V(q_{n+1}).$

We now focus on the gradient modules presented in [30], defined alternating
maps of the form

$\mathcal{G}_{1}(q,p)=\begin{bmatrix}A^{T}\operatorname{diag}(\alpha)\Sigma(Ap+a)+q\\\
p\end{bmatrix}\quad\mathcal{G}_{2}(q,p)=\begin{bmatrix}q\\\
B^{T}\operatorname{diag}(\beta)\Sigma(Bq+b)+p\end{bmatrix}.$

We notice that we can obtain the same map from a time-switching ODE model. We
first introduce the time-dependent Hamiltonian of such a model that is

(17) $H_{s(t)}(z)=\alpha_{s(t)}^{T}\Gamma(A_{s(t)}P_{s(t)}z+b_{s(t)}),\quad
A_{s(t)}\in\mathbb{R}^{n\times n},\,b_{s(t)}\in\mathbb{R}^{n},$

with $s:[0,+\infty)\rightarrow\mathbb{R}^{+}$ being piecewise constant. We can
suppose without loss of generality that $s(t)\in\\{0,1,2,\ldots,K\\}$ and that
$P_{s(t)}$ alternates between $\Pi_{1}=[I_{n},0_{n}]\in\mathbb{R}^{n\times
2n}$ and $\Pi_{2}=[0_{n},I_{n}]\in\mathbb{R}^{n\times 2n}$. Let now
$\Gamma(z)=[\gamma(z_{1}),\ldots,\gamma(z_{n})]$,
$\Sigma(z)=[\sigma(z_{1}),\ldots,\sigma(z_{n})]$, and
$\gamma^{\prime}(s)=\sigma(s)$. We then notice that the Hamiltonian vector
field associated to $H_{s(t)}$ alternates between the following two vector
fields

$X_{H_{1}}=\begin{bmatrix}A_{1}^{T}\operatorname{diag}(\alpha_{1})\Sigma(A_{1}p+b_{1})&0\end{bmatrix}^{T},\quad
X_{H_{2}}=\begin{bmatrix}0&-A_{2}^{T}\operatorname{diag}(\alpha_{2})\Sigma(A_{2}q+b_{2})\end{bmatrix}^{T}.$

We now conclude that if we compute the exact flow of $X_{H_{s(t)}}$ and we
take $s(t)$ to be constant on every interval of length 1, we recover the
gradient module in [30].

Similarly, all the network architectures presented in [19] and related works
are based on defining a neural network $\mathcal{N}(q,p)$ that plays the role
of the Hamiltonian function and then applying a symplectic integrator to the
Hamiltonian system $\dot{z}=\mathbb{J}\nabla\mathcal{N}(z)$. The composition
of discrete flow maps of the a time independent Hamiltonian gives a symplectic
network with shared weights. On the other hand, if the Hamiltonian changes as
in Eq. 17, one gets different weights for different layers and potentially a
more expressive model as presented in [30].

###  4.2 Volume preserving dynamical blocks

Suppose one is interested in defining efficiently invertible and volume
preserving networks. In that case, an approach based on switching systems and
splitting methods can provide a flexible solution. Consider the switching
system defined by $\dot{z}(t)=f_{s(t)}(z(t))$,
$f_{i}\in\mathfrak{X}(\mathbb{R}^{n})$, where for every value of $s(t)$, the
vector field has a specific partitioning that makes it divergence-free. For
example, if we have $n=2m$,

$f_{s(t)}(z)=\begin{bmatrix}u_{s(t)}(z[m:])&v_{s(t)}(z[:m])\end{bmatrix}^{T}$

satisfies such a condition and its flow map will be volume-preserving. We can
numerically integrate such a vector field while preserving this property.
Indeed, we can apply a splitting method based on composing the exact flow maps
of the two-volume preserving vector fields

$f_{s(t)}^{1}(z)=\begin{bmatrix}u_{s(t)}(z[m:])&0\end{bmatrix}^{T},\quad
f_{s(t)}^{2}(z)=\begin{bmatrix}0&v_{s(t)}(z[:m])\end{bmatrix}^{T}.$

This approach gives architectures that are close to the ones of RevNets (see
e.g. [24]). The inverse of the network is efficient to compute in this case,
and this translates into memory efficiency since one does not need to save
intermediate activation values for the backpropagation. A particular class of
these blocks can be obtained with second-order vector fields and, in
particular, with second-order conservative vector fields (hence Hamiltonian):
$\ddot{x}(t)=f_{s(t)}(x(t)),\text{ or }\ddot{x}(t)=-\nabla V_{s(t)}(x(t))$,
where, for example,
$V_{s(t)}(x)=\boldsymbol{\alpha}^{T}\Gamma(A_{s(t)}x+b_{s(t)})$ for some
$\Gamma=[\gamma,\ldots,\gamma]$. With a similar strategy, one can also derive
the volume-preserving neural networks presented in [6].

###  4.3 Mass preserving neural networks

The final property we focus on is mass preservation. By mass preservation, we
refer to the conservation of the sum of the components of a vector (see [9]).
More explicitly, one could be interested in imposing such a structure if the
goal is to approximate a function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$
that is known to satisfy
$T_{n}x=\sum_{i=1}^{n}x_{i}=T_{m}F(x)=\sum_{i=1}^{m}F(x)_{i}$. A simple way to
impose such property is by approximating the target function
$F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ as

$F(x)\approx\frac{\sum_{i=1}^{n}x_{i}}{\sum_{j=1}^{m}\tilde{F}(x)_{j}}\tilde{F}(x)$

where $\tilde{F}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is any sufficiently
expressive neural network. However, this choice might lead to hard training
procedures because of the denominator. Imposing this structure at the level of
network layers is not so intuitive in general. Hence, we rely again on a
suitable ODE formulation. A vector field $X\in\mathfrak{X}(\mathbb{R}^{n})$
whose flow map preserves the sum of the components of the state vector is
simply one having a linear first integral
$g(x)=\boldsymbol{1}^{T}x=\sum_{i=1}^{n}x_{i}$. Thus, we can design vector
fields of the form

(18)
$\dot{y}(t)=(A(y)-A(y)^{T})\boldsymbol{1},\,\,A:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times
n},$

and this property will be then a natural consequence of the exact flow map.
This mass conservation could also be extended to a weighted-mass conservation,
and we would just have to replace $\boldsymbol{1}$ with a vector of weights
$\boldsymbol{\alpha}$. This extension does not however allow to change the
dimensionality from a layer to the next one as easily. To model these vector
fields, we can work with parametric functions like
$\tilde{f}(x)=B^{T}\Sigma(Ax+b)\in\mathbb{R}^{n(n-1)/2}$ and use them to build
the upper triangular matrix-valued function $A$ in Eq. 18. As presented in
[26, Chapter 4], it is also immediate to impose this property at a discrete
level since every Runge–Kutta or multistep method preserves linear first
integrals without time-step restrictions. Thus, a possible strategy to model
mass preserving neural networks is based on combining layers of the following
types:

  1. 1.

Lifting layers: $L:\mathbb{R}^{k}\rightarrow\mathbb{R}^{k+s}$,
$L(x_{1},\ldots,x_{k})=(x_{1},\ldots,x_{k},0,0,\ldots,0)$,

  2. 2.

Projection layers: $P:\mathbb{R}^{k+s}\rightarrow\mathbb{R}^{k}$,
$P(x_{1},\ldots,x_{k},x_{k+1},\ldots,x_{k+s})=(x_{1}+o,\ldots,x_{k}+o)$, with
$o=\sum_{i=1}^{s}x_{k+i}/s,$

  3. 3.

Dynamical blocks: one-step explicit Euler discretisations of Eq. 18.

To test the neural network architecture, we focus on the approximation of the
flow map of the SIR-model

(19)
$\dot{y}=\begin{bmatrix}-y_{1}y_{2}&y_{1}y_{2}-y_{2}&y_{2}\end{bmatrix}^{T}=X(y)^{T}.$

We model the neural network as discussed above. We approximate the $1-$flow
map of Eq. 19 working with pairs of the form
$\\{(x_{i},y_{i}=\Psi^{1}_{X}(x_{i}))\\}_{i=1,\ldots,N}$555Here with
$\Psi^{1}_{X}(x_{i})$ we refer to an accurate approximation of the time-1 flow
map of X applied to $x_{i}$. The plots in Fig. 5 represent the first two
components of the solution for the SIR model. All the line segments connect
the components of the initial conditions with those of the time-1 updates. The
considerable benefit of mass-preservation as a constraint is that it allows
interpretable outputs. . Indeed, in this case the components of $y$ represent
the percentages of three species on the total population, and the network we
train still allows to get this interpretation being mass-preserving.

Figure 5: Plots of the approximation of the time $1-$flow map of the SIR model
Eq. 19 for 10 test initial conditions. We report the first two components of
the solutions. Each point represents either an initial condition (at time 0),
or a time-1 update.

##  5 Conclusion and future directions

In this work, we have introduced a framework to combine the design of ODE
models with the choice of proper numerical methods to, in turn, obtain neural
networks with prescribed properties. After introducing and motivating the
approach, we proved two universal approximation theorems. The first one
relates to sphere preserving and gradient vector fields, while the second one
involves Hamiltonian vector fields. We then obtained Lipschitz constrained
ResNets, focusing mainly on how to introduce layers that are not 1-Lipschitz.
We then applied this construction to get neural networks with adversarial
robustness guarantees. Finally, to show the framework’s flexibility, we
demonstrated how to design dynamical blocks that are symplectic, volume
preserving and mass preserving.

Throughout the manuscript we have focused on explicit numerical methods as
tools to generate neural network architectures. However, many geometric
integrators are implicit (see e.g. [26]); thus, this remains a direction to
pursue in further work so that the framework can be extended to other
properties (see e.g. [5, 38]).

We have adopted the formalism of piecewise-autonomous dynamical systems to
design neural networks without heavily relying on the theory of time-switching
systems. However, switching systems are a well-studied research area (see e.g.
[36, 37, 2]), and it seems natural to study them further and their use to
design neural network architectures.

Finally, we remark that imposing properties on neural networks is a promising
strategy to make them more understandable, interpretable, and reliable. On the
other hand, it is also clear that constraining the architecture can
considerably decrease the network’s expressivity in some cases. Thus, it
remains to understand when it is preferable to replace hard constraints with
soft constraints promoting such properties without imposing them by
construction.

Acknowledgements The authors would like to thank the Isaac Newton Institute
for Mathematical Sciences, Cambridge, for support and hospitality during the
programmes Mathematics of deep learning and Geometry, compatibility and
structure preservation in computational differential equations where work on
this paper was significantly advanced. This work was supported by EPSRC grant
no EP/R014604/1. EC and BO have received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No 860124.

## References

  * [1] A. Agrachev and A. Sarychev, Control on the manifolds of mappings with a view to the deep learning, Journal of Dynamical and Control Systems, (2021), pp. 1–20.
  * [2] T. Alpcan and T. Basar, A stability result for switched systems with multiple equilibria, Dynamics of Continuous, Discrete and Impulsive Systems Series A: Mathematical Analysis, 17 (2010), pp. 949–958.
  * [3] C. Anil, J. Lucas, and R. Grosse, Sorting out Lipschitz function approximation, in International Conference on Machine Learning, PMLR, 2019, pp. 291–301.
  * [4] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb, Solving inverse problems using data-driven models, Acta Numerica, 28 (2019), pp. 1–174.
  * [5] S. Bai, J. Z. Kolter, and V. Koltun, Deep equilibrium models, Advances in Neural Information Processing Systems, 32 (2019).
  * [6] J. Bajārs, Locally-symplectic neural networks for learning volume-preserving dynamics, arXiv preprint arXiv:2109.09151, (2021).
  * [7] P. L. Bartlett, D. J. Foster, and M. J. Telgarsky, Spectrally-normalized margin bounds for neural networks, Advances in Neural Information Processing Systems, 30 (2017).
  * [8] J. Berner, P. Grohs, G. Kutyniok, and P. Petersen, The modern mathematics of deep learning, arXiv preprint arXiv:2105.04026, (2021).
  * [9] S. Blanes, A. Iserles, and S. Macnamara, Positivity-preserving methods for population models, arXiv preprint arXiv:2102.08242, (2021).
  * [10] Y. Brenier and W. Gangbo, $L^{p}$ Approximation of maps by diffeomorphisms, Calculus of Variations and Partial Differential Equations, 16 (2003), pp. 147–164.
  * [11] M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković, Geometric deep learning: Grids, groups, graphs, geodesics, and gauges, arXiv preprint arXiv:2104.13478, (2021).
  * [12] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst, Geometric deep learning: going beyond euclidean data, IEEE Signal Processing Magazine, 34 (2017), pp. 18–42.
  * [13] M. Calvo, D. Hernández-Abreu, J. I. Montijano, and L. Rández, On the preservation of invariants by explicit Runge–Kutta methods, SIAM Journal on Scientific Computing, 28 (2006), pp. 868–885.
  * [14] N. Carlini and D. Wagner, Towards evaluating the robustness of neural networks, in 2017 IEEE Symposium on Security and Privacy, IEEE, 2017, pp. 39–57.
  * [15] E. Celledoni, M. J. Ehrhardt, C. Etmann, R. I. McLachlan, B. Owren, C.-B. Schönlieb, and F. Sherry, Structure-preserving deep learning, European Journal of Applied Mathematics, 32 (2021), pp. 888–936.
  * [16] E. Celledoni, A. Leone, D. Murari, and B. Owren, Learning hamiltonians of constrained mechanical systems, Journal of Computational and Applied Mathematics, 417 (2023), p. 114608.
  * [17] B. Chang, M. Chen, E. Haber, and E. H. Chi, AntisymmetricRNN: A Dynamical System View on Recurrent Neural Networks, in International Conference on Learning Representations, 2019.
  * [18] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, Neural ordinary differential equations, Advances in Neural Information Processing Systems, 31 (2018).
  * [19] Z. Chen, J. Zhang, M. Arjovsky, and L. Bottou, Symplectic Recurrent Neural Networks, in International Conference on Learning Representations, 2020\.
  * [20] G. Dahlquist, Generalized disks of contractivity for explicit and implicit Runge-Kutta methods, tech. report, CM-P00069451, 1979.
  * [21] S. Eidnes, Order theory for discrete gradient methods, BIT Numerical Mathematics, (2022), pp. 1–49.
  * [22] G. França, A. Barp, M. Girolami, and M. I. Jordan, Optimization on manifolds: A symplectic approach, arXiv preprint arXiv:2107.11231, (2021).
  * [23] C. L. Galimberti, L. Furieri, L. Xu, and G. Ferrari-Trecate, Hamiltonian deep neural networks guaranteeing non-vanishing gradients by design, arXiv preprint arXiv:2105.13205, (2021).
  * [24] A. N. Gomez, M. Ren, R. Urtasun, and R. B. Grosse, The reversible residual network: Backpropagation without storing activations, Advances in Neural Information Processing Systems, 30 (2017).
  * [25] I. J. Goodfellow, J. Shlens, and C. Szegedy, Explaining and harnessing adversarial examples, in ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, Y. Bengio and Y. LeCun, eds., 2015.
  * [26] E. Haier, C. Lubich, and G. Wanner, Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations, Springer, 2006.
  * [27] J. Hertrich, S. Neumayer, and G. Steidl, Convolutional proximal neural networks and plug-and-play algorithms, Linear Algebra and its Applications, 631 (2021), pp. 203–234.
  * [28] J.-B. Hiriart-Urruty and C. Lemaréchal, Convex analysis and minimization algorithms I: Fundamentals, vol. 305, Springer Science & Business Media, 2013.
  * [29] K. Huang, Y. Wang, M. Tao, and T. Zhao, Why do deep residual networks generalize better than deep feedforward networks?—a neural tangent kernel perspective, Advances in neural information processing systems, 33 (2020), pp. 2698–2709.
  * [30] P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis, SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems, Neural Networks, 132 (2020), pp. 166–179.
  * [31] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko, et al., Highly accurate protein structure prediction with AlphaFold, Nature, 596 (2021), pp. 583–589.
  * [32] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, Backpropagation applied to handwritten zip code recognition, Neural computation, 1 (1989), pp. 541–551.
  * [33] B. Leimkuhler and S. Reich, Simulating Hamiltonian dynamics, no. 14, Cambridge University Press, 2004.
  * [34] B. Leimkuhler, T. J. Vlaar, T. Pouchon, and A. Storkey, Better Training using Weight-Constrained Stochastic Dynamics, in Proceedings of the 38th International Conference on Machine Learning, 2021, pp. 6200–6211.
  * [35] Q. Li, T. Lin, and Z. Shen, Deep learning via dynamical systems: An approximation perspective, Journal of the European Mathematical Society, (2022).
  * [36] D. Liberzon, Switching in systems and control, vol. 190, Springer, 2003\.
  * [37] D. Liberzon and A. S. Morse, Basic problems in stability and design of switched systems, IEEE Control Systems Magazine, 19 (1999), pp. 59–70.
  * [38] A. Look, S. Doneva, M. Kandemir, R. Gemulla, and J. Peters, Differentiable implicit layers, arXiv preprint arXiv:2010.07078, (2020).
  * [39] Y. Lu, A. Zhong, Q. Li, and B. Dong, Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations, in International Conference on Machine Learning, PMLR, 2018, pp. 3276–3285.
  * [40] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, Towards Deep Learning Models Resistant to Adversarial Attacks, in International Conference on Learning Representations, 2018.
  * [41] S. Massaroli, M. Poli, J. Park, A. Yamashita, and H. Asama, Dissecting neural odes, Advances in Neural Information Processing Systems, 33 (2020), pp. 3952–3963.
  * [42] R. I. McLachlan and G. R. W. Quispel, Splitting methods, Acta Numerica, 11 (2002), pp. 341–434.
  * [43] R. I. McLachlan, G. R. W. Quispel, and N. Robidoux, Geometric integration using discrete gradients, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 357 (1999), pp. 1021–1045.
  * [44] L. Meunier, B. Delattre, A. Araujo, and A. Allauzen, Scalable Lipschitz Residual Networks with Convex Potential Flows, arXiv preprint arXiv:2110.12690, (2021).
  * [45] M. Ozay and T. Okatani, Optimization on submanifolds of convolution kernels in CNNs, arXiv preprint arXiv:1610.07008, (2016).
  * [46] A. Pinkus, Approximation theory of the MLP model in neural networks, Acta Numerica, 8 (1999), pp. 143–195.
  * [47] E. Presnov, Non-local decomposition of vector fields, Chaos, Solitons & Fractals, 14 (2002), pp. 759–764.
  * [48] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics, 378 (2019), pp. 686–707.
  * [49] J. Rauber, W. Brendel, and M. Bethge, Foolbox: A python toolbox to benchmark the robustness of machine learning models, arXiv preprint arXiv:1707.04131, (2017).
  * [50] D. Ruiz-Balet and E. Zuazua, Neural ODE control for classification, approximation and transport, arXiv preprint arXiv:2104.05278, (2021).
  * [51] T. K. Rusch and S. Mishra, UnICORNN: A recurrent model for learning very long time dependencies, in International Conference on Machine Learning, PMLR, 2021, pp. 9168–9178.
  * [52] L. Ruthotto and E. Haber, Deep neural networks motivated by partial differential equations, Journal of Mathematical Imaging and Vision, 62 (2020), pp. 352–364.
  * [53] A. M. Saxe, J. L. McClelland, and S. Ganguli, Exact solutions to the nonlinear dynamics of learning in deep linear neural networks, CoRR, abs/1312.6120 (2014).
  * [54] R. Shwartz-Ziv and N. Tishby, Opening the black box of deep neural networks via information, ArXiv, abs/1703.00810 (2017).
  * [55] P. Simard, Y. LeCun, and J. Denker, Efficient pattern recognition using a new transformation distance, Advances in Neural Information Processing Systems, 5 (1992).
  * [56] B. Smets, J. Portegies, E. J. Bekkers, and R. Duits, PDE-based group equivariant convolutional neural networks, Journal of Mathematical Imaging and Vision, (2022), pp. 1–31.
  * [57] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. J. Goodfellow, and R. Fergus, Intriguing properties of neural networks, CoRR, abs/1312.6199 (2014).
  * [58] T. Teshima, K. Tojo, M. Ikeda, I. Ishikawa, and K. Oono, Universal approximation property of neural ordinary differential equations, arXiv preprint arXiv:2012.02414, (2020).
  * [59] M. Thorpe and Y. van Gennip, Deep limits of residual neural networks, arXiv preprint arXiv:1810.11741, (2018).
  * [60] A. Trockman and J. Z. Kolter, Orthogonalizing Convolutional Layers with the Cayley Transform, in International Conference on Learning Representations, 2021.
  * [61] Y. Tsuzuku, I. Sato, and M. Sugiyama, Lipschitz-margin training: Scalable certification of perturbation invariance for deep neural network, Advances in Neural Information Processing Systems, 31 (2018).
  * [62] J. Wang, Y. Chen, R. Chakraborty, and S. X. Yu, Orthogonal convolutional neural networks, in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 11505–11515.
  * [63] L. Wang, Y. Zhang, and J. Feng, On the Euclidean distance of images, IEEE Transactions on Pattern Analysis and Machine Intelligence, 27 (2005), pp. 1334–1339.
  * [64] E. Weinan, A proposal on machine learning via dynamical systems, Communications in Mathematics and Statistics, 1 (2017), pp. 1–11.
  * [65] L. Xiao, Y. Bahri, J. Sohl-Dickstein, S. Schoenholz, and J. Pennington, Dynamical isometry and a mean field theory of cnns: How to train 10,000-layer vanilla convolutional neural networks, in International Conference on Machine Learning, PMLR, 2018, pp. 5393–5402.
  * [66] M. Zakwan, L. Xu, and G. Ferrari-Trecate, On Robust Classification using Contractive Hamiltonian Neural ODEs, arXiv preprint arXiv:2203.11805, (2022).
  * [67] H. Zhang, X. Gao, J. Unterman, and T. Arodz, Approximation capabilities of neural odes and invertible residual networks, in International Conference on Machine Learning, PMLR, 2020, pp. 11086–11095.
  * [68] A. Zhu, P. Jin, and Y. Tang, Deep hamiltonian networks based on symplectic integrators, arXiv preprint arXiv:2004.13830, (2020).

##  Appendix A Some numerical experiments

We now apply the theoretical background introduced in Section 2 to two
approximation tasks. More precisely, we verify if the introduced architectures
are complicated to train in practice or if they can achieve good performances.
The two tasks of interest are the approximation of a continuous scalar
function and the approximation of a $\mathcal{C}^{1}$ vector field starting
from a set of training trajectories.

We start with the approximation of $f(x)=x^{2}+|x|+\sin{(x)}$,
$x\in\mathbb{R}$, and $g(x,y)=\sqrt{x^{2}+y^{2}}$. The goal here is to
approximate them by composing flow maps of vector fields that are structured
as

$X_{G}(z)=A^{T}\operatorname{diag}(\boldsymbol{\alpha})\Sigma(Az+b)=\nabla\left(\boldsymbol{\alpha}^{T}\Gamma(Az+b)\right),$
$X_{S}(z)=(A(z)-A(z)^{T})z.$

More precisely, following the result in Theorem 2.3, we compose the flow maps
of $X_{G}$ and $X_{S}$, maintaining the same time step for pairs of such flow
maps. We report the results in Fig. 6.

(a) Approximation of $f(x)$

(b) Approximation of $g(x,y)$

Figure 6: Comparison between the trained networks and the true functions.

The second experiment that we report, is the one of approximating a vector
field $X\in\mathfrak{X}(\mathbb{R}^{n})$ starting from a set of training pairs
$\\{(x_{i},y_{i})\\}_{i=1,\ldots,N}$, with $y_{i}=\Psi^{h}_{X}(x_{i})$, for an
accurate approximation $\Psi^{h}_{X}$ of the time-$h$ flow of $X$. We recall
that the universal approximation result based on the Presnov decomposition
allows approximating any vector field $X$ as

$X_{\theta}(z)=A^{T}\operatorname{diag}(\alpha)\Sigma(Az+b)+(B(z)-B(z)^{T})z=X_{G}(z)+X_{S}(z),$

as long as the weights are chosen correctly. In principle, calling $\Psi^{h}$
any numerical method applied to $X_{\theta}$ and said
$\hat{y}_{i}=\Psi^{h}(x_{i})$, one could train the weights of $X_{\theta}$ so
that they minimise

$\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\left\|y_{i}-\hat{y}_{i}\right\|^{2}.$

However, because of the properties of $X_{G}$ and $X_{S}$, we choose to
preserve them with $\Psi^{h}$ and apply a splitting method. In other words, we
apply to $X_{S}$ and $X_{G}$ two integrators $\Psi_{S}^{h}$ and
$\Psi_{G}^{h}$, then compose them to obtain
$\Psi^{h}=\Psi_{G}^{h}\circ\Psi_{S}^{h}$. We choose $\Psi_{S}^{h}$ to be an
explicit method that preserves the conserved quantity $\|x\|^{2}$, while
$\Psi_{G}^{h}$ to be a discrete gradient method (see e.g. [21, 43]) so that it
preserves the dissipative nature of $X_{G}$. As said before, this splitting
strategy is not necessary in principle. However, we propose it as an
alternative inspired by all the works on Hamiltonian neural networks (see e.g.
[19, 16, 21]) where geometric integrators are often utilised. Furthermore, it
would be interesting to understand if this or other splitting strategies give
better approximation results or theoretical guarantees, but this goes beyond
the scope of this work.

We choose $\Psi_{S}^{h}$ as a modified Euler-Heun method, following the
derivation presented in [13], so that it is explicit and it also preserves
$\|x\|^{2}$. For $\Psi_{G}^{h}$ we use instead the Gonzalez discrete gradient
method (see e.g. [43]).

Figure 7: In the two plots, we compare the true and predicted vector fields.

We remark that discrete gradient methods applied to $\dot{x}(t)=-\nabla
V(x(t))$ are of the form

$x^{n+1}=x^{n}-h\overline{\nabla}V(x^{n},x^{n+1})$

and they are hence implicit. However, since we have the trajectories
available, i.e. the $y_{i}$s are known, we do not have to solve a non-linear
system of equations. Indeed, the problem of approximating $X$ amounts to
minimise the following cost function

$\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\left\|y_{i}-\left(\Psi_{S}^{h}(x_{i})-h\overline{\nabla}V(\Psi_{S}^{h}(x_{i}),y_{i})\right)\right\|^{2}.$

In Fig. 7, we report the results obtained for the following two vector fields

$X_{1}(x)=\begin{bmatrix}(\sin{x_{1}})^{3}+x_{1}^{3}\\\
(\sin{x_{2}})^{3}+x_{2}^{3}\\\ (\sin{x_{3}})^{3}+x_{3}^{3}\\\
(\sin{x_{4}})^{3}+x_{4}^{3}\end{bmatrix}\quad X_{2}(x)=\begin{bmatrix}x_{2}\\\
-\sin{x_{1}}\end{bmatrix}.$

##  Appendix B Non-expansive networks with non-Euclidean metric

Figure 8: Three clean images from the CIFAR-10 dataset at the top and a random
perturbation of the red channel in one pixel at the bottom. For humans the
images of the two rows clearly associate to the same object, even with this
pixel perturbation.

In Section 3 we have introduced a way to generate networks that are non-
expansive for the Euclidean metric on the input space. We now propose a
strategy to generalise the reasoning. It is intuitive, and also evident from
Fig. 8, that $\ell^{p}$ norms of $\mathbb{R}^{n}$ are not the ones humans
utilise to compare pictures. There have been many attempts to design
similarity measures between images (see e.g. [63, 55]), but it is still not
evident what should be the preferred choice. To make the model Eq. 12
introduced in Section 3 more general, we hence show how it can be made
contractive for a more generic metric $d$ defined by a symmetric and positive
definite matrix $M\in\mathbb{R}^{n\times n}$. We introduce the notation
$\langle v,w\rangle_{M}=v^{T}Mw$ for any pair of vectors
$v,w\in\mathbb{R}^{n}$. Let again
$\Sigma(z)=[\sigma(z_{1}),\ldots,\sigma(z_{n})]$ where $\sigma$ is an
increasing scalar function. We focus on the autonomous dynamical system

(20) $\dot{z}(t)=-W^{T}\Sigma(MWz+b),\quad W\in\mathbb{R}^{n\times
n},\,b\in\mathbb{R}^{n},$

that is no more a gradient vector field, but it has similar properties to the
one studied above. We suppose $M$ is constant and the same for the other
involved weights. How this reasoning extends to time-switching systems, as
discussed throughout the paper, is quite natural and follows the procedure
seen for $M$ being the identity matrix, compare in particular Section 3. We
now verify the contractivity of the ODE 20 with respect to the metric defined
by $M$:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\frac{1}{2}\|z(t)-y(t)\|^{2}_{M}$
$\displaystyle=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left((z(t)-y(t))^{T}M(z(t)-y(t))\right)$
$\displaystyle=-\langle W^{T}\Sigma(MWz+b)-W^{T}\Sigma(MWy+b),z-y\rangle_{M}$
$\displaystyle=-\langle\Sigma(MWz+b)-\Sigma(MWy+b),MWz-MWy\rangle\leq 0.$

This result implies that all the trajectories of Eq. 20 will converge to a
reference trajectory if the convergence is measured using the metric defined
by $M$. Notice that the scalar product in the last line is the canonical one
of $\mathbb{R}^{n}$. Hence, if we have the $\gamma$ is strongly convex, we can
still combine this dynamics with expansive vector fields.