1. Introduction
Deep learning has emerged as a potent area of research and has enabled a remarkable progress in recent years spanning domains like imaging science [26, 3, 50, 30], biomedical applications [33, 13, 25], satellite imagery, remote sensing [47, 51, 10]
, etc. However, the mathematical foundations of many machine learning architectures are largely lacking
[20, 39, 49, 41, 18]. The current trend of success is largely due to the empirical evidence. Due to the lack of mathematical foundation, it becomes challenging to understand the detailed workings of networks [22, 35].The overarching goal of machine learning algorithms is to learn a function using some known data. Deep Neural Networks (DNN), like Residual Neural Networks (RNN), are a popular family of deep learning architectures which have turned out to be groundbreaking in imaging science. An introductory example of RNN is the ResNet [26] which has been successful for classification problems in imaging science. Compared to the classical DNNs, the innovation of the RNN architecture comes from a simple addition of an identity map between each layer of the network. This ensures a continued flow of information from one layer to another. Despite their success, DNNs are prone to various challenges such as vanishing gradients [8, 20, 48], difficulty in approximating nonsmooth functions, long training time [12], etc.
We remark that recently in [27] the authors have introduced a DenseNet, which is a new approach to prevent the gradient “wash out” by considering dense blocks, in which each layer takes into account all the previous layers (or the memory). They proceed by concatenating the outputs of each dense block which is then fed as an input to the next dense block. Clearly as the number of layers grow, it can become prohibitively expensive for information to propagate through the network. DenseNet can potentially overcome the vanishing gradient issue, but it is only an adhoc method [27, 52]. Some other networks that have attempted to induce multilayer connections are Highway Net [45], AdaNet [16], ResNetPlus [14], etc. All these models, however, largely lack rigorous mathematical frameworks. Furthermore, rigorous approaches to learn nonsmooth functions such as the absolute value function are scarce [28].
There has been a recent push in the scientific community to develop rigorous mathematical models and understanding of the DNNs [18]. One way of doing so is to look at their architecture as dynamical systems. The articles [24, 34, 41, 44, 9] have established that a DNN can be regarded as an optimization problem subject to a discrete ordinary differential equation (ODE) as constraints. The limiting problem in the continuous setting is an ODE constrained optimization problem [41, 44]. Notice that designing the solution algorithms at the continuous level can lead to architecture independence, i.e., the number of iterations remains the same even if the number of layers is increased.
The purpose of this paper is to present a novel fractional deep neural network which allows the network to access historic information of input and gradients across all subsequent layers. This is facilitated via our proposed use of fractional derivative based ODE as constraints. We derive the optimality conditions for this network using the Lagrangian approach. Next, we consider a discretization for this fractional ODE and the resulting DNN is called FractionalDNN. We provide the algorithm and show numerical examples on some standard datasets.
Owing to the fact that fractional time derivatives allow memory effects, in the FractionalDNN all the layers are connected to one another, with an appropriate scaling. In addition, fractional time derivatives can be applied to nonsmooth functions [4]. Thus, we aim to keep the benefits of standard DNN and the ideology of DenseNet, but remove the bottlenecks.
The learning rate in a neural network is an important hyperparameter which influences training [7]. In our numerical experiments, we have observed an improvement in the learning rate via FractionalDNN, which enhances the training capability of the network. Our numerical examples illustrate that, FractionalDNN can potentially solve the vanishing gradient issue (due to memory), and handle nonsmooth data.
The paper is organized as follows. In section 2 we introduce notations and definitions. We introduce our proposed FractionalDNN in section 3. This is followed by section 4 where we discuss its numerical approximation. In section 5, we state our algorithm. The numerical examples given in section 6 show the working and improvements due to the proposed ideas on three different datasets.
2. Preliminaries
The purpose of this section is to introduce some notations and definitions that we will use throughout the paper. We begin with Table 1 where we state the standard notations. In subsection 2.1
we describe the wellknown softmax loss function.
Subsection 2.2 is dedicated to the Caputo fractional time derivative.Symbol  Description 

Number of distinct samples  
Number of sample features  
Number of classes  
Number of network layers (i.e. network depth)  
is the collective feature set of samples.  
are the true class labels of the input data  
Weights  
Linear operator (distinct for each layer)  
Bias (distinct for each layer)  
Lagrange multiplier  
A vector of ones 

Time steplength  
Activation function, acting pointwise  
Order of fractional time derivative  
Derivative w.r.t. the argument  
Trace operator  
Matrix transpose  
Pointwise multiplication  
Max count for randomly selecting a minibatch in training  
Max iteration count for gradientbased optimization solver  
Percentage of training and testing data correctly identified 
2.1. Cross Entropy with Softmax Function
Given collective feature matrix with true labels and the unknown weights , the cross entropy loss function given by
(1) 
measures the discrepancy between the true labels and the predicted labels . Here,
(2) 
is the softmax classifier function, which gives normalized probabilities of samples belonging to the classes.
2.2. Caputo Fractional Derivative
In this section, we define the notion of Caputo fractional derivative and refer [4] and references therein for the following definitions.
Definition 2.1 (Left Caputo Fractional Derivative).
For a fixed real number , and an absolutely continuous function , the left Caputo fractional derivative is defined by:
(3) 
where is the EulerGamma function.
Definition 2.2 (Right Caputo Fractional Derivative).
For a fixed real number , and an absolutely continuous function , the right Caputo fractional derivative is defined by:
(4) 
Notice that, and in definitions Eq. 3 and Eq. 4 exist almost everywhere on , [32, Theorem 2.1], and are represented, respectively, by
Moreover, if and , then one can show that . We note that the fractional derivatives in Eq. 3 and Eq. 4 are nonlocal operators. Indeed, the derivative of at a point depends on all the past and future events, respectively. This behavior is different than the classical case of .
The left and right Caputo fractional derivatives are linked by the fractional integration by parts formula, [5, Lemma 3], which will be stated next. For , let
Lemma 2.3 (Fractional IntegrationbyParts).
For and , the following integrationbyparts formula holds:
(5) 
where and are the left and right RiemannLiouville fractional integrals of order and are given by
3. Continuous Fractional Deep Neural Network
After the above preparations, in this section, we shall introduce the FractionalDNN. First we briefly describe the classical RNN, and then extend it to develop the FractionalDNN. We formulate our problem as a constrained optimization problem. Subsequently, we shall use the Lagrangian approach to derive the optimality conditions.
3.1. Classical RNN
Our goal is to approximate a map . A classical RNN helps approximate , for a known set of inputs and outputs. To construct an RNN, for each layer
, we first consider a lineartransformation of
as,where the pair denotes an unknown linear operator and bias at the layer. When then the network is considered “deep”. Next we introduce non linearity using a nonlinear activation function
(e.g. ReLU or
). The resulting RNN is,(6) 
where is the timestep. Finally, the RNN approximation of is given by,
with as the unknown parameters. In other words, the problem of approximating using classical RNN, intrinsically, is a problem of learning .
Hence, for given datum , the learning problem then reduces to minimizing a loss function , subject to constraint Eq. 6, i.e.,
(7)  
Notice that the system Eq. 6 is the forwardEuler discretization of the following continuous in time ODE, see [26, 23, 41],
(8)  
The continuous learning problem then requires minimizing the loss function at the final time subject to the ODE constraints Eq. 8:
(9)  
s.t.  Eq. 8 
Notice that designing algorithms for the continuous in time problem Eq. 9 instead of the discrete in time problem Eq. 7 has several key advantages. In particular, it will lead to algorithms which are independent of the neural network architecture, i.e., independent of the number of layers. In addition, the approach of Eq. 9 can help us determine the stability of the neural network Eq. 7, see [9, 24]. Moreover, for the neural network Eq. 7, it has been noted that as the information about the input or gradient passes through many layers, it can vanish and “wash out”, or grow and “explode” exponentially [8]. There have been adhoc attempts to address these concerns, see for instance [45, 16, 27], but a satisfactory mathematical explanation and model does not currently exist. One of the main goals of this paper is to introduce such a model.
Notice that Eq. 8, and its discrete version Eq. 6, incorporates many algorithmic processes such as linear solvers, preconditioners, nonlinear solvers, optimization solvers, etc. Furthermore, there are wellestablished numerical algorithms that reuse information from previous iterations to accelerate convergence, e.g. the BFGS method [37], Anderson acceleration [1]
, and variance reduction methods
[40]. These methods account for the history , while choosing . Motivated by these observations we introduce versions of Eq. 6 and Eq. 8 that can account for history (or memory) effects in a rigorous mathematical fashion.3.2. Continuous FractionalDNN
The fractional time derivative in Eq. 3 has a distinct ability to allow a memory effect, for instance in materials with hereditary properties [11]. Fractional time derivative can be derived by using the anomalous random walks where the walker experiences delays between jumps [36]. In contrast, the standard time derivative naturally arises in the case of classical random walks. We use the idea of fractional time derivative to enrich the constraint optimization problem Eq. 9, and subsequently Eq. 7, by replacing the standard time derivative by the fractional time derivative of order . Recall that for , we obtain the classical derivative . Our new continuous in time model, the FractionalDNN, is then given by (cf. Eq. 8),
(10)  
where is the Caputo fractional derivative as defined in Eq. 3. The discrete formulation of FractionalDNN will be discussed in the subsequent section.
The main reason for using the Caputo fractional time derivative over its other counterparts such as the Riemann Liouville fractional derivative is the fact that the Caputo derivative of a constant function is zero and one can impose the initial conditions in a classical manner [42]. Note that is a nonlocal operator in a sense that in order to evaluate the fractional derivative of at a point , we need the cumulative information of over the entire subinterval . This is how the FractionalDNN enables connectivity across all antecedent layers (hence the memory effect). As we shall illustrate with the help of a numerical example in section 6, this feature can help overcome the vanishing gradient issue, as the cumulative effect of the gradient of the precedent layers is less likely to be zero.
Remark 3.1 (Caputo Derivative of Nonsmooth Functions).
Owing to Remark 3.1 we can better account for features, , which are nonsmooth, as a result of which the smoothness requirement on the unknown parameters can be weakened. This, in essence, can help with the exploding gradient issue in DNNs.
3.3. Continuous FractionalDNN and Cross Entropy Loss Functional
Supervised learning problems are a broad class of machine learning problems which use labeled data. These problems are further divided into two types, namely regression problems and classification problems. The specific type of the problem dictates the choice of in Eq. 11. Regression problems often occur in physics informed models, e.g. sample reconstruction inverse problems [3, 25]
. On the other hand, classification problems occur, for instance, in computer vision
[43, 15]. In both the cases, a neural network is used to learn the unknown parameters. In the discussion below we shall focus on classification problems, however, the entire discussion directly applies to regression type problems.Recall that the cross entropy loss functional , defined in Eq. 1, measures the discrepancy between the actual and the predicated classes. Replacing, in Eq. 11 by together with a regularization term , we arrive at
(12)  
s.t. 
Note that, in this case, the unknown parameter , where and are, respectively, the linear operator and bias for each layer, and the weights are a featuretoclass map. Furthermore, is a nonlinear activation function and is the given data, with as the true labels of .
To solve Eq. 12, we rewrite this problem as an unconstrained optimization problem via the Lagrangian functional and derive the optimality conditions. Let denote the Lagrange multiplier, then the Lagrangian functional is given by,
where, is the inner product, and is the Frobenius inner product. Using the fractional integrationbyparts from Eq. 5, we obtain
(13)  
Let denote a stationary point, then the first order necessary optimality conditions are given by the following set of state, adjoint and design equations:

[label=()]

Adjoint Equation. Next, the gradient of with respect to at yields the adjoint equation , equivalently,
(15) where denotes the right Caputo fractional derivative Eq. 4 and is the softmax function defined in Eq. 2. Notice that the adjoint variable in Eq. 15, with its terminal condition, is obtained by marching backward in time. As a result, the equation Eq. 15 is called backward propagation.

Design Equations. Finally, equating , , and to zero, respectively, yields the design equations (with
() as the design variables),
(16) for almost every .
In view of (A)(C), we can use a gradient based solver to find a stationary point to Eq. 12.
Remark 3.2.
(Parametric Kernel ). Throughout our discussion, we have assumed to be some unknown linear operator. We remark that a structure could also be prescribed to , parameterized by a stencil . Then, the kernel is , and the design variables now are . Consequently, can be thought of as a differential operator on the feature space, e.g. discrete Laplacian with a five point stencil. It then remains to compute the sensitivity of the Lagrangian functional w.r.t. to get the design equation. Note that this approach can further reduce the number of unknowns. ∎
Notice that so far the entire discussion has been at the continuous level and it has been independent of the number of network layers. Thus, it is expected that if we discretize (in time) the above optimality system, then the resulting gradient based solver is independent of the number of layers. We shall discretize the above optimality system in the next section.
4. Discrete Fractional Deep Neural Network
We shall adopt the optimizethendiscretize approach. Recall that the first order stationarity conditions for the continuous problem Eq. 12 are given in Eq. 14, Eq. 15, and Item 3. In order to discretize this system of equations, we shall first discuss the approximation of Caputo fractional derivative.
4.1. Approximation of Caputo Derivative
There exist various approaches to discretize the fractional Caputo derivative. We will use the scheme [5, 46] to discretize the left and right Caputo fractional derivative and given in Eq. 3 and Eq. 4, respectively.
Consider the following fractional differential equation involving the left Caputo fractional derivative, for ,
(17) 
We begin by discretizing the time interval uniformly with step size ,
Then using the scheme, the discretization of Eq. 17 is given by
(18) 
where coefficients are given by,
(19) 
Next, let us consider the discretization of the fractional differential equation involving the right Caputo fractional operator, for ,
(20) 
Again using scheme we get the following discretization of Eq. 20:
(21) 
The example below illustrates a numerical implementation of the scheme Eq. 18.
Example 4.1.
Consider the linear differential equation
(22) 
Then, the solution to Eq. 22 is given by, see [42, Section 42], also [38, Section 1.2]
(23) 
where , with , is the Mittag Leffler function defined by
Figure 1 depicts the true solution and the numerical solutions using discretization Eq. 18 for the above example with uniform step size and final time, .
∎
4.2. Discrete Optimality Conditions
Next, we shall discretize the optimality conditions given in Eq. 14 – Item 3. Notice that, each timestep corresponds to one layer of the neural network. It is necessary to do one forward propagation (state solve) and one backward propagation (adjoint solve) to derive an expression of the gradient with respect to the design variables.

[label=()]

Discrete Gradient w.r.t. Design Variables. For , the approximation of the gradient Item 3 with respect to the design variables is given by,
(26)
5. FractionalDNN Algorithm
FractionalDNN is a supervised learning architecture, i.e. it comprises of a training phase and a testing phase. During the training phase, labeled data is passed into the network and the unknown parameters are learnt. Those parameters then define the trained FractionalDNN model for that type of data. Next, a testing dataset, which comprises of data previously unseen by the network, is passed to the trained net, and a prediction of classification is obtained. This stage is known as the testing phase. Here the true classification is not shown to the network when a prediction is being made, but can later be used to compare the network efficiency, as we have done in our numerics. The three important components of the algorithmic structure are forward propagation, backward propagation, and gradient update. The forward and backward propagation structures are given in Algorithms 2 and 1. The gradient update is accomplished in the training phase, discussed in subsection 5.1. Lastly, the testing phase of the algorithm is discussed in subsection 5.2.
5.1. Training Phase
The training phase of FractionalDNN is shown in Algorithm 3.
5.2. Testing Phase
The testing phase of FractionalDNN is shown in Algorithm 4.
6. Numerical Experiments
In this section, we present several numerical experiments where we use our proposed FractionalDNN algorithm from section 5 to solve classification problems for two different datasets. We recall that the goal of classification problems, as the name suggests, is to classify objects into predefined class labels. First we prepare a training dataset and alongwith its classification, pass it to the training phase of FractionalDNN (Algorithm 3). This phase yields the optimal set of parameters learned from the training dataset. They are then used to classify new data points from the testing dataset during the testing phase of FactionalDNN (Algorithm 4). We compare the results of our FractionalDNN with the classical RNN Eq. 9.
The rest of this section is organized as follows: First, we discuss some data preprocessing and implementation details. Then we describe the datasets being used, and finally we present the experimental results.
6.1. Implementation Details

[label=()]

Batch Normalization.
During the training phase, we use the batch normalization (BN) technique
[29]. At each iteration we randomly select a minibatch, which comprises of of the training data. We then normalize the minibatch, to have a zero mean and a standard deviation of one, i.e.
(27) where is the mean and is the standard deviation of the minibatch. The normalized minibatch is then used to train the network in that iteration. At the next iteration, a new minibatch is randomly selected. This process is repeated times. Batch normalization prevents gradient blowup, helps speed up the learning and reduces the variation in parameters being learned.
Since the design variables are learnt on training data processed with BN, we also process the testing data with BN, in which case the minibatch is the whole testing data.

Activation Function. For the experiments we have performed, we have used the hyperbolic tangent function as the activation function, for which,

Regularization. In our experiments, we have used the following regularization:
where is the discrete Laplacian, and
Comments
There are no comments yet.