\section{Introduction}
Steganography is a technique which takes advantage of the redundancy in cover media to conceal secret information. It is one of the main strategies for covert communication. Various types of cover media have been considered in literature for steganography, including image \cite{8789545, 9124671}, text \cite{liu20116security}, audio \cite{8626153}, video \cite{zhang2016steganalytic} or 3D Mesh \cite{zhou2019distortion}. Traditional steganographic schemes focus on designing handcrafted algorithms to modify the cover media subtly for data embedding, which are difficult to achieve a good balance between the undetectability and embedding capacity (termed as the capacity for short). Recently, some deep neural network (DNN)-based steganographic schemes have been proposed \cite{baluja2019hiding, lu2021large, jing2021hinet}, which achieve superior performance to handcrafted ones.

A DNN-based steganographic scheme usually takes advantage of the encoder-decoder structure of DNN for secret data embedding and extraction, which contains two steganographic networks: a DNN-based secret encoder and a DNN-based secret decoder (termed as the secret encoder and secret decoder for short). The secret encoder hides the secret into a cover media to obtain a stego media (i.e., media with hidden data), while the secret decoder recovers the secret from the stego media. Despite the advantage, we find that it is not convenient to deploy the DNN-based steganographic schemes in real world applications. Let's take image steganography as an example, when the receiver receives the stego image, he could only conduct the secret recovery with the possession of the secret decoder, as shown in the upper part of Fig. \ref{fig:illustration}. Unlike the handcrafted steganographic tools, a DNN-based secret decoder is relatively large in size. It would not be secure for the sender and receiver to conduct another transmission for the delivery of the secret decoder.

One possible solution is to treat the secret decoder as ordinary data and hide them into popular cover media such as images or videos using the existing steganographic schemes. Since the size of a secret decoder is relatively large in terms of secret information, using traditional steganographic schemes would result in the need of a large amount of cover media data due to the limited capacity, which would cause a lot of communication burden. Some DNN-based steganographic schemes \cite{baluja2019hiding, weng2019high,jing2021hinet} are able to achieve a high capacity. However, they are designed specifically for hiding images with errors introduced during the secret recovery. It is not appropriate to apply them on the hiding of secret decoder which will be of no use after the recovery.

To deal with the issue of covert communication of the secret decoder, in this paper, we propose a novel scheme for steganography of steganographic networks, where we disguise a secret DNN model (i.e., the secret decoder) into a stego DNN model which performs an ordinary machine learning task (termed as a stego task). As shown in the lower part of Fig. \ref{fig:illustration}, the sender conducts the model disguising using our proposed method to generate a stego DNN model from the secret DNN model and a key. The stego DNN model is then uploaded to public model repositories such as the github, google drive and etc. The receiver downloads the stego DNN model and performs the model recovery using the key to restore the secret decoder. Other users (those without the key) could download the stego DNN model for an ordinary machine learning task. Compared with the existing steganographic schemes, our proposed method is tailored for the covert communication of DNN models, which has the following advantages:
\begin{itemize}
  \item [1)] The communication burden hardly increases before and after steganography, and there is no need to collect a large amount of cover data to accommodate the secret DNN model.
  \item [2)] The behaviour of covert communication is well protected by placing another functional DNN model for an ordinary machine learning task. The stego DNN model could be placed in the model repository that is public available without being noticed.
\end{itemize}

To disguise the secret DNN model, we first select a set of filters according to the corresponding gradients when the model is back-propagated on the datasets for secret and stego tasks. Based on the selected filters, we establish a stego DNN model from the secret DNN model by a partial optimization strategy. The aforementioned process are progressively done to obtain a stego DNN model with good performance on the secret and stego tasks. The main contributions of this paper are summarized below.

\begin{itemize}
    \item [1)] We tackle the problem of the covert communication of the steganographic network by disguising a secret DNN model into a stego DNN model without causing notice.
    \item [2)] We propose a partial optimization scheme for model disguising, where a stego DNN model is trained based on the secret DNN model for both the secret and stego tasks.
    \item [3)] We propose a progressive model disguising strategy for good performance of the stego DNN model.
\end{itemize}

\section{Related Works}
According to whether DNN has been utilized for data embedding or extraction, steganographic schemes could be mainly categorized into traditional steganography and DNN-based steganography. In this section, we give a brief review regarding these two types of steganographic schemes.

\subsection{Traditional Steganography}
Traditional steganography usually involves the alteration of cover data using hand-craft data hiding algorithms. Liao \textit{et al.} \cite{9124671} take a batch of images as the cover data for embedding the secret message, where the image texture features are measured and considered as an indicator to determine the capacity of each image. Zhang \textit{et al.} \cite{zhang2016steganalytic} propose a video steganography scheme by exploring the motion vector for data embedding. Lu \textit{et al.} \cite{8789545} conduct data embedding on a halftoned image, where the histgram of the pixels is slightly modified to accommodate the secret message. Zhou \textit{et al.} \cite{zhou2019distortion} propose to convey secret data through 3D meshes, where the importance of different regions in the 3D meshes is measured for adaptive data embedding. To ensure the undetectability, the capacity of these schemes are limited. Taking the image steganography as an example, each pixel in a grayscale image (8 bits) could accommodate 1 bit of secret message at most for high undetectability.

\subsection{DNN-based Steganography}
Most of the DNN-based steganographic schemes are proposed for image steganography. Some works take advantage of DNN to compute the embedding cost of each pixel \cite{tang2017automatic, tang2020automatic}, which still use handcrafted algorithms for data embedding when the embedding costs are ready.

Recently, more and more DNN-based steganographic schemes are proposed without any handcrafted processes, which usually take advantage of the encoder-decoder structure of DNN for data embedding and extraction. Zhu \textit{et al.} \cite{zhu2018hidden} pioneer the research of such techniques, where the secrets are embedded into a cover image using an end-to-end learnable DNN. Motivated by \cite{zhu2018hidden}, some researchers propose to hide secret images into cover images for high capacity data embedding \cite{baluja2019hiding, lu2021large, jing2021hinet,weng2019high}. The concept of hiding images into images is initially proposed by Baluja \textit{et al.} \cite{baluja2019hiding}, where a secret image is concatenated with the cover image into a set of different channels to learn an image encoding network for generating the stego image. The stego image is then fed into an image decoding network to extract the secret image. The image encoding and decoding networks are jointly learnt for minimized extraction errors and image distortion. Similarly, the works in \cite{jing2021hinet} and \cite{lu2021large} make use of the invertible neural network for image encoding and decoding, where the number of the channels in the secret image branch can be increased to remarkably improve the capacity. Weng \textit{et al.} \cite{weng2019high} propose to hide a secret frame into a cover frame by a temporal residual modeling technique for video steganography. Despite the advantage, the works in \cite{baluja2019hiding, lu2021large, weng2019high, jing2021hinet} are suitable for secrets in terms of images which can not be losslessly recovered in decoding.

Despite the progress in the area of steganography, little attention has been paid to the problem of covert communication of DNN models, which plays an important role in real world applications of the DNN-based steganographic schemes. In this paper, we try to tackle this problem by a network to network steganographic mechanism, which disguises the secret DNN model into a stego DNN model performing an ordinary machine learning task. The secret DNN model can be recovered from the stego DNN model with little performance degradation, while the stego DNN model is able to achieve high fidelity and undetectability.

\section{Problem Formulation}
\label{sec:3_1}
Given a secret DNN model to be disguised, our goal is to select and tune a subset of filters to preserve the function of the secret DNN, where the remaining filters are reactivated such that the whole model works on an ordinary machine learning task (i.e., a stego task) to establish a stego DNN model for model disguising. The secret DNN model could be recovered from the stego DNN model when necessary, as shown in Fig. \ref{fig:PF}.

Let's denote the secret and stego DNN models as ${\Theta}$ and $\tilde{\Theta}$, where $\tilde{\Theta}$ is an optimization of ${\Theta}$ after the model disguising. We use two datasets, a secret dataset and a stego dataset to conduct the model disguising for the secret and stego tasks, which are denoted as $D_{se}=\{\mathbf{x}_e, \mathbf{y}_e\}$ and $D_{st}=\{\mathbf{x}_t, \mathbf{y}_t\}$, with $\mathbf{x}_e$, $\mathbf{x}_t$ being the samples and $\mathbf{y}_e$, $\mathbf{y}_t$ being the labels. We further denote the model disguising process as $\mathcal{P}$ which can be formulated below

\begin{equation}
\begin{aligned} \label{eq1}
& \tilde{\Theta} = \mathcal{P}(\Theta, D_{se}, D_{st}, \mathcal{K}) \\
&\begin{array}{r@{\quad}r@{}l@{\quad}l}
s.t. & \left\{\begin{array}{lc}

    \Theta^\mathrm{S}\in \tilde{\Theta}  \\
    \mathcal{F}(\Theta^\mathrm{S}, \mathbf{x}_{e}) \to \mathbf{y}_{e}  \\
    \mathcal{F}(\tilde{\Theta}, \mathbf{x}_{t}) \to \mathbf{y}_{t}
    \end{array}\right. ,
\end{array}
\end{aligned}
\end{equation}
where $\Theta^\mathrm{S}$ is the selected subset for the secret task and $\mathcal{K}$ is a key. With the stego DNN model available, the receiver can recover the secret DNN model by
\begin{equation}
\Theta^\mathrm{S} = \mathcal{Q}(\tilde{\Theta}, \mathcal{K}),
\label{eq2}
\end{equation}
where $\mathcal{Q}$ refers to the model recovery.

The design of $\mathcal{P}$ and $\mathcal{Q}$ should satisfy the following properties for covert communication.
\begin{itemize}
  \item [1)] \textbf{Recoverability}: the performance of the recovered secret DNN model should be similar to its original version on the secret task.
  \item [2)] \textbf{Fidelity}: the performance of the stego DNN model should be high for the stego task.
  \item [3)] \textbf{Capacity}: the size of the stego DNN model should not expand too much compared with the secret DNN model for efficient transmission.
  \item [4)] \textbf{Undetectability}: it should be difficult for the attackers to identify the existence of the secret DNN model from the stego DNN model.
\end{itemize}

\begin{figure}[ht]
	\centering
	\includegraphics[width=0.9\linewidth]{figs/model_disgusing.pdf}
	\caption{An illustration for the DNN model disguising and recovery.}
	\label{fig:PF}

\end{figure}

\section{The Proposed Method}
\label{sec3}
Our model disguising starts with selecting a subset of filters from the convolutional layers in the secret DNN model, which is important on the secret task but trivial on the stego task. Then, we propose a partial optimization strategy to obtain a stego DNN model from the secret DNN model. The aforementioned process can be progressively performed for optimal performance. We also propose several strategies to deal with the parameters in other layers for real world applications.

\subsection{Filter Selection}
\label{sec:3_2}
Assume the secret DNN model (i.e., ${\Theta}$) contains $L$ convolutional layers, the filters in the $l$-th convolutional layer can be represented as a 4-dimensional tensor: $\mathcal{W}^{l} \in \mathbb{R}^{d^{l} \times c^{l} \times s^{l}_1 \times s^{l}_2 }$, where $d^{l}$ is the number of filters and each filter contains $c^{l}$ 2-dimensional spatial kernels. $s^{l}_1$ and $s^{l}_2$ correspond to the kernel's height and width, respectively. Denote $\mathcal{W}^{l}_{i,:,:,:}$ as the $i$-th filter and $\mathcal{W}^{l}_{:,j,:,:}$ as the $j$-th channel of all the filters in the $l$-th layer. Ignoring fully connected layer, we have: $\Theta = \{\mathcal{W}^{1}, \mathcal{W}^{2}, \cdots, \mathcal{W}^{L}\}$, where $\mathcal{W}^{l} = \{ \mathcal{W}^{l}_{1,:,:,:}, \mathcal{W}^{l}_{2,:,:,:}, \cdots, \mathcal{W}^{l}_{d^{l},:,:,:}\}$ filter-wise or $\mathcal{W}^{l} = \{ \mathcal{W}^{l}_{:,1,:,:}, \mathcal{W}^{l}_{:,2,:,:}, \cdots, \mathcal{W}^{l}_{:,c^{l},:,:}\}$ channel-wise.

The purpose of the filter selection is to choose the filters that are important on the secret task but trivial on the stego task. In our discussions, we term such filters as the important filters on model disguising. For a single machine learning task, researches have indicated that the importance of the filter could be measured according to the gradients of the filter weights \cite{8704878, yang2021robust}. In our case, however, we want to identify the filters which contain weights with large gradients on the secret task and with small gradients on the stego task. For the $i$-th filter in the $l$-th convolutional layer of a secret DNN model $\Theta$, we separately measure its importance on the secret task and stego task by $GoE^{l}_{i}$ and $GoT^{l}_{i}$, where

\begin{equation}
\begin{aligned} \label{eq3}
&\begin{array}{r@{\quad}r@{}l@{\quad}l}
& \left\{\begin{array}{lc}
    GoE^{l}_{i} = g(\mathcal{L}_e, \mathcal{W}^{l}_{i, :, :, :}, \mathbf{x}_{e}, \mathbf{y}_{e}) \\
    \qquad \quad \   + \ g(\mathcal{L}_e, \mathcal{W}^{l+1}_{:, i, :, :}, \mathbf{x}_{e}, \mathbf{y}_{e})  \\
     GoT^{l}_{i} = g(\mathcal{L}_t, \mathcal{W}^{l}_{i, :, :, :}, \mathbf{x}_{t}, \mathbf{y}_{t}) \\
    \qquad \quad \   + \ g(\mathcal{L}_t, \mathcal{W}^{l+1}_{:, i, :, :}, \mathbf{x}_{t}, \mathbf{y}_{t}),
    \end{array}\right.
\end{array}
\end{aligned}
\end{equation}
where
\begin{equation}
g(\mathcal{L}, \mathbf{w}, \mathbf{x}, \mathbf{y}) = \mathcal{AVG} \begin{pmatrix} \begin{vmatrix} \frac{\partial \mathcal{L}(\mathcal{F}(\Theta, \mathbf{x}), \mathbf{y})}{\partial \mathbf{w}} \end{vmatrix} \end{pmatrix},
\label{eq4}
\end{equation}

where $\mathcal{AVG}$ returns the mean of a 3-dimensional tensor, $\mathcal{L}_e$ and $\mathcal{L}_t$ are the loss functions for the secret and stego task, respectively. Here, we take the correlation between two consecutive convolutional layers into consideration. The reason is that the feature map generated by the filter $\mathcal{W}^{l}_{i,:,:,:}$ in the $l$-th layer will be convolved with the channels $\mathcal{W}^{l+1}_{:,i,:,:}$ in the ($l$+1)-th layer. There is a strong correlation between the weights in $\mathcal{W}^{l}_{i,:,:,:}$ and $\mathcal{W}^{l+1}_{:,i,:,:}$ for a certain task. Thus, we consider both the gradients of the weights in the current filter and those in the corresponding channel of the next layer. The importance of the filter on model disguising is then computed as
\begin{equation}
\alpha_i^{l}=GoE^{l}_{i} - \lambda_{g} GoT^{l}_{i},
\label{eq5}
\end{equation}
where $\lambda_{g}$ is a weight for balance. According to $\alpha_i^{l}$, we select the top $N$ most important filters for model disguising from the secret DNN model $\Theta$, which forms the subset $\Theta^\mathrm{S}$ for performing the secret task.

\subsection{Partial Optimization}
The purpose of partial optimization is to make $\Theta^\mathrm{S}$ work on the secret task and $\Theta$ work on the stego task. First of all, we fine-tune the filters in $\Theta^{\mathrm{S}}$ on the secret dataset by
\begin{equation}
    \Theta^{\mathrm{S}} = \Theta^{\mathrm{S}} - \lambda_{e} \nabla_{\Theta^{\mathrm{S}}} \mathcal{L}_{e}(\Theta^{\mathrm{S}}, D_{se}),
\label{eq6}
\end{equation}
where $\lambda_e$ is the learning rate of the secret task. Then, we freeze $\Theta^{\mathrm{S}}$, reinitialize and reactivate the remaining filters in $\Theta$ on the stego dataset $D_{st}$ for the following optimization problem:

\begin{equation}
 \min \limits_{\Theta} \sum \limits_{(\mathbf{x}_t, \mathbf{y}_t) \in D_{st}}  \mathcal{L}_{t}(\mathcal{F}(\Theta ,\mathbf{x}_t), \mathbf{y}_t), \quad
s.t. \  \Theta^{\mathrm{S}} \in \mathbb{C},
\label{eq7}
\end{equation}

where $\mathbb{C}$ is the constant space. To do so, we introduce $\mathrm{M}$, a binary mask with the same size as that of $\Theta$, for partial optimization, which prevents $\Theta^\mathrm{S}$ from updating during the backpropagation. Let $\lambda_t$ be the learning rate of the stego task and $\odot$ denotes the element-wise product, the secret DNN model is optimized below to form a stego DNN model $\tilde{\Theta}$:
\begin{equation}
\begin{aligned} \label{eq8}
& \Theta = \Theta - \lambda_{t} \ \mathrm{M} \odot \ \nabla_{\Theta} \mathcal{L}_{t}(\Theta, D_{st}), \\
\end{aligned}
\end{equation}
where

\begin{equation}
\begin{aligned} \label{eq9}
&\begin{array}{r@{\quad}r@{}l@{\quad}l}
& \mathrm{M}[i] = \left\{\begin{array}{lc}
    0,\quad & \Theta[i] \in \Theta^{\mathrm{S}}  \\
    1,\quad & else \end{array}\right. ,
\end{array}
\end{aligned}
\end{equation}
where $\Theta[i]$ refers to the $i$-th parameter in $\Theta$.

\begin{algorithm}[tb]
\caption{Progressive Model Disguising Algorithm.}
\label{alg:DDPA}
\textbf{Input}: Secret DNN $\Theta$,  Dataset $D_{se}$,  $D_{st}$, Thresholds $\tau_{se}$, $\tau_{st}$ \\
\textbf{Output}: Stego DNN $\tilde{\Theta}$
\begin{algorithmic}[1]
    \STATE $t \gets 1$, $\Theta^{\mathrm{S}}_0 \gets \Theta$, $\tilde{\Theta}_0 \gets \Theta$, $\alpha_{0}^{se} \gets 0$

    \STATE \textbf{while}  $\alpha_{t-1}^{se} < \tau_{se}$ \textbf{do}
    \STATE $\qquad$ Calculate the importance of filters for filter selection from $\Theta^{\mathrm{S}}_{t-1}$ by Eq.~(\textcolor[rgb]{1,0,0}{\ref{eq5}})
    \STATE $\qquad$ $\Theta^{\mathrm{S}}_t \gets \mathrm{Top} \  P_{t}$ important filters in $\Theta^{\mathrm{S}}_{t-1}$
    \STATE $\qquad$ Generate mask $\mathrm{M}$
    \STATE $\qquad$ Reinitialize the remaining filters
    \STATE $\qquad$ Conduct the model disguising according to Eq.~(\textcolor[rgb]{1,0,0}{\ref{eq6}}) and Eq.~(\textcolor[rgb]{1,0,0}{\ref{eq8}}) to get a stego DNN model $\tilde{\Theta}_t$
    \STATE $\qquad$ \textbf{if} $\alpha_t^{st} < \tau_{st}$ \textbf{do}
    \STATE $\qquad$ $\qquad$ Break
    \STATE $\qquad$ \textbf{end if}
    \STATE $\qquad$ $t \gets t+1$

    \STATE \textbf{end while} \\
    \STATE $\tilde{\Theta}=\tilde{\Theta}_t$
\end{algorithmic}
\end{algorithm}

\subsection{Progressive Model Disguising}
To train a disguised stego DNN model from the secret DNN model with good performance for the secret and stego tasks, we propose to progressively select the important filters $\Theta^{\mathrm{S}}$ from the secret DNN model $\Theta$ for model disguising. In particular, we iteratively select a set of important filters from  $\Theta^{\mathrm{S}}$ to train the stego DNN model, which is performed until the stego DNN model achieves satisfactory performance on the secret and stego tasks.

For simplicity, we denote the important filter set as $\Theta_t^{\mathrm{S}}$ and the disguised stego DNN model as $\tilde{\Theta}_t$ in the $t$-th iteration, where we set $\Theta_0^{\mathrm{S}}=\Theta$ and $\tilde{\Theta}_0=\Theta$ for initialization. In the $t$-th iteration, we select the top $P_{t}$ most important filters from $\Theta_{t-1}^{\mathrm{S}}$ to form $\Theta_{t}^{\mathrm{S}}$ for training a stego DNN model $\tilde{\Theta}_{t}$, where
\begin{equation}
  P_{t} = (\lambda_{p})^{t}\times V,
\label{eq10}
\end{equation}
where $\lambda_{p}<1$ is a decay factor and $V$ is the number of filters in $\Theta$. We terminate the iteration until the performance reduction of the stego DNN model is more than $\tau_{se}$ on the secret task, or less than $\tau_{se}$ and $\tau_{st}$ on the secret and stego tasks, respectively. Algorithm \ref{alg:DDPA} gives the pseudo code of the aforementioned process, where $\alpha_t^{se}$ and $\alpha_t^{st}$ are the performance reduction of the stego DNN model in the $t$-th iteration on the secret and stego tasks.

After the progressive model disguising, we record the property of each filter in $\tilde{\Theta}$ into $L$ binary streams $\mathbf{B}=\{b^1, b^2, \cdots, b^L\}$. Each convolutional layer corresponds to a $d^{l}$-bits binary stream $b^l$ with $d^{l}$ being the number of filters in this layer. Each bit corresponds to a filter with $0$ indicating that it belongs to $\Theta^{\mathrm{S}}$ and $1$ otherwise. We consider $\mathbf{B}$ as side information and randomly select a set of parameters from the stego DNN model according to a key $\mathcal{K}$ to host $\mathbf{B}$. In particular, we embed each bit of $\mathbf{B}$ into a selected parameter as follows. We transform the parameter into an integer and flip its least significant bit according to the bit to be hidden.
\begin{figure}[ht]
	\centering
	\includegraphics[width=0.9\linewidth]{figs/output_adaptation.pdf}
	\caption{Output layer adaptation.}
	\label{fig:op}

\end{figure}
The parameters with hidden data are then inversely transformed into the floating numbers to complete the embedding process.

In model recovery, the receiver uses the key $\mathcal{K}$ to recover the side information $\mathbf{B}$ from the stego DNN model, based on which we can identify the important filters $\Theta^{\mathrm{S}}$ to establish a secret DNN model for the secret task.

\subsection{Hiding other layers}

\textbf{Batch Normalization Layer}.~Normalization layers, such as Batch Normalization (BN) \cite{ioffe2015batch}, Group Normalization (GN) \cite{wu2018group}, Layer Normalization (LN) \cite{ba2016layer} and Instance Normalization (IN) \cite{ulyanov2016instance}, play a crucial role in accelerating the convergence of model training and performance boosting. All these normalization layers will normalize the feature maps of the intermediate layers according to some feature statistics. The BN utilizes the moving averaged mean and standard deviation of the feature maps during training, which are completely different between the secret and stego learning tasks. Therefore, when using BN, we store the corresponding feature statistics of the secret DNN model as the side information to be hidden into the stego DNN model with $\mathbf{B}$. Other types of normalization layers, \textit{i.e.} GN, LN, IN, do not introduce additional overhead.

\textbf{Fully Connected Layer}.~
Fully connected layers are very common in DNN models, We treat them as special convolutional layers. For a fully connected layer with $f$ neurons, we represent it as $f$ filters with each being a $c\times1\times1$ tensor, where $c$ is the number of backward neurons that are fully connected to the layer.

\textbf{Output Layer Adaptation}.~In real world applications, the dimension of the output layers are usually different between the secret DNN model and the stego DNN model. For simplicity, we denote the output layers of the secret and stego DNN model as $layer_e$ and $layer_t$, respectively. We further denote the number of neurons in $layer_e$ and $layer_t$ as $O_e$ and $O_t$. When $O_e<O_t$, we upsample $layer_e$ to make it the same dimension as that of $layer_t$ for adaptation. When $O_e>O_t$, however, it is not appropriate to directly downsample $layer_e$ because such a strategy will inevitably distort the function of the secret DNN model. For remedy, we propose to add additional neurons into $layer_e$ to disguise it as the second to last fully connected layer. Then, we newly add a final layer for adaptation, which has the same dimension as that of $layer_t$. Fig. \ref{fig:op} illustrates our output layer adaptation for different cases. All the newly added parameters will not be selected into $\Theta^\mathrm{S}$ to perform the secret task.

\section{Experiments}
\subsection{Setup}
We conduct two types of experiments for our proposed method: 1) steganography of steganographic networks (SSN), where the secret DNN model is a secret decoder, and 2) steganography of general networks (SGN), where the secret DNN model is a DNN model for a general machine learning task.

In SSN, we use the secret decoder of HiDDeN \cite{zhu2018hidden} as the secret DNN model for the secret task, and we randomly select 11000 images from the COCO dataset  \footnote{https://cocodataset.org/\#home} to form the secret dataset which is split into 10000/1000 images for training and testing. We use the bit error rate (BER) between the original and decoded secret information as the performance indicator of the secret decoder. We take the GTSRB \cite{stallkamp2012man} classification task as the stego task for model disguising. We use the GTSRB dataset as the stego dataset, where 80\%/20\% of the images are randomly selected for training/testing. We adopt the classification accuracy (ACC) to be the indicator to evaluate the performance of the stego DNN model.

In SGN, we evaluate our method on two well-known DNN models: ResNet18 \cite{he2016deep} and U-net \cite{ronneberger2015u}. For ResNet18, we assume the Fashion-MNIST \cite{xiao2017fashion} classification task as the secret task and take the CIFAR10 \cite{krizhevsky2009learning} classification task as the stego task for model disguising. For both the secret and stego datasets, we use their default partitions for training and testing with ACC being the performance indicator.

For U-net, we assume the image denoising task as the secret task, where the U-net is trained/tested on a secret dataset with 10000/1000 images randomly selected from ImageNet \cite{deng2009imagenet} after adding Gaussian noise. We compute the peak single noise ratio (PSNR) to evaluate the denoising performance. We take image segmentation task as the stego task for model disguising, where the oxford-pet dataset \cite{parkhi2012cats} is adopted as the stego dataset. We randomly separate the Oxford-Pet dataset into three parts, including 6000 images for training, 1282 images for validating and 100 images for testing. We employ the Mean Intersection over Union (mIOU) as the performance indicator for the stego DNN model.

For both SSN and SGN, we set $\lambda_g$=0.01, $\lambda_e$=$\lambda_t$=0.001, $\lambda_p$=0.9, $\tau_{st}$=0.01. For $\tau_{se}$, we set 0.0001, 0.01, and 0.5 for the secret decoder of HiDDeN, ResNet18 and U-net, respectively. We reinitialize the remaining filters by Kaiming initialization \cite{he2015delving}. All the models take the BN as normalization layer and are optimized using Adam \cite{kingma2014adam} optimizer. All our experiments are conducted on Ubuntu 18.04 system with four NVIDIA
RTX 1080 Ti GPUs.

\begin{table}[t]
\renewcommand{\arraystretch}{1.0}
  \centering

  \resizebox{\linewidth}{!}{
\begin{tabular}{c|c|c|c}
    \hline
    \multirow{2}{*}{Secret} & SSN                        & \multicolumn{2}{c}{SGN} \\
    \cline{2-4}
    \multirow{2}{*}{DNN}    &   HiDDeN-                   & Fashion-MNIST-                    &  ImageNet-    \\
    \multirow{2}{*}{Model} &   Secret Decoder                 & Classification            &  Denoising   \\

                          &   BER \ $(\times10^{-5})$                   &  ACC(\%)                  &   PSNR(dB)            \\
    \hline
    Original          &   6.68                      & 93.99                     &   28.98               \\
    \hline
    Recovered    &   6.68                      & 93.80                     &   28.92               \\
\hline

\end{tabular}}
\caption{The recoverability of the proposed method.}
  \label{tab:Recoverability}
\end{table}

\begin{table}[t]
\renewcommand{\arraystretch}{1.0}
  \centering

  \resizebox{0.91\linewidth}{!}{
\begin{tabular}{c|c|c|c}
    \hline
                             & SSN                        & \multicolumn{2}{c}{SGN} \\
    \cline{2-4}
    \multirow{1}{*}{DNN} &   GTSRB-                  & CIFAR10-                 &  Oxford-Pet-    \\
    \multirow{1}{*}{Model} &   Classification            & Classification            &  Segmentation   \\

                          &   ACC (\%)                  &  ACC (\%)                 &  mIOU (\%)            \\
    \hline
    Cover \quad    &   99.97                      & 93.46                     &   87.68               \\
    \hline
    Stego  \quad     &   99.81                   & 93.02                & 87.05          \\
    \hline
\end{tabular}}
\caption{The fidelity of the proposed method.}
  \label{tab:Fidelity}
\end{table}

\begin{table}[t]
\renewcommand{\arraystretch}{1.0}
  \centering

  \scalebox{0.84}{
  \resizebox{\linewidth}{!}{
\begin{tabular}{c|c|c|c}
    \hline
   \multirow{2}{*}{Expansion}             &   SSN                       & \multicolumn{2}{c}{SGN}    \\
    \cline{2-4}

    \multirow{2}{*}{Rate}  &   HiDDeN-            & \multirow{2}{*}{ResNet18}              & \multirow{2}{*}{U-Net}     \\
                     &   Secret Decoder            &                       &     \\
    \hline
    $e \ (\times10^{-3})$ &   1.65                 & 0.00        & 0.14               \\

\hline
\end{tabular}}}
\caption{The capacity of the proposed method.}

  \label{tab:Capacity}
\end{table}

\subsection{Recoverability} We evaluate the recoverability of our proposed method according to the performance of the secret DNN model before and after the model disguising, i.e., the performance of the original and recovered secret DNN models (from the stego DNN model). Table \ref{tab:Recoverability} gives the results for different steganographic scenarios. It can be seen that, the performance of the secret decoder is not affected at all after model disguising in the SSN case. For SGN, the performance of the secret DNN models slightly decrease, which does not affect the function of the DNN models for the secret tasks.

\subsection{Fidelity}
To evaluate the fidelity of the stego DNN model, we train a cover DNN model which has the same architecture as the stego DNN model for the stego task. The training is conducted on the training set of the stego dataset which has been used for model disguising, where all the parameters in the cover DNN model are tuned to achieve the best performance. Furthermore, we use the same test set to evaluate the performance of the two models.

Table \ref{tab:Fidelity} presents the performance of the stego DNN model and the cover DNN model for both SSN and SGN. For SSN, it can be seen that the ACC of the stego DNN model is reduced with less than $0.2\%$ after the model disguising. The results of the SGN are similar to those of SSN. In particular, after model disguising, the ACC/mIOU of the stego DNN model is slightly degraded when compared with the corresponding cover DNN model on the CIFAR100/Oxford-Pet dataset. These results indicate that the fidelity of our stego DNN model is favorable, which has little performance degradation caused by model disguising on the stego task.

\subsection{Capacity}
The capacity of the proposed method is evaluated in terms of the expansion rate defined as
\begin{equation}
e=\frac{N_{stego} }{N_{sec}}-1,
\end{equation}
where $N_{sec}$ and $N_{stego}$ refer to the number of parameters in the original secret DNN model and the corresponding stego DNN model, respectively. Thus, we have $e \geq 0$, which is closely related to the number of parameters added when adjusting the output layer.  The lower the expansion rate, the higher the capacity, meaning that there is less overhead created for transmission after the model disguising. Table \ref{tab:Capacity} gives the expansion rates $e$ for different steganographic scenarios. It can be seen that, for both SSN and SGN, the expansion rates are extremely small with less that $2 \times 10^{-3}$. In the ResNet18 case, the output layers are with the same dimension for the secret and stego tasks, so the parameters do not expand at all with $e=0$.

\begin{table}[t]
\renewcommand{\arraystretch}{1.0}
  \centering

  \scalebox{0.61}{
  \resizebox{\linewidth}{!}{
\begin{tabular}{c|c}
    \hline
    Classifier & Detection Accuracy (\%) \\
    \hline
     SVM   & 47.50 \\
    \hline
     MLP   & 50.00 \\
\hline
\end{tabular}}}
\caption{The undetectability of the proposed method.}
  \label{tab:Undetectability}
\end{table}

\begin{table*}[t]
\renewcommand{\arraystretch}{1.0}
  \centering

  \resizebox{\linewidth}{!}{
\begin{tabular}{c|c|c|c|c}
\hline

           \multirow{2}{*}{Category}     & \multirow{2}{*}{Method}  &\multirow{2}{*}{Capacity}   &  Recoverability &  Undetectability \\

                &        &   & BER(\%) / PSNR(dB)   &   $P_\mathrm{E}$(\%) / Detection accuracy (\%)\\

           \hline
                              & Tang \textit{et al}.\cite{tang2017automatic}&200.00MB (0.40bpp)    &0 / - &17.40 / -    \\
                              &Yang \textit{et al}. \cite{yang2019embedding}&160.00MB (0.50bpp)    &0 / - &29.71 / -    \\
\multirow{1}{*}{Traditional}   &Liao \textit{et al.} \cite{liao2019new} & 160.00MB (0.50bpp)    &0 / -    & 25.83 / -\\
\multirow{1}{*}{steganography}  &Tang \textit{et al}. \cite{tang2020automatic}& 200.00MB (0.40bpp)   &0 / - &32.19 / -    \\
                                  &Lu \textit{et al.} \cite{8789545}       & 655.35MB (0.015bpp)    &0 / -    & 34.17 / -\\
                                 &Liao \textit{et al.} \cite{9124671}       & 160.00MB (0.50bpp)    &0 / -     & 35.40 / -\\

\hline
                            &Zhu \textit{et al}. \cite{zhu2018hidden}     & 10.00MB (8.00bpp)  &- / 35.70  & - / 76.49\\
\multirow{1}{*}{DNN-based} &Bakuja \textit{et al}. \cite{baluja2019hiding}& 10.00MB (8.00bpp) &- / 34.13 &- / 99.67    \\
\multirow{1}{*}{steganography}   &Weng  \textit{et al}. \cite{weng2019high}      &  10.00MB (8.00bpp) &- / 36.48 &- / 75.03    \\
                                &Jing  \textit{et al}. \cite{jing2021hinet}     & 10.00MB (8.00bpp)  &- / 46.78 &- / 55.86     \\
\hline
                                &\multirow{3}{*}{\textbf{Ours}}       & \multirow{3}{*}{10.00MB ($e$=0.00)}     &   Performance reduction   &  \multirow{2}{*}{- / 47.50 (SVM)} \\
                                &           &  & 0.00\% in BER, 0.19\% in  &\multirow{2}{*}{- / 50.00 (MLP)} \\
                                &           &  & ACC, 0.06dB in PSNR &  \\
\hline
\end{tabular}}
\caption{Performance comparisons among different schemes for hiding DNN models.}
  \label{tab:comparison}
\end{table*}
\subsection{Undetectability}
\label{un}
To evaluate the undetectability, we first establish a model pool with different stego DNN models and cover DNN models trained. We select five benchmark datasets: MNIST \cite{lecun1998gradient}, Fashion-MNIST, CIFAR10, CIFAR100, GTSRB to generate $\mathrm{A}^2_5 = 20$ different secret-stego dataset pairs. Then, we train five well-known DNN models: LeNet-5 \cite{lecun1998gradient}, AlexNet \cite{krizhevsky2012imagenet}, VGG11, VGG19 \cite{simonyan2015vgg} and ResNet34 on the 20 dataset pairs to generate 100 different stego DNN models. Since we can only obtain 25 cover DNN models on the five datasets and five DNN models, we train each cover DNN model four times with different initialized parameters for augmentation. Take the three pairs of cover and stego DNN models (HiDDeN-Secret Decoder/ResNet18/U-net) previously trained into consideration, we obtain a total of 103 stego DNN models and 103 cover DNN models in our model pool.

Next, we compute a 100 bin-histogram from the parameters in each model to get a 100-dimensional feature vector for training a support vector machine (SVM) and a multilayer perceptron (MLP) classifier. Specifically, we randomly select 83 pairs of such features from the cover and stego DNN models for training and the remaining 20 pairs are used for testing. The detection accuracy are shown in Table \ref{tab:Undetectability}, we can see that both the SVM and MLP classifiers are not able to detect the existence of the secret DNN models from the stego DNN models with detection accuracy close to $50\%$.

\subsection{Comparisons}

In this section, we compare our proposed method with the existing state-of-the-art steganographic schemes for hiding DNN models in terms of capacity, recoverability, and undetectability. For fair comparisons among different types of steganographic schemes, we assume the data to be hidden is a secret DNN model with size of 10MB for all the methods. The capacity is evaluated as the size of stego data, the recoverability is reported as the BER or PSNR of the recovered secrets for the existing schemes, the undetectability is indicated as the average detection error $P_{\mathrm{E}}$ (which is the mean of false alarm rate and missed detection rate) or the detection accuracy (50\% means random guess).

Table \ref{tab:comparison} reports the comparisons among different schemes. For the traditional steganographic schemes \cite{liao2019new, 8789545, 9124671, tang2017automatic,yang2019embedding, tang2020automatic}, we duplicate the $P_{\mathrm{E}}$ reported at a certain payload for undetectability from each paper, where the payload is the bit per pixel (bpp), i.e., how many bits can be embedded in an 8-bit gray-scale pixel. For the DNN-based steganographic schemes \cite{baluja2019hiding, zhu2018hidden, weng2019high, jing2021hinet}, we duplicate the PSNR (on imagenet) and detection accuracy at a payload of 8bpp from the results reported in \cite{jing2021hinet}. Here, we assume the 10MB secret data is a natural image with meaningful content for these DNN steganographic schemes.

It can be seen from Table \ref{tab:comparison} that the size of our stego data is less than 1/10 of size of the stego data that are required using the existing traditional steganographic schemes. Compared with the existing DNN-based steganographic schemes \cite{zhu2018hidden, baluja2019hiding, weng2019high, jing2021hinet}, which are not able to losslessly recover the secret and are not suitable for hiding DNN models, our scheme achieves better undetectability (closer to 50\%). In addition, all the existing DNN-based steganographic schemes involve the change of data format between the secret-data (a DNN model) and the stego data (images).

To demonstrate that the existing DNN-based steganographic schemes are not suitable for hiding DNN models, we further adopt a recent DNN-based steganographic scheme
 \cite{jing2021hinet} to hide two DNN models: the HiDDeN secret decoder and the ResNet18.
Specifically, we convert the parameters of a secret DNN model to a set of binary streams to form a set of meaningless images with the size of $256 \times 256 \times 3$. Next, we embed each image into a cover image randomly selected from the COCO dataset using the encoder proposed in \cite{jing2021hinet}, which produces a stego image for decoding. Finally, we recover the images from the stego images using the corresponding decoder, from which we restore the DNN model. Table \ref{tab:NtoI} gives the performance of the recovered DNN models using our proposed method and the method proposed in \cite{jing2021hinet}. It can be seen that, by using the DNN-based steganographic scheme, the recovered DNN models do not work at all.

\begin{table}[t]
\renewcommand{\arraystretch}{1.0}
  \centering

  \scalebox{0.9}{
  \resizebox{\linewidth}{!}{
\begin{tabular}{c|c|c}
    \hline
     \multirow{2}{*}{Steganographic} & HiDDeN- &  \multirow{2}{*}{ResNet18}                        \\

    \multirow{2}{*}{Method} & Secret Decoder &        \\

                             & BER (\%) &   ACC(\%)  \\
    \hline
    Jing  \textit{et al}. \cite{jing2021hinet}  &  50.13 &  10.00  \\
    \hline
     \textbf{Ours}      &        6.68 $\times 10^{-3}$   &  93.80   \\

\hline
\end{tabular}}}
\caption{The performance of the recovered DNN models using different schemes.}
  \label{tab:NtoI}
\end{table}

\section{Conclusion}

In this paper, we propose a novel steganographic scheme to tackle the problem of covert communication of secret DNN models. Our scheme is able to disguise a secret DNN model into a stego DNN model which performs an ordinary machine learning task without being noticed. And the secret DNN model can be recovered from the stego DNN model on the receiver side. In model disguising, we first select a set of filters from the secret DNN model, which are important on the secret task but trivial on the stego task. We fine tune these filters on the secret task and then reactivate the other filters for the stego task according to a partial optimization strategy to train a stego DNN model. These are progressively done to obtain a stego DNN model that works well on both the secret and stego tasks. Various experiments have been conducted to demonstrate the advantage of our proposed method for covert communication of the DNN models, which achieves satisfactory performance in terms of recoverability, fidelity, capacity and undetectability.

\section{Acknowledgments}
This work was supported in part by the National Natural Science Foundation of China under Grants 62072114, U20A20178, U20B2051, U1936214, in part by the Project of Shanghai Science and Technology Commission 21010500200.


\section{Introduction}
Recent advances in deep neural networks have demonstrated remarkable performance in a wide variety of applications, ranging from recommendation systems and improving user experience to natural language processing and speech recognition
\cite{Abiodun2018StateoftheartIA}. Nevertheless, blindly relying on the outcome of these models  can have harmful effects, especially in high-stake domains such as healthcare and autonomous driving, as models can provide inaccurate predictions when queried in out-of-distribution data points
\cite{Amodei2016ConcretePI}.
Consequently, correctly quantifying the uncertainty of models' predictions is an admissible mechanism to
distinguish where a model can or cannot be trusted, and thus, increases the transparency of models about their capabilities and limitations \cite{Abdar2021ARO}. Uncertainty Quantification (UQ) is important for a variety of reasons. For instance, in order to preserve the model’s credibility, it is essential to report and communicate the encountered uncertainties regularly \cite{Volodina2021TheReproducibility}. Additionally, models' predictions are inevitably uncertain in most cases, which has to be addressed to increase their transparency, trustworthiness, and reliability.

In the machine learning literature, uncertainty is usually decomposed into two different types, namely aleatoric uncertainty and epistemic uncertainty \cite{Kiureghian2009AleatoryOE}. \textit{Aleatoric} uncertainty, aka data uncertainty, refers to the inherent uncertainty that stems from the data itself, e.g., noise. On the other hand, \textit{epistemic} uncertainty, also called model uncertainty, is the type of uncertainty that occurs due to the lack of sufficient data. While data uncertainty \textit{cannot} be alleviated, model uncertainty can be addressed by e.g., acquiring more data. Let $\boldsymbol{\sigma}_a^2$ and $\boldsymbol{\sigma}_e^2$ denote the aleatoric and epistemic uncertainties, respectively. Since the distinction between the two is imprecise to some degree \cite{Sullivan2015IntroductionTU}, we focus on the predictive (total) uncertainty, which is defined as the sum of the two
\begin{equation}
\label{eq:unc}
    \boldsymbol{\sigma}_p^2 = \boldsymbol{\sigma}_a^2 + \boldsymbol{\sigma}_e^2.
\end{equation}
Accordingly, the approaches developed for uncertainty estimation can be categorized into three groups: Bayesian UQ methods, ensemble UQ methods, and a combination of both, i.e., Bayesian ensemble UQ  \cite{Abdar2021ARO}. In this paper, we focus on ensemble UQ techniques, either Bayesian or non-Bayesian, as this  group is less explored compared to the solely Bayesian techniques. An ensemble model aggregates the predictions of multiple individual base-learners (or ensemble members), which in our case  are neural networks (NNs), and the empirical variance of their predictions gives an approximate measure of uncertainty. The idea behind this heuristic is highly intuitive: the more the base-learners disagree on the outcome, the more uncertain they are. Therefore, the goal of ensemble members is to have a great level of disagreement (variability) in the areas where little or no data is available, and to have a high level of agreement in regions with abundance of data \cite{Pearce2018UncertaintyIN}.

In this paper, we propose a novel method, called \textit{Random Activation Functions Ensemble (RAFs Ensemble)}, for a more robust uncertainty estimation in (deep) neural networks. RAFs Ensemble is developed on top of Anchored Ensemble technique, proposed by \cite{Pearce2018UncertaintyIN}, however, instead of initializing each NN member in the ensemble with the same activation function, the NNs in RAFs Ensemble are accommodated with different (random) activation functions in the hidden layers. This simple, yet crucial, modification greatly improves the overall diversity of the ensemble, which is one of the most important components in forming a successful ensemble. We empirically show that RAFs Ensemble provides high quality uncertainty estimates compared to five state-of-the-art ensemble methods, that is Deep Ensemble \cite{Lakshminarayanan2017SimpleAS}, Neural Tangent Kernel Gaussian Process Parameter Ensemble \cite{He2020BayesianDE}, Anchored Ensemble \cite{Pearce2018UncertaintyIN}, Bootstrapped Ensemble of NNs Coupled with Random Priors  \cite{Osband2018RandomizedPF}, and Hyperdeep Ensemble \cite{Wenzel2020HyperparameterEF}. The comparisons are performed in a wide range of regression tasks on both synthetic and real-world datasets in terms of negative log-likelihood and root mean squared error.

\section{Related Work}

Uncertainty Quantification (UQ) is an active field of research and various methods have been proposed to efficiently estimate the uncertainty of machine learning models (see \citeauthor{Abdar2021ARO} \citeyear{Abdar2021ARO} for an extensive overview).
While most research focuses on Bayesian deep learning \cite{Srivastava2014DropoutAS, Blundell2015WeightUI, Sensoy2018EvidentialDL, Fan2020BayesianAM, Jrvenp2020BatchSA, Charpentier2020PosteriorNU}, deep ensemble methods, which benefit from the advantages of both deep learning and ensemble learning, have been recently leveraged for empirical uncertainty quantification \cite{Egele2021AutoDEUQAD, Hoffmann2021UncertaintyQB, Brown2020UncertaintyQI, Althoff2021UncertaintyQF}. Although Bayesian UQ methods have solid theoretical foundation, they often require significant changes to the training procedure and are computationally expensive compared to non-Bayesian techniques such as ensembles
\cite{Egele2021AutoDEUQAD, Rahaman2021UncertaintyQA, Lakshminarayanan2017SimpleAS}.

\citeauthor{Lakshminarayanan2017SimpleAS}~\shortcite{Lakshminarayanan2017SimpleAS} are among the first to challenge  Bayesian UQ methods by proposing Deep Ensemble, a simple and scalable technique,
that demonstrates superb empirical performance on a variety of datasets. However, one of the challenges of ensemble techniques when quantifying uncertainty is that they tend to give overconfident predictions \cite{Amodei2016ConcretePI}.
To address this, \citeauthor{Pearce2018UncertaintyIN}~\shortcite{Pearce2018UncertaintyIN} propose to also regularize the model's parameters w.r.t. the initialization values, instead of zero, leading to Anchored Ensembles, which additionally allows for performing Bayesian inference in NNs.

\citeauthor{He2020BayesianDE}~\shortcite{He2020BayesianDE} relate Deep Ensembles to Bayesian inference using neural tangent kernels. Their method, i.e., Neural Tangent Kernel Gaussian Process Parameter Ensemble (NTKGP-param), trains all layers of a finite width NN, obtaining an exact posterior interpretation in the infinite width limit with neural tangent kernel parameterization and squared error loss. They prove that NTKGP-param is always more conservative than Deep Ensemble, yet, its advantages are generally not clear in practice.

A prominent advance to the Bayesian ensemble UQ methods is the bootstrapped ensemble of NNs coupled with random priors, proposed by \cite{Osband2018RandomizedPF}, in which,
the random prior function and neural models share an input and a summed output, but the networks are the only trainable parts, while the random prior remains untrained throughout the whole process.
Furthermore, \citeauthor{Wenzel2020HyperparameterEF}~\shortcite{Wenzel2020HyperparameterEF}
exploit an additional source of randomness in ensembles by designing ensembles not only over weights, but also over hyperparameters. Their method, called Hyperdeep Ensemble, demonstrates high accuracy for a number of different classification tasks. Nevertheless, despite the recent contributions in ensemble UQ methods, the research in this direction still needs further advancement.

\section{Toward Robust Uncertainty Estimation}
\subsection{Preliminaries}
Following the notations of \cite{Lakshminarayanan2017SimpleAS},
let $S_{train}$ be a training dataset  consisting of $n$ independently and identically drawn (i.i.d.) data points, $S_{train} = \{\boldsymbol{x}_i,y_i\}_{i=1}^n$, where $\boldsymbol{x}_i \in \mathbb{R}^d$ denotes  a $d$-dimensional feature vector and $y_i \in \mathbb{R}$ is a scalar output. Similarly, $S_{test}$
indicates the test set.
Subsequently, $X$ represents the design matrix and $\boldsymbol{y}$ indicates the output vector, where $(S_{train}.X,  S_{train}.\boldsymbol{y})$ and $(S_{test}.X,  S_{test}.\boldsymbol{y})$ represent the train and test sets, respectively.
Without the loss of generality, we consider the regression tasks of the form
\begin{equation*}
    \boldsymbol{y} = f(X) + \epsilon,
\end{equation*}
where $\epsilon$ is a normally distributed constant noise, i.e., $\epsilon \sim \mathcal{N}(0, \boldsymbol{\sigma}_a^2)$, and is assumed to be known.
The goal is hence to quantify the predictive uncertainty $\boldsymbol{\sigma}_p^2$ associated with $S_{test}.\boldsymbol{y}$, while optimizing $f$ on the training data.

We adapt the regularized loss function from the Anchored Ensemble technique~\cite{Pearce2018UncertaintyIN}, in which, the regularization of the models' parameters are carried out w.r.t. their initialization values instead of zero. Consequentially, given $\boldsymbol{\theta}_j$ as the parameters of the $j$\textsubscript{th} base-learner, the objective function is as follows
\begin{equation}
\label{anchor}
    \mathcal{L}(\boldsymbol{\theta}_j) = \frac{1}{n}||\boldsymbol{y}-\boldsymbol{\hat{y}}_j||_2^2+\frac{1}{n}||\Gamma^{1/2}(\boldsymbol{\theta}_j - \boldsymbol{\theta}_{0,j})||_2^2,
\end{equation}
where $\boldsymbol{\theta}_{0,j}$ is derived from the prior distribution, $\boldsymbol{\theta}_{0,j} \sim \mathcal{N}(\boldsymbol{\mu}_{0}, \Sigma_{0})$, and $\Gamma$ is the regularization matrix. Furthermore, minimizing this objective allows for performing Bayesian inference in NNs. However, this technique only models the epistemic uncertainty, while aleatoric uncertainty is assumed to be constant \cite{Pearce2018UncertaintyIN}, which is a limitation, as it is not always possible to distinguish the different origins or types of uncertainty in practice (see Equation~\ref{eq:unc}).

Therefore, in this paper, we aim at enhancing the performance of the ensemble toward a more robust uncertainty estimation. The literature suggests that diversifying the ensembles is effective in improving their predictive performance both theoretically and empirically \cite{Zhou2012EnsembleMF, Zhang2012EnsembleML,Hansen1990NeuralNE, Krogh1994NeuralNE}.
Ideally, diversity is achieved when the predictions made by each model in the ensemble are independent and uncorrelated. However, generating diverse ensemble members is not a straightforward task. The main impediment is the fact that each neural network is trained on the same training data to solve the same problem, which usually results in a high correlation among the individual base-learners \cite{Zhou2012EnsembleMF}. In the subsequent section, we introduce a simple technique to efficiently improve the  overall diversity of the ensemble for a more reliable uncertainty quantification.

\subsection{RAFs Ensemble}
In this section, we present Random Activation Functions (RAFs) Ensemble for uncertainty estimation, which can be extended to all ensemble methods in terms of methodological modification.
When a (Bayesian) ensemble is leveraged to estimate the uncertainty of a deep neural network model, we propose to increase the diversity of predictions among the ensemble members using varied activation functions (AFs), in addition to the random initialization of the parameters. To do so, instead of initializing the neural networks with the same activation function, each NN is accommodated with a different (random) activation function. Subsequently, distinct activation functions account for different non-linear properties introduced to each ensemble member, therewith improving the overall diversity of the ensemble.

As mentioned previously, the ensemble diversity is one of the most important building blocks when it comes to creating a successful ensemble \cite{Hansen1990NeuralNE}. Hence, it might be preferable to combine the predictions of top-performing base-learners with the predictions of weaker ones \cite{Zhou2012EnsembleMF}. Otherwise, stacking only strong models will likely result in a poor ensemble as the predictions made by the models will be highly correlated, and thus, the ensemble diversity will be greatly limited. Therefore, the choice of activation functions should be motivated purely by their variability and not their appropriateness for the task at-hand.

Let $\boldsymbol{\mu}_0$ be the prior means, $\Sigma_0$ be the prior covariance, $\boldsymbol{\hat{\sigma}}_a^2$
be an estimate of data noise, $m$ denote the number of base-learners, and $NN_j$ indicate the $j_\text{th}$ member, the entire procedure for both training and prediction is summarized in Algorithm \ref{alg:rafs}. In this algorithm, a regularization matrix is first created and a set of activation functions is defined  (line 1-2). Then, the NNs in the ensemble are trained to minimize the loss function defined in Equation \ref{anchor} with stochastic gradient descent, using arbitrary optimizer and no early stopping (line 3-13). Note that if the size of the ensemble $m$ is smaller or equal to the cardinality of the AFs set $k$, then each NN is trained with a different activation function, and with random functions from the set, otherwise (line 7-11). Consequently, predictions are made with each ensemble member (line 14-16), which are then averaged and an estimate of the predictive uncertainty is computed (line 17-19).

\begin{algorithm}[t]
\caption{RAFs Ensemble}\label{alg:rafs}
\textbf{Input:} $S_{train}, S_{test},$ priors $\boldsymbol{\mu}_0$ and $\Sigma_0$, $m$, $\boldsymbol{\hat{\sigma}}_a^2$\\
\textbf{Output:} Estimate of predictive mean  $\boldsymbol{\hat{y}}$ and variance  $\boldsymbol{\hat{\sigma}}_p^2$
\begin{algorithmic}[1]
\State $\Gamma \leftarrow \boldsymbol{\hat{\sigma}}_a^2 \Sigma_0^{-1}$   \Comment{Regularization matrix}
\State $\mathbb{A}\leftarrow \{a_1,\ldots,a_k\}$ \Comment{Set of $k$ AFs}
\For{$\hspace{1mm} j$ in $1:m$}    \Comment{Train the ensemble}
\State Create $NN_j$ with $\boldsymbol{\theta}_{j,0}\leftarrow \mathcal{N}(\boldsymbol{\mu}_0, \Sigma_0)$
\If{$ j \leq k$}
    \State $\alpha_j=a_j$
\Else
    \State $\alpha_j\leftarrow$ Randomly selected from $\mathbb{A}
\EndIf
\State $NN_j.train(S_{train},\Gamma, \boldsymbol{\theta}_{j,0}, \alpha_j)$ using loss in Eq. \ref{anchor}
\EndFor
\For{$\hspace{1mm} j$ in $1:m$}   \Comment{Predict with the ensemble}
\State $\boldsymbol{\hat{y}}_j = NN_j.predict(S_{test}.X)$
\EndFor
\State $\boldsymbol{\hat{y}} = \frac{1}{m}\displaystyle\sum_{j=1}^{\mathclap{m}} \boldsymbol{\hat{y}}_j $ \Comment{Mean predictions}
\State $\boldsymbol{\hat{\sigma}}_e^2 = \frac{1}{m-1} \displaystyle\sum_{j=1}^{\mathclap{m}} (\boldsymbol{\hat{y}}_j - \boldsymbol{\hat{y}})^2$   \Comment{Epistemic variance}
\State $\boldsymbol{\hat{\sigma}}_p^2 = \boldsymbol{\hat{\sigma}}_e^2 + \boldsymbol{\hat{\sigma}}_a^2$ \Comment{Total variance Eq. \ref{eq:unc}}
\State \Return $\boldsymbol{\hat{y}}, \boldsymbol{\hat{\sigma}}_p^2$
\end{algorithmic}
\end{algorithm}

\section{Empirical Study}

\subsection{Experimental Setups}
In the experiments, the base-learners of RAFs Ensemble are multilayer perceptrons that consist of one hidden layer of 100 neurons. The ensemble size $m$ is set to five. This is standard for the implementations of all methods in this paper, as $m=5$ proved to be empirically sufficient for obtaining predictive uncertainty estimates in the experiments.
In addition, we choose a set of seven activation functions which is comprised of (i) Gaussian Error Linear Unit (GELU) \cite{Hendrycks2016GaussianEL}, (ii) Softsign \cite{Turian2009QuadraticFA}, (iii) Swish \cite{Ramachandran2018SearchingFA}, (iv) Scaled Exponential Linear Unit (SELU) \cite{Klambauer2017SelfNormalizingNN}, (v) hyperbolic tangent (tanh), (vi) error activation function, and (vii) linear (identity) activation function.
Furthermore, the number of testing samples is set to be always larger than the number of training points $n$ to detail the uncertainty. Moreover, to account for epistemic uncertainty, the synthetic testing feature vectors $\boldsymbol{x}\in S_{test}$ range over wider intervals compared to $\boldsymbol{x}\in S_{train}$ and both are sampled uniformly at random.
The code is available at \url{https://github.com/YanasGH/RAFs_code} for reproducibility.

\subsubsection{Baselines.}We include five state-of-the-art methods as baselines for empirical comparison with RAFs Ensemble  as follows. (i) DE \cite{Lakshminarayanan2017SimpleAS}, (ii) AE \cite{Pearce2018UncertaintyIN}, (iii) HDE \cite{Wenzel2020HyperparameterEF}, (iv) RP-param \cite{Osband2018RandomizedPF}, and (v) NTKGP-param \cite{He2020BayesianDE}, on both synthetic and real-world datasets with different dimensionalities (see the Technical appendix for a detailed overview). To ensure fair comparison between the UQ techniques, roughly the same amount of time has been put into hyperparameter tuning for each method.

\subsubsection{Synthetic Data.}
We generate multiple synthetic datasets that fall into four  categories: physical models (PM), many local minima (MLM), trigonometric (T), and others (O). Each set in the PM category is generated from a physical mathematical model, such that all values in $S_{train}$ and $S_{test}$ are achievable in the real world. Generally, the PM datasets have complex modeling dynamics and can be characterized as having predominant epistemic uncertainty due to the considerably wider testing sampling regions by design. Similarly, the MLM data, generated from functions with many local minima, are also designed so that the model uncertainty is higher than the aleatoric one. These datasets are usually hard to approximate due to their inherent high-nonlinearity and multimodality. Another category with higher epistemic uncertainty is trigonometric, such as
data generated by \cite{He2020BayesianDE} and \cite{Forrester2008EngineeringDV}, where the training data is partitioned into two equal-sized clusters in order to detail uncertainty on out-of-distribution data (see Figure \ref{fig:xsinx}).
In contrast, the predominant type of uncertainty in the O category is aleatoric. This category includes datasets generated from various functions such as rational and product integrand functions. It is distinguished from the rest of the categories by its high interaction effects. The dimensionality of all datasets can range from one to ten and we consider two datasets per dimension, thus, the total number of synthetic data is 20. More detail on how the data is created can be found in the Technical appendix.

\subsubsection{Real-world Data.}
Additionally, we use five real-world datasets for evaluation: Boston housing, Abalone shells \cite{Abalone1994}, Naval propulsion plant \cite{Coraddu2013Machine}, Forest fire \cite{Cortez2007ADM}, and Parkinson's disease dataset \cite{Little2007ExploitingNR}. To account for aleatoric uncertainty (some) context factors are disregarded, such that this type of uncertainty is characteristically high (see Technical appendix for more details).

\subsubsection{Evaluation Criteria.}
\label{eval_crit}
We employ two evaluation criteria  to gauge the overall performance of the trained models, namely calibration and robustness to the distribution shift. Both measures are inspired by the practical applications of NNs, as generally there is no theoretical evidence for evaluating uncertainty estimates \cite{Abdar2021ARO}. Calibration is defined as the analytical process of adjusting the inputs with the purpose of making the model to predict the actual observations as precisely as possible \cite{Bijak2021UncertaintyQM}. The quality of calibration can be measured by proper scoring rules such as negative log-likelihood (NLL). NLL is a common choice when it comes to evaluating UQ estimates, as it depends on predictive uncertainty \cite{Lakshminarayanan2017SimpleAS}. Additionally, due to its practical applicability in a wide spectrum of regression tasks, root mean squared error (RMSE) is measured, although it does not depend on the estimated uncertainty \cite{Lakshminarayanan2017SimpleAS}, but serves as a proxy and a secondary assessor of the performance. Moreover, to measure the robustness/generalization of methods to distributional shift, we test the models in out-of-distribution settings, such as the synthetic datasets by \cite{Forrester2008EngineeringDV,He2020BayesianDE}.

\begin{figure}[t]
\begin{subfigure}{.15\textwidth}
  \centering
  \includegraphics[width=\linewidth]{jointxsinxde.pdf}
  \caption{DE}
  \label{fig:xsinxDE}
\end{subfigure}%
\begin{subfigure}{.15\textwidth}
  \centering
  \includegraphics[width=\linewidth]{jointxsinxae.pdf}
  \caption{AE}
  \label{fig:xsinxAE}
\end{subfigure}
\begin{subfigure}{.15\textwidth}
  \centering
  \includegraphics[width=\linewidth]{jointxsinxhde.pdf}
  \caption{HDE}
  \label{fig:xsinxHDE}
\end{subfigure}

\begin{subfigure}{.15\textwidth}
  \centering
  \includegraphics[width=\linewidth]{jointxsinxrp.pdf}
  \caption{RP-param}
  \label{fig:xsinxRP}
\end{subfigure}
\begin{subfigure}{.15\textwidth}
  \centering
  \includegraphics[width=\linewidth]{jointxsinxntkgp.pdf}
  \caption{NTKGP-param}
  \label{fig:xsinxNTKGP}
\end{subfigure}
\begin{subfigure}{.15\textwidth}
  \centering
  \includegraphics[width=\linewidth]{jointxsinxrafs.pdf}
  \caption{RAFs}
  \label{fig:xsinxRAFs}
\end{subfigure}
\caption{Uncertainty quantification of different methods on He et al. dataset. Gaussian process with neural tangent kernel (NTKGP analytic) is included as a reference.}
\label{fig:xsinx}
\end{figure}

\begin{table*}[t]
\centering
{\small
\begin{tabular}{l c c c c c c}
\hline
 &  & & NLL & & & \\
 \hline
 & DE & HDE & AE & NTKGP-p. & RP-p. & RAFs \\
 \hline
 He et al.  1D   & $>$100 $\pm$ 0.18 & 71.31 $\pm$ 0.51 & 38.75 $\pm$ 0.12 & 4.48 $\pm$ 0.18 & 13.05 $\pm$ 0.43 & \textbf{2.21 $\pm$ 0.18}\\

 Forrester et al. 1D   & $>$100 $\pm$ 0.53 & $>$100 $\pm$ 0.51 & 50.82 $\pm$ 0.52 & $>$100 $\pm$ 0.50 & 13.7 $\pm$ 0.58 & \textbf{0.64 $\pm$ 0.74}\\

 Schaffer N.4 2D   & 0.29 $\pm$ 0.01 & -0.71 $\pm$  0.01 & 2.15 $\pm$ 0.01 & -0.55 $\pm$ 0.01 &  -0.36 $\pm$ 0.01 & \textbf{-0.79 $\pm$ 0.01} \\

 Double pendulum 2D   & 2.95 $\pm$ 0.05  & 2.18 $\pm$ 0.84 & -0.36 $\pm$ 0.05 & \textbf{-0.58 $\pm$ 0.05} &  \textbf{-0.47 $\pm$ 0.05} & \textbf{-0.49 $\pm$ 0.04} \\

 Rastrigin 3D & 29.24 $\pm$ 1.30  & \textbf{3.09 $\pm$  1.15} & 35.94 $\pm$ 0.74 & 28.38 $\pm$ 0.64 &  \textbf{4.35 $\pm$ 1.24} & \textbf{3.44 $\pm$ 1.05} \\

 Ishigami 3D   & 6.01 $\pm$ 0.08  & $>$100 $\pm$  0.08 & 8.73 $\pm$ 0.08 & 1.53 $\pm$ 0.08 &  \textbf{-0.01 $\pm$ 0.08} & \textbf{0.06  $\pm$ 0.07} \\

 Environmental 4D & 64.72 $\pm$ 0.23  & 7.84 $\pm$  0.13 & 1.65 $\pm$ 0.20 & 4.5 $\pm$ 0.27 &  3.94 $\pm$ 0.21 & \textbf{0.81 $\pm$ 0.17} \\

 Griewank 4D   & 28.29 $\pm$ 2.43  & \textbf{5.50 $\pm$  1.62} & \textbf{4.64 $\pm$ 3.06} & 10.21 $\pm$ 2.37 &  \textbf{4.29 $\pm$ 2.93} & \textbf{4.79 $\pm$ 2.40} \\

 Roos \& Arnold 5D   & -2.02 $\pm$ 0.01  & \textbf{-2.21 $\pm$  0.00} & -1.89 $\pm$ 0.01 & -1.71 $\pm$ 0.01 &  -1.70 $\pm$ 0.01 & -2.1 $\pm$ 0.01 \\

 Friedman 5D  & 96.94 $\pm$ 0.41  & $>$100 $\pm$  0.51 & 15.04 $\pm$ 0.50 & 41.69 $\pm$ 0.39 &  4.22 $\pm$ 0.44 & \textbf{1.78 $\pm$ 0.39} \\

 Planar arm torque 6D   & 9.58 $\pm$ 0.07  & 4.11 $\pm$  0.08 & 3.07 $\pm$ 0.05 & \textbf{-0.32 $\pm$ 0.08} &  -0.05 $\pm$ 0.07 & -0.16 $\pm$ 0.06 \\

 Sum of powers 6D & $>$100 $\pm$ 0.41  & $>$100 $\pm$  0.62 & 55.03 $\pm$ 0.43 & $>$100 $\pm$ 0.41 &  41.59 $\pm$ 0.40 & \textbf{35.22 $\pm$ 0.35} \\

 Ackley 7D  & 7.11 $\pm$ 0.23  & \textbf{1.38 $\pm$  0.16} & 2.50 $\pm$ 0.36 & 3.11 $\pm$ 0.27 &  2.09 $\pm$ 0.26 & \textbf{1.16 $\pm$ 0.08} \\

 Piston simulation 7D  & \textbf{-2.19 $\pm$ 0.00}  & 14.06 $\pm$  0.00 & 3.50 $\pm$ 2.40 & 2.87 $\pm$ 2.93 &  2.67 $\pm$ 0.42 & 3.63 $\pm$ 0.57 \\

 Robot arm 8D   & 10.71 $\pm$ 0.03  & 6.87 $\pm$  0.01 & 7.11 $\pm$ 0.01 & \textbf{0.27 $\pm$ 0.03} &  0.80 $\pm$ 0.06 & \textbf{0.25 $\pm$ 0.02} \\

 Borehole 8D    & $>$100 $\pm$ 1.01  & $>$100 $\pm$  1.01 & \textbf{4.89 $\pm$ 1.87} & \textbf{5.48 $\pm$ 3.54} &  \textbf{4.06 $\pm$ 1.20} & \textbf{4.36 $\pm$ 1.26} \\

 Styblinski-Tang 9D  & $>$100 $\pm$ 3.05  & $>$100 $\pm$  0.00 & 40.80 $\pm$ 5.33 & $>$100 $\pm$ 3.03 &  \textbf{15.82 $\pm$ 6.31} & 2\textbf{5.23 $\pm$ 4.12} \\

 PUMA560 9D    & 6.59 $\pm$ 0.15  & \textbf{1.62 $\pm$  0.14} & 4.24 $\pm$ 0.14 & 5.93 $\pm$ 0.08 &  6.40 $\pm$ 0.14 & 2.14 $\pm$ 0.13 \\

 Adapted Welch 10D  & $>$100 $\pm$ 0.81  & $>$100 $\pm$  0.75 & $>$100 $\pm$ 0.55 & $>$100 $\pm$ 0.75 &  $>$100 $\pm$ 0.57 & \textbf{78.53 $\pm$ 0.67} \\

 Wing weight 10D  & $>$100 $\pm$ 0.00  & 27.31 $\pm$  4.37 & \textbf{5.46 $\pm$ 4.36} & 67.30 $\pm$ 0.53 &  \textbf{5.54 $\pm$ 4.15} & \textbf{5.39 $\pm$ 1.69} \\
 \hline
 Boston housing & 74.54 $\pm$ 1.06 & $>$100 $\pm$ 1.04 & 71.53 $\pm$ 1.06 & 70.82 $\pm$ 1.06 & $>$100 $\pm$ 1.10 & \textbf{ 40.67 $\pm$ 1.00} \\

 Abalone & $>$100 $\pm$ 0.10 & $>$100 $\pm$ 0.10 & 47.67 $\pm$ 0.10 & $>$100 $\pm$ 0.10 & \textbf{28.37 $\pm$ 0.10} & 28.90 $\pm$ 0.10 \\

 Naval propulsion & \textbf{-2.27 $\pm$ 0.00} & $>$100 $\pm$ 0.00 & 3.92 $\pm$ 0.10 & 2.28 $\pm$ 1.51 & 2.16 $\pm$ 0.16 & 1.91 $\pm$ 0.07 \\

 Forest fire & 15.71  $\pm$ 0.05 & \textbf{3.14 $\pm$ 0.02} & \textbf{2.66 $\pm$ 0.69} &   \textbf{3.10 $\pm$ 1.11} & 4.68 $\pm$ 0.14 & \textbf{2.15 $\pm$ 0.28} \\

 Parkinson's & \textbf{26.74  $\pm$ 0.02} & $>$100 $\pm$ 0.10 & $>$100 $\pm$ 0.16 &   $>$100 $\pm$ 0.03 & $>$100 $\pm$ 0.16 & 45.69 $\pm$ 0.16 \\
\hline
\end{tabular}}
\caption{Performance of methods on all datasets w.r.t. NLL, including 95\% confidence intervals. The best scores are in bold. \label{tab:nll}}
\end{table*}

\subsection{Performance Results}
\subsubsection{Qualitative Comparison.}
Figure \ref{fig:xsinx} shows the performance of different methods compared to a Gaussian process with a neural tangent kernel (NTKGP analytic) as a reference, on a 1D toy dataset generated from $\boldsymbol{y} = \boldsymbol{x}\text{sin}(\boldsymbol{x}) + \boldsymbol{\epsilon}$ (dashed line).
The plots demonstrate that DE, HDE, and AE provide narrow uncertainty bounds in areas where no data has been observed by the model, which translates to high confidence in OOD data. On the contrary, NTKGP-param, RP-param, and RAFs Ensemble bound their uncertainty estimates with wider intervals in areas with no data, accounting for adequate quantification of epistemic uncertainty, while also indicating robustness to OOD data. Among these methods, RAFs Ensemble provides the widest uncertainty which is reasonable considering the amount of data that is available to the methods over each area. Moreover, this observation is quantitatively validated as RAFs Ensemble achieves the lowest NLL compared to the other methods (see Table \ref{tab:nll}).

\subsubsection{Overall Performance.}
We evaluate the overall performance of all methods in terms of both NLL and RMSE. The outcomes of comparing RAFs Ensemble with five baseline methods on twenty synthetic and five real-world datasets are outlined in Table~\ref{tab:nll} and Table~\ref{tab:rmse}. The results illustrate that our approach outperforms the competitors in most scenarios. Furthermore, Table~\ref{tab:rank} summarizes the obtained results in terms of ranking, in which the methods are ranked based on their performance for a particular dataset.
The left integer corresponds to NLL, while the right one points to RMSE, and the bold values indicate the best-performing method.

\subsubsection{Discussion.}
The obtained results in this section illustrate that DE has good uncertainty estimates with respect to NLL for the real-world datasets, and it takes the first place for Naval propulsion and Parkinson's datasets. For the rest of the data categories, when compared to the other methods, DE fails to provide strong performance, usually scoring a very low NLL rank. Therefore, this indicates that Deep Ensemble might have difficulty quantifying epistemic uncertainty in general as displayed by the experiments in this paper, but seemingly manages to capture aleatoric uncertainty well.

Unlike DE, the HDE provides outwardly reliable estimates for datasets with many local minima, despite its unimpressive overall results when compared to the other methods. However, both DE and HDE can produce uncertainty bounds that are unreasonably narrow in areas with unobserved data, as shown in Figure \ref{fig:xsinx} and noted by \cite{Heiss2021NOMUNO}.

\begin{table*}[t]
\centering
{\small
\begin{tabular}{l c c c c c c}
\hline
 &  & & RMSE & & & \\
 \hline
 & DE & HDE & AE & NTKGP-p. & RP-p. & RAFs\\
 \hline
 He et al.  1D   & 3.71 $\pm$ 0.18 & 5.70 $\pm$ 0.51 & \textbf{3.15 $\pm$ 0.12} & 3.64 $\pm$ 0.18 & 5.24 $\pm$ 0.43 & 3.80 $\pm$ 0.18\\

 Forrester et al. 1D   & 5.00 $\pm$ 0.53 & 4.12 $\pm$ 0.51 & 4.09 $\pm$ 0.52 & 6.05 $\pm$ 0.50 & 5.70 $\pm$ 0.58 & \textbf{2.8 $\pm$ 0.74}\\

 Schaffer N.4 2D   & \textbf{0.23 $\pm$ 0.01} & 0.34 $\pm$  0.01 & 0.30 $\pm$ 0.01 & \textbf{0.24 $\pm$ 0.01} &  0.31 $\pm$ 0.01 & 0.27 $\pm$ 0.01 \\

 Double pendulum 2D   & \textbf{0.46 $\pm$ 0.05}  & 2.22 $\pm$ 0.84 & 0.71 $\pm$ 0.05 & \textbf{0.51 $\pm$ 0.05} &  0.74 $\pm$ 0.05 & \textbf{0.58 $\pm$ 0.04} \\

 Rastrigin 3D & 18.41 $\pm$ 1.30  & \textbf{10.96 $\pm$  1.15} & 25.58 $\pm$ 0.74 & 18.10 $\pm$ 0.64 &  \textbf{12.87 $\pm$ 1.24} & \textbf{14.85 $\pm$ 1.05} \\

 Ishigami 3D   & \textbf{0.69 $\pm$ 0.08}  & 1.05 $\pm$  0.08 & 0.88 $\pm$ 0.08 & \textbf{0.69 $\pm$ 0.08} &  \textbf{0.58 $\pm$ 0.08} & \textbf{0.57  $\pm$ 0.07} \\

 Environmental 4D & \textbf{2.04 $\pm$ 0.23}  & 2.51 $\pm$  0.13 & \textbf{1.83 $\pm$ 0.20} & \textbf{2.34 $\pm$ 0.27} &  \textbf{2.03 $\pm$ 0.21} & \textbf{1.68 $\pm$ 0.17} \\

 Griewank 4D   & 83.97 $\pm$ 2.43  & \textbf{45.68 $\pm$  1.62} & \textbf{42.12 $\pm$ 3.06} & 78.47 $\pm$ 2.37 &  \textbf{38.62 $\pm$ 2.93} & 78.79 $\pm$ 2.40 \\

 Roos \& Arnold 5D   & 0.07 $\pm$ 0.01  & \textbf{0.01 $\pm$  0.00} & 0.07 $\pm$ 0.01 & 0.09 $\pm$ 0.01 &  0.08 $\pm$ 0.01 & 0.08 $\pm$ 0.01 \\

 Friedman 5D  & \textbf{3.17 $\pm$ 0.41}  & \textbf{3.63 $\pm$  0.51} & \textbf{2.95 $\pm$ 0.50} & \textbf{3.39 $\pm$ 0.39} &  \textbf{2.74 $\pm$ 0.44} & \textbf{3.1 $\pm$ 0.39} \\

 Planar arm torque 6D   & \textbf{0.65 $\pm$ 0.07}  & \textbf{0.62 $\pm$  0.08} & \textbf{0.71 $\pm$ 0.05} & \textbf{0.71 $\pm$ 0.08} &  1.08 $\pm$ 0.07 & \textbf{0.74 $\pm$ 0.06} \\

 Sum of powers 6D & \textbf{22.81 $\pm$ 0.41}  & \textbf{21.19 $\pm$  0.62} & \textbf{21.87 $\pm$ 0.43} & \textbf{22.79 $\pm$ 0.41} &  \textbf{22.22 $\pm$ 0.40} & \textbf{22.24 $\pm$ 0.35} \\

 Ackley 7D  & 8.92 $\pm$ 0.23  & 2.43 $\pm$  0.16 & 7.28 $\pm$ 0.36 & 8.58 $\pm$ 0.27 &  4.03 $\pm$ 0.26 & \textbf{1.33 $\pm$ 0.08} \\

 Piston simulation 7D  & \textbf{0.02 $\pm$ 0.00}  & 0.04 $\pm$  0.00 & 29.1 $\pm$ 2.40 & $>$100 $\pm$ 2.93 &  5.78 $\pm$ 0.42 & 7.40 $\pm$ 0.57 \\

 Robot arm 8D   & 0.92 $\pm$ 0.03  & \textbf{0.80 $\pm$ 0.01}  & 0.88 $\pm$ 0.01 & 0.93 $\pm$ 0.03 &  1.09 $\pm$ 0.06 & \textbf{0.83 $\pm$ 0.02} \\

 Borehole 8D    & \textbf{32.11 $\pm$ 1.01}  & \textbf{32.12 $\pm$  1.01} & 48.75 $\pm$ 1.87 & $>$100 $\pm$ 3.54 &  38.60 $\pm$ 1.20 & 41.35 $\pm$ 1.26 \\

 Styblinski-Tang 9D  & $>$100 $\pm$ 3.05  & $>$100 $\pm$  0.00 & \textbf{$>$100 $\pm$ 5.33} & $>$100 $\pm$ 3.03 &  \textbf{$>$100 $\pm$ 6.31} & $>$100 $\pm$ 4.12 \\

 PUMA560 9D    & 3.93 $\pm$ 0.15  & \textbf{3.23 $\pm$  0.14} & \textbf{3.40 $\pm$ 0.14} & 3.93 $\pm$ 0.08 &  \textbf{3.24 $\pm$ 0.14} & \textbf{3.4 $\pm$ 0.13} \\

 Adapted Welch 10D  & \textbf{99.51 $\pm$ 0.81}  & \textbf{99.4 $\pm$  0.75} & $>$100 $\pm$ 0.55 & \textbf{99.79 $\pm$ 0.75} &  $>$100 $\pm$ 0.57 & \textbf{100.00 $\pm$ 0.67} \\

 Wing weight 10D  & $>$100 $\pm$ 0.00  & \textbf{58.16 $\pm$  4.37} & \textbf{63.1 $\pm$ 4.36} & $>$100 $\pm$ 0.53 &  \textbf{63.35 $\pm$ 4.15} & $>$100 $\pm$ 1.69 \\
 \hline
 Boston housing & \textbf{11.28 $\pm$ 1.06} & \textbf{11.36 $\pm$ 1.04} & \textbf{11.42 $\pm$ 1.06} & \textbf{11.28 $\pm$ 1.06} & \textbf{11.56 $\pm$ 1.10} & \textbf{11.31 $\pm$ 1.00} \\

 Abalone & \textbf{2.06 $\pm$ 0.10} & \textbf{2.09 $\pm$ 0.10} & \textbf{2.08 $\pm$ 0.10} & \textbf{2.05 $\pm$ 0.10} & \textbf{2.09 $\pm$ 0.10} & \textbf{2.08 $\pm$ 0.10} \\

 Naval propulsion &  \textbf{0.02 $\pm$ 0.00} & \textbf{0.02 $\pm$ 0.00} & 38.86 $\pm$ 0.60 & 62.61 $\pm$ 1.51 & 9.40 $\pm$ 0.16 & 3.45 $\pm$ 0.08 \\

 Forest fire & 1.97  $\pm$ 0.05 & \textbf{1.87 $\pm$ 0.02} & 6.43 $\pm$ 0.69 &   10.43 $\pm$ 1.11 & 2.32 $\pm$ 0.14 & 3.32 $\pm$ 0.28 \\

 Parkinson's & 12.17  $\pm$ 0.02 & 12.40 $\pm$ 0.10 & 12.49 $\pm$ 0.16 &   \textbf{11.97 $\pm$ 0.03} & 12.60 $\pm$ 0.16 & 12.78 $\pm$ 0.16 \\

\hline
\end{tabular}}
\caption{Performance of methods on all datasets w.r.t. RMSE, including 95\% confidence intervals. The best scores are in bold. \label{tab:rmse}}
\end{table*}

\begin{table}[!ht]
\centering
\resizebox{\columnwidth}{!}{
\begin{tabular}{l c c c c c c}
    & \text{DE}        & \text{HDE}       & \text{AE}       & \text{NTKGP-p.}       & \text{RP-p.}    & \text{RAFs} \\
He et al. 1D & 6,2 & 5,3 & 4,1 & 2,2 & 3,3 & \textbf{1},2\\

Forrester et al. 1D & 4,2 & 5,2 & 3,2 & 6,2  & 2,2 &  \textbf{1},1 \\

Schaffer N.4 2D & 5,1 & 2,4 & 6,3 & 3,1 & 4,3 & \textbf{1},2 \\

Double pendulum 2D & 3,1 & 3,3 & 2,2 & \textbf{1},1 & \textbf{1},2 & \textbf{1},1 \\

Rastrigin 3D & 2,2 & \textbf{1},1 & 3,3 & 2,2 & \textbf{1},1 & \textbf{1},1 \\

Ishigami 3D & 3,1 & 5,3 & 4,2 & 2,1 & \textbf{1},1 & \textbf{1},1 \\

Environmental 4D & 6,1 & 5,2 & 2,1 & 4,1 & 3,1 & \textbf{1},1 \\

Griewank 4D & 3,3 & \textbf{1},1 & \textbf{1},1 & 2,2 & \textbf{1},1 & \textbf{1},2 \\

Roos \& Arnold 5D & 3,1 & \textbf{1},1 & 4,3 & 5,1 & 5,1 & 2,2 \\

Friedman 5D & 5,1 & 6,1 & 3,1 & 4,1 & 2,1 & \textbf{1},1 \\

Planar arm torque 6D & 5,1 & 4,1 & 3,1 & \textbf{1},1 & 2,2 & 2,1 \\

Sum of powers 6D & 5,1 & 6,1 & 3,1 & 4,1 & 2,1 & \textbf{1},1 \\

Ackley 7D & 4,5 & \textbf{1},2 & 2,4 & 3,5 & 2,3 & \textbf{1},1 \\

Piston simulation 7D & \textbf{1},1 & 3,2 & 2,5 & 2,6 & 2,3 & 2,4 \\

Robot arm 8D & 5,3 & 3,1 & 4,2 & \textbf{1},3 & 2,4 & \textbf{1},1 \\

Borehole 8D & 2,1 & 3,1 & \textbf{1},4 & \textbf{1},5 & \textbf{1},2 & \textbf{1},3 \\

Styblinski-Tang 9D & 3,2 & 5,5 & 2,1 & 4,3 & \textbf{1},1 & \textbf{1},4 \\

PUMA560 9D & 5,2 & \textbf{1},1 & 3,1 & 4,2 & 4,1 & 2,1 \\

Adapted Welch 10D & 6,1 & 2,1 & 4,3 & 5,1 & 3,2 & \textbf{1},1 \\

Wing weight 10D & 4,4 & 2,1 & \textbf{1},1 & 3,3 & \textbf{1},1 & \textbf{1},2 \\
\hline
Boston housing & 3,1 & 5,1 & 2,1 & 2,1 & 5,1 & \textbf{1},1 \\

Abalone & 5,1 & 6,1 & 3,1 & 4,1 & \textbf{1},1 & 2,1 \\

Naval propulsion plant & \textbf{1},1 & 4,1 & 3,4 & 2,5 & 2,3 & 2,2 \\

Forest fire & 3,2 & \textbf{1},1 & 1,5 & 1,6 & 2,3 & 1,4 \\

Parkinson's & \textbf{1},2 & 6,3 & 4,3 & 5,1 & 3,3 & 2,3
\end{tabular}}
\caption{Rank of the methods corresponding to NLL (left) and RMSE (right). The best overall score is in bold (ties are possible in case of an overlap in confidence intervals). \label{tab:rank}}
\end{table}

Nonetheless, AE demonstrates good performance in the dataset categories that exhibit higher epistemic uncertainty such as the physical models. This is due to the fact that AE is designed for capturing model uncertainty, while aleatoric uncertainty is assumed to be constant. Accordingly, AE achieves inferior performance on the real-world datasets, as those generally have more data uncertainty appropriated.

On the other hand, NTKGP-param achieves its finest performance for datasets in the physical model category, which is normally associated with substantial model uncertainty. A credible rationale to explain this insight is the fact that NTKGP-param tends to be more conservative than Deep Ensemble. However, it is generally unclear in which situations this is beneficial since the ensemble members of NTKGP-param will always be misspecified in practice according to \cite{He2020BayesianDE}.

Furthermore, RP-param manages to rank comparatively high for real-world datasets as well as trigonometric data, that contain vast amounts of aleatoric and epistemic uncertainty, respectively, indicating that it does not quantify either type of uncertainty better than the other. This observation serves as a demonstration that RP-param generalizes well for different types of datasets that exhibit broad characteristics. However, this technique fails to deliver low NLL scores on some occasions, which might be attributed to the fact that RP-param is based on bootstrapping. While bootstrapping can be a successful strategy for inducing diversity, it can sometimes harm the performance when the base-learners have multiple local optima, as is a common case  with NNs \cite{Lakshminarayanan2017SimpleAS}.

Nevertheless, RAFs Ensemble outperforms RP-param, and every other method in the comparisons, for 13 out of 25 datasets. In terms of NLL, our approach does not rank below the second place for any data, which is consistent with the strong results from Table \ref{tab:nll}. Meanwhile, the RMSE scores of this method are altogether satisfactory, although not as prominent compared to the NLL scores. In agreement with the overall outstanding results, RAFs Ensemble holds the highest NLL rank for all data from \textit{MLM} and \textit{T} categories, which can be contemplated as a concluding statement regarding its capabilities to estimate epistemic uncertainty and challenging multimodality. Among all categories, the real-world datasets are least favored by RAFs Ensemble, primarily due to their high level of aleatoric uncertainty. This indicates that RAFs Ensemble captures model uncertainty slightly better than aleatoric uncertainty. Nonetheless, the empirical superiority of this technique is due to the exhaustively exploited added source of randomness via random activation functions, combined with method simplicity and Bayesian behavior, resulted from the anchored loss (Equation \ref{anchor}). This successful combination leads to greatly improved diversity among ensemble members, which can be further confirmed by a direct comparison between RAFs Ensemble and AE.
Note that even though RAFs Ensemble does not provide as prominent results with respect to RMSE in the higher dimensional datasets as it does in datasets of lower dimensions, it still achieves better or on par results compared to the state-of-the-art methods. In addition, RAFs Ensemble can be deployed in both complex and straightforward settings. On a related note, while  DE struggles when dealing with high multimodality and RP-param underperforms when the dataset has interaction effects (from ``others" category), RAFs excels in both such settings.

\subsubsection{Scalability to higher dimensions and larger networks.}
To test the scalability of RAFs Ensemble, we compare it with the strongest baseline, RP-param, on two additional real-world datasets, i.e., a 65-dimensional data with around 20k samples and a 40-dimensional data with almost 40k samples. The former is the superconductivity dataset, where the goal is to predict the critical temperature of superconductors \cite{Hamidieh2018ADS}. The latter summarizes features about articles, where the target is the number of shares in social networks \cite{Fernandes2015API}. Both methods utilize the same neural architecture for their base-learners,
that is two hidden layers of 128 hidden neurons each, which is more complex than the previous experiments. The conclusion of these experiments is conclusive in favor of our approach.
RAFs Ensemble scores NLL of 5.49 and 25.89 on the first and second dataset, respectively, while RP-param scores NLL of over 100 on both datasets.

\begin{figure}[t!]
\centering
\begin{subfigure}{.30\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{FRIEDMAN_ERR_vs_CONF_noleg.png}
  \caption{Friedman (synthetic)}
  \label{fig:conferr_friedman}
\end{subfigure}%
\vspace{1mm}
\begin{subfigure}{.30\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{ABALONE_ERR_vs_CONF.png}
  \caption{Abalone (real-world)}
  \label{fig:conferr_abalone}
\end{subfigure}
\caption{Confidence versus error of estimations.}
\label{fig:caler}
\end{figure}

\subsubsection{Confidence vs. Error.}
We further analyze the relation between the RMSE and the precision thresholds in order to examine the confidence of each method in the prediction task. Figure \ref{fig:caler} displays the confidence versus error plots for one synthetic and one real-world dataset, i.e., Friedman and Abalone (see the Technical Appendix for more detail). In this figure, for each precision threshold $\tau$, the RMSE is plotted for examples where the predicted precision $\boldsymbol{\sigma}_p^{-2}$ is larger than the threshold $\tau$, demonstrating confidence.
In general, reliable estimates are expected to have decreasing error when the confidence is increasing. For Friedman dataset, it is clear that RAFs Ensemble delivers well-calibrated estimates, which is especially in contrast with  DE, NTKGP-param, and HDE (Figure \ref{fig:conferr_friedman}). However, for the Abalone data, RP-param demonstrates the most reliable behavior, although RAFs Ensemble meets its performance at the last precision threshold (Figure \ref{fig:conferr_abalone}). Overall, our approach sustains lower error over most precision thresholds compared to the majority of the other methods, and this contrast in performance is emphasized as the predictions get more confident.

\subsubsection{Ablation.}
We study the effect of number of base-learners in the ensemble on the quality of UQ, which also measures the sensitivity of the results to the cardinality of the set of AFs $k$. We conduct an experiment on two different datasets, one synthetic (PUMA590) and one real-world (Abalone), where the results in terms of NLL are represented in Figure \ref{fig:exprafs}.
Note that Figure \ref{fig:RAFsAbalone} is shown in log-scale for better visibility. According to the theory, in the limit of infinite number of ensemble members, the ensemble error converges to zero \cite{Hansen1990NeuralNE}. However, practically speaking, five NNs in the ensemble provide optimal results regarding the trade-off between empirical performance and computational time \cite{Lakshminarayanan2017SimpleAS}, which is also the case in our experiments. This is further confirmed by the plot in Figure \ref{fig:RAFsPUMA}. In addition, for the PUMA590 dataset, it seems that RAFs Ensemble's performance is not impacted negatively by the number of NNs in the ensemble. Moreover, an interesting observation is the steep through for seven NNs (equal to $k$) in Figure \ref{fig:RAFsAbalone}, which is an indication that there might be a correlation between $k$ and the performance in some cases. A plausible reason for this is the fact that the additional source of randomness is utmostly exploited via a different activation function.

\begin{figure}[t!]
\begin{subfigure}{.2\textwidth}
  \centering
  \includegraphics[width=\linewidth]{RAFs_expPUMA.png}
  \caption{PUMA560}
  \label{fig:RAFsPUMA}
\end{subfigure}~
\begin{subfigure}{.22\textwidth}
  \centering
  \includegraphics[width=1.02\linewidth]{abalon.png}
  \caption{Abalone}
  \label{fig:RAFsAbalone}
\end{subfigure}
\caption{The effect of number of NNs in the ensemble in terms of NLL, including the 95\% confidence interval.}
\label{fig:exprafs}
\end{figure}

\begin{wrapfigure}{r}{0.2\textwidth}
\includegraphics[width=0.2\textwidth]{ablation_supercond.png}
\caption{The effect of cardinality $k$ on NLL for Superconductivity data.}
\label{fig:RAFsSupercond}
\end{wrapfigure}

To further confirm the effectiveness of the random activation functions, we evaluate the performance of RAFs Ensemble (of five NNs) in terms of NLL w.r.t. different cardinalities $k$ of the set of AFs. The dataset used for this experiment is the superconductivity data. As the results in Figure \ref{fig:RAFsSupercond} clearly suggest, by increasing the cardinality $k$, NLL has a decreasing pattern, which shows that having more random AFs significantly improves the performance of the ensemble.

Moreover, we combine our approach with RP-param instead of AE to show that RAFs can be methodologically applied to any ensemble technique. We evaluate the performance of this combination on the Parkinson’s dataset, using the same network architecture for fair comparison.
The obtained results demonstrate that applying RAFs to RP-param leads to reducing the original NLL score of $>100$ to 48.66, which is in line with the results we get when comparing AE with RAFs Ensemble and is a further proof that the methodology indeed increases the performance.

\section{Conclusions}
We introduced a novel method, Random Activation Functions Ensemble, for a more robust uncertainty estimation in approaches based on neural networks, in which, each network in the ensemble is accomodated with a different (random) activation function to increase the diversity of the ensemble. The empirical study illustrates that our approach achieves excellent results in quantifying both epistemic and aleatoric uncertainty compared to five state-of-the-art ensemble uncertainty quantification methods on a series of regression tasks across 25 datasets, which proved there does not have to be a trade-off between simplicity and strong empirical performance. Furthermore, the properties of datasets such as dimensionality or complexity of modeling dynamics do not appear to affect RAFs Ensemble negatively, which also demonstrates robustness in out-of-distribution settings.