Neural Fractional Differential Equations: Optimising the Order of the Fractional Derivative

1. Introduction

Real-world systems are often modelled using integral/differential equations, which are then numerically solved to predict the system behaviour and evolution. This process can be time-consuming, as numerical simulations sometimes take months, and, finding the correct model parameters is often challenging. However, with significant advancements in Neural Networks (NNs) that can learn patterns, real-world systems are increasingly being modelled using a combination of integral/differential models and NNs, or even NNs alone [1,2,3,4].

Neural Ordinary Differential Equations (Neural ODEs) were introduced in 2018 [5] (see also [6,7]) as a continuous version of the discrete Residual Neural Networks, and claimed to offer a continuous modelling solution for real-world systems that incorporate time-dependence, mimicking the dynamics of that system using only discrete data. Once trained, the Neural ODEs result in a hybrid ODE (part analytical, part NN-based) that can be used for making predictions, by numerically solving the resulting ODEs. The numerical solution of these hybrid ODEs is significantly simpler and less time-consuming compared to the numerical solution of complex governing equations, making Neural ODEs an excellent choice for modelling time-dependent real-world systems. However, the simplicity of ODEs sometimes limits their effectiveness in capturing complex behaviours characterised by intricate dynamics, non-linear interactions, and memory. To address this, Neural Fractional Differential Equations (Neural FDEs) were recently proposed [8,9].

Neural FDEs, as described by Equation (1), are a NN architecture designed to fit the solution $\mathbf{h}\left(t\right)$ to given data $\{{\mathbf{x}}_0,{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N\}$ (for example, experimental data), over a specified time range $[t_0,T]$. The Neural FDE combines an analytical part, ${ }_0^C{D}_t^{\alpha }\mathbf{h}\left(t\right)$, with a NN-based part, ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$, leading to the initial value problem,

$$\left\{\begin{array}{c}{ }_0^C{D}_t^{\alpha }\mathbf{h}\left(t\right)={\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right)),\hfill \\ \mathbf{h}\left(t_0\right)={\mathbf{x}}_0,t\in [t_0,T].\hfill \end{array}\right.$$

(1)

Here, ${ }_0^C{D}_t^{\alpha }g\left(t\right)$ denotes the Caputo fractional derivative [10,11], defined for $0<\alpha <1$ (and considering a generic scalar function $g\left(t\right)$) as:

$${ }_0^C{D}_t^{\alpha }g\left(t\right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{(t-s)}^{-\alpha }{g}^{\prime }\left(s\right)ds,$$

(2)

where $\Gamma$ is the Gamma function.

An important feature of Neural FDEs is their ability to learn not only the optimal parameters $\theta$ of the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$, but also the order of the derivative $\alpha$ (when $\alpha =1$ we obtain a Neural ODE). This is achieved using only information from the time-series dataset $\{{\mathbf{x}}_0,{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N\}$, where each ${\mathbf{x}}_i=(x_i^1,x_i^2,\cdots ,x_i^d)\in {\mathbb{R}}^d$, $i=0,\cdots ,N$ is associated with a time instant $t_i$.

In [9] the $\alpha$ value is learned from another NN ${\alpha }_{\varphi }$ with parameters $\varphi$. Therefore, if $\mathcal{L}(\theta ,\varphi )$ represents the loss function, we can train the Neural FDE by solving the minimisation problem (Section 1). The parameters $\theta$ and $\varphi$ are optimised by minimising the error between the predicted $\widehat{\mathbf{h}}\left(t_i\right),i=1,\cdots ,N$ and ground-truth $\{{\mathbf{x}}_0,{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N\}$ values1:

(3)

The popular Mean Squared Error (MSE) loss function was considered in [9] and also in this work. Here, $\mathrm{FDESolve}({\mathbf{f}}_{\theta },{\mathbf{x}}_0,\{t_0,t_1,\cdots ,t_N\})$ refers to any numerical solver used to obtain the numerical solution $\widehat{\mathbf{h}}\left(t\right)$ for each instant $t_i$.

Since Neural FDEs are a recent research topic, there are no studies on the uniqueness of the parameter $\alpha$ and its interaction with the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$. In [9], the authors provide the values of $\alpha$ learned by Neural FDEs for each dataset, however, a closer examination reveals that these values differ significantly from the ground truth values, which were derived from synthetic datasets. The authors attribute this discrepancy to the approximation capabilities of NNs, meaning that, during training, ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$ adapts to any given $\alpha$ (this is a complex interaction since in [9], $\alpha$ is also learned by another NN). Additionally, $\alpha$ must be initialised in the optimisation procedure, and yet no studies have investigated how the initialisation of $\alpha$ affects the learned optimal $\alpha$ and the overall performance of Neural FDEs.

In this work, we address these key open questions about the order of the fractional derivative $\alpha$ in Neural FDEs. We show that Neural FDEs are capable of modelling data dynamics effectively, even when the learned value of $\alpha$ deviates significantly from the true value. Furthermore, we perform a numerical analysis to investigate how the initialisation of $\alpha$ affects the performance of Neural FDEs.

This paper is organised as follows: In Section 2, we provide a brief overview of FDEs and Neural FDEs, highlighting the theoretical results regarding the existence and uniqueness of solutions. We also discuss how the solution depends on the given data. Section 3 presents a series of numerical experiments on the non-uniqueness of the learnt $\alpha$ values. The paper ends with the discussion and conclusions in Section 4.

2. Neural Networks and Theory of Fractional Initial Value Problems

As shown in the introduction, a Neural FDE is composed of two NN:

• a NN with parameters $\theta$, denoted as ${\mathbf{f}}_{\theta }$, that models the right-hand side of a FDE,

  $${ }_0^C{D}_t^{\alpha }\mathbf{h}\left(t\right)={\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right)),$$

  (4)

  where $\mathbf{h}\left(t\right)$ is the state of the system at time-step t;
• a NN with parameters $\varphi$ (or a learnable parameter), referred to as ${\alpha }_{\varphi }$, that models $\alpha$,

  $$\alpha ={\alpha }_{\varphi }.$$

  (5)

As shown in Figure 1 in this work we constructed ${\mathbf{f}}_{\theta }$ with 3 layers: input layer with 1 neuron and hyperbolic tangent (tanh) activation function; hidden layer with 64 neurons and tanh; output layer with 1 neuron. For learning $\alpha$ we considered a NN ${\alpha }_{\varphi }$ with 3 layers: input layer with 1 neuron and hyperbolic tangent (tanh) activation function; hidden layer with 32 neurons and tanh; output layer with 1 neuron and sigmoid activation function (that helps keeping the value of $\alpha$ in the interval $(0,1)$). For ease of understanding we consider $\widehat{h}\left(t\right))$ to be a scalar in Figure 2 and in Equations (6) and (7). We use $\widehat{h}\left(t\right))$ instead of $h\left(t\right)$, because it is assumed that $h\left(t\right)$ is being evaluated during the training of the Neural FDE, where a numerical solution is needed ($\widehat{h}\left(t\right)$ is a numericall approximation of $h\left(t\right)$[9]). After training, we obtain the final/optimal Neural FDE model and use again the notation ${\mathbf{f}}_{\theta }(t,h\left(t\right))$.

In the NN ${\mathbf{f}}_{\theta }(t,\widehat{h}\left(t\right))$ (Figure 1), the values $w_i$ and $w_i^o$ for $i=1,\cdots ,64$ are the weights of the hidden and output layers, respectively, and b is the bias, $\theta =\{w_1,\cdots ,w_{64},w_1^o,\cdots ,w_{64}^o,b\}$. The output of the NN ${\mathbf{f}}_{\theta }(t,\widehat{h}\left(t\right))$ can be written as,

$$w_1^{o}tanh(w_1\widehat{h}\left(t\right)+b)+\cdots w_{64}^{o}tanh(w_{64}\widehat{h}\left(t\right)+b)=\sum _{i=1}^{64}{w}_i^{o}tanh(w_i\widehat{h}\left(t\right)+b).$$

(6)

An important feature is that there is no activation function in the output layer of ${\mathbf{f}}_{\theta }(t,\widehat{h}\left(t\right))$, allowing the NN to approximate any function $f(t,h(t\left)\right)$. If we opt to use, for example, a tanh in the output layer, it would constrain the fractional derivative ${ }_0^C{D}_t^{\alpha }h\left(t\right)$ to vary only from -1 to 1, thus limiting the fitting capabilities of the Neural FDE. This limitation could be mitigated by normalising the given data $\{{\mathbf{x}}_0,{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N\}$.

Remark: In Figure 1 and Equation (6), we observe h as a function of t. However, the NN depicted on the left side of Figure 1 does not use t as an input. Instead, the NN is called at each iteration of the numerical solver that addresses a discretised version of Equation (1). This solver defines all time instants through a mesh over the interval $[t_0,T]$[9], and, consequently, each evaluation of ${\mathbf{f}}_{\theta }$ is always associated with a specific time instant.

Since for a differentiable function $g\left(h\right)$ with continuous derivative on $[h_1,h_2]$ we have that (using the mean value Theorem),

$$|g\left(h_1\right)-g\left(h_2\right)|\le \underset{h\in \left[h_1,h_2\right]}{max}|g^{\prime }\left(h\right)\left|\right|h_1-h_2|$$

and $|{tanh}^{\prime }\left(h\right)|=|1-{tanh}^2\left(h\right)|\le 1$ ($-1\le tanh\left(h\right)\le 1$), we can say that tanh is 1-Lipschitz, that is,

$$|tanh\left(h_1\right)-tanh\left(h_2\right)|\le 1|h_1-h_2|.$$

Define $g\left(h\right)$ as,

$$g\left(h\right)=\sum _{i=1}^{64}{w}_i^{o}tanh(w_{i}h+b),$$

(7)

which is a weighted sum of 64 tanh functions. Then, $g\left(h\right)$ is L-Lipschitz, with $L<64w_{max}^o$ ($w_{max}^o=max\{w_1^o,\cdots ,w_{64}^o\}$). This property will be important to guarantee the uniqueness of the solution of (1).

2.1. Fractional Differential Equations

We now provide some theoretical results that will be fundamental in understanding the expected behaviour of the Neural FDE model. Consider the following fractional initial value problem,

$$\left\{\begin{array}{c}{ }_0^C{D}_t^{\alpha }z\left(t\right)=f(t,z\left(t\right))\hfill \\ z\left(t_0\right)=z_0.\hfill \end{array}\right.$$

(8)

For ease of understanding we restrict the analyses to the case where $z\left(t\right)$ is a scalar function and $0<\alpha <1$.

2.1.1. Existence and Uniqueness

The following Theorem [11,12] provides information on the existence of a continuous solution to problem (8):

Theorem 1. 

Let $0<\alpha$, $x_0\in \mathbb{R}$, $Q>0$ and $T^{*}>0$. Define $G=$$\left[0,T^{*}\right]\times \left[x_0-Q,x_0+Q\right]$, and let the function $f:G\to \mathbb{R}$ be continuous. Furthermore, we define $M:=\underset{(x,y)\in G}{sup}\left|f(x,y)\right|$ and

$$T:=\left\{\begin{array}{cc}{T}^{*},\hfill & \mathrm{if}M=0\hfill \\ min\left\{T^{*},{\left({\displaystyle \frac{Q\Gamma (\alpha +1)}{M}}\right)}^{1/\alpha }\right\}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

Then, there exists a function $z\left(t\right)\in C^0[0,T]$ solving the initial value problem (8).

Note that this continuous solution may be defined in a smaller interval $[0,T]$ compared to the interval $[0,T^{*}]$ where the function $f(t,z(t\left)\right)$ is defined ($T\le T^{*}$). From Theorem (1), we can infer that high values of $|{\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))|$ (see Equation (1)) decrease the interval within which we can guarantee a continuous solution. However, these conclusions should be approached with caution. As will be shown later, we will only access discrete values of the function ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$ in the numerical solution of (1). Note that $\alpha$ also affects the size of the interval. Its contribution can be either positive or negative, depending on the values of Q and M.

The following Theorem establishes the conditions for which we can guarantee the uniqueness of the solution $z\left(t\right)$ [11,12]:

Theorem 2. 

Let $0<\alpha$, $x_0\in \mathbb{R}$, and $T^{*}>0$. Define the set $G=[0,T^{*}]\times [x_0-Q,x_0+Q]$ and let the function $f:G\to \mathbb{R}$ be continuous and satisfy a Lipschitz condition with respect to the second variable:

$$|f(t,z_1)-f(t,z_2)|\le L|z_1-z_2|,$$

where $L>0$ is a constant independent of t, $z_1$, and $z_2$. Then, there exists a uniquely defined function $z\in C^0[0,T]$ solving the initial value problem (8).

As shown above, the function on the right-hand side of the Neural FDE model is Lipschitz (see Equation (7)). Therefore, we can conclude that the solution to Equation (1) is unique.

2.1.2. Analysing the Behaviour of Solutions with Perturbed Data

Other results of interest for Neural FDEs pertain to the dependencies of the solution on $f(t,z(t\left)\right)$ and $\alpha$. In Neural FDEs, both $f(t,z(t\left)\right)$ and $\alpha$ are substituted by NNs that vary with each iteration of the Neural FDE training [9].

Let $u\left(t\right)$ be the solution of the initial value problem,

$$\left\{\begin{array}{c}{ }_0^C{D}_t^{\alpha }u\left(t\right)=\tilde{f}(t,u\left(t\right))\hfill \\ u\left(t_0\right)=z_0.\hfill \end{array}\right.$$

(9)

where $\tilde{f}(t,u\left(t\right))$ is a perturbed version of f, that satisfies the same hypotheses as $f.$

Theorem 3. 

Let

$$\epsilon :=\underset{(x_1,x_2)\in G}{max}\left|f(x_1,x_2)-\tilde{f}(x_1,x_2)\right|.$$

If ε is sufficiently small, there exists some $T>0$ such that both functions z (Equation (8)) and u (Equation (9)) are defined on $[0,T]$, and we have

$$\underset{0\le t\le T}{sup}|z\left(t\right)-u\left(t\right)|=O\left(\underset{(x_1,x_2)\in G}{max}\left|f(x_1,x_2)-\tilde{f}(x_1,x_2)\right|\right),$$

where $z\left(t\right)$ is the solution of (8).

This Theorem provides insight into how the solution of (1) changes in response to variations in the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$. While the variations in both the solution and the function are of the same order for small changes of the function, it is crucial to carefully interpret these results given the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$ as defined by Equation (7).

When training the Neural ODE, one must solve the optimisation problem (Section 1), where the weights and biases are adjusted until an optimal Neural FDE model is obtained (training 1). If a second, independent training (training 2) of the Neural FDE is conducted with the same stopping criterion for the optimisation process, a new Neural FDE model with different weights and biases may be obtained.

The NN learns model parameters based on a set of ordered data, meaning the number of elements in the set will significantly influence the difference between $z\left(t\right)$ and $u\left(t\right)$ as in Theorem (3). This effect is illustrated in Figure 2, where a training dataset of only two observations can be fitted by two distinct functions.

Therefore, Figure 2 tells us that when modelling with Neural FDEs, it is important to have some prior knowledge of the underlying physics of the data. This is crucial because the number of data points available for the training process might be beyond our control. For instance, the data could originate from real experiments where obtaining results is challenging.

Regarding the influence of the order of the derivative on the solution, we have that [11,12]:

Theorem 4. 

Let $u\left(t\right)$ be the solution of the initial value problem,

$$\left\{\begin{array}{c}{ }_0^C{D}_t^{\tilde{\alpha }}u\left(t\right)=f(t,u\left(t\right))\hfill \\ u\left(t_0\right)=z_0.\hfill \end{array}\right.$$

(10)

where $\tilde{\alpha }$ is a perturbed α value. Let $\epsilon :=\tilde{\alpha }-\alpha$. Then, if ε is sufficiently small, there exists some $T>0$ such that both the functions u and z are defined on $[0,T]$, and we have that

$$\underset{0\le t\le T}{sup}|u\left(t\right)-z\left(t\right)|=O(\tilde{\alpha }-\alpha )$$

where $z\left(t\right)$ is the solution of (8).

Once again, for small changes in $\alpha$, the variations in both the solution and the order of the derivative are of the same order. This property will be explored numerically, and with more detail, later in this work when solving (Section 1). It is important to note that the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$ is not fixed in our problem (1) (its structure is fixed, but the weights and bias change along the training optimisation procedure). However, Theorem (4) assumes that the function $f(t,u(t\left)\right)$ is fixed. Therefore, Theorem (4) gives us an idea of the changes in our solution, but does not allow the full understanding of its variation along the optimisation procedure.

2.1.3. Smoothness of the Solution

For the classical case, $\alpha =1$ in Equation (8), we have (under some hypotheses on the interval $[a,b]$) that if $f\in C^{k-1}\left([a,b]\right)$, then $z\left(t\right)$ is k times differentiable.

For the fractional case, even if $f\in C^{\infty }$, it may happen that $z\left(t\right)\notin C^1$. This means that the solutions may not behave well, and solutions with singular derivatives are quite common. See [13] for more results on the smoothness properties of solutions to fractional initial value problems.

These smoothness properties (or lack of smoothness) make it difficult for numerical methods to provide fast and accurate solutions for (1), thus making Neural FDEs more difficult to handle compared to Neural ODEs. It should be highlighted that during the learning process, the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$ will always adjust the Neural FDE model to the data, independent of the amount of error obtained in the numerical solution of the FDE.

3. Neural Fractional Differential Equations and the Learned Order of the Derivative - $\alpha$

In Neural FDEs, the goal is to model systems where the order of the derivative, $\alpha$, is a parameter that needs to be learned along with the parameters of a NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$. A key challenge in this approach is the potential variability and non-uniqueness of the learned $\alpha$ values, especially when the ground-truth $\alpha$ used to generate synthetic datasets is known a priori. This variability arises from the highly flexible nature of NNs, which can adjust their parameters to fit the data well, regardless of the specific value of $\alpha$.

3.1. Numerical Experiments

The results obtained in [9] propose the existence of parameters ${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n$ and corresponding parameter vectors ${\theta }_{\mathbf{1}},{\theta }_{\mathbf{2}},\cdots ,{\theta }_{\mathbf{n}}$, that satisfy the following condition:

$$\mathcal{L}({\theta }_i,{\alpha }_i)\approx 0\mathrm{for}i=1,2,\cdots ,n.$$

(11)

In theory, this implies that for each i, the loss function $\mathcal{L}$ can converge to zero.

There are some results in the literature on the universal approximation capabilities of Neural ODEs and ResNets. Interested readers should consult [14,15] for more information.

Based on these observations, we now conduct numerical experiments to observe in detail the practical outcomes. Specifically, we employ a Neural FDE to model population dynamics with different ground-truth $\alpha$ values. Our goal is to examine how different initialisations of $\alpha$ impact the final learned $\alpha$ values.

We analyse the $\alpha$ values learned across multiple runs and observe how $\alpha$ evolves during training. Additionally, we fix the $\alpha$ value and allow the Neural FDE to learn only the right-hand side of the FDE, comparing the loss for different $\alpha$ values.

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

$${ }_0^C{D}_t^{\alpha }=rP\left(t\right)\left(1-{\displaystyle \frac{P\left(t\right)}{K}}\right),$$

(12)

with initial condition $P\left(0\right)=P\left(t_0\right)=100$, a growth-rate $r=0.1$ and carrying capacity $K=1000$.

To conduct the experiments, three training datasets were generated by numerically solving (12) over the interval $t\in (0,10)$ with a step size of $0.1$ for three different values of $\alpha$, namely $\alpha =0.3,0.5,0.99$. These datasets are denoted as $P_{\alpha =0.3}$, $P_{\alpha =0.5}$, and $P_{\alpha =0.99}$, respectively. For the experiments, Adam optimiser [16] was used with a starting learning rate of 1E-03.

3.1.1. Case Study 1: varying $\alpha$ and fixed number of iterations

The generated datasets were used to train the Neural FDE model (see Equation (1)) with four different initialisations of the derivative order: $\alpha =0.1$, $0.3$, $0.5$, and $0.99$. For each initialisation, three independent runs were conducted considering 200 iterations.

The results obtained for each dataset, $P_{\alpha =0.3}$, $P_{\alpha =0.5}$, and $P_{\alpha =0.99}$, using the various $\alpha$ initialisation values are presented in Table 1, Table 2 and Table 3 and the evolution of the training losses and $\alpha$ values can be visually observed in Figure 3, Figure A1-Figure A2 and Figure 4, Figure A3, Figure A4, respectively.

The results presented in Table 1 through Table 3 indicate that the initial value of $\alpha$ does not significantly impact the final learned value of $\alpha$. While there is some effect, it is challenging to discern a clear pattern, and the observed variations fall within a certain range of randomness. However, the ground-truth value of $\alpha$ in the datasets influences the learnt value, with higher ground-truth values leading to higher learnt values. For example, in Table 3, all learnt $\alpha$ values are above $0.50$, whereas for other datasets, the learnt $\alpha$ values are lower. This is expected since the Neural FDE model should fit the given data (that depends on $\alpha$).

Furthermore, the results demonstrate the approximation capabilities of NNs, which allow them to model the data dynamics with low errors even when the learnt $\alpha$ is far from the ground-truth value. For example, in Table 1, the ground-truth $\alpha$ is $0.3$, but the lowest training loss was achieved with a learnt $\alpha$ of $0.4497$. Similarly, in Table 3, the ground-truth $\alpha$ is $0.99$, but the lowest training loss values were achieved with learnt $\alpha$ of $0.518$ and $0.5876$. It is interesting to note that the values of $\alpha$ are significantly different from 1. The results obtained suggest that the optimisation procedure may be placing greater importance on the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$.

The evolution of the loss and the learned $\alpha$ values along training (Figure 3 and Figure 4, respectively, for the case $P_{\alpha =0.3}$), indicate that the $\alpha$ initialisation and the ground-truth $\alpha$ do not significantly influence the Neural FDE training process. In Figure 4 a pattern is observed. Initially the $\alpha$ values increase drastically, and then slowly decrease approaching a plateau. In some cases, after around 100 iterations, the values begin to increase again. This behaviour results from the complex interplay between the two different NNs. The behaviour is even more irregular for the cases $P_{\alpha =0.5},P_{\alpha =0.99}$, as shown in Appendix A.

3.1.2. Case Study 2: fixed $\alpha$ and fixed number of iterations

To demonstrate the effectiveness of Neural FDEs in minimising the loss function $\mathcal{L}(\theta ,\alpha )$ for various $\alpha$ values, we performed experiments with fixed $\alpha$ values of $0.1$, $0.3$, $0.5$ and $0.99$. In these experiments, the Neural FDEs used exclusively the NN ${\mathbf{f}}_{\theta }$, to fit the solution $\mathbf{h}\left(t\right)$ to the three distinct datasets $P_{\alpha =0.3},P_{\alpha =0.5},P_{\alpha =0.99}$. The experimental conditions remained the same, except that $\alpha$ was not learnable. Note that the stopping criterion was set to a fixed number of iterations, specifically 200.

The final training losses obtained for three runs, for each fixed $\alpha$, are summarised in Table 4, Table 5 and Table 6. The evolution of the training losses can be visualised in Figure 5, Figure 6 and Figure 7.

The results in Table 4, Table 5 and Table 6 show that the final training loss values are generally similar for different values of $\alpha$. Even when the fixed $\alpha$ matches the ground-truth $\alpha$ of the dataset, in general, the final loss values are comparable to other fixed $\alpha$ values.

The training losses evolution in Figure 5, Figure 6 and Figure 7 show similar behaviour for the different fixed $\alpha$ values.

3.1.3. Case Study 3: fixed $\alpha$ and fixed loss threshold

As a final experiment, we aimed to numerically demonstrate that Neural FDEs are capable of achieving the same fit to the data independent of the $\alpha$ order. To unequivocally show this, we modified the stopping criteria of our Neural FDE training from a maximum number of iterations to a loss threshold of 1E-05. This approach ensures that the Neural FDEs are trained until the threshold is achieved, allowing us to demonstrate that they can reach the same loss values regardless of the $\alpha$ value.

In this experiment, we performed one run for each dataset $P_{\alpha =0.3}$, $P_{\alpha =0.5}$, and $P_{\alpha =0.99}$ with fixed $\alpha$ values of 0.1, 0.3, 0.5, and 0.99. The results are organised in Table 7, Table 8 and Table 9

The results presented in Table 7 through Table 9 show that the NN ${\mathbf{f}}_{\theta }$ is capable of optimising the parameters $\theta$ to accurately model the Neural FDE to the data dynamics, regardless of whether the imposed value of $\alpha$ is close to or far from the ground truth.

These results might lead some to believe that changing $\alpha$ does not impact memory. They might also suggest that the Neural FDE can model systems with varying memory levels, independent of $\alpha$. Although, this only happens because the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$, which models the FDE’s right-hand side, adjusts its parameters $\theta$ to fit the FDE for any $\alpha$ value (effectively mimicking the memory effects associated with $\alpha$). Therefore, Neural FDEs do not need a unique $\alpha$ value for each dataset. Instead, they can work with an infinite number of $\alpha$ values to fit the data. This flexibility is beneficial when fitting to given data is required, and the underlying physics is not known.

4. Conclusions

In this work, we present the theory of fractional initial value problems and explore its connection with Neural Fractional Differential Equations (Neural FDEs). We analyse both theoretically and numerically how the solution of a Neural FDE is influenced by two factors: the NN ${\mathbf{f}}_{\theta }$, which models the right-hand side of the FDE, and the NN that learns the order of the derivative.

We also investigate the numerical evolution of the order of the derivative ($\alpha$) and training loss across several iterations, considering different initialisations of the $\alpha$ value. For this experiment, with a fixed number of 200 iterations, we create three synthetic datasets for a fractional initial value problem with ground-truth $\alpha$ values of 0.3, 0.5, and 0.99. We test four different initialisations for $\alpha$: 0.1, 0.3, 0.5, and 0.99. The results indicate that both the initial $\alpha$ value and the ground-truth $\alpha$ have minimal impact on the Neural FDE training process. Initially, the $\alpha$ values increase sharply and then slowly decrease towards a plateau. In some cases, around 100 iterations, the values begin to rise again. This behaviour results from the complex interaction between the two NNs, and it is particularly irregular for $\alpha =0.5$ and $\alpha =0.99$. The loss values achieved are low across all cases.

We then repeated the experiments, fixing the $\alpha$ value, meaning there was no initialisation of $\alpha$, and the only parameters changing in the minimisation problem were those of the NN ${\mathbf{f}}_{\theta }(t,\mathbf{h}\left(t\right))$. The results confirm that the final training loss values are generally similar across different fixed values of $\alpha$. Even when the fixed $\alpha$ matches the ground-truth $\alpha$ of the dataset, the final loss remains comparable to other fixed $\alpha$ values.

In a final experiment, we modified the stopping criteria of our Neural FDE training from a maximum number of iterations to a loss threshold. This ensures that the Neural FDEs are trained until the loss threshold is achieved, demonstrating that they can reach similar loss values regardless of the $\alpha$ value. We conclude that, ${\mathbf{f}}_{\theta }$ struggles to adjust its parameters $\theta$ to fit the FDE to the data for any given derivative order. Consequently, Neural FDEs do not require a unique $\alpha$ value for each dataset. Instead, they can use a wide range of $\alpha$ values to fit the data, suggesting that ${\mathbf{f}}_{\theta }$ is a universal approximator.

If we train the model using data points obtained from an unknown experiment, then the flexibility of the Neural FDE proves to be an effective method for obtaining intermediate information about the system, provided the dataset contains sufficient information. If the physics involved in the given data is known, it is recommended to incorporate this knowledge into the loss function. This additional information helps to improve the extrapolation of results.