\section*{Acknowledgement}
This research is supported by the National Research Foundation, Prime Minister's Office, Singapore under its Corporate Laboratory@University Scheme, National University of Singapore, and Singapore Telecommunications Ltd.


\bibliographystyle{IEEEtran}


\section{Feature Selection via Cross Entropy Method}
\label{sec:cross}

In this section we develop a novel algorithm to solve the feature selection under budget constraints, presented in Eq.~\ref{eq:optim}. As mentioned before, this is a non-convex combinatorial optimization problem. To overcome this computational difficulty, we develop a novel algorithm via a Monte Carlo sampling approach, which has a low computational complexity. Specifically, we use the cross entropy (CE) technique, which is based on the Importance Sampling Technique which is known to be an effective solution in computational approaches to otherwise solving NP-hard problems.
The CE converts the deterministic optimisation problem into a stochastic couterpart, see details in~\cite{asmussen2005heavy}. 
It then approximates the optimal sampling distribution by minimizing the Kullback-Leibler (KL) divergence~\cite{de2005tutorial}. 
Without loss of generality, we may change the optimization problem in Eq. \ref{eq:optim} to a maximization problem, that is: 
\begin{equation}
\mathbf{v}^* = \underset{\mathbf{v} \in\{0,1\}^{m}}{\arg\max} \quad \mathcal{U}(\mathbf{v}), \quad \text{s.t.}\quad \mathbf{c}^T\mathbf{v} \le \lambda
\end{equation}
where $\mathcal{U}(\mathbf{v}) = \frac{1}{R(\mathbf{v})}$.

Instead of searching for the optimal solution directly, like in the case of a brute force search based approach, the CE places a probability distribution on the set of  $m$-dimensional binary  features $(\mathbf{v} \in \{0,1 \}^m)$, thus transforming the problem into estimation. 


In our model, since the selection of a feature is indicated by a binary variable (``select''  or  ``don't select''), 
we use an independent Bernoulli random variable to indicate this choice for each of the $m$ features. The Bernoulli distribution
has a single variable parameter~$p$ , i.e., for the $k$-th feature this is given by $\{p_k\}_{k=1}^m$, and is a member of the Natural Exponential Families (NEF) of distributions \cite{collins2002generalization}. This is a big advantage, since under the NEF, the parameter of the distribution can be estimated analytically in closed form via the Maximum Likelihood Estimator (MLE), making the CE easy to implement. 
At each iteration of the CE method there are three steps:
\begin{enumerate}
    \item[S1] Generate samples: generate $\eta$ independent samples of binary sets given by $\mathbf{w}^{[i]} = \{w^{[i]}_{1},w^{[i]}_{2} \ldots, w^{[i]}_{m}\}$, where $w^{[i]}_{k} \sim \texttt{Bernoulli}(p_{k}) ; 1\le i \le \eta$.

    \item[S2] Selection of elite samples: from those $\eta$ samples, we select only those samples which satisfy the budget constraint  (ie. $\mathbf{c}^T\mathbf{w}^{[i]} \le \lambda$) to obtain a subset of $\eta' \leq \eta$ samples. 
    Next, from this subset we select only the ``elite'' samples---those for which the objective function value exceeds a pre-defined threshold, defined via a quantile value $\rho$. This choice is indicated by 
    $\mathbbm{1}\left(\mathcal{U}(\mathbf{w}^{[i]}) \geq \gamma\right)$,
    where $\gamma$ is the $(1- \rho)^{\text{th}}$-sample quantile of $\mathcal{U}_{1:\eta'}$.
    \item[S3] Estimation of parameters: using only this subset of elite samples, we re-estimate the parameters of the Bernoulli random variable $\{p_k\}_{k=1}^m$ via the MLE, given by
\begin{multline}
\widehat{p}_{k}=  \frac{\sum\limits_{i=1}^{\eta'} \mathbbm{1}\left(\mathbf{w}^{[i]}_{j} =1\right) 
\overbrace{
\mathbbm{1}\left(\mathcal{U}(\mathbf{w}^{[i]}) \geq \gamma\right)}
^{\text{Choose elite samples}}
}{\sum\limits_{i=1}^{\eta'} \mathbbm{1}\left(\mathcal{U}(\mathbf{w}^{[i]}) \geq \gamma\right)}, \\
k=\{1,\cdots,m \} .
\end{multline}
\end{enumerate}
\noindent 
The three steps of the algorithm are then iterated until a stopping rule is met, for example, if the algorithm has converged to a local maxima or has exhausted a pre-defined number of iterations.
The parameter $\rho$ is the quantile value, and is commonly set to $0.9$. 
 The resulting algorithm is presented in Algorithm 1, where in Step 10 we have introduced the parameter $\alpha \in [0,1]$ which controls the learning rate of the algorithm, thus avoiding getting trapped in local maxima, see details in  \cite{de2005tutorial}. 
 
 
It is important to note that $v_{k}$,

where $1 \le k \le m $, denotes whether the $k^\text{th}$ feature is selected by the algorithm. 
We iterate the cross entropy algorithm until the convergence criterion is met, which in this work, translates to reaching the maximum number of iterations ($T_{\max}$). 




\begin{algorithm} 
\label{alg:cross}
\SetAlgoLined
  \KwInput{ $\mathbb{X}^\texttt{train},  \mathbb{X}^\texttt{test}, \eta, \rho \in [0,1), \alpha, \beta, \lambda, \mathbf{c}$} 
 
  
  
  
 
  
  Initialisse vague prior: $\mathbf{\widehat{p}} = \left[\widehat{p}_{1}, \widehat{p}_{2}, \ldots , \widehat{p}_{m}\right]$ such that $p_{k} = 0.5 , \forall k$ \\
  \While{\texttt{stopping criterion not met}} 
    {
    \underline{Step 1:}\\
        Generate $\eta$ independent samples of binary sets given by 
        $\mathbf{w}^{[i]} = \left[w^{[i]}_{1}, w^{[i]}_{2}, \ldots, w^{[i]}_{m}\right]$, where $w^{[i]}_{k} \sim \texttt{Bernoulli}(\widehat{p}_{k}) ; \quad 1\le i \le \eta, 1\le k \le m$.
         
       
       
        \underline{Step 2:}\\
         From $\eta$ samples, remove all samples that do not satisfy the constraint $\mathbf{c}^{T}\mathbf{w}^{[i]} < \lambda$. Let the new indices be $i = {1 , \ldots, \eta'}$.
            
        
            
                
    $\mathcal{X} = \mathcal{M}(\mathbb{X}^\texttt{train}_{\mathbf{w}^{[i]}} , \mathbb{X}^\texttt{test}_{\mathbf{w}^{[i]}})$ 
     
     Compute $\hat{\mathcal{P}}$, from $\mathcal{X}$ as given in Eq. \ref{eq:performance}. 
    
    Compute  $\mathcal{U}_{i} = \mathcal{U}(\mathbf{w}^{[i]}), i = {1 , \ldots, \eta'}$.
            
             
        
        Compute $\gamma$, the $(1- \rho)^{\text{th}}$-sample quantile of $\mathcal{U}_{1:\eta'}$. 

       

\underline{Step 3:}\\
     Re-estimate model parameters via the MLE:  
       $
\widehat{p}_{k}= 
\alpha
\frac{\sum\limits_{i=1}^{\eta'} \mathbbm{1}\left(\mathbf{w}^{[i]}_{j} =1\right) \mathbbm{1}\left(\mathcal{U}(\mathbf{w}^{[i]}) \geq \gamma\right)}{\sum\limits_{i=1}^{\eta'} \mathbbm{1}\left(\mathcal{U}(\mathbf{w}^{[i]}) \geq \gamma\right)} +(1-\alpha)\widehat{p}_{k}$, $\forall  k$

       
        
       
    }

    For each $p_k$, make the final binary decision as follows:
        $v_{k} =
                \left\{\begin{array}{ll}{1,} & {p_{k} \geq \beta} \\ {0,} & {\text { Otherwise }}\end{array}\right. \forall k \in \left[1,\dots, m \right]$
        where $\beta$ is a pre-defined threshold.
        \\

    \Return $\mathbf{v}^* = \left[v_{1}, v_{2}, \ldots , v_{m}\right]. $

\caption{Cross Entropy Based Feature Selection}
\end{algorithm}

\par

\textbf{Computational Complexity:}
The computational complexity of the CE algorithm for feature selection is characterized by the number of independent samples used $\eta$, the feature size $m$ and the stopping criterion. In Algorithm~1, the calculation of $\mathcal{U}(\mathbf{w}^{[i]})$ is the most computationally expensive step. This is because the computation of $\mathcal{U}(\mathbf{w}^{[i]})$ involves using a classifier and matrix multiplication of the misclassification matrix and the loss matrix. Thus, this step becomes the bottleneck for the performance of the  algorithm, which leads us to conclude that the complexity of computation of this step will also be the dominant component of the complexity of the entire algorithm. Now, the complexity of computation of the $\mathcal{U}(\mathbf{w}^{[i]})$ is given by $\mathcal{O}(\eta \times \mathcal{C}(m; \Psi) \times T_{\max})$ where $T_{\max}$ is the number of iterations and $\mathcal{C}(m; \Psi))$ denotes the complexity of the classifier in terms of the feature length $m$ and other parameters $\Psi$ specific to the classifier. For example, if we use the Naive-Bayes classifier, $\mathcal{C}(m; \Psi) = \mathcal{O}(N_{\texttt{train}}m)$, where $N_{\texttt{train}}$ is the number of training examples and $m$ is the number of features. Since the number of samples is a constant for each iteration, the dependent parameters for the complexity of the algorithm may be reduced down to $\eta$ and $T_{\max}$, which are also the only free parameter that can be controlled by us. The values of $\eta$ and $T_{\max}$ are determined based on the trade-off between the computational budget and the required detection performance.










\section{Feature selection via  brute-force and greedy approaches} \label{sec:others}

In this section we develop algorithmic alternatives to the CE technique that will act as comparisons to our proposed approach.
We first present the brute-force approach for solving the optimization problem (Eq.~\ref{eq:optim}). 
Subsequently, we also develop three greedy heuristics that use different strategies for selecting features for device classification. 

\par 


\subsection{Brute Force Method}\label{VB}


The brute force approach is presented in Algorithm~\ref{alg:brute}. The algorithm selects each combination of the features (via the vector $\mathbf v$), and evaluates the risk score if the cost of the features is within the budget constraint ($\lambda$). Of all such feature combinations, one that minimizes the risk score is the optimal set of features for device classification.
Though this method gives the optimal solution over the feature space, the computational complexity of searching over the entire parameter space exhaustively is prohibitively high. Thus, this approach gives us an upper bound on the run-time complexity of the algorithm, and also a lower bound on the risk.

\par \textbf{Computational Complexity:} The exhaustive search in the algorithm means that, for a feature vector of length $m$, the \texttt{while} loop is executed $2^{m}$ times. 
For each of the $2^m$ executions, a classifier has to be trained (and tested) to determine the loss associated with the selected feature vector. Thus, the overall complexity of the brute force search algorithm is given by $\mathcal{O}(2^{m} \times \mathcal{C}(m; \Psi))$.

\begin{algorithm}[t]
\SetAlgoLined

  \KwInput{$\mathbb{X}^\texttt{train}, \mathbb{X}^\texttt{test}$, $\mathbf{c}$, $\lambda$, $\mathbf{f}$, $\mathcal{M}$,$\mathcal{L}$}
  
  $i \leftarrow 1 , R_{\min} \leftarrow \infty$ \Comment{initializations}\\
  $\mathbf{v} \leftarrow \left[0,  \dots, 0\right]_{(1 \times m)}$ \\
    \While{ $i \le 2^{m}$} 
    {
        $\mathbf{v} \leftarrow \texttt{Binary}(i)$  \Comment{$v_{k}$ is $k$-th digit in the Binary representation of $i, \quad 1 \le k \le m$ }
            
        \If{$ \mathbf{c}^{T}\mathbf{v} < \lambda$}
            {
                
    
    $\mathcal{X} \leftarrow \mathcal{M}(\mathbb{X}^\texttt{train}_{\mathbf{v}} , \mathbb{X}^\texttt{test}_{\mathbf{v}})$ 
     
     Compute $\hat{\mathcal{P}}$, from $\mathcal{X}$ as given in Eq. \ref{eq:performance} 
      
    
    Compute $R(\mathbf{v})$ according to Eq. \ref{eq:cyber_risk} 
            
     \If{$R(\mathbf{v}) < R_{\min}$}
     {
     $R_{\min} \leftarrow R(\mathbf{v})$ \\
     $\mathbf{v}_{s} \leftarrow \mathbf{v}$ \\
     }
             }
            $i \leftarrow i+1$ \\
    
    }
    \Return $\mathbf{v}_{s}$
\caption{Brute Force Feature Selection}
\label{alg:brute}
\end{algorithm}










\subsection{Greedy Algorithms}

As the name suggests, greedy algorithms are efficient ways to solve an optimization problem, with the down side being that, they might end up with a local optimum. We define three intuitive greedy approaches in this section; they differ on the parameter (referred to as \textbf{\texttt{key}}) they used to select the next feature for classification. The general steps for these three greedy approaches are given in Algorithm~\ref{alg:greedy}. First, the feature indices are sorted based on an input parameter \textbf{\texttt{key}}. Subsequently, in each iteration within the \texttt{while} loop, the algorithm selects one feature at a time, constrained by the budget, thereby eventually forming the feature vector (indicated by $\mathbf v$) for classification. 



\subsubsection{Cost-based greedy algorithm (CGA)}
   
    
    In this approach, the feature cost ($\mathbf c$) is used as the input parameter \textbf{\texttt{key}} to sort the (indices of) features. Therefore, this simple greedy algorithm selects the minimal cost feature at each iteration. This would reveal the performance of classification, if only feature cost is used as a criterion to select the features. \\
    \textbf{Complexity:} The sorting algorithm takes $\mathcal{O}(m\log m)$ time for $m$ features. The classifier has to be trained only once in this case, after the selection of all the features. Hence, the complexity of this algorithm is $\mathcal{O}(m\log m + \mathcal{C}(m; \Psi))$.
    
    
\subsubsection{Risk-based greedy algorithm (RGA)}
  
   Next, we define a risk based approach. We train a classifier using only one feature at a time, and we do this for all features constituting~$\mathbf f$. For each classifier corresponding to each feature $f_i, 1 \le i \le m$, we compute the risk using Eq.~\ref{eq:cyber_risk}. The vector of risk scores of length $m$ is then provided as input parameter \textbf{\texttt{key}} to Algorithm~\ref{alg:greedy}. Therefore, RGA algorithm attempts to select the features that correspond to the least risk scores, under the assumption that the risk due to (say) two features is equal to the sum of the risk due to the individual features. \\
  
    \textbf{Complexity:} The classifier is run $m$ times to build $m$ models, and obtain their corresponding risks. Therefore, the overall complexity is $\mathcal{O}(m\log m + m \times \mathcal{C}(m; \Psi))$
    
    
\begin{algorithm}[t]
\SetAlgoLined

    \KwInput{$\mathbb{X}^\texttt{train}, \mathbb{X}^\texttt{test}$, $\mathbf{c}$, $\lambda$, $\mathbf{f}$, $\mathcal{M}$, $\mathcal{L}$, \textbf{\texttt{key}}}
  
  $\Phi \leftarrow \arg \texttt{sort}(\mathbf{f}$, \textbf{\texttt{key}}) \Comment{Obtain the sorted indices}\\
    $i \leftarrow 1$\\
    $\mathbf{v} \leftarrow \left[0, \dots, 0\right]_{(1 \times m)}$ \\
    \While{$i \le m$}
    {
    
    $\mathbf{v}_{\left[\Phi(i) \right]} \leftarrow 1$ \Comment{$\mathbf{v}_{\left[j \right]}$ denotes the $j$-th index of $\mathbf{v}$}
    
    \If{$ \mathbf{c}^{T}\mathbf{v} > \lambda$}
    {
    $\mathbf{v}_{\left[\Phi(i) \right]} \leftarrow 0$
    }
    $i \leftarrow i+1$\\
    }
    $\mathcal{X} \leftarrow \mathcal{M}(\mathbb{X}^\texttt{train}_{\mathbf{v}} , \mathbb{X}^\texttt{test}_{\mathbf{v}})$ 
     
     Compute $\hat{\mathcal{P}}$, from $\mathcal{X}$ as given in Eq. \ref{eq:performance} 
      
    
    Compute $R(\mathbf{v})$ according to Eq. \ref{eq:cyber_risk} 

    \Return $\mathbf{v}$
\caption{Greedy Feature Selection}
\label{alg:greedy}
\end{algorithm}


\subsubsection{Value-based greedy algorithm (VGA)}
    Finally, we also consider how valuable a feature is based on both the accuracy and the feature cost. 
   
    For each feature $f_i \in \mathbf{f}$, we train the classifier similar to what we did in the RGA approach. The $\text{F}_1$ score is computed for each of $f_i$ and is denoted by $\texttt{acc}(f_i)$. As defined in Sec. \ref{sec:cost}, $c_i$ denotes the cost of feature $f_i$. Thus, we define the value of a  feature $f_i$ as: \[\texttt{density}_i = \frac{\texttt{acc}(f_i)}{c_i}, \quad 1 \le i \le m\]
    The vector of values is then passed as the parameter \textbf{\texttt{key}} to Algorithm~\ref{alg:greedy}.
    
    
   
   
   
   

  \textbf{Complexity:} Like the previous greedy algorithm, in this algorithm too, the classifier is run $m$ times to train the model for the $m$ features independently. Therefore, the complexity is $\mathcal{O}(m\log m + m \times \mathcal{C}(m; \Psi))$.








































{
\iffalse
\begin{algorithm}
\caption{Feature Selection in IoT networks using Cross Entropy method}\label{feature}
\begin{algorithmic}[1]
\Require{User's query $Q :=\left(\mathbf{x}_{*}, \varepsilon\right), \alpha, w_{h}, w_{l}$ and $\Psi$ }{}
\State Initialization at iteration $t=0 : \operatorname{set} \mathbf{p}_{0}^{H}=\left\{p_{0,1}^{H}, p_{0,2}^{H}, \cdots, p_{0, N_{\mathrm{II}}}^{H}\right\}$ such that $p_{0, j}^{H}=0.5,$ and set
$\mathbf{p}_{0}^{L}=\left\{p_{0,1}^{L}, p_{0,2}^{L}, \cdots, p_{0, N_{\mathrm{L}}}^{L}\right\}$ such that $p_{0, j}^{L}=0.5 .$

\While {stopping criterion} \Do
\State Generate $K$ independent samples of binary sets
$$\Gamma_{i}^{H}=\left\{\gamma_{i, 1}^{H}, \gamma_{i, 2}^{H} \cdots, \gamma_{i, N_{\mathrm{H}}}^{H}\right\},$$ where $\gamma_{i, j}^{H} \sim \operatorname{Ber}\left(p_{t, j}^{H}\right) ; 1 \leq i \leq K$ and 
$$\mathbf{\Gamma}_{i}^{L}=\left\{\gamma_{i, 1}^{L}, \gamma_{i, 2}^{L} \cdots, \gamma_{i, N_{\mathrm{L}}}^{L}\right\},$$ where $\gamma_{i, j}^{L} \sim \operatorname{Ber}\left(p_{t, j}^{L}\right) ; 1 \leq i \leq K$

\State Calculate the MSE values $\sigma_{*}^{2}(k), k=\{1, \ldots, K\},$ which would be obtained by activating the
corresponding sensors to each of the $K$ samples, according to (15).

\State Evaluate for each of the K samples the performance metric $$\mathcal{U}(k)=\left\{\begin{array}{ll}{-\left(w_{h}\left|\mathcal{S}^{H}\right|+w_{l}\left|\mathcal{S}^{L}\right|\right),} & {\sigma_{*}^{2}(k)<\epsilon} \\ {\infty,} & {\text { Otherwise }}\end{array}\right.$$  where   
$$s_{j}^{H}=\left\{\begin{array}{ll}{1,} & {\gamma_{j, 1}^{H}=1} \\ {0,} & {\text { Otherwise }}\end{array} \quad \text { and } \quad s_{j}^{L}=\left\{\begin{array}{l}{1, \quad \gamma_{j, 1}^{L}=1} \\ {0,}  {\text { Otherwise}}\end{array}\quad \right.\right.$$

\State Calculate the $\beta_{t}=(1-\rho)$ quantile level of $\mathcal{U}_{1 : K}$

\State Update $\mathbf{p}^{H}$ as follows:
$$
p_{t, j}^{H}=\alpha \frac{\sum_{i=1}^{K} \mathbbm{1}\left(\Gamma_{i, j}^{H}=1\right) \overbrace{\mathbbm{1}\left(\mathcal{U}(k) \geq \beta_{t}\right)}^{\text {Choose elite samples }}}{\sum_{i=1}^{K} \mathbbm{1}\left(\mathcal{U}(k) \geq \beta_{t}\right)}+(1-\alpha) p_{t-1, j}^{H}
$$

\State Update $\mathbf{p}^{L}$ as follows:
$$
p_{t, j}^{L}=\alpha \frac{\sum_{i=1}^{K} \mathbbm{1}\left(\Gamma_{i, j}^{\mathrm{L}}=1\right) \overbrace{\mathbbm{1}\left(\mathcal{U}(k) \geq \beta_{t}\right)}^{\text {Choose elite samples }}}{\sum_{i=1}^{K} \mathbbm{1}\left(\mathcal{U}(k) \geq \beta_{t}\right)}+(1-\alpha) p_{t-1, j}^{L}
$$

\EndWhile

\State For each element in $\mathbf{p}^{H}$ and $\mathbf{p}^{L}$ make the final binary activation decision as follows:
$$\mathcal{S}_{j}^{H}=\left\{\begin{array}{ll}{1,} & {p_{t, j}^{H} \geq \Psi} \\ {0,} & {\text { Otherwise }}\end{array} \quad \text { and } \quad \mathcal{S}_{j}^{L}=\left\{\begin{array}{ll}{1,} & {p_{t, j}^{L} \geq \Psi} \\ {0,} & {\text { Otherwise }}\end{array}\right.\right.$$
where $\Psi$ is a pre-defined threshold.

\State Evaluate the objective function in (16) without the sensors which do not have enough energy.
If the QoS constraint is met, then no further steps are required; If the QoS constraint is not met,
solve the optimisation problem again, excluding those sensors which were not able to transmit.

\end{algorithmic}
\end{algorithm}


\begin{comment}
\begin{algorithm}[H]
\SetAlgoLined
  \KwInput{Number of important samples $\eta$, $\rho \in (0,1)$ , initial sample parameter $\hat{\boldsymbol{\theta}}_{0}$, feature budget $\lambda$, feature cost vector $f$, objective function $\mathcal{U}(\mathbf{v})$}  
  \KwOutput{$\hat{\mathbf{v}}=\arg \max _{\mathbf{v} \in \{0,1 \}^m} \mathcal{U}(\mathbf{v})$}
 
  $\hat{\gamma_{0}} \leftarrow -\infty, i \leftarrow 1$  \\
 
  \While{$\mathcal{U}(\mathbf{v}) > \hat{\gamma_{0}}$ and $i \le K$} 
    {
        \If{$\sum_{j} c_{j} \times {v}_{j} < \lambda$}
        {
        Generate $\eta$ independent samples $\mathbf{v}^{[1]}, \ldots, \mathbf{v}^{[\eta]}$ from $\mathcal{D}_{\hat{\theta}_{i-1}}$  
        
        Compute $\mathcal{U}(\mathbf{v}^{[1]}) , \dots , \mathcal{U}(\mathbf{v}^{[\eta]})$
        
        $\hat{\gamma}_i \leftarrow $ the $(1- \rho)$-sample quantile of $\mathcal{U}(\mathbf{v}^{[1]}) , \dots , \mathcal{U}(\mathbf{v}^{[\eta]})$. 
        
        $\hat{\theta}_{i} \leftarrow \arg \max_{\boldsymbol{\theta} \in \boldsymbol{\Theta}} \frac{1}{\eta} \sum_{t=1}^{\eta} \mathbbm{1}_{\left\{R\left(\boldsymbol{\mathbf{v}}^{(t)}\right) \ge \hat{\boldsymbol{\gamma}}_{i}\right\}} \log \mathcal{D}_{\boldsymbol{\theta}}\left(\boldsymbol{\mathbf{v}}^{[t]}\right)$ \\
       } 
        $i \leftarrow i+1$

    }
    Set $\hat{\mathbf{v}}$ to be the sample with minimum $R(\hat{\mathbf{v}})$ so far.\\
    \Return $\hat{\mathbf{v}}$

\caption{Cross Entropy Feature Selection}
\end{algorithm}


\begin{algorithm}[H]
\label{alg:cross}
\SetAlgoLined
  \KwInput{Number of independent samples $\eta$, $\rho \in (0,1)$ , initial sample parameter $\boldsymbol{\theta}$, feature budget $\lambda$, feature cost vector $f$, objective function $\mathcal{U}(\mathbf{v})$ , feature vector size = N} 
  \KwOutput{$\mathbf{v}_{O} = \underset{\mathbf{v} \in\{0,1\}^{m}}{\arg \min } \mathcal{U}(\mathbf{v})$}
 
  $\hat{\gamma_{0}} \leftarrow -\infty$  \Comment{Initialize Threshold}
  
  t \leftarrow 0 \Comment{Initialize Iterator}
  
  Set $\mathbf{v}_{0} = {v_{0}^{[1]}, v_{0}^{[2]}, \ldots, v_{0}^{[K]}}$ such that $v_{0}^{[j]} = 0.5 , \forall j$ \\
  \While{$\mathcal{U}(\mathbf{v}) > \hat{\gamma_{0}}$} 
    {
        Generate $\eta$ independent samples of binary sets given by 
        $\Gamma_i = \{\mathbf{w}_{i,1},\mathbf{w}_{i,2} \ldots, \mathbf{w}_{i,K}\}$ where. $\mathbf{w}_{i,j} \sim Ber(v_{t}^{[j]}) ; 1\le i \le \eta$
         
        For each of the $\eta$ samples, compute  $\mathcal{U}(\mathbf{v}^{[m]})$ as given in Eq: (\ref{eq:risk}), where $v_{j}^{[m]} =\left\{\begin{array}{ll}{1,} & {\mathbf{w}_{j,1} = 1} \\ {0,} & {\text { Otherwise }}\end{array}\right.$.\\
        
        $\hat{\gamma}_t \leftarrow $ the $(1- \rho)$-sample quantile of $\mathcal{U}(\mathbf{v}^{[1]}) , \dots , \mathcal{U}(\mathbf{v}^{[\eta]})$. 

        $v_{t}^{[j]}= \beta \frac{\sum_{i=1}^{\eta} \mathbbm{1}_{\left(\Gamma_{i, j} =1\right)} \overbrace{\mathbbm{1}_{\left(\mathcal{U}(\mathbf{v}^{[m]}) \geq \hat{\gamma}_t\right)}}^{\text {Choose elite samples }}}{\sum_{i=1}^{\mathcal{N_S }} \mathbbm{1}_{\left(\mathcal{U}(\mathbf{v}^{[m]}) \geq \hat{\gamma}_t\right)}}$ + $(1-\beta)v_{t-1}^{[j]}$
        
        $t \leftarrow t+1$
    }
    
    
    For each element in $\mathbf{v}$, make the final binary decision as follows:
        $\widehat{v}_{j} =\left\{\begin{array}{ll}{1,} & {v_{t}^{[j]} \geq \Psi} \\ {0,} & {\text { Otherwise }}\end{array}\right.$
        where $\Psi$ is a pre-defined threshold\\
    Set $\mathbf{v}_{O}$ to be the sample with minimum $R(\mathbf{v}_{O})$ so far.\\
    \Return $\mathbf{v}_{O}, R(\mathbf{v}_{O}) $

\caption{Cross Entropy Feature Selection}
\end{algorithm}

\end{comment}
\fi
}


\section*{Appendix: Feature cost}
In this section we demonstrate how feature costs are estimated for a given IoT device.
We defined the cost vector $\mathbf{c}$ to indicate the cost associated with each of the  extracted features. The feature cost is divided into three different components. First, we analyze the memory required for the extraction of the features from the network traffic. Next, we focus on the computational complexity of the algorithm required for extracting the feature. And lastly, we analyze the intrusiveness of a feature into the privacy aspects of a device or user. 
For example, let us consider the feature \texttt{dns\_num\_ans}, which denotes the the number of answers returned for DNS queries. While only a register is required to store the (running) value of this feature, a small buffer has to be maintained to temporarily store (unprocessed) DNS answers extracted from packet payloads (to absorb traffic bursts); therefore, we define the memory cost as \texttt{medium}. Similarly, this feature is computed by simple additions to the counter in the register. Therefore we can assume a \texttt{low} computational complexity. Finally, since the feature is extracted from the DNS responses (i.e., packet payload), which can leak considerable information (e.g., see~\cite{DNS-ads-profit-2011}), we assign a \texttt{high} privacy cost to the feature. 
Table~\ref{tab:feature_costs_details} provides the logic we have used to define the three cost components for extracting any given feature from network traffic. 
To define the total cost of extracting a feature, we take the median of all three costs (memory, compute and privacy). 

\begin{table}[h]
\centering
\caption{Cost logic considered for experiments}
\resizebox{\columnwidth}{!}{%
\begin{tabular}{|c|c|c|c|c|}
\hline
\multirow{2}{*}{\textbf{\#}} & \multicolumn{3}{c|}{\textbf{Parameters}} & \multirow{2}{*}{\textbf{Cost}} \\ \cline{2-4}
 & memory & compute power & privacy &  \\ \hline
1 & \begin{tabular}[c]{@{}c@{}}Single register \\ (constant memory)\end{tabular} & \begin{tabular}[c]{@{}c@{}}Counters \\ (e.g., $\min, \max$, etc.) \end{tabular} & \begin{tabular}[c]{@{}c@{}}Features extracted \\ from packet headers\end{tabular} & \texttt{low} \\ \hline
2 & \begin{tabular}[c]{@{}c@{}}Small buffer \\ (e.g., queue) \end{tabular} & \begin{tabular}[c]{@{}c@{}}Maintenance of \\ hash tables \end{tabular} & \begin{tabular}[c]{@{}c@{}}Extraction of \\ application data (e.g., URL)\end{tabular} & \texttt{medium} \\ \hline
3 & \begin{tabular}[c]{@{}c@{}}Hash table or \\ multiple registers\end{tabular} & \begin{tabular}[c]{@{}c@{}}Pattern matching \\ or sorting \end{tabular} & \begin{tabular}[c]{@{}c@{}}Features extracted from \\ packet payloads \end{tabular} & \texttt{high} \\ \hline
\end{tabular}%
}
\label{tab:feature_costs_details}
\end{table}





\section{Conclusions} \label{Sec:conclusions}

In this work, we motivated, defined and focused on the problem of feature selection for IoT device classification considering the fact that there is a budget constraint for features in practice. For this purpose, we also defined the notion of risk of misclassification. To obtain the optimal solution to this problem, one has to perform a combinatorial search over the solution space. Therefore, we developed a cross entropy based algorithm for solving the optimization problem. We carried out experiments using traffic of real IoT devices; our experiments showed that not only is CE algorithm practical and much faster than the brute force approach, it also obtains low risk score for classification and performs better than than value-based greedy algorithm in most cases. 


\section{Introduction}

The rapid growth of the Internet of Things (IoT) market is also having an impact on the threat landscape of the cyber space~\cite{symantec-iot-attack-report-2019}. 
We now have attacks that compromise large numbers of IoT devices, subsequently using them as bots for launching different kinds of large-scale attacks (e.g., DDoS attacks) that result in significant losses~\cite{Understanding-Mira-2017, hajime-NDSS-2019}. Attackers often explore and exploit vulnerabilities to gain access to machines. In this IoT era, it is now possible to exploit a few vulnerabilities of just one device type, and compromise hundreds of thousands of devices of the same type using the same exploits.

To detect and prevent attacks on and from IoT devices, one of the important tasks is to identify the type of devices connected to a network. 
Device identification is useful in keeping track of various device types in a network, which consequently helps in analyzing and defending against potential vulnerabilities of various IoT devices. In addition, a mitigation solution---such as blocking devices with newly known vulnerabilities---can be implemented quickly, if the network administrator knows the identities (brand, model number, functionality, etc.) of the devices connected to the network. Recent years have seen research proposals on identifying IoT devices, particularly by analyzing their network traffic (e.g., see~\cite{Meidan2017ProfilIoT,miettinen2017iot,fingeprinting-behind-NAT-VPN-2019} and \cite{audi-JSAC-2019}). 
Existing solutions often use machine-learning based approaches to fingerprint IoT devices, by extracting a number of {\em features} from IoT traffic. For instance, in the case of supervised or semi-supervised approaches (e.g.,~\cite{dev-classification-Trans-MC-2018} and~\cite{DEFT-2019}), the features extracted are used to train a classification model to differentiate IoT devices while in operation. 

A common underlying assumption of the device classification works in the literature is that all features can be easily extracted from the traffic and that there is no difference---technical or otherwise---in extracting any two different features. However, we argue that, such an assumption misses out some important points. Extraction of features from traffic incurs a cost, and different features have different costs for extraction. We identify three types of feature-extraction costs: \\

\noindent 1. {\em Computational cost:}  When network traffic characteristics are used for classifying devices, the relevant features are computed by reading and processing traffic at a router, a gateway or a server to which the traffic is mirrored to. That is, there is cost involved in computing a feature from packets, and different features may incur different computational costs. For example, searching for a specific pattern, say, malware signature/hash, in packets is computationally more expensive than counting the number of packets. 
    
\noindent 2. {\em Memory cost:} Feature extraction also involves use of memory to store the (running) value of a feature. For example, storing the number of packets traversing a link requires just one counter, but to count the number of packets of each connection that is active at a link requires maintenance of a hash table or some variants~\cite{HashPipe-SOSR-2017}. And memory is an extremely valuable resource in routers and gateways. Indeed, traditionally ISP routers use sampled NetFlow~\cite{RFC7011} to capture and aggregate meta-information of traffic flows,  due to resource constraints at the router. With NetFlow, information which is often assumed to be easily available (e.g., size of each packet in a connection) is not collected. 
    
\noindent 3. {\em Privacy cost:} The third factor that decides the cost of a feature is privacy. Some features, such as those extracted from payloads of (unencrypted) network traffic, can reveal sensitive or even confidential information of end-users~\cite{smart-home-privacy-2016}. In fact, DNS or even DNS over HTTPS (that is, DNS over encrypted HTTP connection) can reveal sensitive information related to users communication~\cite{DNS-ads-profit-2011, DoH-2019}. Privacy concerns imply that, either some features might be unavailable (i.e., cost is infinite) or there would be some cost incurred to obtain such data, for example by paying for the data with user's consent or by anonymizing the data to minimize privacy leak.
    



Therefore, we argue that, a solution for device identification (or classification) should take the cost of features into consideration. This in turn means, the performance of a device classification solution is  dependent on the cost of features used, since the cost may exclude some important features. 
A solution developer would need to be conscious of the budget available for deploying the solution. 

Previous works also assume no difference in device misclassification, in that they were not concerned by the misclassification class to which an incorrect assignment was made. In practice, not only the fact that a misclassification took place is important, it is also important to consider which class the device was misclassified under as this could have ramifications for actions and costs. Previous works also assume no difference in device misclassification; i.e., they were not concerned of the type a given device is being misclassified into. However, such an approach overlooks important aspects IoT devices. An IoT device, say a camera, might be from vendor X with model name Y; but the same vendor might have multiple models for the product camera. Consider the different misclassifications possible: (i)~[camera, X, Y] gets classified as smart bulb, (ii)~[camera, X, Y] gets classified as [camera, A, B], and (iii)~[camera, X, Y] gets classified as [camera, X, Z]. Clearly, the third misclassification is more acceptable than the first two, since the functionality and brand name were correct; and the first misclassification is the worst result to have among the three. Therefore, we argue the need to consider {\em multiple losses} due to misclassifications that essentially captures the differences in misclassifications. 



In this work, we address the problem of classifying IoT devices, under the realistic assumptions that (i)~features have associated costs, (ii)~a solution developer has a budget constraint, and (iii)~there could be different kinds of losses due misclassifications. Different from existing works, the challenge here is to develop a cost-optimal IoT device classification solution with high accuracy. 

Our contributions are the following:
\begin{enumerate}
    \item We introduce the notion of feature cost and budget constraint in the problem of IoT device classification. Furthermore, going beyond the current way of treating misclassification as a binary value, we consider multiple losses due to misclassification (Sec.~\ref{subsec:loss}), and introduce the notion of {\em risk} (Sec.~\ref{subsec:risk}) to evaluate a classifier's performance. Subsequently, we formulate the feature-selection problem under constrained budget as an optimization problem (Section~\ref{sec:problem}). 
    
    \item We present and develop a cross-entropy (CE) based algorithm (Section~\ref{sec:cross}), that solves the optimization problem efficiently. In comparison to a brute force approach which has to run a classifier exponential number of times, the number of executions (of classifier) that the CE algorithm has to run is linear in its two parameters. 

    \item We conduct extensive experiments using traffic of real IoT devices, analyzing a brute force approach, multiple greedy algorithms and our proposed CE-based algorithm, comprehensively. 
\end{enumerate}


We believe our framework for cost-optimal feature selection is also applicable for other problems, where features have associated costs and solution deployment has a budget constraint. Examples of such problems in the domain of network traffic analysis include the traditional network application classification~\cite{SLIC-2015}, botnet detection~\cite{disclosure-ACSAC-2012}, anomaly detection in traditional~\cite{NADA-2018, GEE-CNS-2019} as well as IoT~\cite{DIoT-ICDCS-2019} networks, etc. 

After defining the system model in Section~\ref{sec:system},   we formulate the optimisation problem in Section \ref{sec:problem}. In Section~\ref{sec:cross}, we develop the CE-based algorithm. Subsequently, in Section~\ref{sec:others}, we present three greedy algorithms for selecting features while respecting the provided budget. We carry out experiments and present results in Section~\ref{sec:simulations}.


\begin{comment}
\subsection{Cost of feature extraction}
Features can be extracted from IP headers of packets or from packet payloads. For example, the ports
accessed by a set of traffic connections can be extracted from the IP header of packets. On the other hand,
to extract the number of DNS queries for a particular domain names requires processing of contents of
DNS packets. Depending on the requirement, features are also represented at different granularities. At
the finest granularity, per-packet features are used for learning differentiating characteristics. At the other
end of the spectrum, features represent sessions. A session is a collection of traffic flows (identified by
the five-tuple of source and destination IP addresses, source and destination ports, and protocol) localized
in time. A session is identified by the time-interval for which the feature vector is extracted. The statistics
of the first n bytes of every packet payload form a per-packet feature vector. Whereas, the mean and
standard deviation of packets in a five-minute time bin is an example of a session-level feature vector.
 
For every feature that can be of use in device classification (or in general, traffic analysis), there
is an associated cost. The costs are due to a few reasons, varying from technical to privacy concerns. Every feature requires some processing of traffic. The number of packets going through a network
interface is one of the simplest and least expensive features that can be extracted from a router. As
the complexity of features increases, more computational and memory resources are required, thereby
increasing the cost of obtaining the feature. Besides, some features might simply be unavailable as the
appliance in operation does not support processing of the corresponding data. For example, extracting
IP header information of each transiting packet is considered impractical in core routers. Even meta
information of packets, such Netflow records {\color{red}ToDo}: , may be supported only with sampling, as unsampled
Netflow requires upgrade of not only existing routers, but also the infrastructure to support large volumes
of data storage and processing. Hence, obtaining such information results in costly upgrades. There
is also costs due to privacy concerns. For example, the destination IP addresses or HTTP URLs that
the clients visited might be considered as private information; and therefore, there may exist a policy
that mandates anonymization of client IP addresses to minimize the possibility of learning behaviors of
particular clients. Anonymization, in a way that is useful for analysis (say, using longest prefix preserving
anonymization {\color{red}ToDo}: cite) adds to the cost of feature extraction. One could also purchase dataset or use
a paid subscription service to enrich the data obtained from network traffic {\color{red}ToDo}: cite geo-location info.

We further motivate the different costs associated with feature extractions using an example scenario
below. {\color{red}ToDo}: is this redundant?

\begin{itemize}
    \item Counting the number of packets in a session defined using predefined static time-interval only
requires maintenance of a counter at the network element of interest (say, an edge router). Let C1
be the associated cost to maintain this counter online. This feature does not require us to read a
packet’s payload
\item Now, to include further information for classification, if we want to segregate the counters based on
the application protocol (HTTP, DNS, SSH, etc.), we need to maintain counters for each protocol;
besides we need to read the port number from the IP header. Let C2 be the associated cost of
extracting and maintaining segregated counters as features for classification.
\item However, port numbers do not always provide reliable information to identify the application protocol
in use {\color{red}ToDo}: cite; besides, there are applications that use dynamic ports {\color{red}ToDo}: . Hence, we may
have to process the payloads, if the solution being developed demands this level of accuracy. Let
C3 be the cost of processing packet header and payload to maintain highly accurate per protocol
counter
\item Due to privacy issues, it is quite likely that the payload information is not available as plain text.
For such cases, there is additional cost to obtain the necessary information from ‘masked’ payloads. Let C4 denote the accumulated costs including the cost due to masking
\end{itemize}
Clearly, in the above example, the cost of feature extraction is such that, $C1 < C2 < C3 < C4$.
Given that different features have different associated costs for extraction, there is a need to analyze a
classification system based on the cost of feature extraction.


\subsection{Risk of misclassification}

A multi-class classification system is usually analyzed using metrics such as recall, precision and
overall accuracy. While these are necessary, we argue that these metrics are insufficient for classifying
IoT devices. There are numerous IoT devices for a single functionality, and even multiple ones from
a particular. The risk of misclassifying an IoT device to a device of the same functionality is different
from misclassifying the device into a completely different functionality. For example, misclassifying a
camera from vendor X to a camera of vendor Y has different (and possibly lower) risk of misclassifying
the camera as a smart socket. Similarly misclassification into the same functionality, and also into the
the same vendor products, but different models, might be considered an ‘acceptable’ error. {\color{red}ToDo}:

\subsection{Objective}

The objective of this work is to analyze a fingerprint classification system based on, both, the cost of
feature extraction and the risk of misclassification. ToDo:

\end{comment}

\section{Performance Evaluation}\label{sec:simulations}

In this section, we describe the experiments performed to evaluate the different algorithms presented in the previous sections. After a briefing on the dataset used, we describe how we estimated costs of features (in Section~\ref{subsec:feature-costs}). We also define the loss matrix used for experiments in this work (Section~\ref{subsec:loss-eval}). The first set of experiments we carry out (in Section~\ref{subsec:compare-classifiers}) are to evaluate a selected set of classification models. Subsequently, we evaluate brute force, CE and greedy algorithms, first under no budget constraints (Section~\ref{subsec:eval-no-budget}), and then later we evaluate CE and greedy algorithms when budget is fixed (Section~\ref{subsec:eval-budget}). In Section~\ref{subsec:eval-budget}, we also perform further analysis to compare CE algorithm with the best greedy algorithm. 





\subsection{Real datasets for experiments}

For the experiments, we used network traffic from 15 different IoT devices as listed in Table I. The devices were connected to the Internet via a gateway, where the network traffic was captured. The traffic was split into intervals of 15 minutes each, and the number of sessions captured for each device is shown in Table~\ref{tab:devices}. From the network data of the IoT devices, 111 features were extracted (see~\cite{DEFT-2019} for the list of features). 




\begin{table}[t]
\centering
\caption{Information on IoT devices used}
\label{tab:devices}
\begin{tabular}{llll}
\hline
Label & Device & Brand  & \begin{tabular}[c]{@{}l@{}}Sessions\\ Captured\end{tabular} \\ \hline
1 & Echo Dot & Amazon  & 490 \\ \hline
2 & Smart Remote & Broadlink & 480 \\ \hline
3 & \begin{tabular}[c]{@{}l@{}}Camera\\ (DCS700L)\end{tabular} & D-Link & 384 \\ \hline
4 & \begin{tabular}[c]{@{}l@{}}Camera\\ (DCS5030L)\end{tabular} & D-Link & 410 \\ \hline
5 & \begin{tabular}[c]{@{}l@{}}Smart Socket\\ (DSPW215)\end{tabular} & D-Link & 672 \\ \hline
6 & Chromecast & Google & 297 \\ \hline
7 & Home Control & Google & 529 \\ \hline
8 & Smart Socket & Oittm & 394 \\ \hline
9 & Hue Light & Phillips & 644 \\ \hline
10 & Smart Things & Samsung & 587 \\ \hline
11 & \begin{tabular}[c]{@{}l@{}}Smart Bulb\\ (LB100)\end{tabular} & TP-Link & 482 \\ \hline
12 & \begin{tabular}[c]{@{}l@{}}Camera\\ (NCS250)\end{tabular} & TP-Link & 587 \\ \hline
13 & \begin{tabular}[c]{@{}l@{}}Camera\\ (NCS450)\end{tabular} & TP-Link & 494 \\ \hline
14 & \begin{tabular}[c]{@{}l@{}}Smart Socket\\ (HS100)\end{tabular} & TP-Link & 452 \\ \hline
15 & \begin{tabular}[c]{@{}l@{}}Smart Socket\\ (HS110)\end{tabular} & TP-Link & 387 \\ \hline
\end{tabular}%
\end{table}



\subsection{Feature costs}
\label{subsec:feature-costs}

As described earlier, the extraction of each feature has a cost. We argue that the cost of feature extraction is a function of three factors: the computational power (or computational complexity) of the extraction process, the memory used, and the confidentiality of the information related to that feature. Therefore, we break the cost into three components, each corresponding to compute power, memory and privacy.
Finding the exact cost due to each component and how to integrate those component costs into a single cost value is outside the scope of this work. Instead, we simplify the cost of each component to be in one of three levels \{\texttt{low, medium, high}\}. We illustrate the concept of assigning costs to feature using examples in the Appendix.

The total cost of a feature $f_i$ can be defined as:
\[\text{Cost}_i = g_i(\text{compute cost}, \text{memory cost}, \text{privacy cost}), \]
where $g_i(.)$ is a function computing the total cost of the feature extraction process. 
For our purposes, we define $g_i, 1 \le i \le m$ to be the median of the three input cost components. While there are other ways of integrating component costs (like, ``cost is always \texttt{high} if privacy cost is \texttt{high}''), we stick to this simple definition for our work here. For example, consider the feature of `{\em connection length, in  number of packets}'. To extract this feature, connection (that is, a 5-tuple flow) identifier has to be hashed and stored in a data structure such as a hash table, and the number of packets needs to be counted~\cite{REX-2017}. 
The computational cost is \texttt{medium} since a hashing is required for every arriving packet, and hash table operations would at worst case be be linear with the number of packets (e.g., for insertion of new 5-tuple flow in the traditional hash table), besides flows have to be regularly removed from the table once they become inactive. 
The memory requirement is assigned as \texttt{high} due to the necessity to maintain a hash table. 
The privacy cost is considered \texttt{low}, as only packet counts of connections are extracted, and no private information (e.g., visited websites) is extracted. Therefore, the cost of extracting {\em connection length} is \texttt{medium}. 







\subsection{Loss matrix, $\mathcal L$}
\label{subsec:loss-eval}

We introduced the concept of misclassificaiton loss in Section~\ref{subsec:loss}. Ideally, the loss matrix is an input provided by the users, based on how they perceive the loss due to misclassification errors. For the purpose of this work, we define the loss matrix as follows.

When the classifier classifies a device $d_i = <t_i, b_i>$ as $d_j = <t_j, b_j>$, the loss value $l_{i,j}$ is calculated as $l_{i,j} = 2*f(t_i, t_j) + f(b_i, b_j),$ where $f(x_i~, x_j) = 1$ if $x_i \neq x_j$ and $0$ if $x_i = x_j$.
Taking every device pair, we follow this definition to derive the values for the loss matrix $\mathcal L$ in the experiments here. It is to be noted that we are giving a higher value of loss if the classifier cannot classify the type of the device even if it gets the right brand. For example, if the classifier classifies a Samsung camera as a Sony camera, the loss value $l_{i,j} = 2 \times 0 + 1 = 1.$ But if it classifies the device as Samsung (Smart things) Hub, then the loss $l_{i,j} = 2*1 + 0 = 2$. 


\vspace{-0.2cm}
\subsection{Comparison of classifiers}
\label{subsec:compare-classifiers}


The cyber risk score $R(\mathbf{v})$ of a feature vector $\mathbf{v}$ depends on the misclassification matrix which is in turn derived from the confusion matrix of the classifier. As such, the choice of classifier plays an important role in evaluating the performance of the system.
We selected four different classification models---Gaussian Naive-Bayes classifier, SVM (support vector machine) with RBF kernel, Decision tree and Random Forest, and compared their performances based on accuracy and the execution time. 
We use $\text{F}_1$-score (the harmonic mean of precision and recall) to represent the accuracy of the classifier. 

All the experiments in this work are carried out on an 8-core Intel Core i7-2600 running at 3.8GHz CPU and equipped with 16GB of RAM.
For comparing the time of execution, we consider both the training and the testing times of the classifiers. Fig.~\ref{fig_time_class} plots the time taken for each classifier for training and predicting on the traffic session of IoT devices. The X-axis here (as well as in all figures in this paper) denotes the number of features considered for an experiment; therefore, each point on the X-axis corresponds to an independent run of the experiment with that many number of features given as input. Since the number of experiments increases with the number of features, we limit the total number of features provided as input in each experiment to 69; beyond 69 features, the experiments did not provide any new insights.

From Fig.~\ref{fig_time_class}, we observe that the Decision tree classifier and the Naive-Bayes classifier take the least time for model building and prediction. 


\begin{comment}
\begin{figure*}
\begin{minipage}{0.32\textwidth}
\centering
\includegraphics[width=1.1\linewidth]{FinalFigures/time_cc_new_100.pdf}
\caption{Execution time of classifiers for varying number of features $m$}
\label{fig_time_class}
\end{minipage}\hspace{0.2cm}
\begin{minipage}{0.32\textwidth}
\centering
\includegraphics[width=1.1\linewidth]{Figures/accuracy_cc_new_100.pdf}
\caption{Classification accuracy as a function of the number of features $m$}
\label{fig_acc_com}
\end{minipage}\hspace{0.2cm}
\begin{minipage}{0.32\textwidth}
\centering
\includegraphics[width=1.1\linewidth]{Figures/timet_updated.pdf}
\caption{Execution time; $\lambda \rightarrow \infty$}
\label{fig:sim_time}
\end{minipage}
\end{figure*}
\end{comment}


\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{FinalFigures/time_cc_new_100.pdf}
\caption{Execution time of classifiers for varying number of features $m$}
\label{fig_time_class}
\vspace{-0.2cm}
\end{figure}






\begin{figure}[t]
\centering
\includegraphics[width=0.9\linewidth]{FinalFigures/accuracy_cc.pdf}
\caption{Classification accuracy as a function of the number of features $m$}
\label{fig_acc_com}
\vspace{-0.2cm}
\end{figure}

We also plot the accuracy of the different classifiers in Fig.~\ref{fig_acc_com}; all of them achieve high $\text{F}_1$-scores. Since Decision tree classifier has much less time of execution, while also having high $\text{F}_1$-score, henceforth we use Decision tree for the rest of the experiments to evaluate cost-aware feature selection algorithms. 















\subsection{Evaluation of algorithms without budget constraint}
\label{subsec:eval-no-budget}

Next, we compare the five different algorithms discussed in sections \ref{sec:cross} and \ref{sec:others}, by varying the number of features $m$ and with no budget constraint, $\lambda \rightarrow \infty$. We evaluate their performance based on (i) the time of execution and (ii) the cyber risk score.
\par 

As mentioned above, for each experiment, we need to provide a set of features; and the algorithms are supposed to find the best features from the given input set of features. To do this, for each evaluation, we can start with a minimum feature set and increase the set by one feature for each experiment. That is, since we have 69 features, if we start with a singleton set and keep adding one feature for every new experiment, we would need 69 experiments. However, this is under the assumption that, we know the sequence in which the features should be added to the set. Consider adding the most discriminatory feature as the last one to the growing set of features; this means no other experiment would use this important feature for evaluation. On the other hand, the ideal case would be to generate all possible sequences --- 69! in our case, and carry out experiments by expanding set from a minimum size to the maximum of 69 features, for each sequence.

Since this is not computationally feasible, unless stated otherwise, we generate sequence of {\em ranked} features. 
To this end, we use the feature ranking mechanism described in our previous work~\cite{desai2019feature} where the features are ranked according to their ability to discriminate between every pair of devices by combining the results from multiple statistical tests.
In order to compare and evaluate the algorithms, we take the full set of ordered feature vector  $\mathbf{f}$ (obtained via ranking),
and select the first nine features from it as $\mathbf{\bar{f}} = \left[f_1, f_2, \dots , f_9 \right]$. For each subsequent iteration, we add the next feature in the ranked order to our selected feature vector and calculate the time of execution and the risk score of all the five algorithms using the new set of selected features.
Thus for the $t^{th}$ iteration, the selected feature vector is represented as $\mathbf{\bar{f}} = \left[f_1, f_2, \dots, f_{t+9} \right]; t = 0, 1, 2 , \dots, 61$, where $|\mathbf{\bar{f}}| = t+9$. Therefore, the algorithm must be iterated 61 times, such that, when $t= 61$, all the 69 features will be explored.
We proceed to the analyses next.






\subsubsection{Execution time analysis}
The comparison of the execution times of the different algorithms help us to get an estimate of their computational complexities in practice. 

In this section, we compare the time of execution for the CE algorithm with respect to the brute force and the greedy algorithms, under no budget constraint ($\lambda \rightarrow \infty$). 





From Fig. \ref{fig:sim_time}, we observe that the time of execution for the brute force algorithm increases exponentially as the number of features in the feature set increases. Thus, though the brute force algorithm has a low time of execution for small feature sets, it is extremely high for larger ones. This is why we limit the experiment on the brute force to 25 features. The greedy algorithms perform the best with respect to the execution time, which gradually increases with increasing $m$. Though the proposed CE algorithm takes more time than the greedy algorithm, its execution time is much lower than the brute force approach.  Also, the growth in run-time with increasing number of features is marginal for the CE algorithm, thereby indicating the ability of the algorithm to scale.

\begin{figure}
\centering
\includegraphics[width = \linewidth]{FinalFigures/time_updated_new.pdf}
\caption{Execution time as a function of number of features $m$ with $\lambda \rightarrow \infty$}
\label{fig:sim_time}
\vspace{-0.3cm}
\end{figure}





\subsubsection{Risk analysis}
The risk score $R(\mathbf{v})$ is another important parameter to evaluate the performance of the feature selection algorithms as we aim to minimize the risk score by the optimal selection of features. 
The graphs for the risk scores of the five different algorithms under consideration without any budget constraint is plotted in Fig. \ref{fig:sim_loss}.
Since this scenario does not constrain the budget, i.e., $\lambda \rightarrow \infty$, note that,  the three greedy algorithms will always choose all the features provided and thus, no difference is expected in the performance (in terms of accuracy and risk score) of the three different greedy approaches. 
We observe that the brute force algorithm gives the best risk score since it performs an exhaustive search over all possible combinations of the provided features. However, it can also be observed that the cross entropy algorithm is close in performance to that of the brute force method, and ultimately converges to the minimal risk score of~1.
Due to this convergence, it may be concluded that any feature more than the first~63 (in the ranked order) is redundant in the classification of the devices when using the cross-entropy algorithm. 
In addition to this, it can also be observed that though the greedy algorithms select all the features, they do not give the best results. This highlights the necessity for the selection of an optimal set of features, and that the selection of a higher number of features does not necessarily imply a better risk score. 

\begin{figure}
\centering
\includegraphics[width = \linewidth]{FinalFigures/risk_updated1.pdf}
\caption{Risk score $R(\mathbf{v})$ as a function of number of features $m$ with $\lambda \rightarrow \infty$}
\label{fig:sim_loss}
\vspace{-0.4cm}
\end{figure}



\subsection{Evaluation of algorithms under budget constraint}
\label{subsec:eval-budget}




\begin{figure}
\centering
\includegraphics[width = \linewidth]{FinalFigures/greedycomp.pdf}
\caption{Risk score $R(\mathbf{v})$ as a function of the number of features $m$ for greedy algorithms with low and high budget $\lambda = \left\{50,200\right\}$}
\label{fig:greedy_budget}
\vspace{-0.3cm}
\end{figure}


\begin{table}[t]
\centering
\caption{Comparison of Greedy Algorithms}
\label{tab:comp_greedy}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
 & \multicolumn{3}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Feature Budget \\ $\lambda = 50$\end{tabular}}} & \multicolumn{3}{c|}{\textbf{\begin{tabular}[c]{@{}c@{}}Feature Budget\\ $\lambda = 200$\end{tabular}}} \\ \hline
 & \textbf{CGA} & \textbf{VGA} & \textbf{RGA} & \textbf{CGA} & \textbf{VGA} & \textbf{RGA} \\ \hline
\textbf{Mean} & 68.28 & 64.82 & 98.03 & 69.79 & 70.04 & 68.11 \\ \hline
\textbf{\begin{tabular}[c]{@{}c@{}}Standard \\ Deviation\end{tabular}} & 41.97 & 36.72 & 51.62 & 18.25 & 18.84 & 20.79 \\ \hline
\end{tabular}
\end{table}


\textbf{Comparison of greedy algorithms:}
We first compare the performance of the three different greedy algorithms under budget constraint to determine the best-performing greedy algorithm. We consider two cases: one with low budget for features, $\lambda = 50$, and another with high feature budget, $\lambda = 200$. We plot the  risk scores of the three greedy algorithms Fig.~\ref{fig:greedy_budget}. 

Similar to previous experiments, each point on the X-axis on Fig.~\ref{fig:greedy_budget} denotes the number of features given as input to the algorithms. However, to limit the number of experiments, we start with nine features and add three features with every new experiment. Therefore the point 30 on the X-axis means, 30 features were given as input to the algorithms in the experiment, and the next point corresponds to 33 features given as input. We observe that, all the three greedy algorithms have similar performance when the budget is high. With high budget, all features given at input can be selected for building the classifier. However, this is not the case for lower budget---under tighter budget constraint, we see that the three algorithms have different performances as they choose different (subsets of) features from the given list of input features. For easier comparison, we also provide the mean and standard deviation of risk scores for each of the algorithms in Table~\ref{tab:comp_greedy}. From these results, it can be seen than VGA performs better than the other two under low budget (while all three have comparable performance under high budget). 

An interesting aspect to note in the low-budget case  is that, the risk scores do not decrease monotonically with increasing features provided at the input, instead we see spikes in the risk score (Fig.~\ref{fig:greedy_budget}). To understand this, let us consider the first spike---the point on the X-axis where the first 24 ranked features are added. Observe, all three algorithms experience a spike in the risk score at this point. We present the analysis of VGA. Compared to the previous point on the X-axis, i.e., the experiment corresponding to the first 21 features, the features \texttt{mdns len mean}, \texttt{mdns duration}, \texttt{mdns num ans}, and \texttt{http time mean} were newly added to the 
classifier model by VGA, but more importantly the features \texttt{dns num qns} , \texttt{dns qry cls},  \texttt{http len mean}, and \texttt{tcp keep alive} were removed by the algorithm. We further used the SHAP method~\cite{SAHP-NIPS-2017} to estimate the contributions of different features in classifying the devices. SHAP  uses the concept of Shapley values from Game theory, to explain the predictions of a machine learning model in a similar way as computing the contributions of players in a game. Using the SHAP method, we note that, the second-most important feature in the list of first 21 features was \texttt{tcp keep alive}, and this was removed by VGA in the next experiment when three more features were added. Besides, none of the MDNS related features that were newly added by VGA were in the top eight most contributing features. We observe a similar removal of important features in both RGA and CGA, thus rendering them ineffective in selecting the most discriminative features to minimize risk.



Based on the above experiments, to compare the performance of the greedy methods with the CE algorithm, we select VGA as the candidate as it gives the best performance among three greedy approaches. \\


\noindent \textbf{Comparison of CE algorithm and VGA, for ranked and random ordering of features:} The risk score calculated is dependent on the ordering of the features $\mathbf{f}$ and hence to test the performance of the algorithms, we carry out the experiments for not only the ranked feature vector, but also for randomly ordered features. The risk score $R(\mathbf{v})$ is a function of the feature budget $\lambda$. Therefore, similar to the previous scenario, we compare the risk score as a function of the number of features $m$ for a low budget ($\lambda = 50$) and high budget ($\lambda = 200$). For CE, we set the number of samples, $\eta = 1000$ and $T_{\max} = 500$; and compare CE with and VGA algorithms. Fig.~\ref{fig:cebudget} plots the results when the input provided was the ranked feature vector. CE is seen to outperform VGA under both low and high budgets. The CE algorithm consistently achieves lower risk score with increasing features, in comparison to VGA. 


Under low budget constraint, we can clearly see that the CE algorithm outperforms the greedy method for all feature lengths $m < 70$.  
We observe a similar trend in Fig.~\ref{fig:cebudget_rand}, when randomly ordered features are provided as input. In this scenario, we evaluate the results for five different randomly ordered feature sequences, and show the standard deviation of the risk scores obtained. We observe that the risk score's standard deviation is high for small feature length, which is quite intuitive as the small number of features (given as input) can be significantly different with every random sample. Thus the risk score also varies with different input features provided. We also note that the standard deviation  of the risk scores of the CE algorithm is always low (and extremely low for the the high-budget scenario), showing the utility of the proposed algorithm.


\begin{figure}
\centering
\includegraphics[width = \linewidth]{FinalFigures/cegreedy_ranked4.pdf}
\caption{Risk score $R(\mathbf{v})$ as a function of the number of features $m$ with discrete $\lambda$ for CE and VGA algorithms for \textbf{ranked features}}
\label{fig:cebudget}
\vspace{-0.3cm}
\end{figure}



\begin{figure}
\centering
\includegraphics[width = \linewidth]{FinalFigures/cegreedy_random.pdf}
\caption{Risk score $R(\mathbf{v})$ as a function of the number of features $m$ with discrete $\lambda$ for CE and VGA algorithms for \textbf{randomly ordered features}}
\label{fig:cebudget_rand}
\vspace{-0.3cm}
\end{figure}








\noindent \textbf{Comparison of CE algorithm and VGA under varying budget:}
Finally, we compare the risk scores for the CE algorithm for varying feature budgets. The features are provided in ranked order.
Since we now have two parameters to vary ($m$ and $\lambda$), for this set of experiments, we keep the numner of samples $\eta = 250$ and the number of iterations $T_{\max} = 50$ for the CE algorithm.

We also conducted a similar experiment for VGA.  In order to have a visual comparative analysis of the performance of the CE algorithm with the greedy algorithm, we take the difference between their risk scores for each $m$ and $\lambda$. The resultant 3D graph for $R_{\texttt{CE}}(\mathbf{v}) - R_{\texttt{VGA}}(\mathbf{v})$ is plotted in Fig.~\ref{fig:diff_budget3d}. The part of the graph which is greater than 0 (above the $z=0$ plane) denotes the conditions where the cross-entropy method has greater risk than the greedy algorithm. %
The risk score achieved by CE is lower than that of VGA for most cases, except when both, the budget is very low, and the number of features is high. 
For large $m$ and a small $\lambda$, the CE algorithm rejects most of the initial samples, prohibiting the algorithm to select the optimal set of features possible, thus resulting in a high risk score. However, on increasing the budget, the risk score is seen to decrease.
Besides, as seen in the previous set of experiments, the performance of CE can be improved by increasing the number of samples and the number of iterations. 



\begin{figure}
\centering
\includegraphics[width=0.95\linewidth]{FinalFigures/diff_budget_new_to235_plane.pdf}
\caption{$R_{\texttt{CE}}(\mathbf{v}) - R_{\texttt{VGA}}(\mathbf{v})$ as a function of the number of features $m$  and the feature budget $\lambda$}
\label{fig:diff_budget3d}
\vspace{-0.3cm}
\end{figure}


\section{Related works}\label{Sec:related}


Fingerprinting or identification of hosts, operating systems, etc. has been a problem of interest for many years now (e.g., see~\cite{shamsi2014hershel}). Similarly, there have been works on fingerprinting devices; for instance 
Gao {\em et al.}~\cite{gao2010passive} proposed to identify wireless access points by applying a wavelet-based approach on frame arrival time differences. With the emergence of IoT devices, research works have looked into inferring different aspects of IoT devices from their network traffic~\cite{maiti2017link, apthorpe2017spying}. For example, the work in~\cite{apthorpe2017spying} reveals private and sensitive information of users at a smart home are leaked by analyzing only the encrypted network traffic. 

Identification of IoT devices is of use to different entities, such as enterprises and ISPs. For an enterprise, identification helps in asset tracking as well as for securing the devices from potential vulnerability exploitation; whereas, for ISPs, knowing the type of different devices connected to its network might help in mitigation of large-scale attacks. For instance, Yang {\em et al.}~\cite{yang2019towards} generated fingerprints of devices using neural algorithms, which were used to discover millions of devices connected to a network. 


Recent works have been exploring machine learning models for IoT device fingerprinting.  A common approach is to employ supervised machine learning, to train a single classifier for each device~\cite{Meidan2017ProfilIoT,miettinen2017iot}; and these set of classifiers are subsequently used for predicting a traffic session (or sequence of first $n$ packets) of a device. If a device traffic is predicted to be of multiple types (or classes), then another metric (such as edit distance~\cite{edit-distance-1964}) is used to break the tie. 
In our previous work~\cite{DEFT-2019}, we combined both supervised and unsupervised approaches for identifying known devices as well as grouping unknown devices of the same type across different networks. We explored large number of features extracted from network traffic, and the empirical study based on 16 IoT devices confirmed with other works that, devices can be classified based on network traffic features with high accuracy. 


There are also works that consider the cost of including different features, although not in the IoT settings (e.g.,~\cite{ma2008penalized,tavenard2016cost}). We briefly discuss a few here. In~\cite{zhang2015multi} the authors used a particle swarm optimization approach for cost-based feature selection in order to get a Pareto front of non-dominated solutions which gives both high accuracy and low cost. Among other works, Min~{\em et. al} formulated  cost-constrained feature selection as a CSP (constraint satisfaction problem)~\cite{min2014feature}; while they proposed a heuristic solution, they assumed decisions as given, and the notion of risk was not considered within the scope of the work. There are also works that included a cost evaluation function as part of existing feature selection approaches; for example, two filter-based feature selection methods were experimented in~\cite{bolon2014framework}.
Another interesting approach was proposed in~\cite{maldonado2017cost}, wherein the cost of feature selection was made part of the SVM classification model for a credit scoring application.




However, to the best of our knowledge, previous works did not consider the cost of different features that are extracted for training and classifying devices. Neither did they consider the {\em risk} of misclassification, in particular when the loss due to wrong classification can be different depending on the classification result. We consider these important aspects in our work here.





\section{System model}
\label{sec:system}




We describe the system model in this section. Below, we first state the assumptions of the system, and then go on to  define misclassification loss as well as our concept of risk.
The commonly referred notations are listed in Table~\ref{tab:def}. 

\begin{table}[]
\centering
\caption{\hlt{Table of Notations}}
\label{tab:def}
\begin{tabular}{|l|l|}
\hline
Term & Definition \\ \hline \hline
$\mathbf{f}$ & \begin{tabular}[c]{@{}l@{}}Feature vector, constituting of the features abstracting \\ network traffic characteristics. $|\mathbf{f}| = m$\end{tabular} \\ \hline
$\mathbf{c}$ & \begin{tabular}[c]{@{}l@{}}Cost vector, constituting of the costs associated\\  with each of the $m$ features. $|\mathbf{c}| = m$\end{tabular} \\ \hline
$\mathbf{v}$ & \begin{tabular}[c]{@{}l@{}}A binary vector denoting the features selected \\ from $\mathbf{f}$. $|\mathbf{v}| = m$\end{tabular} \\ \hline
$\lambda$ & Feature budget \\ \hline
$R(\mathbf{v})$ & \begin{tabular}[c]{@{}l@{}}Risk score, based on the selected\\ feature vector $\mathbf{v}$\end{tabular} \\ \hline
$\mathcal{L}$ & \begin{tabular}[c]{@{}l@{}}Loss Matrix where element $l_{i,j}$ represents\\ misclassification loss. $|\mathcal{L}| = n \times n$\end{tabular} \\ \hline
\end{tabular}
\end{table}


\begin{comment}


\begin{table}[]
\centering
\caption{Table of Notations }
\label{tab:def}
\resizebox{\linewidth}{!}{%
\begin{tabular}{|c|l|}
\hline
 Term & Definition \\ \hline \hline
 
 
 $\mathbf{f}$ & Feature vector, constituting of the features abstracting network traffic characteristics. $|\mathbf{f}| = m$ \\ \hline
 $\mathbf{c}$ & Cost vector, constituting of the costs associated with each of the $m$ features. $|\mathbf{c}| = m$ \\ \hline
 $\mathbf{v}$ & A binary vector denoting the features selected from $\mathbf{f}$. $|\mathbf{v}| = m$ \\ \hline
$\lambda$ & Feature budget \\ \hline

 $R(\mathbf{v})$ & Risk score, based on the selected feature vector $\mathbf{v}$ \\ \hline
 $\mathcal{L}$ & Loss Matrix where element $l_{i,j}$ represents misclassification loss. $|\mathcal{L}| = n \times n$\\ \hline



\end{tabular}%
}
\end{table}


\end{comment}


\subsection{Traffic and features}
 Let $\mathcal D$ denote the set of devices considered for classification in a network; furthermore, let $n = |\mathcal D|$. Corresponding to a device $d \in \mathcal D $, assume $q_d$ units of traffic have been stored and made available for the purpose of device classification. For simplicity, and without loss of generality, we assume $q$ units of traffic are available for every device. Let $\mathcal T^d$ denote the traffic of device $d$. Each unit of traffic (a connection or a session) of $\mathcal T^d$ is processed to extract the {\em feature vector}. We use $\mathbf f$ to represent the $m$ features of interest in our problem: 
  $$ \mathbf f = [f_1, f_2, \dots, f_m ]; \quad \mathbf f \in \mathbb{R}^{m}.$$
 
Denote by $\mathcal E(.)$ the function for feature extraction. For an input traffic data $\mathcal T^d$ of device $d$, $\mathcal E(\mathcal T^d)$ produces a matrix $\mathcal S^d$ of extracted features corresponding to all units of traffic:

\begin{equation}
    \mathcal{S}^d : = \left[\begin{array}{cccc}{s^d_{1,1}} & {s^d_{1,2}} & {\cdots} & {s^d_{1,m}} \\ {s^d_{2,1}} & {s^d_{2,2}} & {\cdots} & {s^d_{2,m}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {s^d_{q,1}} & {s^d_{q,2}} & {\cdots} & {s^d_{q,m}}\end{array} \right],
\end{equation}
\noindent where $s^d_{j, k}$ denotes the $k^\text{th}$ feature extracted from the $j^\text{th}$ traffic unit of device $d$. 
To generalize, we use $\mathbf s$ to denote an extracted feature vector of a device. Obviously, $\mathbf s$ is of length $m$. 
Note that, $\mathcal S^d, \forall d \in \mathcal D$ form the dataset $\mathbb{X}$ used for training and testing the classification model. We use $\mathbb{X}^\texttt{train}$ and  $\mathbb{X}^\texttt{test}$ to denote the partition of the dataset for, training and testing, respectively. 

\subsection{Cost of feature extraction}\label{sec:cost}

Features  characterise different aspects of the network traffic, and are obtained by processing network traffic. Therefore, feature extraction involves cost in terms of resources required for processing and storing the feature, and sometimes for even purchasing the feature. 
With a slight abuse of notation, we use $\mathcal E_{f_k}(.)$ the function to extract feature $f_k, 1 \le k \le m$. Therefore the cost of extracting a feature $f_k$ will be denoted by:
$$ c_k = \texttt{cost}(\mathcal E_{f_k}); \quad \forall k \in [1, 2, \dots, m],$$
 \noindent and the cost vector is defined as:
$$\mathbf c = [c_1, c_2, \dots, c_m]; \quad \mathbf{c} \in  \mathbb{R}_{+}^{m}.$$
  
\subsection{Supervised classification of devices}

 We consider supervised machine learning approaches for IoT classification, wherein traffic data is labelled and provided for training the classifier. The trained model, denoted by $\mathcal M$, is subsequently used in an operational environment to differentiate devices into different types.
 
 
Given the matrices of extracted features $\mathcal S^d, \forall d \in \mathcal D$, the classification model $\mathcal M$ maps $\mathbf s$ to one of the devices types $\{d_1, d_2, \dots, d_n\}$. The output of a machine learning model $\mathcal M$ is usually represented using an $n \times n$ confusion matrix $\mathcal X_{\mathcal M}$ that is processed to obtain relevant performance metrics, such as precision and recall. 
\begin{comment}
\begin{equation}
    \mathcal{X} : = \left[\begin{array}{cccc}{x_{1,1}} & {x_{1,2}} & {\cdots} & {x_{1,n}} \\ {x_{2,1}} & {x_{2,2}} & {\cdots} & {x_{2,n}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {x_{n,1}} & {x_{n,2}} & {\cdots} & {x_{n,n}}\end{array} \right],
\end{equation}
\noindent such that $x_{i,j}$ represents the number of
\end{comment}
The element in the $i^\text{th}$ row and $j^\text{th}$ column of the confusion matrix, $x_{i,j}$, represents the number of data points (or instances) of class $j$ that were predicted to be of class $i$ by the classifier. Given the confusion matrix, we now define $p_{i,j}$, the probability of misclassification, as the probability that an instance of class $j$ is predicted as an instance of class $i$. That is if $\widehat{d}$ is the predicted device and $d$ is the actual device, then we define the probability of classifying $d = j$ as $\widehat{d} = i$ is given as:


\begin{align}
\label{eq:performance}
    {\widehat{p}_{i,j}} = \operatorname{Pr}(\widehat{d}=i | d=j) & 
    =
     \frac{x_{i,j}}{\sum_{k=1}^{n} x_{k,j}}. &
\end{align}

 
For a classification model $\mathcal M$, we define an $n \times n$ misclassification matrix $\hat{\mathcal{P}}_\mathcal{M}$ with elements $\widehat{p}_{i,j}$'s. For simplicity, henceforth, we drop the subscript and denote the misclassification matrix as $\mathcal P$.

\begin{comment}
 \begin{equation}
\hat{\mathcal{P}} :=\left[\begin{array}{cccc}{\widehat{p}_{1,1}} & {\widehat{p}_{1,2}} & {\cdots} & {\widehat{p}_{1,n}} \\ {\widehat{p}_{2,1}} & {\widehat{p}_{2,2}} & {\cdots} & {\widehat{p}_{2,n}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {\widehat{p}_{n,1}} & {\widehat{p}_{n,2}} & {\cdots} & {\widehat{p}_{n,n}}\end{array}\right].
\end{equation}
\end{comment}

\subsection{Misclassification loss} 
\label{subsec:loss}





Assume that the classifier is tasked to predict the class of a device $d_i$. We represent a device $d_{i}$ as as an ordered pair~$<$~\texttt{type}, \texttt{brand}~$>$, where \texttt{type}  indicates the type of the device, e.g., a camera or a speaker, while \texttt{brand} represents the brand of the item, e.g., Sony or Samsung. When the classification model (mis)classifies a device $d_i = <t_i, b_i>$ as $d_j = <t_j, b_j>$, there is a loss incurred; we denote this misclassification loss as $l_{i,j}$, 
and the corresponding loss matrix $\mathcal L$:

\begin{equation}
\mathcal{L} :=\left[\begin{array}{cccc}{0} & {l_{1,2}} & {\cdots} & {l_{1,n}} \\ {l_{2,1}} & {0} & {\cdots} & {l_{2,n}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {l_{n,1}} & {l_{n,2}} & {\cdots} & {0}\end{array}\right].
\end{equation}

This loss matrix is a user-defined parameter and depends on the types of devices considered for classification. While this allows for a generic definition of the misclassification losses as deemed right by the user, we later provide (in Section~\ref{subsec:loss-eval}) the specific definition used in this work for the purpose of illustration as well as for experimentation. 


\subsection{Cyber Risk}
\label{subsec:risk}


Using these definitions, we define the cyber risk of the model.
While $\mathbf f$ is the vector of features of interest, not all features might be available for modeling, given features have cost, and a solution is constrained by the available budget. We use a binary vector $\mathbf{v}$ to indicate whether a feature is selected or not. 
\[\mathbf{v}=\left[v_{1}, \dots, v_{m}\right],\] 
where
\[v_{k}=\left\{\begin{array}{ll}{1,} & {\text {if the } k^{\text {th }} \text { feature selected }} \\ {0,} & {\text {if the } k^{\text {th }} \text { feature is not selected }}\end{array}\right.\]

Next, we define cyber risk. Given a classifier $\mathcal M$, loss matrix $\mathcal L$, feature budget $\lambda$, the extracted data (for the  selected features indicated by $\mathbf v$) $\mathbb X_\mathbf{v}$, and a feature selection vector $\mathbf v$, the cyber risk score $R(\mathbf{v},\mathcal{M}, \mathcal{L}, \lambda, \mathbb X_\mathbf{v})$ is the expected sum of the losses due to misclassifications. Since we are considering the risk score for a fixed classifier model, for a fixed number of devices and traffic data, we may simplify the risk score by denoting it as $R(\mathbf{v})$.
Formally, the \textbf{cyber risk score} or the \textbf{risk score}, based on the set of selected features  (captured by $\mathbf{v}$), is defined as the sum of the product of the probability of misclassification 
and the losses associated with misclassification of each device ($l_{i,j}$), which is given by:
\begin{align}
    R(\mathbf{v}) &= \sum_{i=1}^{n} \sum_{j=1}^{n} \operatorname{Pr}(\widehat{d}=i | d=j; \mathbf{v}) \times l_{i,j} 
    \label{eq:cyber_risk}
\end{align}














\section{Problem Definition - Optimal Feature Selection under Budget Constraint} \label{sec:problem}


This work looks into the problem of minimizing the cyber risk associated with classification of IoT devices, by selecting the optimal set of features from a universal feature set, whilst keeping the feature cost less than the budget available.  \\

\noindent {\bf Problem Definition:} {\em Given the misclassification matrix $\mathcal{P}$ for a given feature vector $\mathbf{f}$, and the loss matrix $\mathcal{L}$, we aim to find the optimal feature vector (via $\mathbf{v}$) which minimizes the cyber risk score $R(\mathbf{v})$, under a budget constraint $\lambda$.} \\

\noindent Thus, the optimization problem may be defined as:
\begin{comment}
\begin{equation}
    \begin{aligned} \widehat{\mathbf{v}}=& \underset{\mathbf{v} \in\{0,1\}^{m}}{\arg \min } R(\mathbf{v}) \\  \text { s.t } & \mathbf{c}^{T}\mathbf{v} \le \lambda \end{aligned}
\end{equation}
\end{comment}


\begin{equation}
\label{eq:optim}
    \begin{aligned} \mathbf{v}_{s}=& \underset{\mathbf{v} \in\{0,1\}^{m}}{\arg \min } R(\mathbf{v}) \\  \text { s.t } & \sum_{j=1}^{m} {c}_{j} v_{j}\le\lambda,
    \end{aligned}
\end{equation}


\noindent where $R(\mathbf{v})$ is as defined in Eq.~\ref{eq:cyber_risk}. 


The optimization problem presented in Eq. \ref{eq:optim} has some unique properties making it a non-trivial problem to solve. Firstly, the cyber risk function $R(\mathbf{v})$ used as the principle objective function, is not an analytic function. Hence, using methods to solve continuous optimization problems, like the gradient descent is not feasible in this scenario. This is mainly because the gradient descent uses the Newton's Method to find the optimal solution which in turn uses a well-defined continuous objective function; and $R(\mathbf{v})$ is not a continuous function. 
Moreover, another important constraint to use such methods requires the objective function, the constraint equations and the domain space to be convex in nature. In our scenario, we do not know whether we have a closed-form analytical expression for the risk score which is the objective function. This is because it is dependent on the classification model and the features selected. Though a probabilistic model for a simple classifier like the Naive-Bayes classifier can be obtained, we cannot say the same for any general classifier, for example the Decision tree or the Random Forest classifiers. Thus, in a general sense, it is not possible to get an analytical expression for the risk score defined in Eq. \ref{eq:cyber_risk}. The best we can do is to obtain a point-wise evaluation of the objective function which is the minimum requirement to attempt a solution to this optimization problem.
We observe that the constraint functions are affine in nature, and since all affine functions are convex, the constraint functions are also convex in nature. This leaves us with the domain space ($m$-dimensional binary subspace) for the optimization problem.
\begin{comment}


\begin{lemma}
A set of $K$-dimensional Binary subspace is not convex in nature
\end{lemma}
\textbf{Proof:} Let us consider a $m$-dimensional binary subspace represented by $\Phi \in \{0,1 \}^m$ is convex in nature. That implies that the line segment between any two points in the set $\Phi$ , i.e., $\lambda_1, \lambda_2 \in \Phi$  and any $\theta$ with $0 \le \theta \le 1$ we will have
\[ \theta \lambda_1 + (1-\theta) \lambda_2 \in \Phi \]

Let us consider the case when we select  the maximum number from the $m$-dimensional binary subspace as one of the points to be considered: 
$ \overbrace{1111\dots111}^{m- \text{digits}}$
The second point is chosen as any point other than the vector constituting of all zeroes. We then take the linear combination of the two points with $\theta \neq 0 \hspace{0.1 cm} \text{or} 1$
Then the resultant point will be in the $(m+1)$-dimension binary subspace which does not belong to $\Phi$, the $m$-dimensional binary space of our problem.\\
Hence, we can conclude that a $m$-dimensional binary subspace $\Phi$ is not a convex set.\\
\end{comment}
\par Since the domain space is not  convex in nature, we can safely conclude that the problem is not a convex optimization problem, and cannot be solved by using convex optimization tools. Therefore, the optimal way to solve such problems is via a combinatorial search over all possible combinations of features by the brute force method which is not a very efficient way.  In classes of problem such as the one posed in this paper,a subsequent brute force enumeration of all possibilities is computationally prohibitive and as such a stochastic search based approach needs to be developed in practice. In the next section, we explore the cross entropy technique for searching the entire subspace in an efficient way to potentially converge to the global optima.


