[1]C. Coelho

1]Centre of Mathematics (CMAT), University of Minho, Braga, 4710-57, Portugal

2]Department of Mechanical Engineering (Section of Mathematics) - FEUP,
University of Porto, Porto, 4200-465, Portugal

# Neural Fractional Differential Equations

<EMAIL_ADDRESS>M. Fernanda P. Costa<EMAIL_ADDRESS>L.L. Ferrás
<EMAIL_ADDRESS>[  [

###### Abstract

Fractional Differential Equations (FDEs) are essential tools for modelling
complex systems in science and engineering. They extend the traditional
concepts of differentiation and integration to non-integer orders, enabling a
more precise representation of processes characterised by non-local and
memory-dependent behaviours. This property is useful in systems where
variables do not respond to changes instantaneously, but instead exhibit a
strong memory of past interactions. Having this in mind, and drawing
inspiration from Neural Ordinary Differential Equations (Neural ODEs), we
propose the Neural FDE, a novel deep neural network architecture that adjusts
a FDE to the dynamics of data. This work provides a comprehensive overview of
the numerical method employed in Neural FDEs and the Neural FDE architecture.
The numerical outcomes suggest that, despite being more computationally
demanding, the Neural FDE may outperform the Neural ODE in modelling systems
with memory or dependencies on past states, and it can effectively be applied
to learn more intricate dynamical systems.

###### keywords:

Neural Fractional Differential Equations, Neural Ordinary Differential
Equations, Neural Networks, Time-Series, Numerical Methods.

##  1 Introduction

Real systems in science and engineering, exhibit complex behaviours, often
characterised by intricate dynamics, non-linear interactions, and emergent
phenomena. These complexities arise in various contexts, such as the complex
interactions of molecules within a cell, the chaotic movement of turbulent
flows, and the challenging task of predicting financial markets.

To effectively predict and understand these complex systems, mathematical
models are employed, allowing to gain insights into the system behaviour
without the need for time-consuming or expensive experiments.

Due to the inherent presence of continuous dynamics in these systems,
Differential Equations (DEs) are commonly employed as mathematical models,
accounting for the continuous evolution of the system’s behaviour and offering
the advantage of enabling predictions throughout the entire time domain and
not only at specific points.

With the emergence of Neural Networks (NNs) and their impressive performance
in fitting mathematical models to data, numerous studies have focused on
modelling real-world systems. However, conventional NNs are designed to model
discrete functions and may not be able to accurately capture the continuous
dynamics observed in several systems. To overcome this limitation, Chen et al.
[1] introduced the Neural Ordinary Differential Equations (Neural ODEs), a NN
architecture that adjusts an Ordinary Differential Equation (ODE) to the
dynamics of a system.

ODEs are favored for their simplicity and effectiveness in describing the
instantaneous rates of change within a system. However, they may not
adequately model systems with more complex memory or strong dependence on past
states beyond the immediate previous states. To address this issue, Fractional
Differential Equations (FDEs) have been introduced [2]. FDEs expand the notion
of derivatives to fractional orders, facilitating the modelling of systems
with memory or non-local interactions. They use all past states to predict the
current and future behaviours of the system and capture fractional-order rates
of variable changes, offering a more thorough representation of the system
dynamics. Consequently, FDEs enable an accurate depiction of systems with
memory or long-range dependencies [3].

Inspired by Neural ODEs and by the inherent memory of FDEs, in a preliminary
work, we briefly proposed a novel deep learning architecture that models a FDE
to the hidden dynamics of the data [4]. To the best of our knowledge, this is
the first time that a NN architecture is proposed to fully fit a FDE to the
dynamics of data, including the order of the fractional derivative. In this
paper, we build upon our prior work [4] by introducing several new and
significant contributions:

  * •

Background: we incorporate a concise background section, providing essential
context to comprehending and understanding the proposed Neural FDE, namely on
Neural ODEs and FDEs;

  * •

Justification: we add a detailed explanation of the mechanisms underlying
Neural FDEs and their theoretical advantages, thereby establishing a robust
justification for their adoption;

  * •

Computational cost: we conduct a thorough analysis of the time and memory
complexity of Neural FDEs, providing insights into the computational
efficiency of our proposed method;

  * •

Numerical experiments: we augment the numerical experiments in [4] by
introducing three additional systems from various fields: relaxation
oscillation process, resistor-capacitor circuit, and a substance decay
process. For each system, we include datasets obtained from the solution of
FDEs considering different orders of the derivative ($\alpha$) totalling five
$\alpha$ values and one ODE ($\alpha=1$). Additionally, we investigate the
model’s robustness and generalisation by testing with higher-frequency
sampling and for larger time interval sets;

  * •

Learnt derivative order analysis: we explore the learnt values of $\alpha$
(the order of the fractional derivate) to discern potential patterns in the
learning process;

  * •

Convergence analysis: we conduct a series of experiments and analysis to
assess the convergence properties and stability of Neural FDEs, providing
valuable insights into the reliability of our proposed method;

  * •

Developed datasets: we meticulously present and describe the datasets used in
our numerical experiments for the four systems, elucidating the processes
behind their generation.

This paper is organised as follows. Section 2 presents a brief review of
essential concepts such as FDEs and Neural ODEs. Section 3 presents the Neural
FDE architecture along with its mathematical formulation and algorithm. In
Section 4 we evaluate the performance of the newly proposed Neural FDE. We use
datasets describing two real-world systems to compare the performance of
Neural FDE with a Neural ODE baseline. The paper ends with the summary of the
findings and future directions in Section 5. Additionally, in Appendix A we
provide an explanation on the impact of refining meshes and in Appendix B we
present numerical results for two additional real-world system datasets.

##  2 Background

Consider a time series, denoted as
$\boldsymbol{X}=(\boldsymbol{x}_{0},\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N-1})$,
where each element $\boldsymbol{x}_{n}\in\mathbb{R}^{d}$ represents the state
of the system at instant $t_{n}$ ($n=0,\dots,N-1$). Similarly, let
$\boldsymbol{Y}=(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\dots,\boldsymbol{y}_{N^{*}})$
denote the associated ground-truth output time series. In this context,
$\boldsymbol{y}_{n}\in\mathbb{R}^{d^{*}}$ is the output at time step $t_{n}$.
Additionally, the predicted output time series is expressed as
$\boldsymbol{\hat{Y}}=(\boldsymbol{\hat{y}}_{1},\boldsymbol{\hat{y}}_{2},\dots,\boldsymbol{\hat{y}}_{N^{*}})$,
where each prediction $\boldsymbol{\hat{y}}_{n}\in\mathbb{R}^{d^{*}}$
corresponds to time step $t_{n}$. Note that $N^{*}$ may differ from $N$.

#### Neural Ordinary Differential Equations

Residual Networks can be seen as a sequence of transformations to a hidden
state, [5, 6, 7, 8, 9, 10]. This sequence of transformations was originally
proposed to be of the form (for a network with $M$ layers):

$\boldsymbol{h}_{j+1}=\boldsymbol{h}_{j}+F(\boldsymbol{h}_{j},\boldsymbol{\theta}_{j}),~{}~{}~{}~{}j=0,...,M-1$
(1)

with $F$ being a function that depends on the state $\boldsymbol{h}_{j}$ and
on $\boldsymbol{\theta}_{j}$ (where $\boldsymbol{\theta}_{j}$ represents all
the parameters of the network determined by the learning process [6] \- the
weight matrix and bias).

The _new_ state, $\boldsymbol{h}_{j+1}$, is obtained by adding
$F(\boldsymbol{h}_{j},\boldsymbol{\theta}_{j})$ to the _old_ state,
$\boldsymbol{h}_{j}$, and, for the residual forward problem presented in (1)
to be stable, it is recommended to control the variation of
$F(\boldsymbol{h}_{j},\boldsymbol{\theta}_{j})$ along the iterations.
Therefore, in [7, 6] (see also [9]), the authors proposed an improved version
of (1) (see (2)) by introducing a constant $\delta$ — a positive constant
enabling control over the variation of $F$ (the smaller the $\delta$ the more
control we have over the variation of $F$). This enhancement results in the
following refined residual network:

$\boldsymbol{h}_{j+1}=\boldsymbol{h}_{j}+\delta
F(\boldsymbol{h}_{j},\boldsymbol{\theta}_{j}),~{}~{}~{}~{}j=0,...,M-1.$ (2)

This equation can be rewritten as,

$\dfrac{\boldsymbol{h}_{j+1}-\boldsymbol{h}_{j}}{\delta}=F(\boldsymbol{h}_{j},\boldsymbol{\theta}_{j}),~{}~{}~{}~{}j=0,...,M-1,$
(3)

which is the explicit Euler discretization (with time-step $\delta$) of the
nonlinear ODE,

$\dfrac{d\boldsymbol{h}(t)}{dt}=F(\boldsymbol{h}(t),\boldsymbol{\theta}(t))$
(4)

over the time interval $[0,T]$ ($t$ is the time). The continuous nature of
this problem is encapsulated in (4), which extends beyond the scope of (2). It
introduces the capability to derive the solution or state for any given moment
within the specified time interval. Several researchers have undertaken
analytical studies on this continuous residual equation (4), which has also
been extended to address more intricate differential equations (dynamical
systems), as exemplified in [10].

Remark 1: To go from (3) to (4), one should have the property that as $\delta$
becomes smaller, the difference between $\boldsymbol{h}_{j+1}$ and
$\boldsymbol{h}_{j}$ should also become smaller. This is in fact true, because
in the limiting case of $\delta=0$, we have from (2) that
$\boldsymbol{h}_{j+1}=\boldsymbol{h}_{j}$.

Remark 2: While it is indeed true that the parameter $\delta$ provides a means
to regulate the stability of the forward problem (2), the issue of stability
is, in fact, a more intricate matter. It is not solely determined by $\delta$
but also influenced by the variations in the weight matrix – refer to [7, 6]
for a comprehensive explanation.

In 2018, the Neural ODE framework was introduced [1]. Building upon (4), the
authors [1] introduced a neural network (NN) architecture specifically
tailored for modeling time-series data characterized by continuous-time
dynamics. Unlike conventional NNs, which are confined to representing discrete
dynamics, the use of Neural ODEs involves adjusting or fine-tuning a function
($\boldsymbol{f_{\theta}}$) to ensure that the solution of the Initial Value
Problem (IVP),

$\dfrac{d\boldsymbol{h}(t)}{dt}=\boldsymbol{f_{\theta}}(\boldsymbol{h}(t),t)\quad\text{with}\quad\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0},$
(5)

accurately represents a dynamic system. Here, $\boldsymbol{f_{\theta}}$
represents the right-hand side of the ODE, given by a NN with parameters
$\boldsymbol{\theta}$, where $t$ denotes time, and
$\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$ is the initial condition of the
IVP, essentially capturing the state at the first instance of the provided
data. This fine-tuning of $\boldsymbol{f_{\theta}}$ allows the solution of (5)
to effectively fit the given data, capturing the behaviour of the dynamic
system. The innovation in the Neural ODE, as compared to (4), lies in the
calculation of $F(\boldsymbol{h}(t),\boldsymbol{\theta}(t))$ using a NN,
$\boldsymbol{f_{\theta}}(\boldsymbol{h}(t),t)$, instead of using an analytical
procedure.

In addition to the NN, the Neural ODEs are made up of a numerical ODE solver
(referred to as $ODESolve$) for the numerical solution of the differential
equation, and each state is then given by,

$\\{\boldsymbol{h}(t)\\}_{t=t_{0},t_{1},\dots,t_{f}}=ODESolve(\boldsymbol{f_{\theta}},\boldsymbol{h}_{0},(t_{0},t_{1},\dots,t_{f})),$
(6)

where $\boldsymbol{h}(t)$ is the vector solution at instant $t$ being $t$ in
the time interval $[t_{0},t_{f}]$.  With analogy to the discrete case,
$\boldsymbol{h}(t_{0})$ is the input layer and $\boldsymbol{h}(t_{f})$ is the
output layer. Since $t\in[t_{0},t_{f}]$, for all values
$t^{*}\in]t_{0},t_{f}[$ we obtain $\boldsymbol{h}(t^{*})$ that can be seen as
a continuous version of hidden layers.

Remark 3: For simplicity we consider that the output of the numerical method
is the output of the Neural ODE, $\boldsymbol{\hat{Y}}(\boldsymbol{\theta})$.
However this may not be always the case and a second NN might be used to
transform the solution of the numerical method into the prediction
$\boldsymbol{\hat{Y}}(\boldsymbol{\theta})$.

The solver makes use of the $\boldsymbol{f_{\theta}}$ function, which takes
the current state of the system $\boldsymbol{h}_{t}$ as input, delivers the
derivative of the state as output ($d\boldsymbol{h}(t)/dt$) and then
integrates the derivative over time to determine the next state
$\boldsymbol{h}_{t+1}$, that is:

$\boldsymbol{h}_{t+1}=\boldsymbol{h}_{t}+\int_{t}^{t+1}\boldsymbol{f_{\theta}}(\boldsymbol{h}(t),t)dt.$
(7)

The Neural ODE [1] allows the use of both adaptive and fixed-step solvers.
When employing a fixed-step solver without explicitly specifying the step size
$\Delta t$ (default case), the time interval discretization occurs solely
based on the time-steps forced by the time-series training data. However, when
the step size is explicitly specified, the discretization takes place for each
sub-interval $(t_{n},t_{n+1})$ with the specified step size yielding
$(t_{n+1}-t_{n})/\Delta t$ time steps for each sub-interval. On the other
hand, when using an adaptive-step step solver, the discretization also occurs
for each sub-interval $(t_{n},t_{n+1})$, but the step size is not
predetermined, instead it is dynamically adjusted based on the local solution
behaviour (local gradient) [11]. The refinement of the mesh inside each sub-
interval allows for the solver to compute the solution with higher accuracy at
an increased computational cost, for a detailed explanation on the
implications of the use of different step sizes see Appendix A.

#### Fractional Differential Equations

FDEs are a class of differential equations that extend the concept of integer-
order derivatives to non-integer (or fractional) orders. There exist different
definitions for fractional derivatives, but for the purpose of this study, our
exclusive focus will be on the Caputo fractional derivative of order $\alpha$.
It is denoted as $\leftindex^{C}_{0}{D}^{\alpha}_{t}g(t)$ [12] (also refer to
[2]), and can be expressed for a general function $g(t)$ as follows:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}g(t)=\dfrac{1}{\Gamma(\lceil\alpha\rceil-\alpha)}\int_{0}^{t}(t-s)^{\lceil\alpha\rceil-\alpha-1}g^{\lceil\alpha\rceil}(s)ds,$
(8)

where $\alpha$ represents the order of the fractional derivative; the symbol
$\lceil\alpha\rceil$ denotes the ceiling function applied to $\alpha$
(producing the smallest integer greater than or equal to $\alpha$).
Additionally, we make use of the Gamma function, denoted by $\Gamma$ [13, 12,
2]. We opt for this specific fractional derivative as it is more suited for
modelling physical processes, since the initial conditions are based on
integer-order derivatives.

In this study, we extend (4) by incorporating the Caputo fractional derivative
(insights into the application of fractional calculus in Neural systems are
discussed in [14], providing a rationale for this modification). This
enhancement allows our dynamic system to encompass fractional rates and
memory, and the resulting outcome is:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}\boldsymbol{h}(t)=\boldsymbol{f}(\boldsymbol{h}(t),t),\,\,\,~{}~{}~{}~{}t>0,\,\,\,\boldsymbol{h}^{(k)}(t_{0})=\boldsymbol{h}^{k}_{0},\,\,\
k=0,\dots,\lceil\alpha\rceil-1,$ (9)

with $(k)$ the derivative of integer-order $k$ and
$\boldsymbol{h}^{(k)}(t_{0})=\boldsymbol{h}^{k}_{0}$ the $k^{th}$ initial
condition of the IVP ($\boldsymbol{h}^{k}_{0}$ is the real value or the vector
of real values associated with the $k^{th}$ inicial condition). Note that
$\boldsymbol{f}(\boldsymbol{h}(t),t)$ represents a function and not a NN. If
this function is continuous and fulfils a Lipschitz condition with respect to
its first argument with Lipschitz constant $L$ on a suitable set $G$ then this
equation has a unique solution on some interval $[0,T]$ [15].

Unlike integer-order derivatives found in ODEs, where the derivative’s value
at a point is solely dependent on the function’s behaviour within an
infinitesimal neighbourhood around that point (given an initial value, the
solution is uniquely determined for any point of the domain [16]), fractional
derivatives, as described by Equation (9), exhibit non-local characteristics.
To compute the state $\boldsymbol{h}_{t}$, it is necessary to take into
account the complete history from the initial state up to the current time, t,
[13]:

$\boldsymbol{h}_{t}=\sum^{\lceil\alpha\rceil-1}_{k=0}\dfrac{t^{k}}{k!}\boldsymbol{h}_{0}^{k}+\dfrac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}\boldsymbol{f}(\boldsymbol{h}(s),s)ds.$
(10)

While the non-locality of FDEs adds complexity and computational expense to
their solution, this characteristic also means that FDEs inherently encompass
memory effects [17]. This property proves beneficial in systems where
variables don’t respond instantaneously to changes but rather exhibit a memory
of past interactions [3].

Remark 4: It’s worth noting that FDEs have been previously proposed for
analytical modelling of NNs, addressing aspects such as stability,
bifurcations, and chaos (refer to [18] and the references therein). It’s
important to clarify that the objective of this work is not the analytical
examination of (9), but rather the development of a Neural FDE.

##  3 Neural Fractional Differential Equations

Having in mind the concepts of Neural ODE [1] and fractional calculus in
Neural systems [14], in [4] we introduced the Neural FDE, a NN architecture
parameterised by $\boldsymbol{\theta}$ (weights and biases) that adjusts the
function $\boldsymbol{f_{\theta}}$ of a FDE of order $\alpha$ (see (11)). The
goal is to ensure that the curve of solutions of this FDE fits the provided
data accurately by changing the weights and biases in
$\boldsymbol{f_{\theta}}$.

Since this model has an extra parameter ($\alpha$) when compared to the Neural
ODE, the order of the FDE, $\alpha$, is also learnt by another NN,
$\alpha_{\boldsymbol{\phi}}$, with parameters $\boldsymbol{\phi}$. This makes
the process of adjusting (11) to the training data completely independent of
any user intervention and gives the Neural FDE the freedom to find the best
possible fit:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}\boldsymbol{h}(t)=\boldsymbol{f_{\theta}}(\boldsymbol{h}(t),t)\text{
with
}\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0},\,\,\,\alpha=\alpha_{\boldsymbol{\phi}}.$
(11)

Remark 5: In this work we consider $\alpha\in(0,1)$ since this is the case in
many applications and simplifies the formulation [17]. This sets $k=0$ in (9),
thus requiring only one initial condition
$\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$ for ensuring the uniqueness of
solution. In case $\alpha>1$, $\lceil\alpha\rceil$ initial conditions are
necessary. Furthermore, a FDE solver is used to solve the IVP with initial
condition $\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$ in time interval
$[t_{0},t_{f}]$:

$\boldsymbol{h}(t)=PCSolve(\alpha,\boldsymbol{f_{\theta}},\boldsymbol{h}_{0},[t_{0},t_{f}]).$
(12)

The name $PCSolve$ comes from the fact that we use a Predictor-Corrector
approach to solve the FDE numerically [17].

The predictions of the Neural ODE are then given by (12), with
$\boldsymbol{h}_{t}=\boldsymbol{\hat{y}}_{t}$. Note, however, that in certain
scenarios, this may not be the case and an additional NN might be employed to
transform the solver output into the desired form.

In summary, the Neural FDE is composed of three main components, as shown in
Figure 1:

  * •

A NN that adjusts the FDE dynamics, $\boldsymbol{f_{\theta}}$: the NN is a
function approximator that learns a function that maps the current state of
the system $\boldsymbol{h}_{t}$ to its derivative,
$\boldsymbol{f}(\boldsymbol{h}(t),t)=\leftindex^{C}_{0}{D}^{\alpha}_{t}\boldsymbol{h}(t)$,
see (9). Any arbitrary NN architecture can be used, from the simpler multi-
layer perceptron to the more complex residual network (ResNet).
$\boldsymbol{f_{\theta}}$ receives the data of the system being modelled as
input and outputs the _results_ to be used by the FDE solver (Equation (12))
to update the state of the system;

  * •

A NN that adjusts the order of the derivative, $\alpha_{\boldsymbol{\phi}}$:
the goal of this NN is to adjust the value of $\alpha\in(0,1)$ in Equation
(11). $\alpha_{\boldsymbol{\phi}}$ receives _a previous/ initial_ $\alpha$
value, so the value of $\alpha$ must be initialised. Similar to
$\boldsymbol{f_{\theta}}$, the choice of architecture for this NN is flexible.
However, given that the objective is to adjust a single value, a
straightforward multi-layer perceptron suffices. It is important to note that,
since $\alpha\in(0,1)$, it is imperative to ensure that the value generated by
$\alpha_{\boldsymbol{\phi}}$ remains within this interval. An approach to
achieve this is to employ a bounded activation function in the output layer.
In the experiments detailed in this paper, a sigmoid activation function was
used;

  * •

A FDE numerical solver: to solve FDEs specialised mathematical techniques and
numerical methods are necessary. Several numerical methods have been
introduced in the literature to solve FDEs, however, during the development of
the Neural FDE, we made the strategic decision to implement the Predictor-
Corrector (PC) algorithm for FDEs, as outlined in [17]. We chose this
algorithm for being general, making it suitable for both linear and nonlinear
problems, as well as for single and multi-term equations. This choice provides
us with greater flexibility in addressing a wide range of scenarios when using
Neural FDEs. Note that any FDE solver can be used with Neural FDEs.

Figure 1: Schematic of the architecture of a Neural FDE. The $PCSolve$ is used
initially during model training to determine the optimal weights and biases,
so that (11) fits the given data. Subsequently, the $PCSolve$ is employed
again during model prediction (in this case, the NN $\boldsymbol{f_{\theta}}$
and $\alpha_{\phi}$ are already known, and their sole purpose is to supply the
essential values required for the numerical solution of (11) to be feasible).
Please refer to Algorithms (1) and (2) for a succinct explanation.

The Neural FDE comprises three components that collectively contribute to
computing the loss function $\mathcal{L}_{\boldsymbol{\theta,\phi}}$ during
the optimisation process (training). This process aims to adjust the
parameters $\boldsymbol{\theta},\boldsymbol{\phi}$ by minimising the loss
function (13) (a measure of the disparity between predicted and ground truth
values, such as, for example, the Mean Squared Error (MSE)); see Algorithm 1.

$\begin{split}\mathcal{L}_{\boldsymbol{\theta},\boldsymbol{\phi}}(\boldsymbol{h}(t_{f}))&=\mathcal{L}_{\boldsymbol{\theta},\boldsymbol{\phi}}\left(\sum^{\lceil\alpha\rceil-1}_{k=0}\dfrac{t_{f}^{k}}{k!}\boldsymbol{h}_{0}^{k}+\dfrac{1}{\Gamma(\alpha)}\int_{0}^{t_{f}}(t_{f}-s)^{\alpha-1}f(\boldsymbol{h}(s),s)ds\right)\\\
&=\mathcal{L}_{\boldsymbol{\theta},\boldsymbol{\phi}}(PCSolve(\alpha,\boldsymbol{f}_{\theta},\boldsymbol{h}_{0},[t_{0},t_{f}])).\end{split}$
(13)

Given that both parameters $\boldsymbol{\theta}$ and $\boldsymbol{\phi}$ are
essential for computing the predictions $\boldsymbol{h}(t)$, we have deduced,
as a first approach, that optimising them jointly would be a good strategy.
Thus, in this work, a single Adam optimiser is used to optimise the full set
of parameters of the Neural FDE. Nevertheless, we recognise that there exist
various possible approaches to train these parameters.

When the training of the Neural FDE has been completed, the result is a FDE,
adjusted by two NNs (Figure 1). To make predictions, a numerical FDE solver is
used to solve the IVP (11) and output the prediction’s vector, Algorithm 2.

Algorithm 1 Neural FDE training process.

Input: start time $t_{0}$, end time $t_{f}$, initial condition
$\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$, maximum number of iterations
_MAXITER_ ;

Choose $Optimiser$;

$\boldsymbol{f_{\theta}}=DynamicsNN()$;

$\alpha_{\boldsymbol{\phi}}=AlphaNN()$;

Initialise $\boldsymbol{\theta},\boldsymbol{\phi}$;

for $k=1:$_MAXITER_ do

$\alpha\leftarrow\alpha_{\boldsymbol{\phi}}$;

${\\{\boldsymbol{h}_{t}\\}}_{t=t_{1},t_{2},\dots,t_{f}}\leftarrow
PCSolve(\alpha,\boldsymbol{f_{\theta}},\boldsymbol{h}_{0},[t_{0},t_{f}])$;

Evaluate loss $\mathcal{L}$;

$\nabla\mathcal{L}\leftarrow\text{Compute gradients of }\mathcal{L}$;

$\boldsymbol{\theta},\boldsymbol{\phi}\leftarrow
Optimiser.Step(\nabla\mathcal{L})$;

end for

Return: $\boldsymbol{\theta},\alpha$;

Algorithm 2 Neural FDE prediction process.

Input: start time $t_{0}$, end time $t_{f}$, initial condition
$\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$, parameters $\boldsymbol{\theta}$,
order $\alpha$;

Load $\boldsymbol{f_{\theta}}$;

${\\{\boldsymbol{h}_{t}\\}}_{t=t_{1},\dots,t_{f}}\leftarrow
PCSolve(\alpha,\boldsymbol{f_{\theta}},\boldsymbol{h}_{0},[t_{0},t_{f}])$;

Return: ${\\{\boldsymbol{h}_{t}\\}}_{t=t_{1},t_{2},\dots,t_{f}}$;

The solver implemented within the Neural FDE allows the use of fixed-step
solvers. As in the Neural ODE, when employing a fixed-step solver without
explicitly specifying the step size $\Delta t$ (default case), the time
interval discretization occurs solely based on the time-steps forced by the
time-series training data. However, when the step size is explicitly
specified, the discretization takes place for each sub-interval
$(t_{n},t_{n+1})$ with the specified step size yielding
$(t_{n+1}-t_{n})/\Delta t$ time steps for each sub-interval. If the mesh
points do not match the training data, an interpolation is used. This is
illustrated in Figure 2.

Remark 5: It’s important to note that, in most cases, solving FDEs of this
nature results in a singularity at the origin. To tackle this issue, certain
solvers incorporate a mesh refinement specifically in the vicinity of the
origin. It’s worth mentioning that, in the future, we plan to implement this
graded mesh approach as a standard practice.

Figure 2: Mesh refinement process employed in the solver, along with an
illustration of linear interpolation. Case 1: the time mesh nodes (dash) used
in the numerical method aligns with the input data series (symbols); Case 2:
the time mesh used in the numerical method is now more refined, but for
specific mesh nodes it aligns with the input data series; Case 3: the time
mesh nodes used in the numerical method can be more refined or not, but they
do not align with the input data series. The numerical solution is obtained
for all mesh nodes considered, but for training purposes only the solutions at
$t_{1}$, $t_{2}$, … are considered, being these obtained by linear
interpolation (this methodology can also be used for irregularly sampled
data).

#### Non-locality and Memory Effects

Neural FDE offer advantages over Neural ODE in modelling time-series data,
primarily because of their inherent ability to capture long-term memory. When
solving an FDE, the current time step $t$ is influenced by the entire past
history of the system under consideration (see Figure 3). This challenges the
traditional concept of causality found in ODEs, and inherent to Neural ODEs.
While ODEs typically relate a future state to its immediate preceding state,
FDEs acknowledge that the current state can be influenced by all past states,
emphasising a non-local and memory-dependent perspective on causality. This
memory effects introduce a broader range of interactions and dependencies,
leading to complex behaviours that cannot be adequately captured by ODE
models.

Figure 3: Memory in FDEs (left) and ODEs (right). Schematic of the computation
of the derivative (fractional or classical) at instant $t$. The case of short
memory (center) is presented for illustrative purposes, aiming to facilitate a
clear understanding of this phenomenon.

#### Computational Cost

The computational cost of the Neural FDE surpasses that of the Neural ODE,
primarily owing to its non-local property. A significant factor contributing
to the computational expenses in Neural FDE is the memory cost associated with
the FDE solver. Unlike ODE solvers in Neural ODE, which rely on information
from the immediate past time step, FDE solvers must retain values of the
function at all preceding time steps to accurately compute subsequent ones.
This storage of historical data substantially increases computational overhead
[17].

The memory cost also plays a crucial role during backpropagation through the
FDE solver. When training Neural FDE models, the necessity to backpropagate
through all preceding time steps further amplifies the computational burden,
resulting in longer training times compared to Neural ODE models.

The characteristic of Neural ODEs makes them a more efficient choice for
modelling systems with simple dynamics and where memory is not a critical
factor. However, while the higher computational cost may be a disadvantage, it
is important to carefully consider the advantages of using Neural FDEs.

##  4 Numerical Experiments

To assess the performance of the architectural framework proposed in this
study, a comprehensive series of experiments was carried out on a variety of
datasets. First, four synthetic toy datasets were generated by numerically
solving their corresponding differential equations. These datasets were
specifically designed to model four distinct and well-established systems
governed by differential equations: Relaxation Oscillation Process (RO),
Population Growth (PG), RC Circuit Capacitor (RC), and Substance Decay (SD).
For each dataset, three experiments were conducted to comprehensively evaluate
the performance of the models:

  * •

Reconstruction: assesses the models’ capability to learn the dynamics of the
training set. The evaluation is carried out using the same set for training
and testing, consisting of $200$ points within the time interval $t=[0,200]$.
This scenario provides insights into how well the models can reproduce the
observed dynamics within the specified temporal range;

  * •

Extrapolation: evaluates the performance of the models in predicting unseen
time horizons to simulate a future prediction scenario. The training set
encompasses $200$ points in the time interval $t=[0,200]$, while the test set
includes $200$ points within the extended time interval $t=[0,300]$. This
experiment represents a challenging task, examining the models’ generalisation
capability and the ability to extrapolate beyond the observed temporal range;

  * •

Completion: evaluates the effectiveness of the models in missing data
imputation, simulating a scenario where data are missing and need to be
estimated. The training set comprises $200$ points in the time interval
$t=[0,200]$, and the test set includes $400$ points within the time interval
$t=[0,200]$.

Furthermore, four different fractional derivative orders, specifically
$0.3,0.5,0.8,0.99$, were used to generate four FDE datasets for each system.
This selection of fractional derivative orders was made to systematically
examine the capabilities and limitations of the Neural FDE in diverse
scenarios. In addition, one ODE dataset was also generated to evaluate the
performance of Neural FDE in modelling ODE generated data. The Neural ODE was
used as a baseline in all experiments and its performance was also tested
using the ODE’s and FDE’s datasets.

To accommodate the randomness of the optimisation process, three independent
runs were conducted for each architecture. Model assessment relied on the
average Mean Squared Error ($MSE_{\text{avg}}$) and its standard deviation
(std).

In all experiments we used the same architecture $\boldsymbol{f_{\theta}}$ for
the Neural ODE and Neural FDE with: an input layer with a neuron and
hyperbolic tangent (tanh) activation function; 2 hidden layers with 64 neurons
and tanh activation function; and an output layer with a neuron. Since the
Neural FDE has another NN, $\alpha_{\boldsymbol{\phi}}$, we use a one neuron
input layer with tanh activation function, a hidden layer with 32 neurons with
tanh activation function, and an output layer with a neuron and a sigmoid
activation function. We always initialise the $\alpha$ value at $0.99$, in the
iterative procedure.

All implementations were done in Pytorch and the Neural ODE was used through
Torchdiffeq [11]. FDE solvers, namely PECE (herein designated PC) [17] were
intentionally developed for this work. These solvers were adapted to ensure
compatibility with NNs and backpropagation through _autograd_ , and they are
CUDA-enabled for efficient GPU utilisation. This tailored implementation
allows for the seamless integration of FDE solvers into the PyTorch framework.
Additionally, the implementation of the PC solver can be used solely for
solving FDEs without using NNs and with CUDA acceleration. To the best of our
knowledge, this is the first publicly available implementation of the PC
solver for FDEs in _Pytorch_ and being CUDA-enabled.

The numerical experiments for the RC and SD systems can be found in Appendix
B.1 and Appendix B.2, respectively.

All computations were performed in a Google Cloud machine type g2-standard-32
with an Intel Cascade Lake 32 vCPU, 128GB RAM and a NVIDIA L4 GPU. The
implementations can be found at [LINK]111available after acceptance.

Remark 6: Please note that the main purpose of this section is not to
establish superiority between the different neural network architectures. Each
architecture has its own distinct advantages and disadvantages, which depend
on the dynamics of the system. Overall, we will see that the Neural FDE
outperforms the Neural ODE, but it’s important to recognise that the mesh
refinement employed during training and testing significantly impacts the
final results, making a fully equitable comparison challenging (see Appendix A
for more details). Consequently, the results presented here should be viewed
as a parametric examination of the behaviour of both Neural ODEs and Neural
FDEs using various datasets. This serves as a means to assess and explore the
methodologies outlined in the preceding sections.

###  4.1 Relaxation Oscillation Process

Relaxation oscillations are a phenomenon observed in various natural and
engineered systems, where the dynamics alternate between periods of slow
relaxation and rapid oscillations. Consider a system governed by the
fractional-order relaxation oscillation equation:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}x(t)+x(t)=1,\quad x(0)=x_{0},$ (14)

where $\leftindex^{C}_{0}{D}^{\alpha}_{t}$ represents a fractional derivative
operator of order $\alpha$ and $x(t)$ represents the variable that undergoes
relaxation oscillations. The initial condition $x_{0}$ sets the initial value
of $x$ at $t=0$, determining the starting point of the oscillatory behaviour
(we consider $x(t_{0})=0.3$). In this fractional-order relaxation oscillation
process, the variable $x(t)$ evolves over time. The system’s behaviour
alternates between periods of slow relaxation and rapid oscillation as it
approaches the equilibrium solution $x(t)=1$.

Considering the well known limitations in the convergence of the numerical
method [17] for highly singular solutions, in this first case study we focused
solely on datasets corresponding to fractional-order values of $\alpha=0.8$,
$0.99$, and $1$ (ODE). This means that the differences between the different
datasets may be smaller when compared to the case studies presented in the
future subsections.

The results are organised in Table 1 and Table 2 for the Neural ODE and the
Neural FDE, respectively. Additionally, the learnt $\alpha$ values by the
Neural FDE, for the three runs, are presented in Table 3. The datasets
$P_{\alpha=0.8}$, $P_{\alpha=0.99}$, and $P_{ODE}$ are constructed with their
respective $\alpha$ values. For $P_{ODE}$, $\alpha$ is set to 1.

Table 1: Performance of Neural ODE when modelling the RO regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.8}$ |  2.69E-2 $\pm$ 3.63E-2 |  6.40E-2 $\pm$ 6.42E-2 |  2.51E-2 $\pm$ 3.42E-2
$P_{\alpha=0.99}$ |  9.52E-2 $\pm$ 1.05E-1 |  8.40E-2 $\pm$ 9.50E-2 |  1.13E-1 $\pm$ 1.19E-1
$P_{\text{ODE}}$ |  3.47E-1 $\pm$ 3.65E-1 |  3.50E-1 $\pm$ 3.68E-1 |  3.47E-1 $\pm$ 3.65E-1
Table 2: Performance of Neural FDE when modelling the RO regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.8}$ |  3.72E-4 $\pm$ 1.31E-4 |  1.43E-2 $\pm$ 6.32E-4 |  3.85E-4 $\pm$ 1.26E-4
$P_{\alpha=0.99}$ |  3.95E-4 $\pm$ 1.44E-5 |  6.36E-4 $\pm$ 6.14E-5 |  1.12E-3 $\pm$ 3.64E-5
$P_{\text{ODE}}$ |  8.80E-5 $\pm$ 1.47E-5 |  2.36E-4 $\pm$ 1.10E-4 |  1.06E-4 $\pm$ 2.90E-5

Table 1 and Table 2 illustrate that the proposed Neural FDE exhibits
significantly superior performance, for all datasets, compared to Neural ODE,
with $MSE_{avg}$ being lower by at least two orders of magnitude. It is
noteworthy that Neural FDE achieved a reduction of at least three orders of
magnitude in error when fitting to the dataset generated by the ODE, in
contrast to Neural ODE. This evidences the performance improvement over Neural
ODEs. Although, it should be mentioned that this results were obtained with no
mesh refinement in both training and testing. This justifies the fact that the
Neural ODE is presenting a higher error when predicting results with a model
trained with its _own dataset_ , $P_{ODE}$. For a more refined mesh, we would
expect for the Neural ODE to perform better than the Neural FDE, when
considering the dataset $P_{ODE}$ (see Appendix A). The disparity in the
errors obtained will be reduced in the upcoming case studies.

Moreover, Figure 4 highlights the notably faster convergence speed of Neural
FDEs compared to Neural ODEs, underscoring the superiority of our proposed
Neural FDE as the preferred choice. It is important to note that the 200
iterations mentioned are unrelated to the 200 data points used in constructing
some datasets.

Table 3: Adjusted $\alpha$ at each run for each dataset of the RO system. | $\alpha=0.8$ | $\alpha=0.99$ | ODE
---|---|---|---
1 | 0.3885 | 0.3451 | 0.3015
2 | 0.5477 | 0.3923 | 0.4431
3 | 0.4472 | 0.4684 | 0.3692

The $\alpha$ values learnt by the Neural FDE, Table 3, show no significant
patterns beyond staying relatively low, _i.e._ lower than $0.5$. The expected
outcome was to achieve $\alpha$ values closer to those used in the training
set. The observed lower values can be attributed to inadequate mesh
refinement, too few training iterations, and the potential minimal impact of
the $\alpha_{\phi}$ NN on the loss function (something that should be
addressed in the future).

Figure 4: Evolution of the loss during training for the three datasets
$P_{\alpha}=0.8,~{}0.99,~{}1$ of the RO system (loss vs. iterations).

###  4.2 Population Growth

Consider a population of organisms that follows a fractional-order logistic
growth. The population size $P(t)$ at time $t$ is governed by the following
fractional-order differential equation:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}P(t)=rP(t)\left(1-\frac{P(t)}{K}\right),\quad
P(t_{0})=100,$ (15)

where $\leftindex^{C}_{0}{D}^{\alpha}_{t}$ represents the fractional
derivative of order $\alpha$, $r=0.1$ is the growth rate, and $K=1000$ is the
carrying capacity of the environment. Consider the initial condition
$P(t_{0})=100$. Please note that when employing fractional derivatives, the
model parameters no longer have classical dimensions, as is the case when
$\alpha=1$. These modified parameters are commonly referred to as
_quasiproperties_.

This equation encapsulates the intricate dynamics of the population’s growth,
incorporating both fractional-order derivatives and the logistic growth model.
The term $rP(t)\left(1-\frac{P(t)}{K}\right)$ characterises the population’s
net reproduction rate as a function of its current size relative to the
environment’s carrying capacity. The fractional derivative
$D_{t}^{\alpha}P(t)$ introduces a nuanced temporal aspect, accounting for non-
integer-order rates of change in population size over time.

The results are organised in Table 4 and Table 5 for the Neural ODE and the
Neural FDE, respectively. Additionally, the learnt $\alpha$ values by the
Neural FDE, for the three runs, are shown in Table 6. To analyse and compare
the convergence of Neural ODE and Neural FDE, the loss values during training
were plotted in Figure 5.

Table 4: Performance of Neural ODE when modelling the PG regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.3}$ |  7.70E-3 $\pm$ 2.15E-5 |  1.04E-2 $\pm$ 1.54E-5 |  7.67E-3 $\pm$ 2.15E-5
$P_{\alpha=0.4}$ |  3.90E-2 $\pm$ 2.45E-2 |  4.17E-2 $\pm$ 1.91E-2 |  3.90E-2 $\pm$ 2.45E-2
$P_{\alpha=0.5}$ |  4.91E-2 $\pm$ 6.26E-3 |  5.57E-2 $\pm$ 5.85E-3 |  4.92E-2 $\pm$ 6.23E-3
$P_{\alpha=0.8}$ |  6.95E-2 $\pm$ 7.11E-3 |  6.74E-2 $\pm$ 1.37E-2 |  6.94E-2 $\pm$ 7.19E-3
$P_{\alpha=0.99}$ |  4.58E-2 $\pm$ 1.50E-3 |  3.45E-2 $\pm$ 4.44E-4 |  4.58E-2 $\pm$ 1.50E-3
$P_{\text{ODE}}$ |  8.33E-2 $\pm$ 5.19E-2 |  7.91E-2 $\pm$ 6.29E-2 |  8.34E-2 $\pm$ 5.19E-2
Table 5: Performance of Neural FDE when modelling the PG regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.3}$ |  1.58E-3 $\pm$ 7.34E-4 |  2.97E-3 $\pm$ 9.06E-4 |  1.96E-3 $\pm$ 7.33E-4
$P_{\alpha=0.4}$ |  4.03E-3 $\pm$ 7.78E-4 |  7.36E-3 $\pm$ 1.25E-3 |  5.01E-3 $\pm$ 6.51E-4
$P_{\alpha=0.5}$ |  1.99E-2 $\pm$ 1.29E-3 |  2.45E-2 $\pm$ 1.61E-3 |  2.06E-2 $\pm$ 1.30E-3
$P_{\alpha=0.8}$ |  1.69E-2 $\pm$ 9.29E-3 |  1.62E-2 $\pm$ 4.83E-3 |  1.77E-2 $\pm$ 9.55E-3
$P_{\alpha=0.99}$ |  1.32E-2 $\pm$ 2.05E-3 |  8.97E-3 $\pm$ 9.89E-4 |  1.45E-2 $\pm$ 1.91E-3
$P_{\text{ODE}}$ |  1.15E-2 $\pm$ 1.49E-3 |  7.84E-3 $\pm$ 9.49E-4 |  1.30E-2 $\pm$ 1.37E-3

Table 4 and Table 5 illustrate that the performance of both the Neural ODE and
the Neural FDE is comparable in the tasks of reconstruction and completion.
However, Neural FDE outperforms Neural ODE in the extrapolation task,
exhibiting a lower $MSE_{avg}$, being similar for $P_{\alpha=0.5}$ and
$P_{\alpha=0.8}$. This superiority underscores the significance of the memory
mechanism employed by Neural FDEs, which allows models to leverage the entire
historical context of the time series for making predictions. Furthermore,
this memory scheme empowers Neural FDEs to capture more intricate
relationships within the data, consequently achieving enhanced performance,
particularly in challenging tasks involving predictions for unseen time
horizons.

From Figure 5, it is evident that Neural FDE exhibits significantly faster
convergence compared to Neural ODE, achieving lower loss values in
considerably less time. This rapid convergence, coupled with its superior
performance, positions Neural FDE as a highly competitive network.

Table 6: Adjusted $\alpha$ at each run for each dataset of the PG regularly-sampled system. | $\alpha=0.3$ | $\alpha=0.4$ | $\alpha=0.5$ | $\alpha=0.8$ | $\alpha=0.99$ | ODE
---|---|---|---|---|---|---
1 | 0.2792 | 0.214 | 0.3582 | 0.4383 | 0.3061 | 0.2824
2 | 0.367 | 0.2369 | 0.3268 | 0.2581 | 0.306 | 0.2847
3 | 0.2849 | 0.2363 | 0.3917 | 0.411 | 0.3631 | 0.3374

Once again, the $\alpha$ values learnt by Neural FDE, Table 6, show no
significant patterns beyond staying relatively low, _i.e._ lower than $0.4$.

Figure 5: Evolution of the loss during training for the four datasets of the
PG system using a regularly sampled dataset (loss vs. iterations).

##  5 Discussion and Conclusion

Considering the concepts of Neural ODE [1] and fractional calculus in Neural
systems [14], we have developed the Neural FDE to deal with time-series data.
This is a FDE featuring a Caputo fractional derivative on the left-hand side
(representing the variation of the state of a dynamical system) and a
_function_ ($\boldsymbol{f_{\theta}}$) of the weights and biases of a NN on
the right-hand side. The goal is to accurately fit the provided data. The
Neural FDE is specifically designed for dynamical systems with strong memory
effects. During the Neural FDE training, we use a specialised numerical solver
for FDEs to compute the solution of the pseudo FDE. Subsequently, a loss
function is employed to compare the numerically predicted outcomes with the
ground truth. Using _autograd_ for backpropagation, we adjust the weights and
biases within the Neural Network $\boldsymbol{f_{\theta}}$ to minimise errors.
It is essential to note that the use of a fractional operator introduces an
additional parameter, namely the order of the fractional derivative denoted as
$\alpha$. Consequently, apart from adjusting the weights and biases in the NN
$\boldsymbol{f_{\theta}}$, we also incorporate a separate NN to determine the
optimal value for $\alpha$. This implies that, during the optimisation
process, adjustments are made to the weights and biases in two distinct NNs.

In the study by Chen et al. [1], the explanation regarding the numerical
methods for solving differential equations was somewhat unclear. Therefore, in
this work, we provide a clear explanation of how mesh refinement can be
performed and how this refinement can be applied to handle irregularly sampled
data.

To assess the performance of the architectural framework proposed in this
study, a comprehensive series of experiments was carried out on a variety of
datasets. First, four synthetic toy datasets were generated by numerically
solving their corresponding differential equations. These datasets were
specifically designed to model four distinct and well-established systems
governed by differential equations: Relaxation Oscillation Process (RO),
Population Growth (PG), RC Circuit Capacitor (RC), and Substance Decay (SD).

The experiments demonstrate that Neural FDE exhibits either significantly
better or similar performance, as measured by $MSE_{avg}$, compared to Neural
ODE. However, even in cases where the produced models demonstrate similar
performance, Neural FDE showcases substantially faster convergence, achieving
lower loss values in significantly fewer iterations than Neural ODE. Although,
it should be mentioned that this results were obtained with no mesh refinement
in both training and testing. This justifies the fact that in some cases the
Neural ODE presents a higher error when predicting results with a model
trained with its _own generated dataset_. For a more refined mesh, we would
expect for the Neural ODE to perform better than the Neural FDE in this
particular case. Note also that in some cases the computational times of the
Neural FDE are prohibitive, and therefore, a mesh refinement study is provided
in Appendix A for a very simple case study.

Although the Neural FDE presented in this work behaves well for the different
case studies considered, there is plenty of room for improvements:

  * •

The use of FDEs needs the computation of the entire history of the dynamical
system at every time step, requiring the storage of all variables.
Additionally, we employed _autograd_ for backpropagation. This results in a
significantly higher computational cost compared to Neural ODEs.

  * •

When computing the loss function, we consider the weights and biases for the
NN $\boldsymbol{f_{\theta}}$, as well as weights and biases for the NN
$\alpha_{\boldsymbol{\phi}}$, optimising the order of the fractional
derivative. Our numerical results indicate that this loss function does not
yield optimal orders for the fractional derivative. Therefore, in future work,
we should explore a different loss function that prioritises minimising the
error while also optimising the order of the fractional derivative.

  * •

While the study by Chen et al. [1] does not explicitly address it, Neural ODEs
exhibit stability issues. These problems are likewise present in Neural FDEs
and stem from significant variations in weight matrices as they evolve in
time. In the future, it is advisable to enhance the model’s stability by
incorporating constraints within the loss function. These constraints would
enable better control over the smoothness of the weight matrix evolution.

  * •

The Predictor-Corrector numerical method, employed to solve the FDE [17],
exhibits a convergence order of $\mathcal{O}(h^{p})$, where $h$ denotes the
time-step and $p=min(2,1+\alpha)$. In simpler terms, this implies that lower
values of $\alpha$ can yield convergence orders of 1. It’s worth noting that
the dependence of the convergence order on $\alpha$ is a common
characteristic. To address this, one potential solution is the use of more
robust numerical methods or the implementation of graded meshes, especially in
the vicinity of the singularity at $t=0$.

Acknowledgements

The authors acknowledge the funding by Fundação para a Ciência e Tecnologia
(Portuguese Foundation for Science and Technology) through CMAT projects
UIDB/00013/2020 and UIDP/00013/2020 and the funding by FCT and Google Cloud
partnership through projects CPCA-IAC/AV/589164/2023 and CPCA-
IAC/AF/589140/2023.

C. Coelho would like to thank FCT for the funding through the scholarship with
reference 2021.05201.BD.

This work is also financially supported by national funds through the
FCT/MCTES (PIDDAC), under the project 2022.06672.PTDC - iMAD - Improving the
Modelling of Anomalous Diffusion and Viscoelasticity: solutions to industrial
problems.

## References

  *   * Chen et al. [2018] Chen, R.T., Rubanova, Y., Bettencourt, J., Duvenaud, D.K.: Neural ordinary differential equations. Advances in neural information processing systems 31 (2018)
  * Podlubny [1999] Podlubny, I.: Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Mathematics in science and engineering 198, 1–340 (1999)
  * Herrmann [2014] Herrmann, R.: Fractional Calculus: An Introduction for Physicists, 2\. ed edn. World Scientific, New Jersey, NJ (2014). https://doi.org/10.1142/9789814551083
  * Coelho et al. [2024] Coelho, C., Costa, M.F.P., Ferrás, L.L.: Tracing footprints: Neural networks meet non-integer order differential equations for modelling systems with memory. In: Tiny Papers @ ICLR (2024). https://openreview.net/forum?id=8518dcW4hc
  * He et al. [2016] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016)
  * Haber et al. [2018] Haber, E., Ruthotto, L., Holtham, E., Jun, S.-H.: Learning across scales—multiscale methods for convolution neural networks. Proceedings of the AAAI Conference on Artificial Intelligence 32(1) (2018) https://doi.org/10.1609/aaai.v32i1.11680
  * Haber and Ruthotto [2017] Haber, E., Ruthotto, L.: Stable architectures for deep neural networks. Inverse Problems 34(1), 014004 (2017) https://doi.org/10.1088/1361-6420/aa9a90
  * E [2017] E, W.: A proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics 5(1), 1–11 (2017) https://doi.org/10.1007/s40304-017-0103-z
  * Lu et al. [2018] Lu, Y., Zhong, A., Li, Q., Dong, B.: Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. In: Dy, J., Krause, A. (eds.) Proceedings of the 35th International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 80, pp. 3276–3285. PMLR, ??? (2018). https://proceedings.mlr.press/v80/lu18d.html
  * Ruthotto and Haber [2019] Ruthotto, L., Haber, E.: Deep neural networks motivated by partial differential equations. Journal of Mathematical Imaging and Vision 62(3), 352–364 (2019) https://doi.org/10.1007/s10851-019-00903-1
  * Chen [2018] Chen, R.T.Q.: torchdiffeq (2018). https://github.com/rtqichen/torchdiffeq
  * Diethelm [2010] Diethelm, K.: The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type. Springer, ??? (2010). https://doi.org/10.1007/978-3-642-14574-2 . http://dx.doi.org/10.1007/978-3-642-14574-2
  * Liu et al. [2018] Liu, Y., Roberts, J., Yan, Y.: A note on finite difference methods for nonlinear fractional differential equations with non-uniform meshes. International Journal of Computer Mathematics 95(6-7), 1151–1169 (2018) https://doi.org/10.1080/00207160.2017.1381691
  * Lundstrom et al. [2008] Lundstrom, B.N., Higgs, M.H., Spain, W.J., Fairhall, A.L.: Fractional differentiation by neocortical pyramidal neurons. Nature Neuroscience 11(11), 1335–1342 (2008) https://doi.org/10.1038/nn.2212
  * Diethelm and Ford [2002] Diethelm, K., Ford, N.J.: Analysis of fractional differential equations. Journal of Mathematical Analysis and Applications 265(2), 229–248 (2002) https://doi.org/10.1006/jmaa.2000.7194
  * Barros et al. [2021] Barros, L.C.d., Lopes, M.M., Pedro, F.S., Esmi, E., Santos, J.P.C.d., Sánchez, D.E.: The memory effect on fractional calculus: an application in the spread of covid-19. Computational and Applied Mathematics 40(3) (2021) https://doi.org/10.1007/s40314-021-01456-z
  * Diethelm et al. [2002] Diethelm, K., Ford, N.J., Freed, A.D.: A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations. Nonlinear Dynamics 29(1), 3–22 (2002)
  * Kaslik and Sivasundaram [2012] Kaslik, E., Sivasundaram, S.: Nonlinear dynamics and chaos in fractional-order neural networks. Neural Networks 32, 245–256 (2012) https://doi.org/10.1016/j.neunet.2012.02.030

##  Appendix A Numerical Method and Mesh Refinement

Consider solving an ODE IVP (Equation (5)) and the FDE IVP (Equation (11)),
with initial condition $\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$ for time-
steps $(t_{0},t_{1},\dots,t_{f})$. When employing a fixed-step solver, the
time interval is discretized according to the time-steps provided (Algorithm
3):

Algorithm 3 Solving an ODE with a fixed-step solver.

Input: discretization time steps $t=(t_{0},t_{1},\dots,t_{f})$, initial
condition $\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$;

${\\{\boldsymbol{h}_{t}\\}}_{t=t_{0},t_{1},\dots,t_{f}}\leftarrow
ODESolve(\boldsymbol{f_{\theta}},\boldsymbol{h}_{0},(t_{0},t_{1},\dots,t_{f}))$;

Return: ${\\{\boldsymbol{h}_{t}\\}}_{t=t_{0},t_{1},\dots,t_{f}}$;

Now, consider solving the same ODE IVP using a fixed-step solver across the
time interval $[t_{0},t_{f}]$ using a step size smaller than $t_{1}-t_{0}$.
This process involves dividing the interval into sub-intervals
$(t_{n},t_{n+1})$ with a uniform step size $\Delta t$, resulting in
$N=(t_{n+1}-t_{n})/\Delta t$ time steps per sub-interval. Subsequently, the
solver proceeds to solve the ODE within each sub-interval (Algorithm 4). If,
however, the selected step size fails to encompass the endpoint $t_{n+1}$, the
state is computed by interpolation (see Figure 2).

Algorithm 4 Solving an ODE with a fixed-step solver with a chosen step size.

Input: time steps $t=(t_{0},t_{1},\dots,t_{f})$, initial condition
$\boldsymbol{h}(t_{0})=\boldsymbol{h}_{0}$, step size $\Delta t$;

for $n=1:f$ do

$N\leftarrow\dfrac{t[n]-t[n-1]}{\Delta t};$

${\\{\boldsymbol{h}_{t}\\}}_{t=t[n-1],t[n]}\leftarrow
ODESolve(\boldsymbol{f_{\theta}},\boldsymbol{h}_{0},(t[n-1],t[n-1]+\Delta
t,\dots,t[n-1]+\Delta t\times N))$;

end for

Return: ${\\{\boldsymbol{h}_{t}\\}}_{t=t_{0},t_{1},\dots,t_{f}}$; (linear
interpolation may be needed)

The choice of the step size is of paramount importance. Smaller $\Delta t$ in
the discretization process using a fixed-step solver with sub-interval
discretization, Algorithm 4, offers a pathway to enhancing the accuracy and
reliability of the approximated solution. By decreasing $\Delta t$, each sub-
interval is divided into more steps, becoming finer. This finer resolution
enables the solver to capture intricate features and rapid changes in the
solution that may occur within shorter time frames. Consequently, the computed
solution becomes more faithful to the true behaviour of the system, especially
in regions of high variability or steep gradients. Also, smaller $\Delta t$
helps mitigate errors inherent to numerical methods. Thus, selecting the
appropriate $\Delta t$ poses a significant challenge due to the trade-off
between accuracy of the solution and computational cost. Moreover, the optimal
$\Delta t$ may vary throughout the solution domain, complicating the task of
selecting a single value that adequately balances accuracy and efficiency
across the entire computation. To remove the need for selecting an appropriate
fixed step size, adaptive-step solvers were introduced (these are only
available in this work for the Neural ODE). Similar to Algorithm 4, adaptive-
step solvers divide $t$ into sub-intervals however, $\Delta t$ is computed and
dynamically adjusted during the integration process, based on the local
behaviour of the solution. This adaptability allows the solver to allocate
smaller steps in regions where the solution changes rapidly or exhibits
intricate features, while larger steps are used in smoother regions where the
solution varies more gradually. Thus, when employing Algorithm 3, if the step
size inherent to $t$, given by $t_{1}-t_{0}$, is substantial then the
solutions approximation may be inaccurate and most likely lack the accuracy of
the approximations done with Algorithm 4 and an adaptive-step solver.
Furthermore, when the step size of a fixed-step solver is sufficiently small,
the approximation accuracy can reach, and in some cases exceed, that achieved
by adaptive-step solvers.

To demonstrate how the choice of step size influences the fit of Neural ODEs
and Neural FDEs to the system under study, we conducted some simple
experiments. We used the PG system to generate four datasets, each covering a
smaller time domain $t=[0,10]$ with $50$ data points. These datasets were
created in two variations: two were derived by solving an ODE with step sizes
of $1$ and $0.1$, while the other two were obtained by solving an FDE with
$\alpha=0.8$ using step sizes of $1$ and $0.1$. Subsequently, we employed both
a Neural ODE to learn the ODE generated datasets and a Neural FDE to fit the
FDE generated datasets. We then performed an analysis of the plots depicting
the ground truth and predicted curves, as illustrated in Figures 7-10.

From Figures 7-10, it can be seen that smaller $\Delta t$ results in a more
accurate approximation of the system’s behaviour. This finer granularity also
improves the fitting by both the Neural ODE and Neural FDE models.

Figure 6: ODE dataset with $\Delta t=1$.

Figure 7: ODE dataset with $\Delta t=0.1$.

Figure 8: FDE dataset $\alpha=0.8$ with $\Delta t=1$.

Figure 9: FDE dataset $\alpha=0.8$ with $\Delta t=0.1$.

Figure 10: Modelling PG datasets with Neural FDE.

##  Appendix B Additional Numerical Experiments

###  B.1 RC Circuit Capacitor

Consider an RC circuit with a fractional-order capacitor. The fractional-order
equation describing the circuit is given by:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}V(t)=\frac{1}{RC}V(t),\quad V(t_{0})=10,$
(16)

where $\alpha$ is the fractional order between 0 and 1, $V(t)$ is the voltage
across the capacitor, $R=100$ is the resistance, and $C=0.2$ is the fractional
capacitor. The initial condition for the voltage is $V(t_{0})=10$.

This equation captures the intricate dynamics of the circuit’s behaviour,
combining the effects of fractional-order elements with the fundamental
principles of RC circuits. The term $\frac{1}{RC}V(t)$ represents the rate of
change of voltage across the fractional-order capacitor, illustrating how the
circuit responds to the interplay between resistance ($R$), fractional
capacitance ($C$), and the time-dependent voltage ($V(t)$).

The results are organised in Table 7 and Table 8 for the Neural ODE and Neural
FDE, respectively. Additionally, the learnt $\alpha$ values by Neural FDE, for
the three runs, are shown in Table 9. To analyse and compare the convergence
of Neural ODE and Neural FDE, the loss values during training were plotted in
Figure 11.

Table 7: Performance of Neural ODE when modelling the RC regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.3}$ |  5.84E-2 $\pm$ 7.88E-2 |  5.64E-2 $\pm$ 7.44E-2 |  5.85E-2 $\pm$ 7.89E-2
$P_{\alpha=0.4}$ |  9.56E-3 $\pm$ 1.34E-4 |  1.50E-2 $\pm$ 1.01E-4 |  9.56E-3 $\pm$ 1.34E-4
$P_{\alpha=0.5}$ |  2.64E-2 $\pm$ 6.06E-5 |  4.23E-2 $\pm$ 1.27E-4 |  2.64E-2 $\pm$ 6.09E-5
$P_{\alpha=0.8}$ |  6.25E-2 $\pm$ 9.16E-5 |  5.63E-2 $\pm$ 7.09E-4 |  6.22E-2 $\pm$ 8.63E-5
$P_{\alpha=0.99}$ |  5.35E-2 $\pm$ 1.62E-2 |  4.86E-2 $\pm$ 2.23E-2 |  5.33E-2 $\pm$ 1.63E-2
$P_{\text{ODE}}$ |  4.11E-2 $\pm$ 1.19E-5 |  3.14E-2 $\pm$ 9.43E-07 |  4.08E-2 $\pm$ 1.15E-5
Table 8: Performance of Neural FDE when modelling the RC regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.3}$ |  1.04E-3 $\pm$ 5.66E-4 |  1.70E-3 $\pm$ 8.13E-4 |  1.21E-3 $\pm$ 5.06E-4
$P_{\alpha=0.4}$ |  2.83E-3 $\pm$ 1.03E-3 |  6.39E-3 $\pm$ 1.50E-3 |  3.32E-3 $\pm$ 1.03E-3
$P_{\alpha=0.5}$ |  8.08E-3 $\pm$ 2.11E-3 |  2.11E-2 $\pm$ 3.40E-3 |  9.22E-3 $\pm$ 2.01E-3
$P_{\alpha=0.8}$ |  4.31E-2 $\pm$ 5.17E-3 |  4.54E-2 $\pm$ 2.41E-3 |  4.47E-2 $\pm$ 5.19E-3
$P_{\alpha=0.99}$ |  3.77E-2 $\pm$ 2.28E-3 |  2.97E-2 $\pm$ 1.18E-3 |  3.79E-2 $\pm$ 1.97E-3
$P_{\text{ODE}}$ |  3.35E-2 $\pm$ 9.63E-4 |  2.73E-2 $\pm$ 2.67E-4 |  3.42E-2 $\pm$ 1.19E-3

Table 7 and Table 8 demonstrate that, overall, the Neural ODE and the Neural
FDE exhibit comparable performance in modelling this system. Notably, Neural
FDE displays lower $MSE_{avg}$ values for lower derivative orders
($P_{\alpha=0.3}$, $P_{\alpha=0.4}$, and $P_{\alpha=0.5}$). However, as
depicted in Figure 11, while both Neural ODE and Neural FDE yield similarly
performant models, Neural FDE demonstrates faster convergence of the loss
function, requiring less training time.

Table 9: Adjusted $\alpha$ at each run for each dataset of the RC regularly-sampled system. | $\alpha=0.3$ | $\alpha=0.4$ | $\alpha=0.5$ | $\alpha=0.8$ | $\alpha=0.99$ | ODE
---|---|---|---|---|---|---
1 | 0.4298 | 0.2956 | 0.2136 | 0.2196 | 0.4585 | 0.2427
2 | 0.4741 | 0.3883 | 0.2287 | 0.2163 | 0.391 | 0.2511
3 | 0.3068 | 0.2229 | 0.3092 | 0.3802 | 0.3226 | 0.3397

Again, the $\alpha$ values learnt by Neural FDE, Table 9, show no apparent
pattern, being always lower than $0.5$.

Figure 11: Evolution of the loss during training for the four datasets of the
RC system (loss vs. iterations).

###  B.2 Substance Decay

Consider a substance undergoing fractional-order decay with the following
equation:

$\leftindex^{C}_{0}{D}^{\alpha}_{t}N(t)=-\lambda N(t),\quad N(t_{0})=100,$
(17)

where $\alpha\in(0,1)$ represents the fractional order, $N(t)$ denotes the
quantity of the decaying substance, and $\lambda=0.5$ is the decay constant.
The initial condition is $N(t_{0})=100$.

This equation captures the intricate dynamics of substance decay,
incorporating both fractional-order derivatives and exponential decay
principles. The term $-\lambda N(t)$ signifies the decay rate proportional to
the current quantity of the substance. The fractional derivative
$\leftindex^{C}_{0}{D}^{\alpha}_{t}$ introduces a nuanced temporal aspect,
considering non-integer-order rates of decay over time.

Interpreting and studying this fractional-order decay model could provide
valuable insights into the behaviour of decaying substances, potentially
offering more accurate predictions and informed decay management strategies.

The results are organised in Table 10 and Table 11 for the Neural ODE and the
Neural FDE, respectively. Additionally, the learnt $\alpha$ values by Neural
FDE, for the three runs, are shown in Table 12. To analyse and compare the
convergence of the Neural ODE and the Neural FDE, the loss values during
training were plotted, Figure 12.

Table 10: Performance of Neural ODE when modelling the SD regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.3}$ |  1.43E-2 $\pm$ 6.40E-3 |  3.59E-2 $\pm$ 6.66E-3 |  2.86E-2 $\pm$ 1.17E-2
$P_{\alpha=0.4}$ |  3.01E-2 $\pm$ 2.84E-2 |  7.35E-2 $\pm$ 6.32E-2 |  2.42E-2 $\pm$ 1.02E-2
$P_{\alpha=0.5}$ |  1.83E-2 $\pm$ 5.68E-3 |  2.32E-2 $\pm$ 8.19E-3 |  1.48E-2 $\pm$ 2.38E-3
$P_{\alpha=0.8}$ |  2.00E-3 $\pm$ 2.00E-4 |  2.57E-3 $\pm$ 1.95E-4 |  1.47E-3 $\pm$ 1.67E-4
$P_{\alpha=0.99}$ |  9.17E-3 $\pm$ 1.24E-2 |  9.47E-3 $\pm$ 1.25E-2 |  9.06E-3 $\pm$ 1.24E-2
$P_{\text{ODE}}$ |  6.53E-3 $\pm$ 8.24E-3 |  6.55E-3 $\pm$ 8.29E-3 |  6.53E-3 $\pm$ 8.24E-3
Table 11: Performance of Neural ODE when modelling the SD regularly-sampled system ($MSE_{avg}$ $\pm$ std). DATASET | RECONSTRUCTION | EXTRAPOLATION | COMPLETION
---|---|---|---
$P_{\alpha=0.3}$ |  8.78E-4 $\pm$ 4.10E-4 |  4.73E-2 $\pm$ 2.73E-3 |  4.55E-3 $\pm$ 2.94E-4
$P_{\alpha=0.4}$ |  2.17E-3 $\pm$ 2.01E-3 |  3.44E-2 $\pm$ 5.25E-3 |  7.90E-3 $\pm$ 1.75E-3
$P_{\alpha=0.5}$ |  1.47E-3 $\pm$ 4.56E-4 |  3.77E-3 $\pm$ 3.22E-4 |  1.42E-3 $\pm$ 1.54E-4
$P_{\alpha=0.8}$ |  3.35E-3 $\pm$ 3.27E-3 |  4.21E-3 $\pm$ 3.88E-3 |  2.66E-3 $\pm$ 2.77E-3
$P_{\alpha=0.99}$ |  8.07E-3 $\pm$ 4.96E-3 |  9.00E-3 $\pm$ 5.16E-3 |  7.06E-3 $\pm$ 4.77E-3
$P_{\text{ODE}}$ |  1.35E-3 $\pm$ 5.95E-4 |  1.86E-3 $\pm$ 6.62E-4 |  9.62E-4 $\pm$ 4.75E-4

Table 10 and Table 11 indicate that at the reconstruction task, Neural FDE
exhibits superior performance for datasets generated with FDEs with
$\alpha\leq 0.5$, while showing similar performance for others. Overall,
Neural FDE also demonstrates better performance at the completion task
compared to Neural ODE. Both Neural ODE and Neural FDE demonstrate similar
performance in extrapolation. Although the performance of Neural FDE in
modelling this system may not show a substantial improvement in terms of
$MSE_{avg}$, Figure 12 reveals that once again, the convergence rate of Neural
FDE is significantly higher, leading to the need of fewer iterations when
compared to Neural ODE, while achieving comparable performance.

Table 12: Adjusted $\alpha$ at each run for each dataset of the SD regularly-sampled system. | $\alpha=0.3$ | $\alpha=0.4$ | $\alpha=0.5$ | $\alpha=0.8$ | $\alpha=0.99$ | ODE
---|---|---|---|---|---|---
1 | 0.2535 | 0.4945 | 0.4228 | 0.6092 | 0.4132 | 0.759
2 | 0.3598 | 0.508 | 0.3512 | 0.4379 | 0.5958 | 0.7184
3 | 0.4019 | 0.435 | 0.3446 | 0.6627 | 0.5444 | 0.6529

Unlike for the other systems, the $\alpha$ values learned by Neural FDE are
generally higher, but no discernible pattern emerges, as shown in Table 12. As
said before, the potential minimal impact of the $\alpha_{\phi}$ NN on the
loss function may be influencing the final $\alpha$ values obtained.

Figure 12: Evolution of the loss during training for the four datasets of the
SD system (loss vs. iterations).