Unification of popular artificial neural network activation functions
We present a unified representation of the most popular neural network
activation functions. Adopting Mittag-Leffler functions of fractional calculus,
we propose a flexible and compact functional form that is able to interpolate
between various activation functions and mitigate common problems in training
neural networks such as vanishing and exploding gradients. The presented gated
representation extends the scope of fixed-shape activation functions to their
adaptive counterparts whose shape can be learnt from the training data. The
derivatives of the proposed functional form can also be expressed in terms of
Mittag-Leffler functions making it a suitable candidate for gradient-based
backpropagation algorithms. By training multiple neural networks of different
complexities on various datasets with different sizes, we demonstrate that
adopting a unified gated representation of activation functions offers a
promising and affordable alternative to individual built-in implementations of
activation functions in conventional machine learning frameworks.
Mohammad Mostafanejad
# Unification of popular artificial neural network activation functions

###### Abstract

We present a unified representation of the most popular neural network activation functions. Adopting Mittag-Leffler functions of fractional calculus, we propose a flexible and compact functional form that is able to interpolate between various activation functions and mitigate common problems in training neural networks such as vanishing and exploding gradients. The presented gated representation extends the scope of fixed-shape activation functions to their adaptive counterparts whose shape can be learnt from the training data. The derivatives of the proposed functional form can also be expressed in terms of Mittag-Leffler functions making it a suitable candidate for gradient-based backpropagation algorithms. By training LeNet-5 neural network on MNIST and CIFAR-10 datasets, we demonstrate that adopting a unified gated representation of activation functions offers a promising and affordable alternative to individual built-in implementations of activation functions in conventional machine learning frameworks.

## I Introduction

Activation functions are one of the key building blocks in artificial neural networks (ANNs) that control the richness of the neural response and determine the accuracy, efficiency and performance[1] of multilayer neural networks as universal approximators.[2] Due to their biological links[3, 4] and optimization performance, saturating activation functions[5] such as logistic sigmoid and hyperbolic tangent[6] were commonly adopted in early neural networks. Nevertheless, both activation functions suffered from the vanishing gradient problem.[7] Later studies on image classification using restricted Boltzmann machines[8] and deep neural networks[9] demonstrated that rectified linear units (ReLUs) can mitigate the vanishing gradient problem and improve the performance of neural networks. Furthermore, the sparse coding produced by ReLUs not only creates a more robust and disentangled feature representation but also accelerates the learning process.[9]

The computational benefits and the current popularity of ReLUs should be taken with a grain of salt due to their notable disadvantages such as bias shift,[1] ill-conditioned parameter scaling[9] and dying ReLU.[10] Furthermore, the unbounded nature of ReLUs for positive inputs, while potentially helpful for training deep neural networks, can aggravate the exploding gradient problem in recurrent neural networks.[11, 12] In order to address the dying ReLU and the vanishing/exploding gradient problems, a multitude of ReLU variants has been proposed[13, 14, 15, 16, 17] but none has managed to consistently outperform the vanilla ReLUs in a wide range of experiments.[18] Alternative activation functions such as exponential linear units (ELUs)[1] and scaled exponential linear units (SELUs)[19] have also been proposed to build upon the benefits of ReLU and its variants and provide more robustness and resistance towards the input noise. Yet, among the existing slew of activation functions in the literature,[20, 21, 22] no activation function seems to offer global superiority across all modalities and application domains.

Trainable activation functions,[23, 24, 25, 26] whose functional form is learnt from the training data, offer a more flexible option than their fixed-shape counterparts. In order be able to fine-tune the shape of activation functions during backpropagation,[27] partial derivatives of activation functions with respect to unknown learning parameters are required. It is important to note that some trainable activation functions can also be replaced by simpler multilayer feed-forward subnetworks with
constrained parameters and classical fixed-shape activation functions [21]. The ability to replace a trainable activation function with a simpler sub-neural network highlights a deep connection between the choice of activation functions and performance of neural networks. As such, pre-setting the best possible trainable activation function parameters or fine-tuning the experimental settings [28] such as data preprocessing methods, gradient and weight clipping [12], optimizers [29, 30, 31, 32, 33], regularization methods such as \(L_{1}\), \(L_{2}\) and drop out [34], batch normalization [35], learning rate scheduling [28, 36], (mini-)batch size, or network design [37] variables such as depth (number of layers) and width (number of neurons per layer) of the neural network as well as weight initialization methods [38, 39, 16] becomes an important but challenging task. Several strategies such as neural architecture search [40] and network design space design [41] have been proposed to assist the automation of the network design [37] process but they have to deal with an insurmountable computational cost barrier for practical applications.

In this manuscript, we take a theoretical neuroscientific standpoint [3] towards activation functions by emphasizing the existing connections among them from a mathematical perspective. As such, we resort to the expressive power of rational functions as well as higher transcendental special functions of fractional calculus to propose a unified gated representation of activation functions. The presented functional form is conformant with the outcome of a semi-automated search, performed by Ramachandran _et al._[42], in order to find the optimal functional form of activation functions over a pre-selected set of functions. The unification of activation functions offers several significant benefits: The unified form requires fewer lines of code to be implemented and leads to less confusion in dealing with a wide variety of empirical guidelines on activation functions because individual activation functions correspond to special parameter sets in a single functional representation. The derivatives of the proposed functional form can also be expressed in terms of its constituent special functions, making it a suitable choice for an efficient implementation of backpropagation algorithms for training ANNs. Finally, the proposed functional form can be adopted as a fixed-shape or trainable activation function or both when training neural networks. In other words, one can access different activation functions or interpolate between them by fixing or varying a set of parameters in the gated functional, respectively.

The manuscript is organized as follows: In Sec. II, we introduce Mittag-Leffler functions of one- and two-parameters and discuss their important analytical and numerical properties. Next, we use Mittag-Leffler functions to create a gated representation that can unify a set of most commonly used activation functions. Section III delineates the computational details of our experiments presented in Sec. IV, where we provide numerical evidence for the efficiency and accuracy of the proposed functional form. Concluding remarks and future directions are presented in Sec. V.

## II Theory

### Mittag-Leffler functions of one- and two-parameters

Mittag-Leffler functions, sometimes referred to as "_the queen of functions in fractional calculus_" [43, 44] are one of the most important higher transcendental functions that play a fundamental role in fractional calculus [45, 46, 47]. The interested reader is referred to Refs. [48, 49, 46] for a survey of scientific and engineering applications. The one-parameter Mittag-Leffler function is defined as [50, 49]

\[E_{\alpha}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(\alpha k+1)},\qquad \alpha\in\mathbb{C}, \tag{1}\]
where \(\Gamma(\cdot)\) stands for Euler Gamma function [51, 52] and \(\mathbb{C}\) denotes the set of complex numbers. For all values of \(\mathrm{Re}(\alpha)>0\), the series in Eq. 1 converges everywhere in the complex plane and the one-parameter Mittag-Leffler function becomes an entire function of the complex variable \(z\).[49] However, when \(\mathrm{Re}(\alpha)<0\), the series in Eq. 1 diverges everywhere on \(\mathbb{C}\setminus\{0\}\). As \(\alpha\to 0^{+}\), the Mittag-Leffler function can be expressed as [49]

\[E_{0}(\pm z)=\frac{1}{1\mp z},\qquad|z|<1. \tag{2}\]

Although Mittag-Leffler series has a finite radius of convergence, the restriction in Eq. 2 can be lifted and the asymptotic geometric series form can be adopted as a part of the definition of Mittag-Leffler function for \(\alpha=0\).[53] The aforementioned definition seems to be matching the implementation of Mittag-Leffler function in Mathematica 13.2.[54] Note that for \(x>0\) and \(0\leq\alpha\leq 1\), the one-parameter Mittag-Leffler function with negative arguments, \(E_{\alpha}(-x)\), is a completely monotonic [55] function with no real zeros.[49] The two-parameter Mittag-Leffler function can be similarly defined as

\[E_{\alpha,\beta}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(\alpha k+\beta)}, \qquad\mathrm{where}\quad\mathrm{Re}(\alpha)>0,\quad\mathrm{and}\quad\beta \in\mathbb{C}. \tag{3}\]

The exponential form of Mittag-Leffler functions, \(E_{1}(z)=E_{1,1}(z)\), has no zeros in the complex plane. Nonetheless, for all \(m\in\mathbb{N}\), where \(\mathbb{N}\) is the set of natural numbers, \(E_{1,-m}\) has its only \(m+1\)-order zero located at \(z=0\). All zeros of \(E_{2}(z)\) are simple and can be found on the negative real semi-axis. For a more detailed discussion on the distribution of zeros and the asymptotic properties of Mittag-Leffler functions, see Ref. [49].

Parallel to the study of analytic properties, the realization of accurate and efficient numerical methods for calculating Mittag-Leffler functions is still an open and active area of research.[49, 56] In particular, the existence of free, open-source and accessible software for computing Mittag-Leffler functions is key to their usability in practical applications. We must note that the code base and programmatic details of recent updates to the implementation of Mittag-Leffler functions in Mathematica are not publicly available for further analysis in this manuscript. Nevertheless, several open-source numerical modules for numerical computation of Mittag-Leffler functions are available in the public domain. Gorenflo _et al.[57]_ have proposed an algorithm for computing two-parameter Mittag-Leffler functions that is suitable for use in Mathematica. Podlubny's algorithm is implemented in MATLAB which allows the computation of Mittag-Leffler functions with arbitrary accuracy.[58] Garrappa has proposed an efficient method for calculating one- and two-parameter Mittag-Leffler functions using hyperbolic path integral transform and quadrature.[59] Both MATLAB [60] and Python[61] implementations of Garrappa's algorithm are also available in the public domain. Zeng and Chen have also constructed global Pade approximations for special cases of parameters \(0<\alpha\leq 1\) and \(\beta\geq\alpha\), based on the complete monoticity of \(E_{\alpha,\beta}(-x)\).[62] Another powerful feature of Mittag-Leffler functions is its relation to other higher transcendental special functions such as hypergeometric, Wright, Meijer \(G\) and Fox \(H\)-functions.[49, 51, 52, 63, 64] which allows for more general analytic manipulations and efficient numerical computations. For example, Mathematica automatically simplifies the one-parameter Mittag-Leffler functions with non-negative (half-)integer \(\alpha\) to (sum of) generalized hypergeometric functions.[49, 65] In addition to the algorithm complexity and implementation specifics, the total number of activation functions in a neural network can strongly affect its runtime on computing accelerators such as graphics processing units (GPUs). The neural network architecture[37] is also a major factor in determining the computational cost.[41] We will consider the impact of these factors in our numerical experiments in Sec. IV.
Gated representation of activation functions

In order to unify the most common classical fixed-shape activation functions, listed in a recent survey,[21] we propose the following functional form

\[x\Phi\bigg{[}x\bigg{|}\gamma\begin{array}{c}\alpha_{1}&\beta_{1}&f\\ \alpha_{2}&\beta_{2}&g\end{array}\bigg{]}:=x\bigg{\{}x^{\gamma-1}\left(\frac{E _{\alpha_{1},\beta_{1}}\left[f(x)\right]}{E_{\alpha_{2},\beta_{2}}\left[g(x) \right]}\right)\bigg{\}}, \tag{4}\]

where the gate function, \(\Phi[f(x),g(x)]\), is a rational binary composition of two "well-behaved" functional mappings \(f,g:\mathbb{R}\rightarrow\mathbb{R}\) and is responsible for generating a (non-)linear neural response. Here, \(\mathbb{R}\) denotes the set of real numbers. The gated representation in Eq. 4 incorporates the functional form obtained from an automated search over a set of a pre-selected functions [42] and is consistent with the functional form of popular activation functions such as ReLU and Swish. Throughout this manuscript, we restrict ourselves to \(\gamma\geq 0\), \(\text{Re}(\alpha)>0\) and \(\beta\in\mathbb{R}\).

Table 1 presents a shortlist of popular fixed-shape classical activation functions that are accessible to the proposed gated representation as special cases via different sets of parameters.

The ReLU activation function is commonly represented in a piecewise functional form as \(\max(0,x)\). In order to mimic this behavior, the gate function in Eq. 4 should reduce to identity for \(x>0\) and zero otherwise. The former condition is satisfied when \(\gamma=1\) and \(E_{\alpha_{1},\beta_{1}}\left[f(x)\right]/E_{\alpha_{2},\beta_{2}}\left[g(x) \right]=1\), for which \(\Phi\bigg{[}x\bigg{|}1\begin{smallmatrix}\alpha&\beta&f\\ \alpha&\beta&f\end{smallmatrix}\bigg{]}=1\) is a trivial case. Plots of ReLU activation function and its gated representation are shown in Fig. 1(a). Note that the y-axis label, \(a(x)\), collectively refers to activation functions regardless of their functional form. The gate functional for sigmoid activation function, \(\sigma(x)\), takes \(f(x)=-e^{-x}\) and \(g(x)=0\) to yield

\[x\Phi\bigg{[}x\bigg{|}0\begin{array}{ccc}0&1&-e^{-x}\\ 1&1&0\end{array}\bigg{]}=\frac{1}{1+e^{-x}}=\sigma(x). \tag{5}\]

Plots of sigmoid activation function and its gated representation are shown in Fig. 1(b). The sigmoid gate function in Eq. 5 can be morped into that of Swish by setting \(\gamma=1\) and \(f(x)=e^{-cx}\) to obtain

\[x\Phi\bigg{[}x\bigg{|}1\begin{array}{ccc}0&1&-e^{-cx}\\ 1&1&0\end{array}\bigg{]}=x\,\sigma(cx), \tag{6}\]

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline Activation Function & Argument & \(\gamma\) & \(\alpha_{1}\) & \(\alpha_{2}\) & \(\beta_{1}\) & \(\beta_{2}\) & \(f(x)\) & \(g(x)\) \\ \hline ReLU\({}^{a}\) & \(x\) & 1 & 1 & 1 & 1 & 0 & 0 \\ Sigmoid & \(x\) & 0 & 0 & 1 & 1 & 1 & \(-e^{-x}\) & 0 \\ Swish & \(x\) & 1 & 0 & 1 & 1 & 1 & \(-e^{-cx}\) & 0 \\ Softsign & \(x\) & 1 & 0 & 1 & 1 & 1 & \(-|x|\) & 0 \\ Hyperbolic Tangent & \(x\) & 1 & 2 & 2 & 2 & 1 & \(x^{2}\) & \(x^{2}\) \\ Mish & \(\log(1+e^{x})\) & 2 & 2 & 2 & 2 & 1 & \(x^{2}\) & \(x^{2}\) \\ Bipolar Sigmoid & \(x\) & 0 & 0 & 0 & 1 & \(-e^{-x}\) & \(e^{-x}\) \\  & \(x/2\) & 1 & 2 & 2 & 2 & 1 & \(x^{2}\) & \(x^{2}\) \\ GELU & \(x\) & 1 & 1/2 & 1 & 1 & 1 & \(x/\sqrt{2}\) & \(x^{2}/2\) \\ \hline \hline \end{tabular} \({}^{a}\) The set of parameters are pertinent to \(x\geq 0\), otherwise \(\Phi=0\).

\({}^{b}\) At least two representations exist for the bipolar sigmoid.

\end{table}
Table 1: Special cases of gate function \(\Phi\) in Eq. 4
Figure 1: Plots of built-in and gated representation of various activation functions
where \(c\) is a trainable parameter. For \(c=1\), the resulting activation function in Eq. 6 is referred to as Swish-1.[42] Plots of Swish-1 activation function and its gated variant are shown in Fig. 1(c). Setting \(f(x)=-|x|\) in the gate function, one can convert Swish into Softsign, defined as

\[x\Phi\bigg{[}x\bigg{|}1\begin{array}{ccc}0&1&-|x|\\ 1&1&0\end{array}\bigg{]}=\frac{x}{1+|x|}. \tag{7}\]

Plots of Softsign activation function and its gated representation are illustrated in Fig. 1(d). The gate functional for the hyperbolic tangent activation function takes \(f(x)=g(x)=x^{2}\) to yield

\[x\Phi\bigg{[}x\bigg{|}1\begin{array}{ccc}2&2&x^{2}\\ 2&1&x^{2}\end{array}\bigg{]}=\tanh(x). \tag{8}\]

Plots of hyperbolic tangent activation function and its gated representation are presented in Fig. 1(f). As mentioned in Sec. I, the gated functional form in Eq. 4 arms us with significant variational flexibility. In addition to accessing a set of fixed-shape activation functions via setting the gate function parameters, we can also interpolate between different functional forms by varying those parameters over a finite domain. Figure 2 illustrates an example where by fixing all parameters in the gated representation of hyperbolic tangent except \(\beta_{2}\), one can smoothly interpolate between linear (\(\beta_{2}=2\)) and hyperbolic tangent (\(\beta_{2}=1\)) activation functions. Thus, it is possible to tune the saturation behavior of gated representation of saturating functions such as hyperbolic tangent and mitigate their vanishing/exploding gradient problem in a controlled fashion.[8; 9; 38] Furthermore, one can turn \(\beta_{2}\) (or in principle, any other parameter) into a trainable parameter and allow the hosting neural network to learn its optimal value from the training data.

Our unification strategy can go beyond the aforementioned list of fixed-shape or trainable activation functions. For instance, Mish[66] can be obtained by passing Softplus, \(\log(1+e^{x})\), to hyperbolic tangent gate function as an argument and setting \(\gamma=2\) to get

\[x\Phi\bigg{[}\log(1+e^{x})\bigg{|}2\begin{array}{ccc}2&2&x^{2}\\ 2&1&x^{2}\end{array}\bigg{]}=x\,\tanh\big{[}\log(1+e^{x})\big{]}. \tag{9}\]
Plots of Mish and its gated representation are illustrated in Fig. 1(e). The bipolar sigmoid function can also be expressed by using the hyperbolic tangent gate function and passing a scaled linear function as an argument to get

\[x\Phi\bigg{[}\frac{x}{2}\bigg{|}1\begin{array}{ccc}2&2&x^{2}\\ 2&1&x^{2}\end{array}\bigg{]}=\tanh\left(\frac{x}{2}\right). \tag{10}\]

Equivalently, one can also express the bipolar sigmoid function with a different set of paramenters and arguments in the gate function as

\[x\Phi\bigg{[}x\bigg{|}0\begin{array}{ccc}0&1&-e^{-x}\\ 0&1&e^{-x}\end{array}\bigg{]}=\frac{1-e^{-x}}{1+e^{-x}}. \tag{11}\]

Setting \(f(x)=\frac{x}{\sqrt{2}}\) and \(g(x)=\frac{x^{2}}{2}\), Gaussian error linear units (GELUs) can also be written as

\[\frac{x}{2}\Phi\bigg{[}x\bigg{|}1\begin{array}{ccc}\frac{1}{2}&1&\frac{x}{ \sqrt{2}}\\ 1&1&\frac{x^{2}}{2}\end{array}\bigg{]}=\frac{x}{2}\bigg{[}1+\mathrm{erf}\left( \frac{x}{\sqrt{2}}\right)\bigg{]}, \tag{12}\]

where \(\mathrm{erf}(\cdot)\) is the error function.[51, 52]

## III Computational details

In order to illustrate the efficiency and performance of the neural networks armed with gated activation functions, we train the classical LeNet-5 (Fig. 3) neural network[67] on Modified National Institute of Standards and Technology (MNIST)[68] and CIFAR-10[69] image classification datasets. The MNIST dataset consists of 60,000 training and 10,000 testing grayscale images of hand-written digits \((0,1,2,\ldots,9)\) that are normalized and centered to a 28\(\times\)28 fixed size. CIFAR-10 dataset contains 50,000 training and 10,000 test images from 10 object classes (airplane, automobile, bird, cat, deer, dog, frog, horse, ship, truck). Each data point in CIFAR-10 dataset is a 32\(\times\)32 RGB image. Both datasets are pulled from Wolfram Data Repository[70, 71] and all training experiments are performed using Wolfram Mathematica 13.2.[54] Individual training sessions (excluding that of the baseline) involve replacing all three element-wise ReLU activation layers in LeNet-5 architecture with their (gated) counterparts from Table I. Bipolar sigmoid and GELU are excluded from our study as their hosting neural networks fail to converge without deviation from the selected default settings in Mathematica and further modifications to avoid the divergence. The adaptive momentum (ADAM) stochastic gradient descent optimizer is used for training with the stability parameter, the first and the second moment exponential decay rates set to \(\epsilon=10^{-5}\), \(\beta_{1}=0.9\) and \(\beta_{2}=0.999\), respectively. All network parameters are randomly initialized using "1234" as seed in order to ensure each individual training session in the ensemble starts with the same set of initial parameters. We ran each experiment corresponding to MNIST and CIFAR-10 datasets for 10 and 20 epochs, respectively with a fixed batch size of 64. Each training session is repeated twenty times in single-precision (32-bit floating-point numbers) and the averaged results are rounded to four significant digits as shown in Tables II and III.

Two computational frameworks have been deployed in this study: A single personal laptop armed with a NVIDIA GeForce GTX 1650 GPU and a data center unit equipped with an NVIDIA A100 80GB PCIe GPU. Results from the former setup can be found in the Supporting Information. The trajectory of all training instances and a simple Mathematica code snippet for reproducing the results can be found in Ref. [72].
## IV Results and discussion

We analyze the performance of LeNet-5 neural network architecture by replacing the element-wise layers of ReLU activation functions with their counterparts from Table 1 as well as their gated representations. In order to quantify the effect of various activation functions on the performance of LeNet-5 classifier, we mainly focus on two performance metrics: Loss and accuracy. The training loss is measured using multi-class cross-entropy, which is defined as[73, 74]

\[\mathscr{L}=-\sum_{i=1}^{N}\sum_{k=1}^{K}y_{i,k}\ln(\hat{y}_{i,k}). \tag{13}\]

For each data point \(i\) in a dataset of size \(N\), cross-entropy can measure how well the estimated probabilities, \(\hat{y}_{i,k}\) for each class \(k\), where \(k\in\{1,2,3,\ldots,K\}\), match those of the target class labels, \(y_{i,k}\). Compared with other loss functions, cross-entropy can also improve the convergence rate of the optimization process in our case by more aggressively penalizing the incorrect predictions and generating larger gradients.[74] Accuracy is defined as fraction of number of times that the classifier is correct in its predictions.[73] Since accuracy is not an appropriate metric for imbalanced datasets,[74] we also include the macro-average values of precision, recall, and F1 score in the Supporting Information despite our knowledge about the balanced nature of class distributions in MNIST and CIFAR-10 datasets. Table 2 shows the training and validation performance results where those of LeNet-5 neural network with ReLU element-wise layers are taken as baseline. For each activation function, there are two entries: the first entry refers to the results of built-in activation functions in Mathematica 13.2 and the second one corresponds to those of gated activation functions.

Table 2 reveals that the ensemble average validation accuracy of LeNet-5 classifier on MNIST test set is not significantly sensitive towards the choice of activation functions in the element-wise layers. In particular, the validation accuracy of LeNet-5 classifier armed with sigmoid activation functions is slightly smaller than that of the baseline neural network with ReLUs. Furthermore, choosing other activation functions such as Softsign and Mish further deteriorates the

Figure 3: LeNet-5 neural network architecture
corresponding average validation accuracies compared with that of ReLUs in the baseline LeNet-5 architecture. Plots of training/validation loss and accuracy versus epochs can be found in Supporting Information.

Our main interest in Table 2 is in the average values of total wall-clock time spent on training LeNet-5 on MNIST data set using a NVIDIA A100 80GB PCIe GPU. The average timings reveal that the added cost of calculating one- or two-parameter Mittag-Leffler functions in the gated representation of activation functions is small compared with their built-in variants implemented in Mathematica 13.2 program package.[75] Specifically, the largest measured time gap is observed between built-in and gated representations of Mish which mainly stems from the overhead of calculating Softplus and passing it as an argument to two-parameter Mittag-Leffler functions for each neural response in the element-wise activation layers. On the other hand, the computational time gap between built-in and gated representation of Softsign is very small.

Table 3 presents the performance results for training/validation of LeNet-5 network on CIFAR-10 dataset. All results correspond to the average of 20 individual training sessions, each running for 20 epochs. A comparison of the average accuracy values in Table 3 with those in Table 2 reflects the more intricate nature of the CIFAR-10 dataset which requires deeper neural networks or more advanced architectural design and training strategies. The interested reader is referred to a study on "Convolutional Deep Belief Neural Networks on CIFAR-10" [76] and ImageNet competition for a chronological survey of efforts on this topic.[77] Table 3 reveals that the validation accuracy of LeNet-5 neural network can be improved by replacing ReLUs with any other activation function considered in this study with the exception of sigmoid and hyperbolic tangent. In particular, replacing ReLUs with Swish-1 or Mish yields the largest improvements in validation accuracy of \(\approx\) 1.2 %.

The average timings for training LeNet-5 neural network on CIFAR-10 dataset show similar trends to those presented in Table 2. Specifically, the time difference between computations pertinent to built-in and gated representation of Mish activation function peaks at about 30 seconds. Note that the timings in Table 3 are roughly two times larger in magnitude than their counter

\begin{table}
\begin{tabular}{l c c c} \hline \hline Activation Function\({}^{c}\) & Loss & Accuracy (\%) & Wall-Clock Time (s) \\ \hline Sigmoid & 0.0270 (0.0290) & 99.13 (99.13) & 21.93 \\  & 0.0270 (0.0290) & 99.13 (99.13) & 30.39 \\ Swish-1 & 0.0411 (0.0295) & 98.73 (98.98) & 28.38 \\  & 0.0422 (0.0295) & 98.69 (98.99) & 30.60 \\ Softsign & 0.0138 (0.0324) & 99.56 (99.00) & 28.39 \\  & 0.0137 (0.0315) & 99.56 (99.03) & 28.61 \\ tanh & 0.0146 (0.0329) & 99.52 (98.95) & 20.77 \\  & 0.0126 (0.0325) & 99.58 (99.01) & 29.72 \\ Mish & 0.0279 (0.0299) & 99.14 (99.07) & 24.38 \\  & 0.0339 (0.0304) & 98.96 (99.03) & 40.73 \\ \hline ReLU\({}^{d}\) & 0.0144 (0.0292) & 99.54 (99.16) & 23.78 \\ \hline \hline \end{tabular}

* All results are ensemble averages over 20 independent training and testing experiments.
* All calculations are performed using a NVIDIA A100 80GB PCIe GPU.
* The test results are given in parentheses. The first and second rows in each activation function entry correspond to the built-in and gated representations, respectively.
* The LeNet-5 neural network architecture with ReLU activation functions is taken as the base architecture.

\end{table}
Table 2: The best performance metrics and timings pertinent to training and testing of LeNet neural network on MNIST dataset with various activation functions\({}^{a,b}\)
parts in Table II due to the adopted number of training epochs (20 for CIFAR-10 compared with 10 for MNIST). The average timings reported in both tables suggest that a unified implementation of the most popular activation functions using Mittag-Leffler functions is possible only at a small and affordable additional computational cost compared with their built-in individual implementations. The gap between the built-in implementations of activation functions and their gated representations can be further reduced as more efficient algorithms and implementations of special functions such as Mittag-Leffler function become available. Comparing the aforementioned timings in II and III obtained using a NVIDIA A100 80GB PCIe GPU with those in Supporting Information computed by a NVIDIA GeForce GTX 1650 GPU demonstrates that the availability of more powerful computing accelerators is yet another asset to benefit from when training ANNs with a large number of gated activation functions.

## V Conclusion and future work

In this manuscript, we have presented a unified representation of some of the most popular neural network activation functions. The proposed functional form not only sheds light on direct analytical connections between several well-established activation functions in the literature, but also allows for interpolating between different functional forms through varying the gate function parameters. Furthermore, the gated functional form and its derivatives can both be expressed in terms of Mittag-Leffler functions which makes them a suitable candidate for training neural networks using backpropagation algorithms. A unified representation of activation functions can also lead to large savings in programming efforts compared with what is otherwise required for individual implementations of activation functions in popular machine learning frameworks via inheritance and/or customized classes. Through training the classic LeNet-5 neural network on standard benchmark datasets such as MNIST and CIFAR-10, we have established the possibility of implementing an efficient unified representation of activation functions without sacrificing the accu

\begin{table}
\begin{tabular}{l c c c} \hline \hline Activation Function\({}^{c}\) & Loss & Accuracy (\%) & Wall-Clock Time (s) \\ \hline \multirow{2}{*}{Sigmoid} & 0.7353 (1.0151) & 74.79 (65.39) & 44.26 \\  & 0.7353 (1.0152) & 74.78 (65.38) & 54.97 \\ \multirow{2}{*}{Swish-1} & 0.7323 (0.9016) & 74.55 (69.39) & 54.44 \\  & 0.7322 (0.9018) & 74.55 (69.40) & 60.77 \\ \multirow{2}{*}{Softsign} & 0.6380 (0.9297) & 77.89 (68.64) & 53.09 \\  & 0.6280 (0.9313) & 78.26 (68.56) & 54.52 \\ \multirow{2}{*}{tanh} & 0.8129 (0.9766) & 71.63 (66.79) & 45.58 \\  & 0.8015 (0.9738) & 72.05 (66.90) & 56.09 \\ \multirow{2}{*}{Mish} & 0.6842 (0.8914) & 76.18 (69.40) & 48.37 \\  & 0.6841 (0.8913) & 76.18 (69.44) & 78.09 \\ \hline ReLU\({}^{d}\) & 0.7299 (0.9279) & 74.57 (68.24) & 45.08 \\ \hline \hline \end{tabular}

* All results are ensemble averages over 20 independent training and testing experiments.
* All calculations are performed using a NVIDIA A100 80GB PCIe GPU.
* The test results are given in parentheses. The first and second rows in each activation function entry correspond to the built-in and gated representations, respectively.
* The LeNet-5 neural network architecture with ReLU activation functions is taken as the base architecture.

\end{table}
Table 3: The best performance metrics and timings pertinent to training and testing of the LeNet-5 neural network on CIFAR-10 dataset with various activation functions\({}^{a,b}\)
racy. The analytic properties of one- and two-parameter Mittag-Leffler functions and their relation to other generalized and special functions [52] such as hypergeometric and Wright functions [44] also opens a door to a pristine area of research in fractional ANNs [78] and backpropagation algorithms [79] which is currently under investigation by us.

## 1 General formula for derivatives of Mittag-Leffler function

The differentials of the one- and two-parameter Mittag-Leffler functions can be expressed in terms of Mittag-Leffler function itself. The aforementioned closeness property is computationally beneficial for an efficient implementation of the gradient descent-based backpropagation algorithms for training ANNs. For a more in-depth discussion on differential and recurrence relations of Mittag-Leffler functions of one-, two- and three-parameter(s), see Refs. [49] and [80].

Let \(p\in\mathbb{N}\), where \(\mathbb{N}\) denotes the set of natural numbers. Then, the general derivatives of one-parameter Mittag-Leffler function can be given as

\[\begin{split}\frac{d^{p}}{dz^{p}}E_{p}(z^{p})&=E_{ p}(z^{p}),\qquad\text{and}\\ \frac{d^{p}}{dz^{p}}E_{p/q}(z^{p/q})&=E_{p/q}(z^{p/ q})+\sum_{k=1}^{q-1}\frac{z^{-kp/q}}{\Gamma(1-kp/q)}\qquad q=2,3,\dots.\end{split} \tag{14}\]

Assuming \(\alpha>0\) and \(\beta\in\mathbb{R}\), the first-derivative of the two-parameter Mittag-Leffler functions can be written as a sum of two instances of two-parameter Mittag-Leffler functions as [80]

\[\frac{d}{dz}E_{\alpha,\beta}(z)=\frac{E_{\alpha,\alpha+\beta-1}(z)+(1-\beta)E _{\alpha,\alpha+\beta}(z)}{\alpha}. \tag{15}\]

In general, one can write

\[\frac{d^{m}}{dz^{m}}E_{\alpha,\beta}(z)=\frac{1}{\alpha^{m}}\sum_{k=0}^{m}c_{ k}^{(m)}E_{\alpha,\alpha m+\beta-k}(z),\qquad m\in\mathbb{N}, \tag{16}\]

where \(c_{0}^{(0)}=1\) and the remaining coefficients for \(k=1,2,\dots\) can be computed using the following recurrence relation

\[c_{k}^{(m)}=\begin{cases}\left[1-\beta-\alpha(m-1)\right]c_{0}^{(m-1)},&k=0, \\ c_{k-1}^{(m-1)}+\left[1-\beta-\alpha(m-1)+k\right]c_{k}^{(m-1)},&1\leq k\leq m -1,\\ 1,&k=m.\end{cases} \tag{17}\]

## Acknowledgements

The author would like to thank Dr. Reza Hemmati for proofreading the manuscript, NVIDIA Corporation for the generous Academic Hardware Grant and Virginia Tech for providing an institutional license to Mathematica 13.2.