Bayesian optimization for sparse neural networks with trainable
  activation functions
In the literature on deep neural networks, there is considerable interest in
developing activation functions that can enhance neural network performance. In
recent years, there has been renewed scientific interest in proposing
activation functions that can be trained throughout the learning process, as
they appear to improve network performance, especially by reducing overfitting.
In this paper, we propose a trainable activation function whose parameters need
to be estimated. A fully Bayesian model is developed to automatically estimate
from the learning data both the model weights and activation function
parameters. An MCMC-based optimization scheme is developed to build the
inference. The proposed method aims to solve the aforementioned problems and
improve convergence time by using an efficient sampling scheme that guarantees
convergence to the global maximum. The proposed scheme is tested on three
datasets with three different CNNs. Promising results demonstrate the
usefulness of our proposed approach in improving model accuracy due to the
proposed activation function and Bayesian estimation of the parameters.
Mohamed Fakhfakh, Lotfi Chaari
# Bayesian optimization for sparse neural networks

###### Abstract

In the literature on deep neural networks, there is considerable interest in developing activation functions that can enhance neural network performance. In recent years, there has been renewed scientific interest in proposing activation functions that can be trained throughout the learning process, as they appear to improve network performance, especially by reducing overfitting. In this paper, we propose a trainable activation function whose parameters need to be estimated. A fully Bayesian model is developed to automatically estimate from the learning data both the model weights and activation function parameters. An MCMC-based optimization scheme is developed to build the inference. The proposed method aims to solve the aforementioned problems and improve convergence time by using an efficient sampling scheme that guarantees convergence to the global maximum. The proposed scheme is tested on three datasets with three different CNNs. Promising results demonstrate the usefulness of our proposed approach in improving model accuracy due to the proposed activation function and Bayesian estimation of the parameters.

 Activation function, Deep neural networks, Optimization, MCMC, Hamiltonian dynamics

## I Introduction

Classification is a machine-learning task that identifies which objects are present in an image or video. It is critical for various applications such as computer vision [1, 2, 3], medical diagnostics [4], signal processing [5], and others. The process involves learning significant or nontrivial relationships from a set of training data and extending these relationships to interpret new test data [6]. The task involves categorizing elements into one of the finite set of classes by comparing the measured attributes of a given object with the known properties of objects to determine whether the object belongs to a specific category.

Convolutional Neural Networks (CNNs) [7, 8, 9, 10, 11] have become the industry standard in numerous applications over the past two decades, as they can process complex high-dimensional input data into simple low-dimensional concepts through a series of nonlinear transformations. Each feature layer in CNNs comprises features from the layer below, creating a hierarchical organization of ever-more-abstract concepts. They are particularly effective at capturing high-level abstractions in real-world observations, making them a popular choice for image classification tasks. CNNs can learn features from raw image pixels, reducing the need for manual feature engineering. This makes them well-suited for tasks such as object recognition, where the goal is to identify the presence and location of specific objects in an image. In this context, the activation function plays a critical role in learning representative features. The Rectified Linear Unit (ReLU) [12] is currently the most widely used activation function for neural networks. When provided with a positive argument, the ReLU activation function is either zero or the identity. ReLUs have the additional benefit of alleviating the vanishing gradient problem, in addition to providing sparse codes [13]. According to the state of the art, activation functions can be fixed or trainable during a learning phase [14]. In most cases, gradient descent is the most widely used method in the literature for parameter estimation.

From another side, Bayesian methods have advanced significantly in many domains over time and have numerous useful applications. The primary idea is to represent all uncertainties in the model using probabilities. Bayesian techniques are distinctive in that they treat the problem as an inference problem [15]. One of the most significant advantages is the ability to incorporate prior information about the model parameters and hyperparameters. Recent advancements in Markov Chain Monte Carlo (MCMC) methods [16, 17, 18, 19] make it easier to use Bayesian analyses in complex datasets with missing observations and to handle multidimensional outcomes. Recent studies have demonstrated that using a Bayesian framework in CNNs for the optimization process leads to more promising performances than standard gradient descent [20, 21].

In this study, we introduce a new trainable activation function. The parameters of the proposed function are automatically estimated from the data. For doing so, a Bayesian framework is used where these parameters as well as the network weights are assumed to be realizations of random variables. With an adequate likelihood, a hierarchical Bayesian model is built with priors and hyperpriors. An MCMC-based inference is then use, specifically a Gibbs sampler, to derive estimators from the target distributions. Our method is an extension of our previous work described in [22], where we employed non-smooth Hamiltonian methods to fit sparse artificial neural networks. Our main objective with the Bayesian scheme is to minimize the target cost function of the learning model. The use of non-smooth Hamiltonian techniques enables us to perform efficient and fast sampling, even when dealing with non-differentiable energy functions that arise due to the use of sparse regularization functions.

The contribution of this paper is therefore twofold: _i)_ proposing a new trainable activation function, and _ii)_ a more general Bayesian formulation than in [22] to integrate the estimation of all parameters from the data.
The rest of this paper is organized as follows. The state of the art is the focus of Section II. Then, the Problem statement is in Section III. In section IV we detail the adopted hierarchical Bayesian model. The proposed Bayesian inference scheme is developed in Section V and validated in Section VI. Finally, The conclusion and future work are drawn in Section VII.

## II Related work

Finding the best activation function to integrate into an architecture is a challenging task. This section proposes a taxonomy of different activation functions described in the literature. A primary established classification is based on the ability to modify the shape of the activation function during the training phase, resulting in two major categories that can be distinguished [14].

### _Fixed-shape activation functions_

This category pertains to the use of activation functions in neural network research, including sigmoid [23], hyperbolic tangent (tanh) [24], and ReLU, all of which have a defined form. The introduction of rectified functions, particularly ReLU, has led to a marked enhancement in neural network performance and heightened scientific interest. Consequently, this category can be subdivided into subcategories based on their distinctive characteristics.

**Classic activation functions:**

The findings presented in [25] demonstrated that a feed-forward and shallow network can efficiently manage any continuous function defined on a compact subset. For many years, bounded activation functions like Sigmoid and tanh were the preferred choices for neural networks, with researchers demonstrating their efficacy, particularly in shallow network architectures [26]. Although Sigmoid, Bipolar sigmoid [27], Hyperbolic tangent, Absolute value [28], and other bounded activation functions are commonly used, their efficacy is limited when training multi-layer neural networks due to the vanishing gradient problem [29].

**Rectifier-based activation functions:**

The primary advantage of using rectified activation functions is to alleviate the problem of vanishing gradient. The success of ReLU [30] has inspired the development of many new activation functions during the last years [31, 32]. While ReLU has numerous benefits such as solving the vanishing gradient issue and making sparse coding easier [33], it is not without flaws. The "dying" ReLU problem [34] and non-differentiability at zero are the main concerns.

To address these issues, several variations of ReLU have been developed, such as Leaky ReLU (LReLU) [34], Truncated rectified [35], softplus [36], Exponential linear unit (ELU) [37], E-swish [38], and Flatten-T Swish [39].

### _Trainable activation functions_

The concept of using trainable activation functions is not new in the field of neural network research, with many studies published on this topic as early as the 1990s [40, 41, 42]. However, the growing interest in neural networks in recent years has led researchers to reconsider the potential benefits of trainable activation functions in improving network performance.

In this section, we discuss the main strategies proposed in the literature for learning activation functions from data. Based on their primary characteristics, these strategies can be classified into three families.

**Parameterized standard activation functions:**

In [40], a generalized hyperbolic tangent function was proposed by introducing two trainable parameters \(\alpha\) and \(\beta\), which adjust the saturation level and slope, respectively. Similarly, a sigmoid function with two trainable parameters was used in [41] to modify the activation function's shape. These parameters are learned along with the network weights using the backpropagation algorithm. More recently, the work in [43] aimed to avoid manually setting the parameter of the ELU unit by proposing an alternative based on two trainable parameters. The proposed activation function, called PELU, is defined as follows:

\[PELU(x)=\begin{cases}\frac{\beta}{\gamma}x&\text{if }x\geq 0\\ \beta\times(exp(\frac{x}{\gamma})-1)&\text{otherwise,}\end{cases} \tag{1}\]

where \(\beta\) is a trainable parameter.

A flexible ReLU function has been proposed in [44]:

\[frelu(x)=ReLU(x+\alpha)+\beta, \tag{2}\]

where \(\alpha\) and \(\beta\) are parameters learned from data. This is done to capture negative information that is lost in the classic ReLU function [12].

The activation function introduced in [45] is another type of ReLU function that partially learns its shape from the training set. Indeed, it can modify the negative part of the data via the parameter \(\alpha\). This function is called Parametric ReLU (PReLU) and can be defined as follows:

\[PReLU(x)=\begin{cases}x&\text{if }x>0\\ \alpha\times x&\text{otherwise}\end{cases} \tag{3}\]

where the parameter \(\alpha\) is learned jointly with the model using a gradient method.

**Functions based on ensemble methods:**

Functions based on ensemble methods involve combining multiple basic activation functions to form a more complex function. In [46], a method for investigating activation functions built as compositions of various basic activation functions are proposed. Similarly, [47] introduces a similar method using a genetic algorithm composed of two new activation functions: Exponential Linear Sigmoid SquasHing (EliSH) and HardELiSH. EliSH is defined as follows:

\[EliSH(x)=\begin{cases}x/(1+e^{-x})&\text{if }x\geq 0\\ (e^{x}-1)/(1+e^{-x})&\text{otherwise.}\end{cases} \tag{4}\]
The negative part of EliSH is a multiplication of two functions, ELU and Sigmoid, while the positive part is shared with Swish. HardELiSH is defined as a multiplication of HardSigmoid and ELU in the negative portion and HardSigmoid and Linear in the positive part. In [48], another interesting activation function is proposed, called the Mexican Hat Linear Unit (MeLU). This activation function solves the problems of unstable learning related to trainable parameters. Unstable learning can cause a decrease in accuracy and an increase in generalization error when the model's performance varies significantly in response to slight changes in the data or parameters. Mexican hat-type functions have a smoother curve than ReLU, which prevents saturation and allows for optimal performance.

Let \(f\) be the function defined by

\[f_{\gamma,\lambda}(x)=max(\lambda-|x-\gamma|,0), \tag{5}\]

where \(\lambda\), \(\gamma\) are real numbers. This function returns zero when \(|x-\gamma|>\lambda\). Moreover, it increases with a derivative of 1 between \(\gamma-\lambda\) and \(\gamma\), then decreases using a derivative of -1 between \(\gamma\) and \(\gamma+\lambda\). MeLU is defined for each layer as

\[MeLU(x)=PReLU(x)+\sum_{j=1}^{k-1}c_{j}f_{\gamma_{j},\lambda_{j}}(x), \tag{6}\]

where \(k\) represents the number of parameters that can be learned in each neuron. The \(c_{j}\) parameters are learnable real numbers, while \(\gamma_{j}\), \(\lambda_{j}\) are fixed parameters.

In [49], the authors introduced a new function called Adaptive Blending Units (ABU) as a trainable linear combination of a set of activation functions

\[ABU(x)=\sum_{i=1}^{k}\alpha_{i}\times f_{i}(x), \tag{7}\]

where (\(\alpha_{1}\), \(\alpha_{2}\),..., \(\alpha_{k}\)) are parameters to be learned, and {\(f_{1}\)(.), \(f_{1}\)(.),..., \(f_{k}\)(.)} is a set of activation functions. The parameters \(\alpha_{i}\) are all initialized to the value \(\frac{1}{k}\) and are trained using the gradient descent method.

Likewise, similar methods were introduced recently like Kernel-based activation function [50, 51] and Trained activation function [52].

**Activation functions based on other techniques:**

In [42], a method using polynomial functions with adjustable coefficients is proposed. Similarly, in the context of fuzzy activation functions, a neural unit based on Type-2 fuzzy logic [53] was developed. Other works have proposed functions using interpolation and spline approaches [54]. However, depending on the chosen technique, these strategies may require additional input.

### _Comparison and analysis_

Learning activation functions is a popular topic in the field of machine learning because the performance of learning architectures can be improved by using more suitable activation functions. Various approaches have been explored to enhance these performances. Most trainable activation functions are variants of standard activation functions, whose shape is adjusted using trainable parameters. Several studies on trainable activation functions have reported substantial performance improvements compared to neural network architectures equipped with classic fixed activation functions such as ReLU or sigmoid.

Other trainable activation functions can be expressed as sub-networks nested within the main architectures or those based on different approaches than classical activation functions, such as fuzzy logic.

While these activation functions have significant potential to improve neural network model performance, their implementation can be more complex and require more time and resources for learning.

Despite encouraging results, it is still challenging to identify a strategy for automatic learning of an activation function that would solve different problems and significantly improve performance. Most trainable activation functions use the gradient descent method for hyperparameter estimation. The main limitations of these techniques lie in the computation time and gradient vanishing. This process prevents the network from learning deep features and may even lead to excessive processing capacity during training.

Indeed, neural networks can get stuck in local minima [55], which can harm the model's performance. This phenomenon is partly due to the gradient vanishing that occurs as derivatives decrease as the model deepens. This gradient decrease makes it harder to optimize the model, leading to a decrease in performance.

## III Problem statement

In the previous Section, we have examined different activation functions presented in the literature. In this paper, we introduce a modification of the MeLU activation function [56], while integrating the estimation of its parameters into a global Bayesian optimization framework. The choice of this function is mainly justified by its promising performance, as well as its form which promotes non-linearity and sparsity. However, the limits of the MeLU function are mainly related to memory requirements. This limit indicates that using this function may require a larger amount of memory compared to other activation functions, which can affect computational efficiency.

As a first contribution of this paper, the Modified Mexican ReLU (MMeLU) activation function is proposed to solve the complexity problem and improve model performance. The specificity of this new activation function is that it requires fewer parameters to estimate than MeLU. The second contribution is related to integrating the parameters estimation of the MMeLU function into a Bayesian optimizer [57], rather than using a standard optimization procedure such as the ADAM optimizer.

To define MMeLU, let

\[f_{\gamma,b}(x)=max(b-|x-\gamma|,0), \tag{8}\]
where \(\gamma\) and \(b\) are real numbers.

Using \(f_{\gamma,b}\), the proposed activation function MMeLU can be defined as

\[\textit{MMeLU}(x)=c\times f_{\gamma,b}(x)+(1-c)\times\textit{ReLU}(x) \tag{9}\]

where \(c\) is a real number belonging to the interval \([0,1]\). For the proposed MMeLU function, \(c\), \(b\), and \(\gamma\) are the parameters to be estimated.

The shape of the proposed MMeLU function, as well as those of competing functions ReLU, FReLU, PReLU, and MeLU, are illustrated in Figure 1. The activation function proposed in this paper is made up of a mixture of the ReLU and the Mexican hat functions. The curves in Figure 1((a)-(d)) clearly show the flexibility and non-linearity of our MMeLU function with different configurations of the parameters \(\gamma\), \(b\), and \(c\).

The Mexican hat function is often used as an activation function in neural networks due to its advantages. Mexican hat functions are continuous, which means that changes in inputs produce continuous changes in outputs. This is important in neural networks because it allows a gradual update of weights. Additionally, they have high representation capacity, which means that they can accurately model complex functions. This allows neural networks to learn non-linear relationships between inputs and outputs.

The Mexican hat function looks like a bell but with a peak in the center that makes it more pronounced for small values, as shown by the MMeLU curves in Figures 1 and 1. This means that for small input values, the Mexican hat function will have a stronger response than other activation functions, such as ReLU or FReLU. This behavior can be useful in certain situations, for example when the input data has a restricted value range or when the neural network's response needs to be more sensitive to small variations in the input data.

Let us now consider the convolutional neural network in Figure 2. The MMeLU activation function is applied after each convolutional layer, replacing the _ReLU_ function. This strategy has the potential to increase the non-linearity of the model and improve its ability to represent features in images. It is worth noting that when \(c\) is estimated close to 0 (\(c\sim 0\)), the MMeLU function tends to the ReLU behavior.

Regarding the model fitting, let us assume that the estimated label (or numerical value) is obtained by applying the proposed activation function MMeLU(x,\(W\)), where **x** is the input data and \(W\in\mathbb{R}^{N}\) denotes the weights vector. The model parameters (weight vector and parameters of the activation function) can be determined during the training phase using a generic error function \(D\) (Euclidean, Minkowski, etc.). Akin to [21], the target CNN is assumed to be sparse. We also use the same Bayesian formulation of the optimization problem. For the \(M\) input data, we can write

\[\widehat{W}=\arg\min_{W}\mathcal{L}(W)\] \[=\arg\min_{W,b,\gamma,c}\sum_{m=1}^{M}D(MMeLU(x^{m};W)-y^{m})+ \sum_{l=1}^{L}\lambda_{l}\|W^{l}\|_{1}, \tag{10}\]

where \(y^{m}\) is the ground truth for input data \(x^{m}\), \(L\) is the number of layers in the network, and \(\lambda_{l}\) is a regularization parameter to be estimated for layer \(l\) that balances the solution between the data attachment term and the \(\ell_{1}\) sparse regularization terms.

In the following Section, we formulate the adopted hierarchical Bayesian model to conduct the inference and fit the model weights and the activation function parameters.

## IV Hierarchical Bayesian Model

The problem of estimating the parameters of the MMeLU activation function is formulated within a Bayesian framework. In this sense, all parameters and hyperparameters are supposed to follow probability distributions. A likelihood is defined to model the relationship between the target weight vector, the activation function parameters, and the data. A prior distribution is defined to model the prior knowledge about the target weights and all activation function parameters.

### _Likelihood_

Following the principle of minimizing the error between the reference vector **y** (labels or continuous values) and its estimate \(\widehat{\textbf{y}}\), we define the likelihood distribution as

Fig. 1: Illustrations of ReLU, FReLU, PReLU, MeLU, and MMeLU curves with different configurations.
\[f(\mathbf{y},\mathbf{x};W,c,\gamma,b)\propto\] \[\prod_{m=1}^{M}\exp\left[-D;(MMeLU(x^{m};W,c,\gamma,b)-y^{m})\right]. \tag{11}\]

It is worth noting that when a Euclidean distance is used for \(D\), the adopted likelihood is nothing but a Gaussian distribution.

### _Priors_

In our model, the unknown parameters are grouped in the unknown vector \(\theta\) = {\(W\), \(c\), \(\gamma\), \(b\), \(\lambda\)}, where \(\lambda=\{\lambda_{1},\ldots,\lambda_{L}\}\).

**Prior for \(W\):**

To promote sparsity in the neural network, we use a Laplace distribution for the weight vector \(W\) akin to [21]:

\[f(W;\lambda)\propto\prod_{l=1}^{L}\prod_{k=1}^{K_{l}}\left[\frac{1}{\lambda_ {l}}\exp\left(-\frac{|W_{k}^{l}|}{\lambda_{l}}\right)\right], \tag{12}\]

where \(K_{l}\) is the number of weights in layer \(l\) of the network and \(\lambda_{l}\) is a parameter to be estimated.

**Prior for \(\lambda_{l}\):**

Since \(\lambda_{l}\in\mathbb{R}_{+}\), we chose to use an inverse gamma (IG) distribution:

\[f(\lambda_{l};\delta,\mu)=IG(\lambda_{l};\delta,\mu)\propto(\lambda_{l})^{-1- \delta}\exp\left(-\frac{\mu}{\lambda_{l}}\right), \tag{13}\]

where \(\delta\) and \(\mu\) are positive parameters that were fixed at \(10^{-3}\) to have a non-informative prior.

**Prior for \(c\):**

Regarding the parameter \(c\), we consider a uniform distribution over the interval \([0,1]\), denoted as

\[c\sim U_{[0,1]}(c). \tag{14}\]

**Prior for \(\gamma\) :**

Since \(\gamma\) is a real value, a Gaussian distribution is used as follows

\[f(\gamma;\sigma^{2})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{\gamma^{2 }}{2\sigma^{2}}\right), \tag{15}\]

where \(\sigma^{2}\) is a hyperparameter to be estimated.

**Prior for \(b\):**

Since \(b\) is a positive real number, an exponential distribution is used as follows:

\[f(b;\lambda_{b})\propto\begin{cases}\frac{1}{\lambda_{b}}\exp\left(-\frac{b}{ \lambda_{b}}\right);\text{if};b\geq 0\\ 0;\text{otherwise}.\end{cases} \tag{16}\]

where \(\lambda_{b}\) is a hyperparameter to be estimated. This prior penalizes large values of \(b\).

### _Hyperpriors_

Since \(\lambda_{b}\) and \(\sigma^{2}\) are positive real numbers, an inverse gamma (IG) distribution was used as a hyper-_a priori_:

\[f(\lambda_{b};\delta,\mu)=IG(\lambda_{b};\delta,\mu)\propto(\lambda_{b})^{-1- \delta}\exp\left(-\frac{\mu}{\lambda_{b}}\right) \tag{17}\]

and

\[f(\sigma^{2};\delta,\mu)=IG(\sigma^{2};\delta,\mu)\propto(\sigma^{2})^{-1- \delta}\exp\left(-\frac{\mu}{\sigma^{2}}\right), \tag{18}\]

where \(\delta\) and \(\mu\) are positive parameters that were fixed at \(10^{-3}\).

Fig. 2: A general diagram of a convolutional neural network with the MMeLU activation function.
## V Inference scheme

By adopting a Maximum _a Posteriori_ (MAP) approach, we first need to express the posterior distribution. Let \(\Phi_{e}\) be the hyperparameters to be estimated, represented by \(\Phi_{e}\) = \(\{\sigma^{2},\)\(\lambda_{b}\}\), and \(\Phi_{m}\) be the hyperparameters to be fixed, \(\Phi_{m}\) = \(\{\delta,\)\(\mu\}\). Using the likelihood, the prior distributions, and the defined hyperpriors, we can write the posterior distribution as:

\[f(\theta,\Phi_{e};y,\Phi_{m})\propto f(y;\theta)f(\theta;\Phi_{ e})f(\Phi_{e};\Phi_{m}) \tag{19}\]

which can be reformulated in a detailed version as

\[f(\theta,\Phi_{e};\mathbf{y},\mathbf{x},\Phi_{m})\propto\] \[\prod_{m=1}^{M}\exp\left[-D\ (MMeLU(x^{m};W,c,\gamma,b)-y^{m})\right]\times\] \[\prod_{l=1}^{L}\frac{1}{\lambda_{l}^{K_{l}}}\prod_{k=1}^{K_{l}} \left[\exp\left(-\frac{|W_{k}^{l}|}{\lambda_{l}}\right)\right]\times(\lambda_ {l})^{-1-\delta}\exp\left(-\frac{\mu}{\lambda_{l}}\right)\times\] \[\exp\left(-\frac{\gamma^{2}}{2\sigma^{2}}\right)\times\frac{1}{ \lambda_{b}}\exp\left(-\frac{b}{\lambda_{b}}\right)1_{\mathbb{R}_{+}}(b) \times 1_{[0,1]}(c)\times\] \[(\lambda_{b})^{-1-\delta}\exp\left(-\frac{\mu}{\lambda_{b}} \right)\times(\sigma^{2})^{-1-\delta}\exp\left(-\frac{\mu}{\sigma^{2}}\right).\]

One can clearly notice that the posterior in (IV-B) is complicated to deal with in order to derive close-form estimators. We, therefore, resort to numerical approximations using a Markov Chain Monte Carlo technique (MCMC) [17, 58]. Specifically, we use a Gibbs sampler to sequentially sample according to the conditional posteriors. To calculate the conditional distributions associated with each parameter of the model, one needs to integrate the joint posterior distribution in (IV-B) with respect to all the other parameters.

Regarding the parameter \(W\), calculations based on (IV-B) lead to the following form:

\[f(W;c,\gamma,b,\lambda)\propto\exp\left[-\sum_{l=1}^{L}\sum_{k=1 }^{K_{l}}\frac{|W_{k}^{l}|}{\lambda_{l}}\right]\times\] \[\exp\left[-\sum_{m=1}^{M}\left(D;(MMeLU(x^{m};W,c,\gamma,b)-y^{m })\right)\right]. \tag{21}\]

The conditional distribution for the parameter \(c\) is given by:

\[f(c;W,b,\gamma)\propto 1_{[0,1]}(c)\times\] \[\exp\left[-\sum_{m=1}^{M}\left(D;(MMeLU(x^{m};W,c,\gamma,b)-y^{m })\right)\right]. \tag{22}\]

For the parameter \(b\), the condition distribution is given by:

\[f(b;W,c,\gamma,\lambda_{b})\propto\exp\left(-\frac{b}{\lambda_{ b}}\right)\times\] \[\exp\left[-\sum_{m=1}^{M}\left(D;(MMeLU(x^{m};W,c,\gamma,b)-y^{m })\right)\right]. \tag{23}\]

As regards \(\gamma\), the conditional distribution writes:

\[f(\gamma;W,b,c,\sigma^{2})\propto\exp\left(-\frac{\gamma^{2}}{ 2\sigma^{2}}\right)\times\] \[\exp\left[-\sum_{m=1}^{M}\left(D(MMeLU(x^{m};W,c,\gamma,b)-y^{m })\right)\right]. \tag{24}\]

The conditional distribution for the parameter \(\lambda_{l}\) is given by:

\[f(\lambda_{l};\delta,\mu) \propto\lambda_{l}^{-1-(\delta+K_{l})}\exp\left(-\frac{\mu}{ \lambda_{l}}\right)\] \[\propto IG(\delta+K_{l},\mu). \tag{25}\]

For the hyperparameter vector \(\Phi_{e}\), it is necessary to calculate the conditional distributions from which it is possible to sample based on the likelihood and adopted priors.

The conditional distribution for the hyperparameter \(\lambda_{b}\) is given by:

\[f(\lambda_{b};,b,\mu,\delta) \propto\lambda_{b}^{-2-\delta}\exp\left(-\frac{b+\mu}{\lambda_{b}}\right)\] \[\propto IG(\delta+1,b+\mu) \tag{26}\]

The conditional distribution for the hyperparameter \(\sigma^{2}\) is given by:

\[f(\sigma^{2};\mu,\gamma,\delta) \propto(\sigma^{2})^{-1-\delta}\exp\left(-\frac{\gamma^{2}+2\mu}{ 2\sigma^{2}}\right)\] \[\propto IG(\delta,\gamma+2\mu). \tag{27}\]

The sampling scheme is summarized in Algorithm 1, where the model weights \(W\) and the parameters of the proposed MMeLU function are sampled.

```
Fix the hyperparameters \(\Phi_{m}\) ; for\(r=1,\ldots,S\)do  * Sample \(c\) according to \(f(c;W,b,\gamma)\) ;  * Sample \(\gamma\) according to \(f(\gamma;W,b,c,\sigma^{2})\) ;  * Sample \(b\) according to \(f(b;W,c,\gamma,\lambda_{b})\) ;  * Sample \(\sigma^{2}\) according to \(f(\sigma^{2};\mu,\gamma,\delta)\) ;  * Sample \(\lambda_{b}\) according to \(f(\lambda_{b};b,\mu,\delta)\) ;  * Sample \(\lambda_{l}\) according to \(f(\lambda_{l};\delta,\mu)\) \(\forall\ l\in\{1,\ldots,L\}\) ;  * Sample \(W\) as described in [22] ;  end for
```

**Algorithm 1**Gibbs sampler for the proposed method.

In Algorithm 1, \(S\) denotes the number of MCMC sampling iterations. After the burn-in period, the sampled coefficients are used to calculate the estimators \(\widehat{W}\), \(\widehat{c}\), \(\widehat{b}\), \(\widehat{\gamma}\), in addition to \(\widehat{\sigma^{2}}\), \(\widehat{\lambda}\) and \(\widehat{\lambda_{b}}\).

## VI Experimental validation

In order to validate the proposed method, three image classification experiments are conducted using different datasets: COVID-19 dataset including Computed tomography (CT) images for challenging classification [59], and two standard
datasets, namely Fashion-MNIST [60] and CIFAR-10 [61]. Table I illustrates the setting details of the different datasets.

To compare the proposed method with the state of the art, two types of optimizers are used: _i)_ our previous Bayesian optimizer that uses non-smooth Hamiltonian methods described in [22] with the standard ReLU activation function, and _ii)_ the standard Adam optimizer (with a learning rate of \(10^{-3}\)) with eight of the most well-known activation functions: ReLU, LReLU, ELU, PReLU, SeLU [62], swish, FReLU, and MeLU. As regards coding, we used python programming language with Keras and Tensorflow libraries on an Intel(R) Core(TM) i7-2720QM CPU 2.20GHZ architecture with 16 Go memory.

### _ConvNet models_

In this work, three CNN architectures are utilized. Similar to the LeNet model [63], the first one has two fully-connected and three convolutional (Conv-32, Conv-64, and Conv-128) layers (FC-64 and FC-softmax). The second architecture has nine convolutional (3XConv-32, 3XConv-64, and 3XConv-128) and three FC layers (FC-128, FC-64 and FC-softmax) which are organized similarly to VGG-Net [64]. These architectures are shown in Table II. The third model is a deeper CNN with 25 convolutional layers and 4 FC layers (for more information, see section VI-F). Each one uses convolutional layers with a stride size of 1 and \(3\times 3\) filters in addition to \(2\times 2\) max-pooling.

Deep neural networks are expanded with three regularizing strategies since they can easily overfit when trained on small datasets: Batch Normalization [65], \(\ell_{1}\) Regularization [66] and Dropout [67].

We chose to use different architecture depths in the experiments, mainly to test the ability of our proposed method to achieve better performance in the shallow model. In this sense, training with large and complex data can be expensive.

### _Sampling Results_

After using the proposed Bayesian optimization method to train the CNN models detailed above for the classification of Covid-19 CT images, we analyzed the convergence behavior. Figure 3 presents the sampling chains for the \(\gamma\), \(b\), and \(c\) parameters of the proposed MMeLU function (a-c), as well as the histograms of the corresponding samples (d-f). The sampling chains and histograms of the sampled coefficients confirm the good convergence properties of the designed Gibbs sampler. After a burn-in period of 350 iterations, the algorithm achieves stable convergence and exhibits a good mixing rate of the sampled chains.

Fig. 3: Sampling of parameters \(c\) (a,b), \(b\) (c,d), and \(\gamma\) (e,f): chains and histograms.
### _Experiment 1: COVID-19 classification using CT images_

This section examines the effectiveness of our approach in classifying Covid-19 infections from other pneumonia in CT data. The COVID-CT dataset1 includes 397 images that are negative for COVID-19 and 349 images that are positive for COVID-19, and belong to 216 patients. We used 566 images for the train and 180 images for the test.

Footnote 1: [https://www.kaggle.com/luisblanche/covidct](https://www.kaggle.com/luisblanche/covidct)

Table III presents compelling evidence that the proposed Bayesian method outperforms all other activation functions for both the \(CNN_{1}\) and \(CNN_{2}\) architectures. Moreover, except with respect to "ReLU + [22]", the computational time for convergence is shorter that all the other activation functions. With respect to our previous work ("ReLU + [22]"), it is worth noting that the proposed method performs better in terms of accuracy and loss value, which confirms the usefulness of the used trainable activation function. However, due to the use of additional parameters, the proposed method is approximately 15\(\%\) slower. These conclusions hold for both \(CNN_{1}\) and \(CNN_{2}\).

Notably, when regularization is employed, a considerable drop in performance is observed for all competing activation functions, which is largely attributable to the inherent difficulty of classifying CT images due to their content richness and the similarities between Covid-19 infection and other types of pneumonia.

Learning and test curves (accuracy and loss) illustrated in Figures 4 and 5 confirm the good behavior of the proposed method, which is not necessarily the case of the other competing models where a marked difference between the precision and loss curves can be noticed. For example, while the LReLU function introduces a negative bias that suppresses excessive activations, an inappropriate bias value can lead to underfitting. Similarly, although the ELU function allows negative activation, its exponential form can lead to an explosion of the activation value for large values of \(x\).

The Swish function is known to accelerate learning convergence, but it can also lead to overfitting by being more sensitive to outliers. Likewise, although the FReLU function can capture complex data patterns, it can also suffer from overfitting if the parameters are not well chosen.

These remarkable differences confirm the interest and efficiency of our MReLU function, which outperforms all competing activation functions in terms of accuracy and robustness to regularization.

### _Experiment 2: Fashion-MNIST image classification_

The learning performance of the competing activation function algorithms is assessed in this scenario using the standard _Fashion-MNIST_ dataset. A training set of 60,000 images is employed, with a test set of 10,000 images. Each example is a \(28\) grayscale image paired with a label from one of ten classes, with 7,000 images in each class. We used 48,000 images for the train set and 12,000 for the test set for model training.

Table IV reports the results obtained for the _Fashion-MNIST_ dataset. Our proposed method outperformed all other competing activation functions, achieving a minimum accuracy of 93% for both CNN models. Additionally, Table IV shows that the processing time for competing methods is greater than 150 minutes for \(CNN_{1}\), while only almost 112 minutes are required for our method. Similar conclusions hold with \(CNN_{2}\).

Moreover, our proposed MMeLU function showed similar global performance compared to "ReLU + [22]" for the two CNN models. One more advantage of MMeLU is its fully automatic estimation of parameters, which is not the case for ReLU + [22], where some parameters need to be manually set. Furthermore, our proposed method has demonstrated its ability to achieve remarkable performances in challenging cases, as shown in the first experiment.

### _Experiment 3: CIFAR-10 image classification_

In this experiment, the learning performance of the proposed model (made up of the trainable activation function and the Bayesian sparse optimization) is evaluated using the standard _CIFAR-10_ dataset. The CIFAR-10 dataset contains 60000 \(32\times 32\) color images divided into 10 classes with 6000 images per class. There are 50,000 training and 10,000 test images.

Table V presents the classification results for the _CIFAR-10_ dataset. It can be observed that the proposed Bayesian model performed well overall, even when multiple classes were used. In contrast, all competing activation functions (with Adam optimizer) had similar accuracy, around 88%, with a
loss rate almost double that of the proposed model for both architectures.

Moreover, the MMeLU activation function takes significantly less time for training than other activation functions + Adam, which often take a long time before convergence. For the competing functions when used with Adam, we noticed similar performance between the trainable and standard activation functions, with a slight superiority observed with the trainable function ELU.

The same conclusions can be drawn by examining the results of "ReLU + [21]" on this dataset. These findings suggest that the proposed MMeLU function provides better flexibility for the activation task, even for the multi-class case.

### _Experiment 4: Comparison on Deep CNN_

This section investigates the performance of our algorithm on the standard Fashion-MNIST dataset using a deep CNN. This CNN is claimed to be deeper than \(CNN1\) and \(CNN2\). It has 25 convolutional layers (5 X Conv3x3-32, 5 X Conv3x3-64, 5 X Conv3x3-128, 5 X Conv3x3-256 and 5 X Conv3x3-512) and four layers (FC-512, FC-256, Fc-128 and FC-softmax). All of them use convolutional layers with \(3\times 3\)

Fig. 4: Experiment 1: Train and test curves using \(CNN_{1}\) for all competing activation functions.
kernel filters as well as \(2\times 2\) max-pooling and stride size of 1.

The results demonstrate that our proposed MMeLU method maintains good performance on this deep architecture, as shown in table VI. The same conclusions can be drawn as from previous experiments regarding the superiority of MMeLU over competing activation functions, including our previous work in [22]. Interestingly, in addition to better scores in comparison to [22], the proposed method provides convergence time which is only 2 \(\%\) slower. It turns out that when the depth of the network increases (similarly the number of parameters), the additional cost to estimate the activation function parameters is no longer significant while reaching better performance.

## VII Conclusion

In this paper, we proposed a Bayesian approach for training sparse deep neural networks with trainable activation functions. Our method learns the weights and parameters of the proposed activation function directly from the data without any user configuration. Compared to competing algorithms, our

Fig. 5: Experiment 1: Train and test curves using \(CNN_{2}\) for all competing activation functions.
approach achieves promising results with better classification accuracy and high generalization performance, while also being computationally efficient, particularly in challenging cases. Our experiments on standard datasets such as Fashion-MNIST and CIFAR-10 demonstrate the precision and stability of our approach, as well as its general applicability.

Future work will focus on parallelizing the algorithm to enable GPU calculations and further reduce computational time.
