# On divergence-free neural ODE for classification

Zakaria Jarraya, Lucas Drumetz, Simon Benaïchouche, Douraied Ben Salem, François Rousseau

► To cite this version:

Zakaria Jarraya, Lucas Drumetz, Simon Benaïchouche, Douraied Ben Salem, François Rousseau. On divergence-free neural ODE for classification. ICPR 2024: 27th International Conference on Pattern Recognition, Dec 2024, Kolkata, India. pp.254-265, 10.1007/978-3-031-78186-5_17. hal-05120232

HAL Id: hal-05120232

https://imt-atlantique.hal.science/hal-05120232v1

Submitted on 19 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On divergence-free neural ODE for classification

Zakaria Jarraya<fnref n="1" />, Lucas Drumetz<fnref n="2" />, Simon Benaïchouche<fnref n="2" />, Douraied Ben Salem<fnref n="3" />, and François Rousseau<fnref n="1" />

**Abstract.** Inspired by the principles of fluid dynamics, divergence-free neural ODEs represent a promising class of neural networks that can improve classification performance and robustness. In this work, we propose methods for constructing divergence-free vector fields in high-dimensional latent spaces and thus extending their application to real-world scenarios. Our experimental results show that divergence-free neural ODEs can outperform their conventional counterparts in three classification tasks. In particular, divergence-free models show remarkable resilience to adversarial perturbations. Furthermore, divergence-free models are more effective at preserving the underlying structure of the data. This research provides a physics-based perspective on machine learning applications, in particular classification tasks, by interpreting neural networks as dynamical systems.

**Keywords:** Classification · neural ODE · divergence-free · adversarial robustness.

## 1 Introduction

Residual networks (ResNets)[13], a prominent type of neural network architecture, have attracted attention due to their remarkable performance improvements compared to traditional feedforward networks.

Their connection to ordinary differential equations (ODEs) [28, 26, 25] has opened a new avenue for integrating machine learning with dynamical systems research. This link has spurred the development of normalizing flows [19], a technique closely related to diffusion models, which are currently at the forefront of generative modelling.

Increasing the number of layers of a neural network generally leads to better performance. In the infinite depth limit, residual networks tend to a continuous ordinary differential equation representation:

$$ \frac{d\mathbf{z}}{dt} = f_{\theta}(\mathbf{z}, t) \qquad (1) $$

where $t$ is often refered as 'time', $f_{\theta}: \mathbb{R}^n \times [0, +\infty[ \to \mathbb{R}^n$ denotes a vector field parameterized by $\theta$ and $\mathbf{z}: [0, +\infty[ \to \mathbb{R}^n$ represents the $n$-dimensional state of the ODE. This continuous formulation is known as neural ODE [3].

<footnote n="1">¹ IMT Atlantique, LaTIM U1101 INSERM, Brest, France</footnote>

<footnote n="2">² IMT Atlantique, Lab-STICC, UMR CNRS 6285, Brest, France</footnote>

<footnote n="3">³ Univ Brest, CHU Brest, UMR 1101, F 29200 Brest, France</footnote>
Considering $z(0)$ as the input of neural network, the output $z(T)$ can be computed by using the following equation:

$$z(T) = \int_{0}^{T} f_{\theta}(z, t)dt$$

where $f_{\theta}$ is a neural network corresponding to the velocity field of the dynamical system.

There are two main ways of training neural ODEs: the first kind uses a continuous analog of backpropagation called the adjoint sensitivity method, leading to a second differential equation that is then solved numerically [3] with off-the-shelf methods. The second kind, on which we focus in this paper, first discretizes the ODE (using a differentiable implementation of e.g. Euler or Runge-Kutta methods) and optimizes a cost function on $z(T)$ via classical automatic differentiation [18]. A discretization of eq. (1) with Euler's forward integration scheme corresponds to a standard ResNet.

This theoretical framework allows deep networks to be viewed as dynamical systems, facilitating the smooth transformation of inputs into desired output configurations. From a probabilistic point of view, the underlying data distribution at time $T$ can be seen as the pushforward distribution of the input data distribution by the flow of the vector field $f_{\theta}$. If $f_{\theta}$ is regular enough, the data distribution is equivalently transported via a (weak) continuity equation [27] (Theorem 4.4) on the probability measure $\mu_t$ of $z(t)$:

$$\frac{\partial \mu_t}{\partial t} + \nabla \cdot (\mu_t f_{\theta}) = 0,$$

where $\nabla \cdot \mu$ is the divergence operator. Physically speaking, the data is transported in a continuous fashion without spontaneous appearance or disappearance of probability mass. The associated probability density function can be derived using the change of variables formula, incorporating the Jacobian matrix of the mapping, which is used among other applications in normalizing flows. In the context of classification tasks, this framework translates into considering the pushforward as an intermediate step before applying a linear classifier.

The study of the robustness of a classification algorithm, in particular with respect to adversarial attacks, becomes in this context a study of the stability of the dynamical system modelled by the neural ODE approach. This leads to an analysis of the properties of the underlying velocity field $f$. Empirical studies have demonstrated that neural ODEs exhibit superior robustness against adversarial attacks compared to convolutional neural networks (CNNs) [12]. This enhanced robustness can be attributed to the inherent property of continuous flows in neural ODEs, where the trajectories of data points do not intersect. Following this line of research, in order to improve the robustness and convergence of these networks, the neural ODE framework has been enriched by frameworks relying on Lyapunov functions [16, 24].

Recent research has explored embedding physical principles into deep neural networks to gain deeper insights into how they work. In particular, the divergence of the vector fields is a fundamental concept in fluid dynamics. Recent
works have focused on constructing neural ODE networks whose flow is derived
from divergence-free vector fields [23, 21]. By incorporating such constraints,
neural networks gain a physical basis that enhances their interpretability and
trustworthiness. This approach complements traditional machine learning tech-
niques by providing a deeper understanding of the network's behaviour and its
relationship to the underlying data.

This article focuses on neural ODE framework using divergence-free vector
fields for classification purpose. Our main focus is on investigating such networks
in conventional machine learning settings with high-dimensional data spaces.
Our main contributions include the development of scalable techniques for the
construction of divergence-free vectors fields, and conducting comprehensive ex-
perimental evaluations to assess the classification performance and robustness
of these approaches.

## 2 Methods

A common neural ODE-based network architecture for classification is depicted
in Figure 1. The input $\mathbf{x} \in X$ (i.e. a signal or an image) is first passed through
a feature extractor (also called encoder or embedding) to obtain a new repre-
sentation of the data: a latent variable $\mathbf{z}$. A neural ODE layer than acts as a
continuous nonlinear mapping in the feature space. Finally, a linear classifier $\mathbf{V}$
is applied to make a prediction $\hat{\mathbf{y}}$ of the true label $\mathbf{y}$ of the input $\mathbf{x}$.

In this work, we focus on the neural ODE part and investigate the impact
of divergence-free constraints for classification. To this end, we consider two
concurrent approaches: one that relies on penalizing the divergence term of vector
field in the training loss, and one that satisfies the divergence free condition by
design.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Neural ODE-based network architecture for classification. In this work, we focus on the neural ODE part and the divergence of the vector fields $f_\theta$.</figcaption></figure>
## 2.1 Divergence-free functions by penalization approach

A way to softly enforce a constraint on the solutions of a ML model is to penalize the squared norm of the equation defining the constraint [22]. The divergence-free property can be approximately satisfied in a neural network by penalizing the squared divergence at several randomly sampled points in the domain of $z$. For computational efficiency, these samples are often chosen as the training points. The neural network will have to find a balance between the classification task and satisfying the divergence-free property. The cost function of such a model writes:

$$ \arg \min_{\phi, \theta, V, \mathbf{b}} \sum_{i=1}^{N} \mathcal{L}(\sigma(\mathbf{V}\mathbf{z}_i(T) + \mathbf{b}), \mathbf{y}_i) + \lambda(\nabla \cdot f_{\theta}(\mathbf{z}_i(0)))^2 \quad (2) $$

where $N$ is the number of training samples, $\mathcal{L}$ is the classification loss, $\sigma$ is the softmax function, and $\lambda$ is a weight parameter that dictates how precisely the divergence-free property is enforced, at the expense of the classification loss. Here, the divergence is only penalized at initial time for all training points to limit the computational burden.

One advantage of this method of construction is its generality because it can applied to any existing network structure and adds the divergence-free property on top of it. The drawbacks include satisfying the divergence-free only approximately, which may be problematic in particular in physical applications [23], and the added regularizer in the training loss, which may harm the classification performance without a careful tuning of $\lambda$, in addition to increasing the computational burden since the optimizer has to differentiate through the divergence calculation, leading to a quadratic scaling in the dimension [2].

To alleviate this computational burden, a statistical estimator is used to approximate the divergence term. The divergence of vector field being trace of its Jacobian, the problem of computing the divergence of $f_{\theta}(\mathbf{z}(t))$ becomes that of computing efficiently the trace of its Jacobian. When Jacobian-vector products can be obtained efficiently, as is the case when using automatic differentiation, one may resort to statistical estimators that give an unbiased estimation of the trace in a way that scales linearly with the dimension. Without delving too much into details, most of those estimators are built around improving the convergence rate of the original Girard-Hutchinson estimator. It relies on the following property of the trace of a matrix $\mathbf{M} \in \mathbb{R}^{n \times n}$:

$$ \mathrm{Tr}(\mathbf{M}) = \mathbb{E}_{\varepsilon}[\varepsilon^T \mathbf{M} \varepsilon] \quad (3) $$

where $\varepsilon$ follows a centered and isotropic ($\mathbb{E}[\varepsilon\varepsilon^T] = \mathbf{I}_n$) distribution, typically the standard normal distribution. Then, for $m$ i.i.d. samples $\varepsilon_i$, $i = 1, ..., m$ of $\varepsilon$, an unbiased estimator of the trace is given by:

$$ \hat{\mathrm{Tr}}(\mathbf{M}) = \frac{1}{m} \sum_{i=1}^{m} \varepsilon_i^T \mathbf{M} \varepsilon_i \quad (4) $$
where *m* is the number of samples drawn.

## 2.2 Divergence-free functions by design

The other line of work for the construction of divergence-free functions for arbitrary high dimensions is by constructing functions that have this property baked-in. Following the notations of Richter-Powell et al. [23], let $\mathbf{A}: \mathbb{R}^n \to \mathbb{R}^{n \times n}$ be a smooth skew-symmetric matrix-valued function, then the vector field

$$v = \begin{bmatrix} \nabla \cdot \mathbf{A}_1 \\ \vdots \\ \nabla \cdot \mathbf{A}_n \end{bmatrix} \qquad (5)$$

is divergence-free. This is equivalent to proposition 2 in Morel et al. [21]. One key element of such an approach is that we have access to any smooth divergence-free vector fields (i.e. universality property). However, this universal method requires parameterizing $O(n^2)$ functions. As noted in [23, 21], this implies that these models are difficult to scale due to extensive use of automatic differentiation.

To handle high dimensional latent space in the Neural ODE part, we propose a scalable approach to construct divergence-free functions. Let $H$ be a smooth scalar function: $\mathbb{R}^n \to \mathbb{R}$ and $\mathbf{M} \in \mathbb{R}^{n \times n}$ a skew-symmetric matrix. Considering $\mathbf{A}(\mathbf{x}) = H(\mathbf{x})\mathbf{M}$, the vector field $v$ computed using Eq. (5) is divergence-free.

Our construction is a particular case of the universal method proposed in [23, 21]. As such, this construction is not universal but scales in $O(n)$. It is then possible to consider applications such as the classification of high-dimensional data, rather than being restricted to low-dimensional physical applications. Moreover, given an arbitrary skew-symmetric matrix $\mathbf{M}$, this approach leads to the straightforward computation of a divergence-free vector field given by $\mathbf{M}\nabla H$.

This research investigates methods employing penalization and the proposed approach for constructing divergence-free functions within the context of classification applications.

# 3 Experimental Setup

## 3.1 Data

In this section, we introduce the three datasets on which we test the different neural ODE constructions. We use first a 2D toy dataset consisting in two coiling spirals, each associated to a single class. This dataset has the advantage of being easy to visualise and will serve to build intuition about the different models, though of course one should be cautious about extending intuitions gained from a low dimensional case in higher dimensions. For this dataset, 1000 samples are generated equally distributed between the two classes. This dataset is split into a training and a testing dataset following a 70%-30% rule respectively.
Then, to get closer to relevant classification problems, in terms of dimension-
ality namely, we consider two standard image datasets: MNIST and CIFAR-10.
MNIST contains 70k samples in total of 28 × 28 gray-level images of digits be-
tween 0 and 9, 60k of which are used for training, the remaining 10k being used
for testing. CIFAR-10, on the other hand, contains low resolution 32 × 32 color
images of 10 different classes. There are 60k samples, 50k for training, and 10k
for testing.

Each dataset is normalized prior to the training using a standard scaler.

## 3.2 Architecture

The general architecture of the networks tested in the experiments comprises a feature extractor which is used to reduce the dimension of the input while at the same time capturing the most important features for the classification task.

Note that the feature extractor varies depending on the input data type and dimension. For the spiral data, there is no actual feature extractor in this case since the data are low dimensional. For the image datasets, we use opensource residual convolutional encoders as feature extractors: a ResNet18 [20] for MNIST and a WideResNet [20] for CIFAR10. The latent representation is vectorized into a 64-dimensional space.

Then, the Neural ODE layer transports the representation into its final state.
In all cases, the neural ODE is first discretized using a Dormand-Price (DOPRI)
5 solver [7]. *T* is set to 1 without loss of generality. After the integration of the
ODE, a trainable linear classifier outputs the prediction.

In the vanilla and penalization cases, the vector field $f_{\theta}$ is parameterized using 3 fully connected layers and Tanh activations. For the built-in divergence free model, the vector field $f_{\theta}$ is obtained analytically from the scalar function $H$ as described in Section 2.2. For fair comparison, $H$ is also parametrized using 3 fully connected layers, but the final layer outputs only a single neuron to obtain a scalar field.

## 3.3 Training

The training loss used is the Cross-Entropy while the optimization algorithm is Adadelta (in its standard Pytorch implementation) with a learning rate of 1 in all cases and a batch size of 50. In each experiment, the training procedure lasts for 100 epochs and each experiment is repeated 10 times to provide standard deviations for each evaluation metric in addition to average values. For the penalization method, the regularization parameter $\lambda$ is set to 0.01 for image datasets and 10 in the case of the coiling spirals. Those values were chosen empirically to maximize robustness (see next section for the associated metrics).

## 3.4 Evaluation metrics

In this section, we present the different metrics used to compare all competing methods. We use two types of metrics for each dataset, related to classification performance, and representation quality.
**Classification performance and Robustness** The first kind is related to classification performance, in a classical setting and in the case where the inputs are adversarial attacks generated using the FGSM algorithm [9] with a radius $\epsilon = 0.1$ for the coiling spirals and CIFAR10 datasets and $\epsilon = 0.3$ in case of MNIST for the $L_\infty$ norm. In both cases, we use the overall accuracy.

**Representation quality** The second kind are descriptive metrics that are proxies for robustness and describe the learned representation. We use the average intra class and inter class variances as descriptors of the average compactness of the representation of classes and separability between classes, respectively. Intuitively, the ratio between the intra and interclass variance should not become too large. However, we expect the scale of these two metrics to be smaller in the divergence free case, as the volume preserving property should lead to less expansion of the learned representations compared to unconstrained models.

We also use the overall energy of the learned vector field as a measure of the total displacement of the particles [17]:

$$E(f_{\theta}) = \int_{X} \int_{0}^{T} \|f_{\theta}(x(t), t)\|_{2}^{2} dt dx \quad (6)$$

where the integral is computed over the trajectories of the particles only for simplicity. Other discretizations in space and time of this metric are possible, for example by gridding the domain, but this gets intractable even for moderate dimensions, though it arguably provides a more robust measure regarding generalization outside the domain of the training data.

We also report the $L_2$ distortion metric between an initial latent vector $\mathbf{z}(0)$ and the final latent representation $\mathbf{z}(T)$, which indicates by how much does the neural ODE deviates from an isometry, using the usual Euclidean distances:

$$D_{f_{\theta}}(\mathbf{z}_1, \mathbf{z}_2) = \frac{\|\Phi_f^T(\mathbf{z}_1) - \Phi_f^T(\mathbf{z}_2)\|_2}{\|\mathbf{z}_1 - \mathbf{z}_2\|_2} \quad (7)$$

where $\Phi_f^T(\mathbf{z})$ is the flow of the vector field $f_{\theta}$ at $t$, starting from a given latent vector $\Phi_f^0(\mathbf{z}) = \mathbf{z}$. This quantity is averaged over 200 randomly sampled pairs of (encoded) training samples and thus gives an idea of the average local Lipschitz constant of each model over the support of the data, which is known to be related to robustness [30].

# 4 Resultats

All the metrics described in the previous section are reported in Table 1 for the three datasets (2D spirals, MNIST and CIFAR10) and the three considered methods (vanilla neural ODE, penalization-based neural ODE and divergence-free neural ODE).

First of all, all three methods, regardless of the dataset used, exhibit comparable classification accuracy. However, when subjected to FGSM-type adversarial
<table><caption>Table 1: Quantitative evaluation of the three methods (vanilla neural ODE, penalization-based approach and divergence-free model) on the three datasets (2D spirals, MNIST, CIFAR10) using six metrics (accuracy, accuracy under adversarial attacks, energy, intra and inter-class variance and distortion).</caption><thead><tr><th>Data</th><th>Method</th><th>Accuracy ↑</th><th>FGSM ↑</th><th>Energy ↓</th></tr></thead><tbody><tr><td rowspan="3">Spirals</td><td>Vanilla</td><td>95.25 ±5.3%</td><td>89.74 ±8.79%</td><td>102.44 ±54.81</td></tr><tr><td>Penalization</td><td>95.85 ±5.45%</td><td>90.3 ±8.3%</td><td>18.57 ±10.07</td></tr><tr><td>DivFree</td><td>97.76 ±3.21%</td><td>93.34 ±4.13%</td><td>32.80 ±19.43</td></tr><tr><td rowspan="3">MNIST</td><td>Vanilla</td><td>99.37 ±0.03%</td><td>76.11 ±2.28%</td><td>5.44 ±0.22</td></tr><tr><td>Penalization</td><td>99.16 ± 0.49%</td><td>81.89 ± 1.62%</td><td>3.21 ± 0.71</td></tr><tr><td>DivFree</td><td>98.62 ± 1.17%</td><td>79.44 ±1.71%</td><td>1.82 ± 0.42</td></tr><tr><td rowspan="3">CIFAR10</td><td>Vanilla</td><td>95.54 ± 0.03%</td><td>33.85 ± 0.11%</td><td>2.10 ± 0.01</td></tr><tr><td>Penalization</td><td>93.31 ± 0.02%</td><td>39.23 ± 0.90%</td><td>1.75 ± 0.36</td></tr><tr><td>DivFree</td><td>95.86 ± 0.00%</td><td>47.43 ± 0.66%</td><td>0.89 ± 0.25</td></tr><tr><th>Data</th><th>Method</th><th>Interclass variance ↓</th><th>Intraclass variance ↓</th><th>Distortion ↓</th></tr><tr><td rowspan="3">Spirals</td><td>Vanilla</td><td>7.66 ±3.06</td><td>19.24 ±10.81</td><td>6.27 ±3.54</td></tr><tr><td>Penalization</td><td>2.88 ±0.77</td><td>4.61 ±2.89</td><td>2.48 ±0.91</td></tr><tr><td>DivFree</td><td>4.31 ±1.26</td><td>8.93 ±4.29</td><td>3.53 ±1.86</td></tr><tr><td rowspan="3">MNIST</td><td>Vanilla</td><td>101.03 ±3.10</td><td>0.12 ±0.00</td><td>1.38 ±0.01</td></tr><tr><td>Penalization</td><td>86.61 ± 1.79</td><td>0.07 ± 0.05</td><td>1.20 ± 0.08</td></tr><tr><td>DivFree</td><td>59.93 ±4.15</td><td>0.05 ± 0.00</td><td>1.03 ± 0.01</td></tr><tr><td rowspan="3">CIFAR10</td><td>Vanilla</td><td>57.85 ± 0.00</td><td>9.29 ± 0.19</td><td>1.18 ± 0.01</td></tr><tr><td>Penalization</td><td>44.58 ± 0.39</td><td>6.57 ± 0.52</td><td>1.01 ± 0.07</td></tr><tr><td>DivFree</td><td>41.97 ± 0.82</td><td>7.38 ± 0.11</td><td>1.00 ± 0.00</td></tr></tbody></table>

attacks, a distinct pattern emerges. The introduction of the divergence-free constraint noticeably enhances the network's robustness compared to the vanilla neural ODE. This effect is particularly pronounced for the CIFAR10 dataset.

Examining the metrics indirectly related to robustness, it becomes apparent that incorporating the divergence-free constraint consistently results in a reduction in energy. Kinetic energy serves as a direct indicator of the distance traversed by the data points as they evolve from the initial state to the final output. The underlying assumption is that transformations characterized by shorter paths are less susceptible to perturbations, as minor alterations in the input data will not result in substantial changes in the output. From a physical perspective, this phenomenon aligns with the principle of least action, a fundamental concept in physics that has been applied to neural network training [17].

This trend extends to metrics associated with variance and distortion, further highlighting the robustness-enhancing effect of this constraint. Divergence-free neural networks exhibit a significantly lower ratio of intraclass to interclass variance compared to the vanilla method. This suggests that divergence-free networks are better at preserving the inherent structure of the data, leading to more
distinct class separation. The distortion metric tells us how much of the initial
information do we recover, the closer the value to one is, the more information
we recover. Here once again the divergence-free methods beat the vanilla one
when it comes to how much information they restore.

These quantitative results provide compelling evidence supporting the effec-
tiveness of the divergence-free constraint in enhancing the classification perfor-
mance of Neural ODE-based models. This analysis is further enriched by a visual
representation of the results for the 2D spirals dataset, presented in Figure 2.
To begin, we present an average of the decision boundaries obtained from ten
experiments, accompanied by the test set points overlaid. As evident from the
Table 1, all three methods effectively classify the data, achieving an accuracy
exceeding 95%. Notably, the divergence-free method inherently generates higher
uncertainty beyond the data support. This characteristic appears to be beneficial
for preventing overgeneralization beyond the data’s true distribution. In other
words, the divergence-free method seems to be more cautious about assigning
classifications to data points that lie far away from the training data. This vi-
sual interpretation is in line with the quantitative study of robustness described
previously.

Visualization of the vector fields shows a more structured pattern for the
vanilla method compared to the divergence-free method. In addition, the vector
fields differ in amplitude, reflecting significant variations in energy. The visual-
alization of the divergence field for the vanilla method reveals a pronounced
concentration of high-divergence values at the centers of the spirals, indicating
a "source" for the corresponding vector fields at the origin. As expected, the
penalization-based method occupies an intermediate position between the other
two methods, striking a balance between classification accuracy and adherence
to the divergence-free constraint.

# 5 Conclusion

In this study, we adopted a dynamical systems perspective to investigate neural
networks in the context of classification tasks. Our focus was on the Neural
ODE architecture, a continuous variant of residual networks. Building on recent
advances [23, 21], we proposed two strategies for incorporating divergence-free
constraint: a penalization method and a direct construction method.

The penalization method involves the introduction of a penalty term in the
loss function, which discourages the network from generating high divergence
vector fields. In contrast, the direct construction method obtains divergence-
free vector fields by computing the gradient of a smooth scalar function. This
approach allows the use of divergence-free functions in high-dimensional spaces.

We subjected these methods to rigorous experimental evaluation on three dif-
ferent datasets. Using six comprehensive evaluation metrics, we analysed their
performance and demonstrated their remarkable robustness against adversarial
attacks. Further investigations revealed that the introduction of the divergence-
free constraint leads to lower-energy solutions, consistent with fundamental phys-
Decision boundaries

<figure><img src="image_3.png" /></figure>

Vector field

<figure><img src="image_4.png" /></figure>

Divergence

<figure><img src="image_5.png" /><figcaption>Fig. 2: Qualitative evaluation on 2D spirals of the vanilla neural ODE (first column), penalization-based approach (second column) and divergence-free model (third column). First row shows the decision boundaries averaged over 10 runs, overlayed with testing points. The second row shows the vector fields displayed with their amplitude ($L_2$ norm) of one run. The third row displays the divergence (in log scale) of one run.</figcaption></figure>
ical principles such as the principle of least action. Furthermore, divergence-free
models are more effective at preserving the underlying structure of the data. This
research offers a physics-based perspective on machine learning applications,
particularly classification tasks, by interpreting neural networks as dynamical
systems.

**Acknowledgements** Please place your acknowledgments at the end of the paper, preceded by an unnumbered run-in heading (i.e. 3rd-level heading).

References

[1] Omri Azencot, N. Benjamin Erichson, Mirela Ben-Chen, and Michael W. Mahoney. "A Differential Geometry Perspective on Orthogonal Recurrent Models". ArXiv. 2021.

[2] Christopher Bishop and Hugh Bishop. *Deep Learning: Foundations and Concepts*. New York: Springer, 2023.

[3] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duve-
naud. "Neural ordinary differential equations". *Advances in neural infor-
mation processing systems*. 2018.

[4] Leena Chennuru Vankadara and Ulrike von Luxburg. "Measures of dis-
tortion for machine learning". *Advances in Neural Information Processing
Systems*. 2018.

[5] Moustapha Cisse, Piotr Bojanowski, Edouard Grave, Yann Dauphin, and Nicolas Usunier. "Parseval networks: Improving robustness to adversarial examples". *International conference on machine learning*. 2017.

[6] Francesco Croce, Maksym Andriushchenko, Vikash Sehwag, Edoardo Debenedetti, Nicolas Flammarion, Mung Chiang, Prateek Mittal, and Matthias Hein. "RobustBench: a standardized adversarial robustness benchmark". *Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks*. 2021.

[7] Gerhard Wanner Ernst Hairer and Syvert P. Nørsett. *Solving ordinary differential equations I*. Springer, 1993.

[8] Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam Ober-
man. "How to train your neural ODE: the world of Jacobian and kinetic
regularization". *International conference on machine learning*. 2020.

[9] Ian Goodfellow, Jonathon Shlens, and Christian Szegedy. "Explaining and Harnessing Adversarial Examples". *International Conference on Learning Representations*. 2015.

[10] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy. "Explaining and Harnessing Adversarial Examples". stat 1050 (2015), p. 20.

[11] Eldad Haber and Lars Ruthotto. "Stable architectures for deep neural networks". *Inverse problems*. Vol. 34. 1. 2017, p. 014004.

[12] Yan Hanshu, Du Jiawei, Tan Vincent, and Feng Jiashi. "On Robustness of Neural Ordinary Differential Equations". *International Conference on Learning Representations*. 2020.
[13] Kaiming He, X. Zhang, Shaoqing Ren, and Jian Sun. "Identity Mappings in Deep Residual Networks". *European Conference on Computer Vision*. 2016.

[14] Judy Hoffman, Daniel A Roberts, and Sho Yaida. "Robust Learning with Jacobian Regularization". *stat. Vol*. 1050. 2019, p. 7.

[15] Yifei Huang, Yaodong Yu, Hongyang Zhang, Yi Ma, and Yuan Yao. "Adversarial robustness of stabilized neural ode might be from obfuscated gradients". *Proceedings of Machine Learning Research*. 2021.

[16] Qiyu Kang, Yang Song, Qinxu Ding, and Wee Peng Tay. "Stable neural ode with lyapunov-stable equilibrium points for defending against adversarial attacks". *Advances in Neural Information Processing Systems*. 2021.

[17] Skander Karkar, Ibrahim Ayed, Emmanuel de Bézenac, and Patrick Galinari. "A principle of least action for the training of neural networks". *Machine Learning and Knowledge Discovery in Databases: European Conference*. 2020.

[18] P Kidger. "On neural differential equations". PhD thesis. University of Oxford, 2021.

[19] Ivan Kobyzev, Simon JD Prince, and Marcus A Brubaker. "Normalizing flows: An introduction and review of current methods". *IEEE transactions on pattern analysis and machine intelligence*. 2020.

[20] Kuang liu. *Github pytorch-cifar repository*. URL: https://github.com/kuangliu/pytorch-cifar.

[21] Guillaume Morel, Lucas Drumetz, Simon Benaichouche, Nicolas Courty, and François Rousseau. "Turning Normalizing Flows into Monge Maps with Geodesic Gaussian Preserving Flows". *Transactions on Machine Learning Research*. 2023.

[22] Maziar Raissi, Paris Perdikaris, and George 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*. 2019.

[23] Jack Richter-Powell, Yaron Lipman, and Ricky TQ Chen. "Neural conservation laws: A divergence-free perspective". *Advances in Neural Information Processing Systems*. 2022.

[24] Ivan Dario Jimenez Rodriguez, Aaron Ames, and Yisong Yue. "LyaNet: A Lyapunov framework for training neural ODEs". *International Conference on Machine Learning*. 2022.

[25] François Rousseau, Lucas Drumetz, and Ronan Fablet. "Residual networks as flows of diffeomorphisms". *Journal of Mathematical Imaging and Vision*. Vol. 62. 2020, pp. 365–375.

[26] Lars Ruthotto and Eldad Haber. "Deep neural networks motivated by partial differential equations". *Journal of Mathematical Imaging and Vision*. Vol. 378. 2020, pp. 686–707.

[27] Filippo Santambrogio. "Optimal transport for applied mathematicians". *Birkhäuser, NY*. Vol. 55. 58-63. 2015, p. 94.
[28] E Weinan. "A proposal on machine learning via dynamical systems". *Com-munications in Mathematics and Statistics*. Vol. 1. 5. 2017, pp. 1-11.

[29] Yandong Wen, Kaipeng Zhang, Zhifeng Li, and Yu Qiao. "A Discrimi-
native Feature Learning Approach for Deep Face Recognition". *European
Conference on Computer Vision*. 2016.

[30] Yao-Yuan Yang, Cyrus Rashtchian, Hongyang Zhang, Russ R Salakhutdi-
nov, and Kamalika Chaudhuri. "A closer look at accuracy vs. robustness".
*Advances in neural information processing systems*. 2020.

[31] Muhammad Zakwan, Liang Xu, and Giancarlo Ferrari-Trecate. "Robust classification using contractive Hamiltonian neural ODEs". *IEEE Control Systems Letters*. 2022.